G <- 1 #Gender (G=1 for women)
R <- 0 #Type of risk (R=0 for mortality, R=1 for longevity)
Y <- 2023 #Year

r <- 0.03 #interest

v <- function(t) exp(-r*t)

#contract info:
age0 <- x0 <- 50 #Initial age
T <- 20 #contract time
E1 <- 500000 #survival benefit
E2 <- 2000000 #death benefit
n <- 1000 #simulations
ages <- life.gen(age0,n) #n random life-times for age0-year old
cost <- c(1:n) #we save cost of each insurance outcome

for(i in 1:n){
  cost[i] <- (E1*v(T)*1*(ages[i]>=(age0+T))+v(ages[i]-age0)*E2*1*(ages[i]<(age0+T)))
}

mean(cost)
# ((mean(cost)*exp(r*T)-E1)*p.K13(0,T))/(1-p.K13(0,T))
# should be close to what we need to cover death benefits

par(mfrow=c(2,2))
hist(cost, col="red", xlab="PV payments", main="Histogram")
boxplot(cost, col="red", ylab="PV payments", main="Boxplot")

#costs that are not "survival"
#par(mfrow=c(1,2))
hist(cost[cost>E1*v(T)], col="red", xlab="PV death benefits", main="Histogram")
boxplot(cost[cost>E1*v(T)], col="red", ylab="PV death benefits", main="Boxplot")

#proportion of deaths should be approx:
1-p.K13(0,T)

##Theoretical premium (not proved in the lecture yet)
integrand <- function(s) v(s)*p.K13(0,s)*my.kol(x0+s,Y+s)

prem <- v(T)*E1*p.K13(0,T)+E2*(1/12)*sum(sapply(seq(0,T,by=1/12),integrand))