#Example 9.1 (Term Insurance Unit-linked Black-Scholes)
require(plot3D)

#Mortality used in the example
#Gompertz-Makeham law of mortality parameters
a <- -9.13275
b <- 0.0809438
c <- 0.000011018

#Gompertz-Makeham law
mort <- function(u){
  return(exp(a+b*u-c*u^2))
}

#Survival probability: s to t
surv_prob <- function(s,t){
  mu <- b/(2*c)
  sigma <- sqrt(1/(2*c))
  #if(s>t){ return("s>t") }
  val <- exp(-sigma*exp(a+(b*b)/(4*c))*sqrt(2*pi)*(pnorm((t-mu)/sigma,0,1)-pnorm((s-mu)/sigma,0,1)))
  return(val)
}

#interest
r <- 0.03

#Black-Scholes specifications
mu <- 0.07
sig <- 0.2

#contract specifications
x0 <- 30 #age
T <- 10
G <- 1 #guarantee

#Black-Scholes function

BS <- function(t,T,x,K){
  d1 <- (log(x/K)+(r+0.5*sig*sig)*(T-t))/(sig*sqrt(T-t))
  d2 <- (log(x/K)+(r-0.5*sig*sig)*(T-t))/(sig*sqrt(T-t))
  
  return(x*pnorm(d1)-K*exp(-r*(T-t))*pnorm(d2))
}

BSE <- function(t,T,x,K){ return(BS(t,T,x,K)+K*exp(-r*(T-t)))}


#surface for pure endowment
h <- 1/12
time <- seq(0,T,by=h)
n.t <- length(time)-1
fund <- seq(0,2,by=0.1)
surf <- matrix(rep(0,(length(time))*(length(fund))), nrow=length(time))

for(i in 1:length(time)){
  for(j in 1:length(fund)){
    int1 <- function(s){
      val <- surv_prob(x0+time[i],x0+s)*mort(x0+s)*BSE(time[i],s,fund[j],G)
      return(val)
    }
    surf[i,j] <- as.numeric(integrate(int1,time[i],T)[1])
  }
}

persp3D(time,fund,surf, theta=320, phi=30, expand=0.5, ticktype="detailed", clab = c("Value"), xlab="time (years)", ylab="fund", zlab="", axes=TRUE)

#Simulating stochastic reserves
#We know that W_t-W_s is indep of W_s, s<t.
noise <- c(0,cumsum(rnorm(n.t,0,h)))
stock <- S0*exp((mu-0.5*sig*sig)*time+sig*noise)

stoch.V <- stock

for(i in 1:length(time)){
    int1 <- function(s){
      val <- surv_prob(x0+time[i],x0+s)*mort(x0+s)*BSE(time[i],s,stock[i],G)
      return(val)
    }
    stoch.V[i] <- as.numeric(integrate(int1,time[i],T)[1])
}

plot(time,stoch.V, type="l", col="blue")