Gradient code bug

[Mark Liik]

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:

Additionally, some common markings:

Two intermediate variables, related to the likelihood and gradient:

In addition, one partial derivative, used in corresponding gradient:

And here is the gradient itself:

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.

[Phil Bromiley]

You'll increase your chances of a useful answer by following the FAQ on asking questions - provide Stata code in code delimiters, readable Stata output, and sample data using dataex. It also helps if you post the simplest example that represents your problem. Your post assumes we'll take the time to puzzle through all the equations, etc.

[Maarten Buis]

I don't know if this is the problem you run into, but

(1-normal(-`z'))

Is a classic example where "computer math" does not always agree with "human math". There is a nice discussion of why that is in this Stata Journal article: https://www.stata-journal.com/articl...article=pr0025 You can avoid the problem by coding that as

normal(`z')

In "human math" the two are equivalent; the standard normal distribution is symmetric around 0, so the area under the curve left of a number x is equal to the area under the curve right of -x. But computers find this way of computing much easier.

[Mark Liik]

This is a good tip, thank you! I checked the model and it helped a bit, but unfortunately not with that problem.