Maximum Likelihood Estimation of a Dependent Variable Model with ml

[Sebastian Kaumanns]

Dear all,

I am currently trying to estimate the parameters of a limited dependent variable model, as introduced by Lesmond et al. (1999) (http://rfs.oxfordjournals.org/conten.../1113.abstract). I attached the model and the log-likelihood function (as jpgs) below.

I have the data for R_jt (R) and R_mt (MR) and would like to estimate alpha₁, alpha₂, beta, and sigma. I am using Stata 13 on Windows 7 and tried the following:

program myprog
version 13.1
args lnfj a1 a2 b sigma

quietly replace `lnfj' = ln((2*c(pi)*`sigma')^(-0.5)))-((R+`a1'-`b'*MR)^2)/(2*(`sigma')^2) if R<0

quietly replace `lnfj' = ln((2*c(pi)*`sigma')^(-0.5)))-((R+`a2'-`b'*MR)^2)/(2*(`sigma')^2) if R>0

quietly replace `lnfj' = ln(normalden((`a2'-`b'*MR)/`sigma')-normalden((`a1'-`b'*MR)/`sigma') if R=0
end

ml model lf myprog (`a1'=`b'*MR R) (`a2'=`b'*MR R) (`b'=R`a1' MR) ()

Please feel free to laugh heartily, as it is my first attempt at ml or any Stata programming at all. Right now, I am getting the following error:

. ml check

Test 1:  Calling myprog to check if it computes log likelihood and
         does not alter coefficient vector...
         FAILED; myprog returned error 132.

Here is a trace of its execution:
------------------------------------------------------------------------------
-> myprog __00000E __00000F __00000G __00000H __00000I
        - `begin'
        = capture noisily version 13.1: myprog __00000E __00000F __00000G __00000H __00000I
          ---------------------------------------------------------------------------------------------------------------------------------------------------------------- begin myprog ---
          - version 13.1
          - args lnfj a1 a2 b sigma
          - quietly replace `lnfj' = ln((2*c(pi)*`sigma')^(-0.5)))-((R+`a1'-`b'*MR)^2)/(2*`sigma') if R<0
          = quietly replace __00000E = ln((2*c(pi)*__00000I)^(-0.5)))-((R+__00000F-__00000H*MR)^2)/(2*__00000I) if R<0
too many ')' or ']'
          ------------------------------------------------------------------------------------------------------------------------------------------------------------------ end myprog ---
        - `end'
        = set trace off
------------------------------------------------------------------------------
Fix myprog.

I feel rather allright about my translation of the log-likelihood function into the program (but are happy about corrections). But for the model, I feel like I have no idea what I am doing. I read that the model needs to contain a separate equation (as identified by brackets) for every paramater that needs to be estimated, so that is what I tried to do. Any help is greatly appreciated!

[Scott Merryman]

You have a syntax error in the third equation (it is R==0 not R=0). Also, you have too many parentheses in the first and second equation and too few in the third.

The code and example below appears to work

clear
set seed 93921
set obs 100
gen MR =runiform()
gen Rstar = 2*MR + rnormal()
gen R = Rstar -1 if Rstar<1
replace R = 0 if  Rstar>1 & Rstar < 5
replace R = Rstar - 5 if Rstar > 5


capture program drop myprog
program myprog
version 13.1
    args lnfj a1 a2 b sigma
    quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R+`a1'-`b'*MR)^2)/(2*(`sigma'^2)^2) if R<0
    quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R+`a2'-`b'*MR)^2)/(2*(`sigma'^2)^2) if R>0
    quietly replace `lnfj' = ln(normalden((`a2'-`b'*MR)/`sigma')-normalden((`a1'-`b'*MR)/`sigma')) if R==0
end

ml model lf myprog (`a1'=`b'*MR R) (`a2'=`b'*MR R) (`b'=R`a1' MR) ()
ml max, difficult

[Sebastian Kaumanns]

Dear Scott,

sorry that I am a bit late. Your code worked great for me! Thank you very much!

[Sebastian Kaumanns]

Dear all:
I still struggle with my ML problem. I adjusted the code recommended by Scott Merryman for my empirical setup, i.e., the estimation of total trading costs based on a limited dependent variable model as introduced by Lesmond et al. (1999, http://rfs.oxfordjournals.org/conten.../1113.abstract).
I used the following final stata code:

    capture program drop myprog
    program myprog
    version 13.0
        args lnfj a1 a2 b sigma
        quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R1+`a1'-`b'*MR1)^2)/(2*(`sigma'^2)^2) if R1<0
        quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R1+`a2'-`b'*MR1)^2)/(2*(`sigma'^2)^2) if R1>0
        quietly replace `lnfj' = ln(normalden((`a2'-`b'*MR1)/`sigma')-normalden((`a1'-`b'*MR1)/`sigma')) if R1==0
    end
  
    set more off
    qui ml model lf myprog (R1 = `a1' `b'*MR1   ) (R1 = `a2'`b'*MR1 ) () ()  
    qui ml max ,  difficult   iterate(10)  nowarning  nooutput

Based on simulated return data, my programming seems to work very well. However, when applying it to my "real world" data, i.e., daily firm return and market return data for different markets and time periods between 2001 and 2014, I often receive the following error message: "could not find feasible values" r(491).

I read related topics where fellow Stata-users had the same error message and it seems that it is often caused by errors in the likelihood function. Does anyone recognize any flaws in my likelihood function? In a related vein, is it possible that data structure might (also) cause such error messages? So far I did not come across any strange data structure in those instances of the error message.

Best
Sebastian

For completeness, you will find my complete simulation programming below as well:

*******************************************************************************
******************** TOTAL TRADING COSTS (Oct. 2015) *************
*******************************************************************************


***Step (1): Generate Dataset: 25 return times series (varying across different a1 and a2 values)



forvalues i=1(1)25 {

clear
set seed 18973
set obs 500
gen MR =runiform()
gen Rstar = 2*MR + rnormal()

set more off
macro dir
macro drop Aplha1
macro drop Aplha2
global Aplha1 "`i'/10"
global Aplha2 "-`i'/10"

***

gen help1 = 1 if Rstar>$Aplha1
replace help1 = 1 if  Rstar<$Aplha2
replace help1 = 0 if help1==.

gen help2=help1*Rstar
gen R= help2-$Aplha1 if Rstar>$Aplha1
replace R= help2+$Aplha2 if Rstar<$Aplha2
replace R=0 if R==.
drop help1 help2

gen ZeroR=1 if R==0
replace ZeroR=0 if ZeroR==.

egen meanZeroR=mean(ZeroR)

gen Round=`i'


save data`i', replace


di `i'
}

***

forvalues i=2(1)25 {

set more off
use data1, clear
append using data`i'
save data1, replace

di `i'
}


******************************************************************************


***Step (2): Alpha_1 estimation: LOOP ESTIMATION

use data1, clear

set more off
capture drop a1 a2
egen group=group(Round)
gen a1=.

capture program drop myprog
program myprog
version 13.0
    args lnfj a1 a2 b sigma
    quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R1+`a1'-`b'*MR1)^2)/(2*(`sigma'^2)^2) if R1<0
    quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R1+`a2'-`b'*MR1)^2)/(2*(`sigma'^2)^2) if R1>0
    quietly replace `lnfj' = ln(normalden((`a2'-`b'*MR1)/`sigma')-normalden((`a1'-`b'*MR1)/`sigma')) if R1==0
end

***

forvalues i=1(1)25 {


gen R1=R
replace R1=. if group!=`i'

gen MR1=MR
replace MR1=. if group!=`i'


***
set more off
ml model lf myprog (R1 = `a1' `b'*MR1  ) (R1 = `a2'`b'*MR1 ) () ()  
ml max ,  difficult   iterate(100)  nowarning  nooutput

***

qui replace a1=_b[_cons] if group==`i'

drop R1 MR1

di `i'
}


save dataALPHA1, replace


******************************************************************************


***Step (3): Alpha_2 estimation: LOOP ESTIMATION

use data1, clear

set more off
capture drop a1 a2
egen group=group(Round)
gen a2=.

capture program drop myprog
program myprog
version 13.0
    args lnfj a2 a1  b sigma
    quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R1+`a2'-`b'*MR1)^2)/(2*(`sigma'^2)^2) if R1>0
    quietly replace `lnfj' = ln((2*c(pi)*`sigma'^2)^(-0.5))-((R1+`a1'-`b'*MR1)^2)/(2*(`sigma'^2)^2) if R1<0
    quietly replace `lnfj' = ln(normalden((`a2'-`b'*MR1)/`sigma')-normalden((`a1'-`b'*MR1)/`sigma')) if R1==0
end

***

forvalues i=1(1)25 {


gen R1=R
replace R1=. if group!=`i'

gen MR1=MR
replace MR1=. if group!=`i'


***
set more off
ml model lf myprog  (R1 = `a2'`b'*MR1 ) (R1 = `a1' `b'*MR1  ) () ()  
ml max ,  difficult   iterate(100)  nowarning  nooutput

***

qui replace a2=_b[_cons] if group==`i'

drop R1 MR1

di `i'
}


save dataALPHA2, replace

******************************************************************************


***Step (4): Merge Alpha_1 and Alpha_2 data and estimate total trading costs (TTC), and generate correlation matrix for TTC, proportion of zero return days (meanZeroR), and theoretical alpha values (ALPHA)


use dataALPHA1, clear
merge using dataALPHA2

***

gen TTC=a1-a2


gen ALPHA=.
forvalues i=1(1)25 {
set more off
replace ALPHA=  `i'/10 if group == `i'
di `i'
}

***

spearman TTC meanZeroR a1 a2 ALPHA, stats( rho p obs)