disableprob=function(theta,psi,zeta,l,K,M)
{
# Transition probabailities are computed for a disability scheme
# with rehabilitation. Dying, disability and rehabilitation
# are based on Gomperz-Makeham intensities with parameter vectors
# "theta", "psi" and "zeta" respectively. The inital age is "l". "M" is
# the technical parameter for computing introduced in Section 12.5
#
# Output are the transition probabilities given age "l" up to time "K"
# stored as the three-way array "kp" where the first index is current state,
# the second the state left and the third one is time.
#Initializing
K1=K*M
P=array(0,dim=c(3,3,K1+1))
ind=1+0:K*M
P[,,1]=diag(c(1,1,1),3)
#Loop over k
for(k in 1:K1){
lk=l+(k-1)/M
p=matrix(0,3,3)
p[2,1]=psi[1]+psi[2]*exp(psi[3]*lk)
p[3,1]=theta[1]+theta[2]*exp(theta[3]*lk)
p[1,2]= zeta[1]+zeta[2]*exp(zeta[3]*lk)
p[3,2]=p[3,1]
p=p/M
diag(p)=1-apply(p,2,sum)
P[,,k+1]=P[,,k]%*%p}
#Loop finished
list(kp=P[,,ind])
}