Dear forum visitors!

I got a problem, working with panel data and max likelihood model. Specifically, I am looking for help in gradient coding. I am using d1 method and some of gradients are fine, but some are not. I would take those non-working one by one, now only one. So, d1debug mode shows a big difference between numerical and analytical gradient. I’ll show here a formula for gradient and then Stata 14 code, which is given only for related part, to keep it short. I expect, that the formula of the gradient is correct, I consulted about this. Also, d0 method works. The code that needs to be checked is a partial derivative of the likelihood w.r.t. the variance of the mean of tehnical inefficiency.

The following are the intermediate and help variables, also the parameters to be estimated. As I use panel data, two indexes are involved: i – group variable, t – time variable. Capital letter I means total number of groups (panels). Estimable parameters, included in gradient:

Array
Additionally, some common markings:
Array
Two intermediate variables, related to the likelihood and gradient:
Array
In addition, one partial derivative, used in corresponding gradient:
Array
And here is the gradient itself:
Array

Here is the code for that gradient, with related variables. Here `id' is a panel ID variable.

Related and intermediate variables

quie gen double `z' = `mu'/`sigmau'

quie gen double `mustar' = (`sigmav2'*`mu' - `sigmau2'*`Sres') /`sig'

quie gen double `sigstar' = sqrt(`sigmau2'*`sigmav2') / sqrt(`sig')

quie gen double `dstar' = `mustar' / `sigstar'

quie gen double `sig' = `T'*`sigmau2'+`sigmav2'

quie by `id': gen double `Bb' = (-`sigmav2'*(`mu'*`T' + `Sres' ) /`sigstar'*(`sig')^2)

by `id': replace `Bb' = `Bb'[_N]

here `Sres' means the sum of residual and a code for this

quie gen double `res' = ($ML_y1 -`bx') /*residual*/

quie by `id': gen double `Sres' = sum(`res')


quie by `id': gen double `B' = `Bb'-(`mustar'*`sigmav'^3) / (2*`sigstar'^2*`sigmau'*`sig'^1.5)

quie gen double `Sdmu' = normalden(-`z') / (1-normal(-`z'))

quie by `id': gen double `dB' = (`dstar'*`B')

quie by `id': gen double `musig' = ((`mu'^2) / 2*`sigmau2'^2)

Here is the code for the gradient

mlvecsum `lnf' `d3' = -0.5*`T' / `sig' + (`Sdmu'*`mu') / (2*`sigmau2'^1.5) + `Sdb'*`B' / 2 + `musig'*`dB' if `last'==1, eq(3)

I suspect something is wrong somewhere in the summation. Main suspects are the variable B and the last member of the gradient. Pls, do not look at the order of the commands. If anyone could see that bug, I would be very grateful.