Hi,

I am running into an issue with using _max_CovPar = 3. It is supposed to return the QML covariance matrix but I am almost certain it returns (QML covariance matrix) * (number of observations) instead.

The proper calculation is inv(hessian)*(gradient*gradient'/N)*inv(hessian) since the gradient is K by N. However, I believe _max_CovPar = 3 returns inv(hessian)*(gradient*gradient')*inv(hessian).

Can someone please clarify for me?

Thanks!

## 8 Answers

```
h0 = hessp(&lpr,b,z);
g0 = gradp(&lpr,b,z);
A = -h0;
B = g0'g0;
qml1 = invpd(A)*(B)*invpd(A);
qml2 = invpd(A)*(B/N)*invpd(A);
print "ols " stderr';
print "Maxlik qml " sqrt(diag(cov))';
print "qml1 " sqrt(diag(qml1))';
print "qml2 " sqrt(diag(qml2))';
```

ols 0.0099537875 0.010130433 0.010065362 Maxlik qml 0.0099536075 0.010143001 0.010164126 qml1 0.0099547791 0.010146236 0.010170286 qml2 9.9547791e-05 0.00010146236 0.00010170286

Here is the full command file I used:

```
library maxlik;
#include maxlik.ext
maxset;
proc lpr(b,z);
local s,m;
m = z[.,1] - b[1] - z[.,2:3]*b[2:3];
s = m'm/rows(m);
retp(lnpdfmvn(m,s));
endp;
N = 10000;
rndseed 4562456;
b0 = { 1, 1, 1 };
x = rndn(N,2);
y = b0[1] + x*b0[2:3] + rndn(N,1);
z = y~x;
_max_covpar = 3;
{b,f,g,cov,ret} = maxlik(z,0,&lpr,b0);
{ vnam, m, b, stb, vc, stderr, sigma, cx, rsq, resid, dwstat } = ols("",y,x);
h0 = hessp(&lpr,b,z);
g0 = gradp(&lpr,b,z);
A = -h0;
B = g0'g0;
qml1 = invpd(A)*(B)*invpd(A);
qml2 = invpd(A)*(B/N)*invpd(A);
print "ols " stderr';
print "Maxlik qml " sqrt(diag(cov))';
print "qml1 " sqrt(diag(qml1))';
print "qml2 " sqrt(diag(qml2))';
;
```

It might be worth noting that Maxlik minimizes - logL / N.

Thanks,

I agree that this works when you use the built in gradient function. That function returns a Kx1 vector, correct?

However, when you use a user defined analytic gradient that returns a KxN matrix rather than a Kx1 vector the correct formula is your qml2 formula. If _max_CovPar = 3 is returning qml1 when using an analytic gradient then there is a problem. Can you verify this for me?

Thanks,

Casey

When the function returns a vector of log-probabilities, Maxlik computes an NxK matrix of derivatives, which is the case in my example. You can see the results. Your suggestion dividing B by N produces a clearly wrong answer, and the version where B is not divided by N is clearly a correct answer. So if there's a problem, please describe it in more detail.

I'm sorry. You're saying that g0 is NxK?

I hate to bug you but what are the dimensions of your g0?

Hey Crozowski,

You can copy and paste the code that Ron posted into your GAUSS editor and run it. If you do that you can look at all the variables and will see that 'g0' is 10,000x3, which is NxK for this example.

Thanks! I will when I get the chance.