Multiple State Transition Probability calculated using R

80 views
edited Sep 12, 2017

This question refers to the Contingencies text book: Actuarial Mathematics Life Contingent Risks, Second Edition, David C.M. Dickson, et al. Section 8.6 - Question 8.6.

Question 8.6 uses the disability income multiple state model with specific transition rates and "numerical integration" in order to solve for the premiums.

I however was more interested in trying to make a function in R which calculates the transition probabilities: $$_tp_x^{00}$$ and $$_tp_x^{01}$$. This can also be used to solve the question but that is more of a side note.

The transition probabilities were given and using the in build integration function in R, I was able to make a function which calculates $$_tp_x^{\overline{00}}$$.

The main problem is that in the disability income model, transitions are possible between both state 0 and state 1 so $$_tp_x^{00}$$ and $$_tp_x^{01}$$ are functions of each other so I did not see how I would be able to solve for them using this simple 1 parameter integration function in R. Is there a way to solve $$_tp_x^{00}$$ and $$_tp_x^{01}$$ using integration in R or simplifying it by hand without the use of approximations? See below my R code which has the transition rate functions and the $$_tp_x^{\overline{00}}$$ function.

a1=4*10^-4

a2=5*10^-4

b1=3.4674*10^-6

b2=7.5858*10^-5

c1=0.138155

c2=0.087498

x=60

mu01 = function(x)

a1+b1*exp(c1*x)

mu10 = function(x)

0.1*mu01(x)

mu02 = function(x)

a2+b2*exp(c2*x)

mu12 = function(x)

mu02(x)

mu0102 = function(x)

mu01(x)+mu02(x)

tpx00_bar = function(t)

exp(-integrate(mu0102,x,x+t)\$value)

tpx00_bar(60)

The authors have given there the Kolmogorov forward equations which $$_tp_{60}^{00}$$ and $$_tp_{60}^{01}$$ satisfy. The authors also discuss that one would need an Euler discretization method in order to obtain the required probabilities.
Ofcourse, it would have been ideal if we could write down an expression for $$_tp_{60}^{00}$$ as an integral, so that we could try solve this either analytically or at the very least, numerically. However, this is not possible here.