Maximum likelihood estimator for multivariate normal distribution

[Alex Stead]

Hi all,

What I'm trying to do is to create a maximum likelihood estimator for a multivariate normal model, which can estimate a system of equations. After having tried and failed to do this in the usual way, I've come to the conclusion that this is best done in Mata, since the lnmvnormal function I am using takes matrices as arguments.

Unfortunately Mata doesn't seem to have as much documentation as MLE in Stata, and my experience with it is extremely limited. Based on what I could find, I've made an attempt. However it seems I've gone wrong somewhere. I've tried getting this to work for a very simple example using simulated data to start (no cross-equation restrictions, no correlation between error terms).

capture drop e1 e2 e3 y1 y2 y3 x1 x2 x3
capture gen double e1= invnorm(uniform())
capture gen double e2= invnorm(uniform())
capture gen double e3= invnorm(uniform())
capture gen double x1= invnorm(uniform())
capture gen double x2= invnorm(uniform())
capture gen double x3= invnorm(uniform())
capture gen double y1= x1 + e1
capture gen double y2= x2 + e2
capture gen double y3= x3 + e3

mata:

function mynormal_d0(transmorphic M, real scalar todo, real rowvector b, fv, g, H){
    y1 = moptimize_calc_depvar(M, 1)
    y2 = moptimize_calc_depvar(M, 2)
    y3 = moptimize_calc_depvar(M, 3)
    xb1 = moptimize_calc_xb(M, b, 1)
    xb2 = moptimize_calc_xb(M, b, 2)
    xb3 = moptimize_calc_xb(M, b, 3)
    lnsigma1 = moptimize_calc_xb(M, b, 4)
    lnsigma2 = moptimize_calc_xb(M, b, 5)
    lnsigma3 = moptimize_calc_xb(M, b, 6)
    y = (y1, y2, y3):
    mu = (xb1, xb2, xb3):
    sigmasq1 = exp(lnsigma1)^2:
    sigmasq2 = exp(lnsigma2)^2:
    sigmasq3 = exp(lnsigma3)^2:
    sigma = (sigmasq1, 0, 0 \ 0, sigmasq2, 0 \ 0, 0, sigmasq3):
    fv = moptimize_util_sum(M, lnmvnormalden(mu, sigma, y):)
}

end

: M = moptimize_init()
: moptimize_init_evaluator(M, &mynormal_d0())
: moptimize_init_evaluatortype(M, "d0")
: moptimize_init_depvar(M, 1, "y1")
: moptimize_init_depvar(M, 2, "y1")
: moptimize_init_depvar(M, 3, "y1")
: moptimize_init_eq_indepvars(M, 1, "x1")
: moptimize_init_eq_indepvars(M, 2, "x2")
: moptimize_init_eq_indepvars(M, 3, "x3")
: moptimize(M)
: moptimize_result_display(M)

What I get when I try to run this is:

mata:
------------------------------------------------- mata (type end to exit) -----------------------------
: 
: function mynormal_d0(transmorphic M, real scalar todo, real rowvector b, fv, g, H){
>         y1 = moptimize_calc_depvar(M, 1)
>         y2 = moptimize_calc_depvar(M, 2)
>         y3 = moptimize_calc_depvar(M, 3)
>         xb1 = moptimize_calc_xb(M, b, 1)
>         xb2 = moptimize_calc_xb(M, b, 2)
>         xb3 = moptimize_calc_xb(M, b, 3)
>         lnsigma1 = moptimize_calc_xb(M, b, 4)
>         lnsigma2 = moptimize_calc_xb(M, b, 5)
>         lnsigma3 = moptimize_calc_xb(M, b, 6)
>         y = (y1, y2, y3):
>         mu = (xb1, xb2, xb3):
>         sigmasq1 = exp(lnsigma1)^2:
>         sigmasq2 = exp(lnsigma2)^2:
>         sigmasq3 = exp(lnsigma3)^2:
>         sigma = (sigmasq1, 0, 0 \ 0, sigmasq2, 0 \ 0, 0, sigmasq3):
>         fv = moptimize_util_sum(M, lnmvnormalden(mu, sigma, y):)
illegal arglist
(2 lines skipped)
-------------------------------------------------------------------------------------------------------
r(3000);

I've done a lot of guesswork and interpretation here, and I'm not sure if I'm even along the right tracks. Any help appreciated.