###############
## Lecture 6 ##
###############
library(ggplot2)
source("K2013.r")
source("probabilityK2013.r")

##################
#### Pension #####
##################

#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(t) exp(-r*t)

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

#Single premium
int <- Vectorize(function(ss) v(ss)*p(0,ss) )

pi0 <- P*integrate(int,T0,T)$value
pi0

#Yearly premium
pi <- pi0/integrate(int,0,T0)$value
pi

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

#PV premiums
h <- 1/12 #monthly
tt <- seq(0,T,by=h)

V.res <- V.prem <- V.pen <- tt

for(t in 1:length(tt)){
  V.prem[t] <- -pi*(1/v((t-1)*h))*(integrate(int,(t-1)*h,T0)$value)*1*(0<=((t-1)*h))*(((t-1)*h)<T0)
  V.pen[t] <- P*(1/v((t-1)*h))*(integrate(int,T0,T)$value)*1*(0<=((t-1)*h))*(((t-1)*h)<T0)+P*(1/v((t-1)*h))*(integrate(int,(t-1)*h,T)$value)*1*(T0<=((t-1)*h))*(((t-1)*h)<T)
}

pV <- data.frame(tt,V.res/1000000,V.prem/1000000,V.pen/1000000)
ggplot(data = pV, aes(x = tt)) +
  geom_line(aes(y = V.prem/1000000, colour = "PV premiums")) +
  geom_line(aes(y = V.pen/1000000, colour = "PV pensions"),size=0.9) +
  geom_line(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")

####################
#### Endowment #####
####################

#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(t) exp(-r*t)

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

#Single premium
int1 <- Vectorize(function(ss) v(ss)*p(0,ss)*mu(ss) )
int2 <- Vectorize(function(ss) v(ss)*p(0,ss) )
pi0 <- v(T)*E1*p(0,T) + E2*integrate(int1,0,T)$value
pi0
pi <- pi0/integrate(int2,0,T)$value
pi

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

#PV premiums
h <- 1/12 #monthly
tt <- seq(0,T,by=h)

V.res <- V.prem <- V.E1 <- V.E2 <- tt

for(t in 1:length(tt)){
  V.prem[t] <- -pi*(1/v((t-1)*h))*(integrate(int2,(t-1)*h,T)$value)*1*(0<=((t-1)*h))*(((t-1)*h)<T)
  V.E1[t] <- (v(T)/v((t-1)*h))*E1*p((t-1)*h,T)
  V.E2[t] <- E2*(1/v((t-1)*h))*(integrate(int1,(t-1)*h,T)$value)*1*(0<=((t-1)*h))*(((t-1)*h)<T)
}

V.res <- V.prem+V.E1+V.E2

pV <- data.frame(tt,V.res/1000000,V.prem/1000000,V.E1/1000000,V.E2/1000000)
ggplot(data = pV, aes(x = tt)) +
  geom_line(aes(y = V.prem/1000000, colour = "PV premiums")) +
  geom_line(aes(y = V.E1/1000000, colour = "PV survival benefit")) +
  geom_line(aes(y = V.E2/1000000, colour = "PV death benefit")) +
  geom_line(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")