MaxLik QML Covariance Matrix

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.