Projecting a Thiele differential equation with Black Scholes returns

Question

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?

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.