#Example 9.1 (Pure Endowment Unit-linked Black-Scholes)
require(plot3D)
require(MASS)

#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
S0 <- 0.5
mu <- 0.07
sig <- 0.2

#contract specifications
x0 <-30
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((surv_prob(x0+t,x0+T))*(BS(t,T,x,K)+K*exp(-r*(T-t))))}

#surface for pure endowment
h <- 1/12
time <- seq(0,T,by=1/12)
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)){
    surf[i,j] <- BSE(time[i],T,fund[j],G)
  }
}


persp3D(time,fund,surf, theta=300, phi=25, 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)
plot(time,stock, type="l", col="red")
plot(time,BSE(time,T,stock,G), type="l", col="blue")

#many paths
noise <- c(0,cumsum(rnorm(n.t,0,h)))
stock <- S0*exp((mu-0.5*sig*sig)*time+sig*noise)
plot(time,BSE(time,T,stock,G), type="l", col="blue")
for(iter in 1:10){
  noise <- c(0,cumsum(rnorm(n.t,0,h)))
  stock <- S0*exp((mu-0.5*sig*sig)*time+sig*noise)
  lines(time,BSE(time,T,stock,G), type="l", col=iter)
}