Statalist The Stata Listserver

[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: timing of events model

From   "marcel spijkerman" <>
Subject   st: timing of events model
Date   Thu, 20 Jul 2006 13:22:21 +0000

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,


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'



* For searches and help try:

© Copyright 1996–2015 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index