Using matrices in user-written maximum-likelihood estimation programs

[Jonathan Old]

I am trying to estimate a (structural) model with maximum-likelihood. The function I try to estimate has a matrix inverse, which I calculate as intermediate step in the ML program. However, I get an error when I use the parameter I try to estimate inside the matrix inverse function. The error is a r(111) (not found) error.

In the MWE below, I first generate a (random) symmetric matrix axx. I then generate an example dataset. I estimate the parameter beta via maximum likelihood with the mlestimation-function written next. One step within that function is to generate the inverse of a matrix that is a function of beta, and the not-found error happens here. If I replace beta with any scalar, there is no error.

/* [> Generate symmetric matrix of relations, axx <] */
clear
set obs 100
gen x = _n
expand 100
gsort x
gen z=_n-(x-1)*100
gen rel = rbinomial(1,0.2)

reshape wide rel  , i(x) j(z)

mkmat rel*, mat(aplus)
mata aplus = st_matrix("aplus")
mata a = makesymmetric(aplus)
mata st_matrix("axx",a)


/* [> Generate dataset <] */
drop _all
set obs 100
gen id = _n
gen rel = rbinomial(10,0.3)
gen rand1 = rnormal(0,1)
gen y = rand1 + rel + rel*rand1

/* [> Define ML-function <] */
cap program drop mlestimation
program mlestimation
   args lnf beta sigma
   tempvar GAMMA LAMBDA CVEC

   /* [> Generate temporary values/variables <] */
   qui gen `GAMMA' = 1/(1+`beta'*rel)
   mkmat `GAMMA', mat(G)
   qui egen `LAMBDA' = total(`GAMMA')
   /* [> Error happens here, as beta is "not found": r(111) <] */
   matrix C = invsym(I(100)+`beta'*axx)*G
   svmat C, name(`CVEC')

   /* [> Likelihood function <] */
   qui replace `lnf' = ln(normalden($ML_y1 - `LAMBDA'*(1-`LAMBDA')*`CVEC')/`sigma')-ln(`sigma')
end


/* [> Estimate ML <] */
ml model lf mlestimation ( y = ) /beta /sigma
ml maximize, difficult

[Joseph Coveney]

For one thing, your model has three equations, y:_cons, /beta and /sigma, which, combined with lnf, means that ml model will feed your evaluator program a total of four arguments. In its args line, you've set your evaluator program mlestimation to see only three of of the four that it's being given.

[Joseph Coveney]

For another thing, `beta' is a temporary variable name that ml model places in the dataset. invsym() is expecting only matrix names, and not variable names, for its arguments. Try something like

tempname Beta
mkmat `beta', matrix(`Beta')
matrix define C = invsym(I(100) + `Beta' * axx) * G

so that the matrix function gets a matrix as an argument instead of a variable. (After taking care of that, you'll still need to assure that all of the matrixes and vectors are conformable for the operations that you're doing.)

Check to make sure that you don't have the same problem in reverse a couple of lines later when you to refer to `CVEC' inside of the density function.

In general, because you're doing a fair amount of variable-to-matrix-to-variable-to-matrix operations (which are inherently clunky), I suspect that your code would look cleaner and be less liable to problems if you wrote the bulk of your likelihood evaluator (if not the entire thing) in Mata.