*
* Example 22.10, page 800
* Duration models
*
open data hazard.txt
data(format=prn,org=columns) 1 62 dur prod
*
* Exponential
*
nonlin lambda
frml expduration = log(lambda)-lambda*dur
sstats(mean) / dur>>mean
compute lambda=1.0/mean
maximize(method=bfgs) expduration
*
* Weibull
*
nonlin lambda p
frml weibduration = log(p)+(p-1)*log(dur)+p*log(lambda)-(dur*lambda)**p
compute p=1.0
maximize(method=bfgs) weibduration
*
* Log logistic
*
nonlin lambda p
frml loglogduration = log(lambda)+log(p)+(p-1)*log(lambda*dur)-2*log(1+(lambda*dur)**p)
maximize(method=bfgs) loglogduration
*
* Log normal
*
compute lambda=.045,p=1
frml lognormduration = log(p/dur)+log(%density(p*log(lambda*dur)))+log(%cdf(-p*log(lambda*dur)))
maximize(method=bfgs) lognormduration
*
* Weibull with covariate
*
nonlin b0 b1 p
frml zfrml = b0+b1*prod
frml weibduration = (z=zfrml),(w=-exp(z)*dur**p),z+w+log(p)+(p-1)*log(dur)
compute p=1.0
maximize(method=bfgs) weibduration
*
* Proportional hazard
*
* RISKSET is a vector of vectors. Each individual has a vector and each vector has a slot
* for each individual. For individual i, slot j will be 1 if j’s exit time is greater than
* or equal to i’s.
*
compute nobs=62
dec vect[vect] riskset(nobs)
dec vect hazards(nobs)
do i=1,nobs
   dim riskset(i)(nobs)
   ewise riskset(i)(j)=dur(j)>=dur(i)
end do i
*
* The function HazardCalc gets computed at the start of each function evaluation to
* recalculate the relative hazard rates. These are put into the vector HAZARDS. The
* probability that an individual is the one to exit at her exit time is the ratio of her
* relative hazard to the sum of those hazards across the risk set for her exit time.
*
nonlin b1
compute b1=0.0
frml hazard = exp(b1*prod)
*
function HazardCalc
do i=1,nobs
   compute hazards(i)=hazard(i)
end do i
end
*
frml logl   = log(hazards(t))-log(%dot(riskset(t),hazards))
maximize(method=bfgs,start=HazardCalc()) logl