stata ml estimation

[Wong Wan Man]

Hello everyone,

I am trying to create a clogit equivalent. Below is my program.

Note that :
1. x1 x2 x3 id choice are harded-coded variables. They have the same name as the one I loaded into the dataset.
2. I am trying to estimate beta1 to beta3, which have nonlinear effect




program myconditional_logit
args todo beta1 beta2 beta3 beta4 lnL
version 11
tempvar den p last xb

gen double `xb' = `beta1' * x1^`beta1' + `beta2' * x2^`beta2' +`beta3' * x3 ^`beta3'
local y choice
local by1 id
sort `by1'
quietly{
by `by1': egen double `den' = sum(exp(`xb'))
gen double `p' = exp(`xb')/`den'
mlsum `lnL' = `y' *log(`p') if `y'==1
if (`todo'==0 | `lnL' > =.) exit
}
end

My function call is
ml model d0 myconditional_logit () () () ()



However, when I try to run the program. It issues the following error:

myconditional_logit 0 __000009 __00000A __00000B __00000C
- `begin'
= capture noisily version 13: myconditional_logit 0 __000009 __00000A __0
> 0000B __00000C
----------------------------------------- begin myconditional_logit ---
- args todo beta1 beta2 beta3 beta4 lnL
- version 11
- tempvar den p last xb
- gen double `xb' = `beta1' * x1 + `beta2' * x2 +`beta3' * x3 +`beta4'
> * x4
= gen double __00000G = __000009 * x1 + __00000A * x2 +__00000B * x3 +_
> _00000C * x4
matrix operators that return matrices not allowed in this context
------------------------------------------- end myconditional_logit ---
- `end'


I am thinking this is because stata treat x1 as a full vector. Is there any way I can make the x1 to be observation specific?

Thanks everyone!

[Christophe Kolodziejczyk]

your syntax is for a lf evaluator

program myconditional_logit

args todo b lnf
version 11


tempvar beta1 beta2 beta3 beta4 xb den p last xb

mleval `beta1' = `b' , eq(1)
mleval `beta2' = `b' , eq(2)
mleval `beta3' = `b' , eq(3)
mleval `beta4' = `b' , eq(4)

// gen double `xb' = `beta1' * x1^`beta1' + `beta2' * x2^`beta2' +`beta3' * x3 ^`beta3'

local y choice
local by1 id
sort `by1'

by `by1': egen double `den' = sum(exp(`xb'))
gen double `p' = exp(`xb')/`den'
mlsum `lnL' = `y' *log(`p') if `y'==1
if (`todo'==0 | `lnL' > =.) exit

end

try to avoid egen. It will slow down your code since it has to call ado-files which are unnecessary.

tempvar last
by `by1' : gen `last' = _n==_N
by `by1': gen double `den' = sum(exp(`xb'))
gen double `p' = sum(cond( `y'==1,exp(`xb'),0)
replace `p' = `p'/`den'
mlsum `lnL' = `y' *log(`p') if `last'
if (`todo'==0 | `lnL' > =.) exit

[Wong Wan Man]

your syntax is for a lf evaluator

program myconditional_logit

args todo b lnf
version 11


tempvar beta1 beta2 beta3 beta4 xb den p last xb

mleval `beta1' = `b' , eq(1)
mleval `beta2' = `b' , eq(2)
mleval `beta3' = `b' , eq(3)
mleval `beta4' = `b' , eq(4)

// gen double `xb' = `beta1' * x1^`beta1' + `beta2' * x2^`beta2' +`beta3' * x3 ^`beta3'

local y choice
local by1 id
sort `by1'

by `by1': egen double `den' = sum(exp(`xb'))
gen double `p' = exp(`xb')/`den'
mlsum `lnL' = `y' *log(`p') if `y'==1
if (`todo'==0 | `lnL' > =.) exit

end

try to avoid egen. It will slow down your code since it has to call ado-files which are unnecessary.

tempvar last
by `by1' : gen `last' = _n==_N
by `by1': gen double `den' = sum(exp(`xb'))
gen double `p' = sum(cond( `y'==1,exp(`xb'),0)
replace `p' = `p'/`den'
mlsum `lnL' = `y' *log(`p') if `last'
if (`todo'==0 | `lnL' > =.) exit

Thanks Christophe.

However, I test the code by running and it seems the 'xb' variable is missing.

Notice that in the original code I set 'xb' to be the sum of of a nonlinear transformation of the choice-specific attribute x1 to x3.

How should I handle such a case? My target is to specify 'xb' as above.\


Attached is the error:

myconditional_logit 0 __000009 __00000A __00000B __00000C
- `begin'
= capture noisily version 13: myconditional_logit 0 __000009 __00000A __00000B __00000C
------------------------------------------------------------------------------ begin myconditional_logit ---
- args todo b lnf
- version 11
- tempvar beta1 beta2 beta3 beta4 xb den p last xb
- mleval `beta1' = `b' , eq(1)
= mleval __00000D = __000009 , eq(1)
- mleval `beta2' = `b' , eq(2)
= mleval __00000E = __000009 , eq(2)
- mleval `beta3' = `b' , eq(3)
= mleval __00000F = __000009 , eq(3)
- mleval `beta4' = `b' , eq(4)
= mleval __00000G = __000009 , eq(4)
- local y choice
- local by1 id
- sort `by1'
= sort id
- tempvar last
- by `by1' : gen `last' = _n==_N
= by id : gen __00000M = _n==_N
- by `by1': gen double `den' = sum(exp(`xb'))
= by id: gen double __00000I = sum(exp(__00000L))
__00000L not found

[Christophe Kolodziejczyk]

try to uncomment the line with gen double `xb'

[Wong Wan Man]

try to uncomment the line with gen double `xb'

wont solve the problem and you know it as well.

anyway I code in python and it solved the problem.

forget about stata. just use python.

just in case someone who have difficulty dealing with nonlinear model specification come across this in the future, here's my message:

you can use package called "ad" in python to derive automatically the gradient and hessian of a likelihood function. you only need to write down the likelihood and it would compute the gradient and hessian automatically for you.

Then simply write a newton raphson to optimize the likelihood. simple and neat!
no need worry about stata syntax, d1 d0 d2 d3 d4 lf01234 etc etc.