st: What's wrong with my differentiation?

["Mak, Timothy" <timothy.mak07@imperial.ac.uk>]

Dear list, 

Sorry if this is not much of a Stata question, but it has been bugging me for hours. 

Let's say we input the following: 


input yesno freq
1 4
0 6
end

logit yesno [fw=freq]
matrix list e(V)

We get a variance of 0.416667, suggesting the hessian = -1/0.416667 = -2.4

However, for some reason, when I tried to get this hessian analytically, I fail. I know it must be to do with my calculations, but I haven't been able to spot the mistake for 5 hours, and I checked it with Wolframalpha, and it's still wrong. So here's my derivation: 

y=4
z=6

loglikelihood = lik = y*log(p) + z*log(1-p) + constant
p = 1/(1+exp(-x)) = exp(x)/(1+exp(x))
dlik/dx = dlik/dp * dp/dx (chain rule)

dlik/dp = y/p - z/(1-p)
dp/dx = exp(x)/(1+exp(x))^2 = exp(-x)/(1+exp(-x))^2

Therefore, dlik/dx = (y/p - z/(1-p)) * exp(x)/(1+exp(x))^2

Denote: 
A=exp(x)/(1+exp(x))^2 = dp/dx
B=(y/p - z/(1-p))

********************
dlik/dx = A*B
********************

d2lik/dx2 = d/dx (dlik/dx) = A dB/dx + B dA/dx (product rule)
dB/dx = dB/dp * dp/dx = dB/dp * A (chain rule)
dB/dp = -yp^(-2) + z(1-p)^(-2) = z/(1-p)^2 - y/p^2 

********************************
dB/dx = A* (z/(1-p)^2 - y/p^2 )
********************************

dA/dx = exp(x) d/dx (1+exp(x))^(-2) + (1+exp(x))^(-2) * d/dx (exp(x)) (product rule)
d/dx(1+exp(x))^(-2) = -2exp(x)(1+exp(x))^(-3) (chain rule)
d/dx(exp(x)) = exp(x)

Therefore, 
*****************************************************************************
dA/dx = -2 * exp(x) * exp(x) * (1+exp(x))^(-3)  + (1+exp(x))^(-2) * exp(x) 
	= A * (-2exp(x)/(1+exp(x)) + 1) = A * (1 - 2/(1+exp(-x)))
*****************************************************************************

Therefore, 

********************************************************************
d2lik/dx2 = A*A*( z/(1-p)^2 - y/p^2 ) + B * A * (1 - 2/(1+exp(-x)))
********************************************************************

Now, for the above toy example, the ML of p = 0.4, which means on the logit scale, x=logit(0.4)=-0.405

Plugging x=-0.405 into the above equations: 

A=0.24
B=0
d2lik/dx2 = 0.24^2 * (6/0.36 - 4/0.16) = -0.48

Hessian = -0.48 != -2.4, which is what the answer should be! 

Thank you very much for your help if you've read so far! As I mentioned before I tried putting the above into Wolframalpha, but it ends up with such a complicated formula that after 1 or 2 hours, I still haven't been able to work out what's gone wrong. Following Wolframalpha's formula, however, I did manage to get -2.4, so obviously my formula above contains a mistake. 

(The reason I need the analytic formula is so that I can do write an estimation procedure for something, for which the above represents a special case. What I'm writing is a small extension to -logit-, although I can't exactly use -logit- in my engine, hence my wish to write down the analytic formula.) 

Thanks again, 

Tim

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/