Problem of singular covariance matrix in a maximum likelihood estimation

[Khalid MAMAN WAZIRI]

Dear Statalist members,

Thanks to the helpful advice of some members, I was able to write a MATA program to maximize my own likelihood function which contains a tri-variate Normal CDF (using GHK2 mata-built-in function). However, I keep having error messages that the variance used by the GHK2 function is not definite-positive. Please find below component that the function GHK2 is calculating.


In fact, the variance matrix is V and it contains 4 different parameters to be estimated. The program seems to work well, except this issue with the variance V.
Can someone please advise me on what I need to do to overcome such an issue or get around it ?

Thank you very much,

Khalid WAZIRI.

[Sergio Correia]

Hi Khalid,

I can't see the picture you are referring to. Also, it might be better to the command & result as a [CODE] snippet (see the buttons at the top of the text box

[Khalid MAMAN WAZIRI]

Hi Khalid,

I can't see the picture you are referring to. Also, it might be better to the command & result as a [CODE] snippet (see the buttons at the top of the text box

Thank you Sergio, I will correct that.
Please see below the code, and the error message.

mata
void myds(transmorphic scalar ML, real rowvector b, real colvector lnf)

 
 {
 real colvector y, xb, wr1, wr2
 real colvector sig_v2, sig_u2, rho_v1, rho_v2

y  = moptimize_util_depvar(ML, 1)
wr1 = moptimize_util_depvar(ML, 2)
wr2 = moptimize_util_depvar(ML, 3)

xb = moptimize_util_xb(ML, b, 1)
sig_u2 = moptimize_util_xb(ML, b, 2)
sig_v2 = moptimize_util_xb(ML, b, 3)
rho_v1 = moptimize_util_xb(ML, b, 4)
rho_v2 = moptimize_util_xb(ML, b, 5)


real colvector rhov1, rhov2, sigv2, sigv, sigu2, sigu, sigma2, sigma, d1, d2, d3, d4, d5, d6
real matrix xe1, xe2, xe3, xg, v1, v2, v3, v1t, v2t, v3t, va
real scalar i

rhov1 = (exp(2 :* rho_v1) :- 1) :/ (exp(2 :* rho_v1) :+ 1)
rhov2 = (exp(2 :* rho_v2) :- 1) :/ (exp(2 :* rho_v2) :+ 1)


sigv2 = exp(sig_v2)
sigv = exp(0.5 :* sig_v2)

sigu2 = exp(sig_u2)
sigu = exp(0.5 :* sig_u2)

sigma2 = sigv2 :+ sigu2
sigma  = sqrt(sigma2)
epsilon = y :- xb

tvar1 = (1 :/ sigma) :* normalden((epsilon :+ 0):/ sigma)
tvar2 = normal(0 :/ sigu)


tmp1 = epsilon :+ 0


gam1 = rhov1 :* sigv :* tmp1 :/ sigma2
gam2 = rhov2 :* sigv :* tmp1 :/ sigma2
gam3 = -1 :* sigu2 :* tmp1 :/ sigma2


kap1 = -1 :* wr1
kap2 = -1 :* wr2
kap3 =  0

xe1 = gam1 :- kap1
xe2 = gam2 :- kap2
xe3 = gam3 :- kap3

xg = (xe1, xe2, xe3)


d1 = ((1 :- rhov1 :* rhov1) :* sigv2 :+ sigu2) :/ sigma2
d2 = ((0 :- rhov1 :* rhov2) :* sigv2 :+ 0 :* sigu2) :/ sigma2
d3 = rhov1 :* sigv :* sigu2 :/ sigma2
d4 = ((1 :- rhov2 :* rhov2) :* sigv2 :+ sigu2) :/ sigma2
d5 = rhov2 :* sigv :* sigu2 :/ sigma2
d6 = sigv2 :* sigu2 :/ sigma2

v1=(d1,d2,d3)
v1t=v1'
v2=(d2,d4,d5)
v2t=v2'
v3=(d3,d5,d6)
v3t=v3'

nobs=rows(y)

va=J(3,3,.)
tvar3=J(nobs,1,.)


p2 = ghk2setup(10000, 5, 3, "halton")
anti = 0
start = .

for (i=1; i<=nobs; i++) {

va=(v1t[.,i],v2t[.,i],v3t[.,i])

tvar3[i,1]=ghk2(p2, xg[i,.], va, anti, start)
 
va=J(3,3,.)

}


lnf = ln(tvar1) :- ln(tvar2) :+ ln(tvar3)
  }

 end
 
 
  ml model lf myds() (frontier: $yvar $WW1 $WW2= $xvar) (sigu2:) (sigv2: ) (trhov1: ) (trhov2: ) if $Indi ==1

  mat coef_lai = .15044475  , .02353871 , -.00039912 ,  .62243346  , -1.631312  ,-1.5822784   ,0, 0
  
  ml init coef_lai, copy
  ml check

Therefore, there are four parameters to be estimated in the matrix va : rhov1 , rhov2 , sigv2 , sigu2 (sigma2 is just the sum of sigv2 and sigu2).
and the error message after ml check is :

.   ml check

Test 1:  Calling myds() to check if it computes log likelihood and
         does not alter coefficient vector...
         Passed.

Test 2:  Calling myds() again to check if the same log likelihood value is
         returned...
         Passed.

Test 3:  Calling myds() to check if 1st derivatives are computed...
         test not relevant for type lf evaluators.

Test 4:  Calling myds() again to check if the same 1st derivatives are
         returned...
         test not relevant for type lf evaluators.

Test 5:  Calling myds() to check if 2nd derivatives are computed...
         test not relevant for type lf evaluators.

Test 6:  Calling myds() again to check if the same 2nd derivatives are
         returned...
         test not relevant for type lf evaluators.

------------------------------------------------------------------------------
Searching for alternate values for the coefficient vector to verify that
myds() returns different results when fed a different coefficient vector:

Searching...
initial:       log likelihood =     -<inf>  (could not be evaluated)
searching for feasible values ghk2: covariance matrix is not positive-definite.
r(3352);

I re-attach the picture where you can see the variance matrix composition (in the code it is "va" ).

Thanks a lot !
Khalid.

[Sergio Correia]

You might want to print the matrix to screen every time the fn is called, to see what is exactly the problem. Maybe the matrix is not symmetric, maybe one of the rhos have a typo so it's not positive def, etc. But to debug that out you need to carefully check the input values and the resulting V.

[Khalid MAMAN WAZIRI]

You might want to print the matrix to screen every time the fn is called, to see what is exactly the problem. Maybe the matrix is not symmetric, maybe one of the rhos have a typo so it's not positive def, etc. But to debug that out you need to carefully check the input values and the resulting V.

Thank you Sergio, actually I did that and the matrix is symmetric but it could become singular at some iterations. It depends on the parameter values passed to the likelihood function.
Thank you for your answer!