#Example 9.5 (Explicit PDE 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 <- 1
mu <- 0.07
sig <- 0.2

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


########################################################
#step size
dt <- 1/12
dx <- 0.3
x.max <- 6 #maximum value of fund
#grid
time <- seq(0,T,by=dt)
space <- seq(0,x.max,by=dx)
Nt <- length(time)-1
Nx <- length(space)-1

#Solving the PDE:

#boundary conditions
init.x <- function(x) max(G,x)
boundary.down <- function(t) surv_prob(x0+t,x0+T)*G*exp(-r*(T-t))
boundary.up <- function(t) surv_prob(x0+t,x0+T)*x.max

v <- matrix(rep(0,(Nt+1)*(Nx+1)), nrow=Nt+1)

v[Nt+1,] <- sapply(space,init.x)
v[,1] <- sapply(time,boundary.down)
v[,Nx+1] <- sapply(time,boundary.up)

#with varying interest
#r.curve <- 0.03 + runif(Nt+1,-0.03,0.03)
#plot(time, r.curve, type="l",col="blue")

r.curve <- rep(r,Nt+1)

for(i in Nt:1){
  print(i)
  for(j in Nx:2){
    v[i,j] <- v[i+1,j]-dt*(r.curve[i+1]*v[i+1,j]+mort(x0+time[i+1])*v[i+1,j]-( r.curve[i+1]*space[j]*((v[i+1,j+1]-v[i+1,j])/dx) )-0.5*(space[j]^2)*(sig^2)*((v[i+1,j+1]-2*v[i+1,j]+v[i+1,j-1])/(dx^2) ) )
  }
}

persp3D(time,space,v, theta=300, phi=25, expand=0.5, ticktype="detailed", clab = c("Value"), xlab="years", ylab="fund", zlab="", axes=TRUE)



##############################################################
#Using the direct formula
#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))) }

#Explicit reserve:
Reserve <- function(t,x){ return( surv_prob(x0+t,x0+T)*BSE(t,T,x,G) ) }

surf.res <- matrix(rep(0,(Nt+1)*(Nx+1)), nrow=Nt+1)
for(i in 1:(Nt+1)){
  for(j in 1:(Nx+1)){
    surf.res[i,j] <- Reserve(time[i],space[j])
  }
}
persp3D(time,space,surf.res, theta=300, phi=25, expand=0.5, ticktype="detailed", clab = c("Value"), xlab="years", ylab="fund", zlab="", axes=TRUE)
##############################################################

persp3D(time,space,v-surf.res, theta=300, phi=25, expand=0.5, ticktype="detailed", clab = c("Value"), xlab="years", ylab="fund", zlab="", axes=TRUE)