 |
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]
st: problem with MLE
I am trying to estimate a fertility model with birth history data using maximum likelihood. The data comes from a survey (DHS) and I am using Stata 8.1.
When I run *ml check*, I get the message that all tests have been passed. Then, once the optimization starts, the program spits out the error message that "numerical derivatives could not be calculated" and stops. Any suggestions on how to go about finding the error, correcting it and making the program work?
In Stata language, my likelihood evaluator is of the type *lf*; and so I have not programmed the derivative or the hessian. Below I give the likelihood evaluator program and also the log file containing the output.
Any suggestions to get this working will be most welcome.
Deepankar
--------- Start: Likelihood function -------------
clear
set more 1
capture program drop step1mle_lf
program step1mle_lf
version 8.1
#delimit ;
args lnf theta1 theta2 theta3 theta4 theta5 theta6 theta7 theta8 theta9
theta10 theta11 theta12 theta13 theta14;
/* theta4 = generates probability of male birth (will be estimated as a parameter) */
tempvar q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 x1 x2 x3 x4 x5 x6
p21 p31 p32 p41 p42 p43 p51 p52 p53 p54 p61 p62 p63 p64 p65;
#delimit cr
quietly { /*--------------- BEGINNING OF QUIET COMPUTATION ---------------*/
/* ------------ GENERATING THE CONDITIONAL PROBABILITIES ------------------- */
gen double `q1' = 1/(1+exp(`theta5'))
gen double `q2' = exp(`theta5')/(1+exp(`theta5'))
gen double `q3' = 1/(1+exp(`theta6')+exp(`theta7'))
gen double `q4' = exp(`theta6')/(1+exp(`theta6')+exp(`theta7'))
gen double `q5' = exp(`theta7')/(1+exp(`theta6')+exp(`theta7'))
gen double `q6' = 1/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))
gen double `q7' = exp(`theta8')/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))
gen double `q8' = exp(`theta9')/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))
gen double `q9' = exp(`theta10')/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))
gen double `q10' = 1/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))
gen double `q11' = exp(`theta11')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))
gen double `q12' = exp(`theta12')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))
gen double `q13' = exp(`theta13')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))
gen double `q14' = exp(`theta14')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))
/* -------------------------------------------------------------------------------- */
/* ----------------------GENERATING THE PROBABILITY TABLE ------------------------ */
/* p`i'`k' IS THE CONDITIONAL PROBABILITY OF OBSERVING THE BIRTH SEQUENCE */
/* GIVEN THAT THE DESIRED FAMILY SIZE IS `i' AND THE DESIRED TARGET FOR BOYS IS `k' */
/* STEP 1 */
forval i=2/6 {
local m = `i'-1
forval j=1/`m' {
gen double `p`i'`j'' = 0
}
}
/* STEP 2 */
tempvar cboys c1 c2 c3 c4 c5 c6
gen byte `cboys' = 0
gen `c1' = ($ML_y5==1)
replace `cboys' = `cboys' + `c1'
replace `c1' = `c1'*`cboys'
gen `c2' = ($ML_y6==1)
replace `cboys' = `cboys' + `c2'
replace `c2' = `c2'*`cboys'
gen `c3' = ($ML_y7==1)
replace `cboys' = `cboys' + `c3'
replace `c3' = `c3'*`cboys'
gen `c4' = ($ML_y8==1)
replace `cboys' = `cboys' + `c4'
replace `c4' = `c4'*`cboys'
gen `c5' = ($ML_y9==1)
replace `cboys' = `cboys' + `c5'
replace `c5' = `c5'*`cboys'
gen `c6' = ($ML_y10==1)
replace `cboys' = `cboys' + `c6'
replace `c6' = `c6'*`cboys'
/*---- PROBABILITY OF BIRTH OF A MALE OFFSPRING ----*/
tempvar pmale
gen `pmale' = 1/(1+exp(`theta4'))
forval i=1/6 {
forval j=1/6 {
forval l=1/6 {
local m=`i'-1
forval k=1/`m' {
replace `p`i'`k''=((`pmale')^`j')*((1-`pmale')^`l') if ($ML_y1==`i' & $ML_y4==`j' & $ML_y11==`l')
}
}
}
}
/* STEP 3 */
forval i=1/6 {
forval j=1/6 {
forval l=1/6 {
local m = `i'-1
forval k=1/`m' {
replace `p`i'`k''=0 if ($ML_y1==`i' & $ML_y2==`j' & $ML_y1==$ML_y2 & `c`j''!=`k' & $ML_y4>=`k')
replace `p`i'`k''=0 if ($ML_y1==`i' & $ML_y2==`j' & $ML_y1>$ML_y2 & `c`j''!=`k')
replace `p`i'`k''=((`pmale')^`j')*((1-`pmale')^`l') if ($ML_y1==`i' & $ML_y4==`j' & $ML_y11==`l' &
$ML_y1==$ML_y2 & $ML_y2==0)
}
}
}
}
/* ----------------- END OF THE PART THAT GENERATES THE PROBABILITY TABLE --------------------------- */
gen double `x1' = (`pmale'^($ML_y4))*((1-`pmale')^($ML_y11))
gen double `x2' = exp(lnfact($ML_y1))
gen double `x3' = `x1'*`x2'
gen double `x4' = ((exp(`theta1'))^($ML_y1))/exp(exp(`theta1'))
gen double `x5' = norm(-`theta3')
gen double `x6' = ((exp(`theta1'+`theta2'))^($ML_y1))/exp(exp(`theta1'+`theta2'))
} /* ----------------- END OF QUIET COMPUTATION ---------------- */
#delimit ;
quietly replace `lnf' = ln( ($ML_y3*(`x1')*(1/`x2')*`x4'*`x5') + ((1-`x5')*`x6'*(1/`x2')*( `p21' + `p31'*`q1'
+ `p32'*`q2' + `p41'*`q3' + `p42'*`q4' + `p43'*`q5' + `p51'*`q6' + `p52'*`q7'
+ `p53'*`q8' + `p54'*`q9' + `p61'*`q10' + `p62'*`q11' + `p63'*`q12'
+ `p64'*`q13' + `p65'*`q14')));
#delimit cr
end
---------- End: Likelihood function -------------
---------- Start: Log file ----------------------
log: c:\deepankar\gender-data\step1\india-mle-step1-results-1998.smcl
log type: smcl
opened on: 2 Dec 2007, 16:55:03
.
. set mem 100m
Current memory allocation
current memory usage
settable value description (1M = 1024k)
--------------------------------------------------------------------
set maxvar 5000 max. variables allowed 1.733M
set memory 100M max. data space 100.000M
set matsize 400 max. RHS vars in models 1.254M
-----------
102.987M
. use india-step1-mle-data-1998.dta /* DATSET: Contains observations with N=2,3,4,5,6 */
.
.
. /*---------- DETAILS OF SURVEY DESIGN -------------*/
. gen smpwt = v005/1000000
. svyset [pweight=smpwt]
pweight is smpwt
.
.
. /* ------------------ UNRESTRICTED MODEL ---------------------------*/
. #delimit ;
delimiter now ;
. ml model lf step1mle_lf (fsize: dfsize = edu age edu rur work middle poor hindu, nocons) (alpha:)
> (target: alive dfsize1 nboy cord1 cord2 cord3 cord4 cord5 cord6 ngirl = age edu rur work middle poor contra jfam hindu
> )
> /four /five /six /seven /eight /nine /ten /eleven /twelve /thirteen /fourteen, svy robust tech(bhhh nr);
note: edu dropped due to collinearity
. #delimit cr
delimiter now cr
.
. ml check
Test 1: Calling step1mle_lf to check if it computes log pseudolikelihood and
does not alter coefficient vector...
Passed.
Test 2: Calling step1mle_lf again to check if the same log pseudolikelihood value
is returned...
Passed.
Test 3: Calling step1mle_lf to check if 1st derivatives are computed...
test not relevant for method lf.
Test 4: Calling step1mle_lf again to check if the same 1st derivatives are
returned...
test not relevant for method lf.
Test 5: Calling step1mle_lf to check if 2nd derivatives are computed...
test not relevant for method lf.
Test 6: Calling step1mle_lf again to check if the same 2nd derivatives are
returned...
test not relevant for method lf.
------------------------------------------------------------------------------
Searching for alternate values for the coefficient vector to verify that
step1mle_lf returns different results when fed a different coefficient vector:
Searching...
initial: log pseudolikelihood = -<inf> (could not be evaluated)
searching for feasible values +
feasible: log pseudolikelihood = -29212.563
improving initial values ..........
improve: log pseudolikelihood = -29212.563
continuing with tests...
------------------------------------------------------------------------------
Test 7: Calling step1mle_lf to check log pseudolikelihood at the new values...
Passed.
Test 8: Calling step1mle_lf requesting 1st derivatives at the new values...
test not relevant for method lf.
Test 9: Calling step1mle_lf requesting 2nd derivatives at the new values...
test not relevant for method lf.
------------------------------------------------------------------------------
step1mle_lf HAS PASSED ALL TESTS
------------------------------------------------------------------------------
Test 10: Does step1mle_lf produce unanticipated output?
This is a minor issue. Stata has been running step1mle_lf with all
output suppressed. This time Stata will not suppress the output.
If you see any unanticipated output, you need to place quietly in
front of some of the commands in step1mle_lf.
-------------------------------------------------------------- begin execution
---------------------------------------------------------------- end execution
.
. ml search
initial: log pseudolikelihood = -29212.563
improve: log pseudolikelihood = -26594.685
rescale: log pseudolikelihood = -17651.076
rescale eq: log pseudolikelihood = -15128.187
.
. ml max, difficult
initial: log pseudolikelihood = -15128.187
rescale: log pseudolikelihood = -15128.187
rescale eq: log pseudolikelihood = -15128.187
(setting optimization to BHHH)
could not calculate numerical derivatives
missing values encountered
r(430);
end of do-file
r(430);
. log close
log: c:\deepankar\gender-data\step1\india-mle-step1-results-1998.smcl
log type: smcl
closed on: 2 Dec 2007, 17:04:00
-------------- End: Log file ----------------
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/