# Projecting a Thiele differential equation with Black Scholes returns

I am trying to solve the equation

$\frac{d}{dt}V(t)=r(t)V(t)+\pi-\mu(x+t)(b_d-V(t))$

numerically using the R function 'ode'. This is a Thiele differential equation for a life insurance reserve with premium rate $\pi$, mortality intensity $\mu$ (for an $x+t$ year old), death benefit $b_d$ and interest rate process $r$.

I want to do this in a so called unit-linked setting, where the returns on the policy are generated by investments in stocks, hence I have assumed a Black Scholes model for simplicity.

I generate a Geometric Brownian Motion. As seen the time interval is 40 years and the step size of the simulation is $40/100000$. For each simulated point I calculate the returns as

$r[i]=\frac{S[i+1]-S[i]}{S[i]}$

Such that I have $100000$ return values. Plotting these two for two very different scenarios yields

My problem is, at I would expect the reserve process to vary a lot more, much like the simulated Geometric Brownian Motion. In essence, it is too smooth, i think. Also if I generate GBM trajectories which end in, say, the value 500 and 10 respectively, the difference in the final value of the reserve varies very little.

Does anyone know why this is, am I doing something wrong? The R-code is attached.

maturity <- 40
simulation.length <- 100001
dt <-  maturity/(simulation.length-1)

timeline <- seq(0,maturity, dt)

BM <- GBM <- EV <- rep(0, times=simulation.length)
EV[1] <- GBM[1] <- S0

for(i in 2:simulation.length){
BM[i] <- BM[i-1]+sqrt(dt)*rnorm(1)
GBM[i] <- GBM[1]*exp((mu-(sigma^2)/2)*(i-1)*dt+sigma*BM[i])
EV[i] <- EV[1]*exp(mu*(i-1)*dt)
}

return <- rep(0,length(GBM))
returns[1] <- 0
for (i in 2:length(GBM)-1)
{
returns[i] <- (GBM[i+1] - GBM[i]) / GBM[i]
}

dV <- function(t,V, parms)
{
list(returns[t/0.0004+1] * V + premiumRate - mortalityIntensity(t+25) *
(deathBenefit - V))
}

out <- ode(y = 0, times = timeline, func = dV, parms = NULL)

interpolatedReserve <- approxfun(times,out[,2], method="linear")


The problem is not that the approxfun function is not good enough, because I used it on the GBM to plot the green trajectory.

I use 60000 as the premium rate (around 1000 dollars per month paid to pension) and a death benefit of 1.000.000 (realistic numbers in Danish Kroner). Is it because the returns $r\cdot V$ are too small to be noticed compared to the yearly premium rate? I would just suspect the reserve plot to follow the movements of the GBM roughly?

• You are looking for a rate of return $r(t)$, computing $(S_{t+\Delta t}-S_{t})/S_{t}$ gives you the discrete return $r(t) \Delta t$ not $r(t)$. Take the degenerate situation where the premium rate and mortality rate are zero, + zero volatility, to check your results. Starting from the same initial value, you should get that the risky asset path is the same as $V$ (exact if you are working with log-returns vs. arithmetic returns). Commented Oct 5, 2017 at 9:24
• Thank you for your answer, I thought the Idea was quite good, so I did as told. The result is as expected, the degenerate case does not equal the path of the GBM, but the shape is resembled. I am aware that I use discerete return, but are there other options? I have a strong feeling that I should be able to project the differential equation using simulated returns. The plots that sould be equal are found in the link below. It is seen that the shape is equal but the value not. Is there an alternative to using discrete returns? imgur.com/HPQEOEL Commented Oct 5, 2017 at 15:26
• The plots do not seem to start from the same values. If I understand well, you simulated a single realisation of some stochastic process $(S)_{t \geq 0}$ (GBM$(\mu, \sigma)$) which you are basically recasting as a deterministic time evolution $dS(t)/S(t) = r(t) dt$ where $r(t)$ now incorporates the Brownian shocks. By construction it should thus give the same result as $dV(t) = r(t)V(t) dt$ if you start from the same initial values. Also in that degenerate case there exists a close form solution. Commented Oct 5, 2017 at 15:40

$dV(t)=\alpha(t,V(t))dt+\sigma(t,V(t))dW(t)\\=(\alpha V(t)+\pi-\mu(x+t)(b_d-V(t)))dt+\sigma V(t)dW(t)$
• Glad it works for you! Since you were working with an ODE I assumed you were not familiar with SDEs but isn't your original problem more like: $dV(t) = (r(t)V(t)+ \pi - \mu(x+t)(b_d-V(t)))dt$ with $dr(t) = a dt + b dW(t)$? Commented Oct 6, 2017 at 7:18