st: Non-linear random effect model: write my own ml program

["Yuyan Shi" <shiyuyan@gmail.com>]

Dear All,

I want to estimate a user-specified non-linear random effect model, which
does not have routine in Stata.

The function is:
*************************************
Pij = 1-normal[(Tij-mu)/sigma]+ui+eij
*************************************
Where i denotes individual, ij denotes "j"th observation for individual i.
This is a panel data.
Pij and Tij are dependent variables, mu and sigma are functions of x that I
want to estimate, ui is individual random effect, and eij is random error
term. Assume ui and eij is normally distributed with mean 0 and unknown
standard error.

I know I should integrate j observation for individual i, then integrate all
individuals to get summation of log likelihood. However, my experience in
non-linear panel data programming is limited. Can someone outline the steps
I should do for a ml estimation?

This is what I have so far based on linear random-effect model(not
converging). Not sure I am in the right direction. Really appreciate your
help.

----------------------------------------------------------------------------
---
program define myre
  args todo b lnf
  tempvar xb mu lgmu sigma lgsigma sigmau lgsigmau z T S_z2 Sz_2 S_temp a
first
  tempname sigmae lgsigmae
  mleval `lgmu' = `b', eq(1)
  qui gen double `mu' = exp(`lgmu')
  mleval `lgsigma' = `b', eq(2)
  qui gen double `sigma' = exp(`lgsigma')
  qui gen double `xb' = normal(($ML_y1-exp(`lgmu'))/exp(`lgsigma'))
  mleval `lgsigmau' = `b', eq(3)
  mleval `lgsigmae' = `b', eq(4) scalar
  qui gen double `sigmau' = exp(`lgsigmau')
  scalar `sigmae' = exp(`lgsigmae')

  sort id
  qui {
    gen double `z' = $ML_y2 - 1 + `xb'
    *number of observations for each individual i
    by id: gen `T' = _N
    gen double `a' = (`sigmau'^2) / (`T'*(`sigmau'^2) + `sigmae'^2)
    by id: egen double `S_z2' = sum(`z'^2)
    by id: egen double `S_temp' = sum(`z')
    by id: gen double `Sz_2'    = `S_temp'^2
    by id: gen `first' = (_n == 1)
    mlsum `lnf' = -.5 *                                ///
          (   (`S_z2' - `a'*`Sz_2')/(`sigmae'^2)  +   ///
              log(`T'*`sigmau'^2/`sigmae'^2 + 1) +   ///
              `T'*log(2*_pi*`sigmae'^2)            ///
          )  if `first' == 1
  }
end
ml model d0 myre (lgmu:y1 y2=$xvar) (lgsigma:$xvar) (lgsigmau:$zvar)
(lgsigmae:)

- Yuyan

*
*   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/