* DISCLAIMER: The following code is provided as is and should be used at your own risk. 
* No user support is available.

* Nigel Rice, November 2006


* This program relates to the material described in Chapter 7 of the book 
* Jones, A.M., Rice, N., Bago d'Uva, T. & Balia, S. (2007)
* "Applied Health Economics", London: Routledge {ISBN: 9780415397728}  


/***********************************/
/* DURATION MODELS: AN APPLICATION */
/* TO HEALTH AND RETIREMENT        */
/***********************************/

clear 
capture log close 
log using "H&R application_W.log", replace 

set trace off 
set more off 
set mem 700m 
program drop _all 

use H&R_application.dta, clear 

/* Analysis for females only */

drop if male == 1 

/* Summarise Labour Market Transitions */

forvalues j=1(1)12 {
 table jbstat if wavenum == `j' 
 }

/* summarise data by retirement status */

local summvars "retired hlltyes sah4ex sah4gvg sah4fair sah4vpp shlltyes ssah4ex ssah4gvg ssah4fair ssah4vpp m2lnhinc HseOwn HseMort HseRent HseAuthAss marcoup deghdeg hndalev ocse everppenr everemppr privcomp0 civlocgov jbsecto0 lspjb"

summ `summvars' 
summ `summvars' if retired == 0 
summ `summvars' if retired == 1

/* Estimate health stock */
/* Ordered Probit to generate latent self-assessed health variable */

oprobit hlstatc4 hlparms hlpsee hlphear hlpskin hlpchest hlpheart hlpstom hlpdiab hlpanx hlpalch hlpepil hlpmigr hlpothr 
predict sahlat if e(sample), xb

sort pid wavenum 
by pid: gen lsahlat = sahlat[_n-1] 
by pid: gen lhlltyes = hlltyes[_n-1]
by pid: gen sahlat0 = sahlat[1] 
by pid: gen hlltyes0 = hlltyes[1]

/* Empirical approach to duration modelling */
/* Descriptive statistics */

sort pid wavenum 
iis pid 
tis wavenum 
xtdes, patterns(20) 

/* Set failure time data */

sort pid wavenum 
stset wavenum , id(pid) failure(retired == 1) origin(wavenum == 1) 
tab _t _d 
tab _d lhlltyes 

/* Life tables for health limitations */

sort pid wavenum 
ltable _t (_d) , by(hlltyes0) test tvid(pid) failure 

sts graph, by(hlltyes0) l1title("Proportion not retired") clpattern("1" "_-") noshow 


/* JENKIN'S .ado file for discrete time proportional hazards models */
/* NB: pgmhaz command only works for values of deadvar = 0 or 1  */
/* Need to remove missing data */

sort pid wavenum 

drop if _d == . 
iis pid 
tis wavenum 
xtdes , patterns(20) 

/* Create time dummies for complementary log-log baseline hazard  */

forvalues j = 1/11 { 
     gen t`j' = 0 
     recode t`j' 0 = 1 if _t == `j' 
     } 

/* Use health limitations */
/* Set local list of variables */

local survars "lhlltyes hlltyes0 age5054 age5559 age6064 m2lnhinc lHseMort lHseAuthAss lHseRent marcoup deghdeg hndalev ocse everppenr everemppr privcomp0 civlocgov0 jbsecto0 lspjb shlltyes NorthW NorthE SouthW London Midland Scot Wales"

/* Analysis using complementary log-log function and glm */

cloglog _d `survars' t2-t11, nolog

glm _d `survars' t2-t11 , family(binomial) link(cloglog) 

pgmhaz8 `survars' t2-t11 , i(pid) d(_d) s(_t) nobeta0 nolog

/* Use latent SAH */
/* Set local list of variables */

local survars2 "lsahlat sahlat0 age5054 age5559 age6064 m2lnhinc lHseMort lHseAuthAss lHseRent marcoup deghdeg hndalev ocse everppenr everemppr privcomp0 civlocgov0 jbsecto0 lspjb slatsah NorthW NorthE SouthW London Midland Scot Wales"

pgmhaz8 `survars2' t2-t11 , i(pid) d(_d) s(_t) nobeta0 nolog


end