st: Problem with ml maximize: numerical derivatives are approximate

[Amir Fekrazad <amir.fekrazad@gmail.com>]

Dear statalisters,

I wrote a program to estimate the elements of the inverse
variance-covariance matrix (Q=inverse(V)) of a joint student t
distribution.  I assumed that Q is symmetric, toeplitz and 3-banded,
ie. q_ij=f( |i-j| ) and q_ij=0 if |i-j|>3
(I know that my assumptions are not true for non-independent data but
I tested the program only with independent data)

With these assumptions, I only have to estimate parameters q0
(diagonal element), and q1, q2 and q3 (qn = elements with distant n
from the diagonal)

I checked and traced my program and I am quite sure it has no problem.
I was careful to use double variables wherever needed and I gave it
reasonable initial values (or let "ml search" search for them itself).
however, it never converges and I get the following message with
different datasets and data sizes:

numerical derivatives are approximate
flat or discontinuous region encountered

I'm stuck and don't know what to do... I searched statalist archive
but couldn't find any helpful thread. I attach my program below. It is
somehow unconventional and different from other ml programs. Any
advice is really appreciated.

======================================================================

program myprog0
version 9.1
args todo b lnf

tempname q0 q1 q2 q3 q

local nu = 30
local T = _N

mleval `q0' = `b' , eq(1) scalar
mleval `q1' = `b' , eq(2) scalar
mleval `q2' = `b' , eq(3) scalar
mleval `q3' = `b' , eq(4) scalar

mat define `q' = `q0' * I(`T')

// start: here we should specify elements of q matrix
forval i=1/`T' {

if(`i'-1 > 0){
mat `q'[`i',`i'-1]=`q1'
}
if(`i'+1<=`T'){
mat `q'[`i',`i'+1]=`q1'
}

if(`i'-2 > 0){
mat `q'[`i',`i'-2]=`q2'
}
if(`i'+2<=`T'){
mat `q'[`i',`i'+2]=`q2'
}


if(`i'-3 > 0){
mat `q'[`i',`i'-3]=`q3'
}
if(`i'+3<=`T'){
mat `q'[`i',`i'+3]=`q3'
}

}
// end: here we should specify elements of q

tempname p1 p2 p3 p4 p5 mm nn cc
tempvar mu ymu ymumat

egen double `mu' = mean($ML_y1)
gen double  `ymu' = $ML_y1 - `mu'
mkmat `ymu', mat(`ymumat')
scalar `p1' = -0.5*`T'*ln(c(pi)*`nu')
scalar `p2' = lngamma( 0.5*(`nu'+`T') )
scalar `p3' = -lngamma(0.5*`nu')
scalar `mm' = det(`q')
if(`mm'<=0){
scalar `p4' = (-1+`mm')* (10^8)
}
else{
scalar `p4' = 0.5*ln(`mm')
}
matrix define `cc' = (`ymumat')'*`q'*`ymumat'
scalar `nn' = 1+ (  (1/`nu')*det(`cc')  )
if(`nn'<=0){
scalar `p5' = (-1+`nn')*(10^8)
}
else{
scalar `p5' = -0.5*(`nu'+`T')*ln(`nn')
}

scalar `lnf' = `p1' + `p2'+`p3'+`p4' + `p5'

end

////////////


ml model d0 myprog0 (q0: y=) /q1 /q2 /q3
ml search
ml check
ml max

======================================================================
output:
initial:       log likelihood = -112.82474
rescale:       log likelihood = -102.29964
rescale eq:    log likelihood = -75.125482
Iteration 0:   log likelihood = -75.125482  (not concave)
Iteration 1:   log likelihood =  93.746918  (not concave)
Iteration 2:   log likelihood =  151.28476  (not concave)
Iteration 3:   log likelihood =  188.49715  (not concave)
Iteration 4:   log likelihood =   481.5772  (not concave)
Iteration 5:   log likelihood =  634.17554  (not concave)
Iteration 6:   log likelihood =  820.38195  (not concave)
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 7:   log likelihood =  901.99555  (not concave)
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
============================================================


Best regards,
Amir

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/