###############################
## Lecture 7 (discrete time) ##
###############################
library(ggplot2)
source("K2013.r")
source("probabilityK2013.r")

##################################
#### Endowment discrete time #####
##################################

#Individual information
x0 <- 50
G <- 1
R <- 0
Y <- 2023

#Contract information
T <- 20
E1 <- 500000
E2 <- 2000000

#Interest
r <- 0.03 #interest
v <- function(n) exp(-r*n)

#Transition probabilities
mu <- function(s) my.kol(x0+s,Y+s)
p <- function(a,b) exp(-integrate(mu,a,b)$value)

#Single premium
sum1 <- function(n) v(n)*p(0,n)*(1-p(n,n+1))
pi0 <- v(T)*p(0,T)*E1+E2*sum(sapply(0:(T-1),sum1))
pi0

#Yearly premium
sum2 <- function(n) v(n)*p(0,n)
pi <- pi0/sum(sapply(0:(T-1),sum2))
pi

###########################################
## Present value and prospective reserve ##
###########################################

#PV premiums
tt <- 0:T
V.res <- V.pi <- V.E1 <- V.E2 <- tt
V.E1[T+1] <- E1

for(t in 0:(T-1)){
  aux1 <- 0
  for(n in t:(T-1)) aux1 <- aux1 + v(n)*p(t,n)
  V.pi[t+1] <- -pi*(1/v(t))*aux1
  V.E1[t+1] <- (1/v(t))*v(T)*p(t,T)*E1
  aux2 <- 0
  for(n in t:(T-1)) aux2 <- aux2 + v(n+1)*p(t,n)*(1-p(n,n+1))
  V.E2[t+1] <- E2*(1/v(t))*aux2
}
V.res <- V.pi+V.E1+V.E2

pV <- data.frame(tt,V.res/1000000,V.pi/1000000,V.E1/1000000,V.E2/1000000)
ggplot(data = pV, aes(x = tt)) +
  geom_point(aes(y = V.pi/1000000, colour = "PV premiums")) +
  geom_point(aes(y = V.E1/1000000, colour = "PV survival benefit")) +
  geom_point(aes(y = V.E2/1000000, colour = "PV death benefit")) +
  geom_point(aes(y = V.res/1000000, colour = "Total reserve")) +
  geom_hline(yintercept = 0, linetype="dotted",color = "black")+
  scale_colour_manual("", 
                      breaks = c("PV premiums","PV survival benefit","PV death benefit", "Total reserve"),
                      values = c("PV premiums"="green","PV survival benefit"="blue","PV death benefit"="red", 
                                 "Total reserve"="black")) +
  xlab("Policy time") +
  scale_x_continuous(breaks=seq(0,100, 20))+
  scale_y_continuous("Value (MNOK)") + 
  labs(title="Expected prospective values for endowment")

################################
#### Pension discrete time #####
################################

#Individual information
x0 <- 30
G <- 0
R <- 1
Y <- 2023

#Contract information
T0 <- 40
T <- 90
P <- 200000

#Interest
r <- 0.03 #interest
v <- function(n) exp(-r*n)

#Transition probabilities
mu <- function(s) my.kol(x0+s,Y+s)
p <- function(a,b) exp(-integrate(mu,a,b)$value)

#Single premium
sum1 <- function(n) v(n)*p(0,n)

pi0 <- P*sum(sapply(T0:(T-1),sum1))
pi0

#Yearly premium
pi <- pi0/sum(sapply((0:(T0-1)),sum1))
pi

###########################################
## Present value and prospective reserve ##
###########################################

tt <- 0:T
V.res <- V.prem <- V.pen <- 0*tt

for(t in 0:(T-1)){
  if(t<=(T0-1)){
    vpfun <- function(t,n) v(n)*p(t,n)
    V.prem[t+1] <- -pi*(1/v(t))*sum(sapply(t:(T0-1),vpfun, t=t))
    V.pen[t+1] <- P*(1/v(t))*sum(sapply(T0:(T-1),vpfun, t=t))
  }
  if(t>=T0){
    vpfun <- function(t,n) v(n)*p(t,n)
    V.pen[t+1] <- P*(1/v(t))*sum(sapply(t:(T-1), vpfun, t=t) )
  }
}

V.res <- V.prem + V.pen

pV <- data.frame(tt,V.res/1000000,V.prem/1000000,V.pen/1000000)
ggplot(data = pV, aes(x = tt)) +
  geom_point(aes(y = V.prem/1000000, colour = "PV premiums"),size=0.9) +
  geom_point(aes(y = V.pen/1000000, colour = "PV pensions"),size=0.9) +
  geom_point(aes(y = V.res/1000000, colour = "Total reserve"),size=0.3) +
  geom_vline(xintercept = 40, linetype="dotted",color = "black")+
  scale_colour_manual("", 
                      breaks = c("PV premiums","PV pensions", "Total reserve"),
                      values = c("PV premiums"="green","PV pensions"="red", 
                                 "Total reserve"="blue")) +
  xlab("Policy time") +
  scale_x_continuous(breaks=seq(0,100, 20))+
  scale_y_continuous("Value (MNOK)") + 
  labs(title="Expected prospective values for pension")

###################################
#### Disability discrete time #####
###################################

source("Lecture2_Euler.r")
x0 <- 30
T <- 20
D <- 100000
P0 <- diag(3)
h <- 1/12
N <- T/h

Euler_sol <- Euler(x0,P0,h,x0+T)

p.act.act <- function(t,s) Euler(x0+t,P0,h,x0+s)[1,1,((s-t)/h)+1]
p.act.dis <- function(t,s) Euler(x0+t,P0,h,x0+s)[1,2,((s-t)/h)+1]
p.act.dead <- function(t,s) 1-p.act.act(t,s)-p.act.dis(t,s)
p.dis.act <- function(t,s) Euler(x0+t,P0,h,x0+s)[2,1,((s-t)/h)+1]
p.dis.dis <- function(t,s) Euler(x0+t,P0,h,x0+s)[2,2,((s-t)/h)+1]
p.dis.dead <- function(t,s) 1-p.dis.act(t,s)-p.dis.dis(t,s)

#Interest
r <- 0.03 #interest
v <- function(n) exp(-r*n)

dis.fun <- function(n) v(n)*p.act.dis(0,n)
pi0 <- D*sum(sapply(0:(T-1),dis.fun))
pi0
act.fun <- function(n) v(n)*p.act.act(0,n)
pi <- pi0/sum(sapply(0:(T-1),act.fun))
pi
###########################################
## Present value and prospective reserve ##
###########################################
tt <- 0:T
V.res <- V.prem <- V.act <- V.dis <- 0*tt

T0 <- 10 #premium waiver, set equal to T for no premium waiver.
for(t in 0:(T-1)){
    vactpremfun <- function(t,n) v(n)*p.act.act(t,n)
    if(t<T0) V.prem[t+1] <- -pi*(1/v(t))*sum(sapply(t:(T0-1),vactpremfun, t=t))
    vdisfun <- function(t,n) v(n)*p.act.dis(t,n)
    V.act[t+1] <- D*(1/v(t))*sum(sapply(t:(T-1),vdisfun, t=t))
    vdisfun2 <- function(t,n) v(n)*p.dis.dis(t,n)
    V.dis[t+1] <- D*(1/v(t))*sum(sapply(t:(T-1),vdisfun2, t=t))
}

V.res <- V.prem + V.act


pV <- data.frame(tt,V.res,V.prem,V.act)
ggplot(data = pV, aes(x = tt)) +
  geom_point(aes(y = V.prem, colour = "PV premiums")) +
  geom_point(aes(y = V.act, colour = "PV disability")) +
  geom_point(aes(y = V.res, colour = "Total reserve")) +
  geom_hline(yintercept = 0, linetype="dotted",color = "black")+
  scale_colour_manual("", 
                      breaks = c("PV premiums","PV disability", "Total reserve"),
                      values = c("PV premiums"="green","PV disability"="red", 
                                 "Total reserve"="blue")) +
  xlab("Policy time") +
  scale_x_continuous(breaks=seq(0,100, 10))+
  scale_y_continuous("Value (NOK)") + 
  labs(title="Expected prospective values for pension")