I am trying to solve the equation


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


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

Scenario 1

Scenario 2

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?

  • $\begingroup$ 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). $\endgroup$
    – Quantuple
    Oct 5, 2017 at 9:24
  • $\begingroup$ 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 $\endgroup$ Oct 5, 2017 at 15:26
  • $\begingroup$ 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. $\endgroup$
    – Quantuple
    Oct 5, 2017 at 15:40

1 Answer 1


I found a solution, but I am not exactly sure why what I did was wrong.

Never the less, the solution was, using the R package Sim.DiffProc to simulate the Stochastic Diffusion Process

$ 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) $

Thanks to Quantuple for leading me in the right direction.

  • $\begingroup$ 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)$? $\endgroup$
    – Quantuple
    Oct 6, 2017 at 7:18
  • $\begingroup$ I am doing my masters thesis in actuarial mathematics, so I know plenty of advanced probability and continuous time finance. I just programmed it in a deterministic rate scenario because I could then check my programme was correct, and then wanted to extend it to stochastic returns, but I realized late that I transitioned from an ODE to an SDE and hence could not use non-Ito calculus anymore. $\endgroup$ Oct 6, 2017 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.