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)