#Example 9.1 (Pension 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
S0 <- 1
mu <- 0.07
sig <- 0.2

#contract specifications
x0 <- 30 #age
T0 <- 40
T <- 90
P <- 1 #pension

#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)))}

#yearly premium
int1 <- function(s){
  return( surv_prob(x0,x0+s)*BSE(0,s,S0,P) )
}
int2 <- function(s){
  return( exp(-r*s)*surv_prob(x0,x0+s) )
}

pi0 <- as.numeric(integrate(int1,T0,T)[1])/as.numeric(integrate(int2,0,T0)[1])

#surface for pension
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)){
    int3 <- function(s){
      val <- exp(-r*s)*surv_prob(x0+time[i],x0+s)
      return(val)
    }
    int4 <-function(s){
      val <- surv_prob(x0+time[i],x0+s)*BSE(time[i],s,fund[j],G)
    }
    surf[i,j] <- -pi0*exp(r*time[i])*as.numeric(integrate(int3,time[i],max(T0,time[i]))[1])+as.numeric(integrate(int4,max(time[i],T0),T)[1])
  }
}

persp3D(time,fund,surf, theta=320, phi=20, expand=0.4, 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)){
  int3 <- function(s){
    val <- exp(-r*s)*surv_prob(x0+time[i],x0+s)
    return(val)
  }
  int4 <-function(s){
    val <- surv_prob(x0+time[i],x0+s)*BSE(time[i],s,stock[i],G)
  }
  stoch.V[i] <- -pi0*as.numeric(integrate(int3,time[i],max(T0,time[i]))[1])+as.numeric(integrate(int4,max(time[i],T0),T)[1])
}
plot(time,stoch.V, type="l", col="red")
abline(v=40)