st: timing of events model

["marcel spijkerman" <marcel_spijkerman@hotmail.com>]

Dear all,

I want to estimate treatment effects with duration data. An often used model 
for this is the so called timing-of- events model that allows for 
selectivity through the correlation of unobserved variabels between the 
proces of entering treatment and transition from one state to the other.
To program this I tried to made use of HSHAZ (discrete duration with 
heterogeneity) written for STATA by Stephen Jenkins. My program is the one 
below.

The model simultaneously estimates two duration processes (treatment and 
unemployment), indicated by $ML_y1 and $ML_y2. Because it is possible  
entering treatment before the transition to employment I added a third 
variable $ML_y3 that ensures sums are taken over the correct period.

Although I think this should work, several experiments to estimate the model 
did not produce any results. Usually maximisation breaks down.

My questions are these:

to those who are familiar with this model (timing of events with 
Heckman-Singer unobserved heterogeneity, is this a correct way to model it;

are there perhaps other ways to program this model (i am not hooked up to a 
discrete version)

Help is kindly appreciated,

Marcel

ps: I specify all four mass points because I do not estimate constants.


program define hshazllb4

	version 8.2

	args todo b lnf


		tempvar I1 I2 lp2 lp3 lp4 m1 m2 m3 m4 h1a h1b h2a h2b sum1a sum1b sum2a 
sum2b ST1a ST1b ST2a ST2b last lnfi
		tempname p1 p2 p3 p4


		mleval `I1' = `b' , eq(1)
		mleval `I2' = `b' , eq(2)
               mleval `m1' = `b', scalar eq(3)
		mleval `m2' = `b', scalar eq(4)
		mleval `m3' = `b', scalar eq(5)
		mleval `m4' = `b', scalar eq(6)
		mleval `lp2' = `b', scalar eq(7)
		mleval `lp3' = `b', scalar eq(8)
		mleval `lp4' = `b', scalar eq(9)


		quietly {

			scalar `p1' = 1 / (1 + exp(`lp2') + exp(`lp3') + exp(`lp4'))
			scalar `p2' = exp(`lp2') / (1 + exp(`lp2') + exp(`lp3') + exp(`lp4'))
			scalar `p3' = exp(`lp3') / (1 + exp(`lp2') + exp(`lp3') + exp(`lp4'))
			scalar `p4' = exp(`lp4') / (1 + exp(`lp2') + exp(`lp3') + exp(`lp4'))



			by id: gen double `h1a' = 1-exp(-exp(`I1' + `m1'))
			by id: gen double `h1b' = 1-exp(-exp(`I1' + `m2'))
			by id: gen double `h2a' = 1-exp(-exp(`I2' + `m3'))
			by id: gen double `h2b' = 1-exp(-exp(`I2' + `m4'))

			by id: gen double `sum1a' = sum(exp(`I1' + `m1'))
			by id: gen double `sum1b' = sum(exp(`I1' + `m2'))
			by id: gen double `sum2a' = sum((exp(`I2' + `m3'))*$ML_y3)
			by id: gen double `sum2b' = sum(exp((`I2' + `m4'))*$ML_y3)

			by id: gen double `ST1a' = exp(-`sum1a'[_N]) if _n==_N
			by id: gen double `ST1b' = exp(-`sum1b'[_N]) if _n==_N
			by id: gen double `ST2a' = exp(-`sum2a'[_N]) if _n==_N
			by id: gen double `ST2b' = exp(-`sum2b'[_N]) if _n==_N

			by id: gen byte `last' = (_n==_N)

			gen double `lnfi' = cond(!`last',0,	///
				ln( (`p1'*(((`ST1a' * (`h1a'/(1-`h1a'))^$ML_y1) * 
(`ST2a'*(`h2a'/(1-`h2a'))^$ML_y2)))) + 	///
                           (`p2'*(((`ST1a' * (`h1a'/(1-`h1a'))^$ML_y1) * 
(`ST2b'*(`h2b'/(1-`h2b'))^$ML_y2))))             +  ///
                           (`p3'*(((`ST1b' * (`h1b'/(1-`h1b'))^$ML_y1) * 
(`ST2a'*(`h2a'/(1-`h2a'))^$ML_y2)))) 		   +	///
                           (`p4'*(((`ST1b' * (`h1b'/(1-`h1b'))^$ML_y1) * 
(`ST2b'*(`h2b'/(1-`h2b'))^$ML_y2)))) 		))

			mlsum `lnf' = `lnfi'

                        }

end


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/