options(encoding = 'UTF-8')
install.packages("distr")
install.packages("distrEx")
install.packages("ggplot2")

rm(list=ls())

library("distrEx")
library("distr")
library("ggplot2")

#n=100
n=100 #Number of insured
lmax=n*0.25*1000+n*0.5*2000+n*0.25*3000 #Maximum cost

k<-seq(1000,lmax,by=1000) #All cost values

bin1 <- 1000*Binom(n*0.25,0.03) #First category
bin2 <- 2000*Binom(n*0.5,0.025) #Second category
bin3 <- 3000*Binom(n*0.25,0.015) #Third category

bin <- bin1 + bin2 + bin3 #Portfolio distribution

plot(k,p(bin)(k),type="s",main="Distribution function for n=100",
     xlab="Cost",ylab="Cumulative probability", col = "blue")

#Removing extrem values
eps <- 10^(-6)/2
low <- q(bin)(eps)
up <- q(bin)(1-eps)

k <- seq(low,up,1000)

plot(k,d(bin)(k),type="l",main="Density function for n=100",
     xlab="Cost",ylab="Probability")

#Stop-Loss function
f <- function(x, k) {
  ifelse(x-k > 0, x-k, 0)
}

stoploss <- c()

for (i in 1:(up/1000-low/1000+1)) {
  stoploss[i] <- sum(f(k,(i-1)*1000)*d(bin)(k))
}

plot(k,stoploss,type="b",main="Stop-Loss for n=100",
     xlab="K",ylab="Pure Premium")

#Participation

bin1m <- 1000*Binom(floor((n-1)*0.25),0.03)
binm1 <- bin1m + bin2 + bin3

x1 <- c()

for (i in 1:((up-low)/1000+1)) {
  x1[i] <- (1000*d(binm1)((i-1)*1000-1000)*0.03)/d(bin)((i-1)*1000)
}

bin2m <- 2000*Binom(floor((n-1)*0.5),0.025)
binm2 <- bin1 + bin2m + bin3

x2 <- c()

for (i in 1:((up-low)/1000+1)) {
  x2[i] <- (2000*d(binm2)((i-1)*1000-2000)*0.025)/d(bin)((i-1)*1000)
}

bin3m <- 3000*Binom(floor((n-1)*0.25),0.015)
binm3 <- bin1m + bin2 + bin3m

x3 <- c()

for (i in 1:((up-low)/1000+1)) {
  x3[i] <- (3000*d(binm3)((i-1)*1000-3000)*0.015)/d(bin)((i-1)*1000)
}

plot(k,x1,type="l",main="Conditional mean risk sharing for n=100",
     xlab="Total Cost",ylab="Contribution", col = "black",ylim=c(0,200))
lines(k, x2, col = "blue", lty = 1)
lines(k,x3, col = "red", lty = 1)
lines(k,d(bin)(k)*1000,type="h",col="grey")
legend("topleft", legend = c("p=1000", "p=2000", "p=3000"), col = c("black", "blue", "red"), lty=1)

#n=1000
n=1000
lmax=n*0.25*1000+n*0.5*2000+n*0.25*3000

k<-seq(1000,lmax,by=1000)

bin1 <- 1000*Binom(n*0.25,0.03)
bin2 <- 2000*Binom(n*0.5,0.025)
bin3 <- 3000*Binom(n*0.25,0.015)

bin <- bin1 + bin2 + bin3
plot(k,p(bin)(k),type="s",main="Distribution function for n=1000",
     xlab="Cost",ylab="Cumulative probability", col="blue")

#Removing extrem values
eps <- 10^(-6)/2
low <- q(bin)(eps)
up <- q(bin)(1-eps)

k <- seq(low,up,1000)

plot(k,d(bin)(k),type="l",main="Density function for n=1000",
     xlab="Cost",ylab="Probability")

stoploss <- c()

for (i in 1:(up/1000-low/1000+1)) {
  stoploss[i] <- sum(f(k,(i-1)*1000+low)*d(bin)(k))
}

plot(k,stoploss,type="b",main="Stop-Loss for n=1000",
     xlab="K",ylab="Pure Premium")

bin1m <- 1000*Binom(floor((n-1)*0.25),0.03)
binm1 <- bin1m + bin2 + bin3

x1 <- c()

for (i in 1:((up-low)/1000+1)) {
  x1[i] <- (1000*d(binm1)((i-1)*1000+low-1000)*0.03)/d(bin)((i-1)*1000+low)
}

bin2m <- 2000*Binom(floor((n-1)*0.5),0.025)
binm2 <- bin1 + bin2m + bin3

x2 <- c()

for (i in 1:((up-low)/1000+1)) {
  x2[i] <- (2000*d(binm2)((i-1)*1000+low-2000)*0.025)/d(bin)((i-1)*1000+low)
}

bin3m <- 3000*Binom(floor((n-1)*0.25),0.015)
binm3 <- bin1 + bin2 + bin3m

x3 <- c()

for (i in 1:((up-low)/1000+1)) {
  x3[i] <- (3000*d(binm3)((i-1)*1000+low-3000)*0.015)/d(bin)((i-1)*1000+low)
}

plot(k,x1,type="l",main="Conditional mean risk sharing for n=1000",
     xlab="Total Cost",ylab="Contribution", col = "black",ylim=c(0,120))
lines(k, x2, col = "blue", lty = 1)
lines(k,x3, col = "red", lty = 1)
lines(k,d(bin)(k)*1000,type="h",col="grey")
legend("topleft", legend = c("p=1000", "p=2000", "p=3000"), col = c("black", "blue", "red"), lty=1)

w1000 <- q(bin)(0.8)

w11000 <- (1000*d(binm1)(w1000-1000)*0.03)/d(bin)(w1000)
w21000 <- (2000*d(binm2)(w1000-2000)*0.025)/d(bin)(w1000)
w31000 <- (3000*d(binm3)(w1000-3000)*0.015)/d(bin)(w1000)

SL11000 <- 1.5*sum(f(x1,w11000)*d(bin)(k))
SL21000 <- 1.5*sum(f(x2,w21000)*d(bin)(k))
SL31000 <- 1.5*sum(f(x3,w31000)*d(bin)(k))

purepremium1000 <- c(w11000+SL11000,w21000+SL21000,w31000+SL31000)

CB11000 <- sum(f(w11000,x1)*d(bin)(k))
CB21000 <- sum(f(w21000,x2)*d(bin)(k))
CB31000 <- sum(f(w31000,x3)*d(bin)(k))

#Second method

w <- seq(low,up,10)

stoploss2 <- c()

for (i in 1:(up/10-low/10+1)) {
  stoploss2[i] <- sum(f(w,(i-1)*10+low)*d(bin)(w))
}

Totalprem1000 <- w + 1.5*stoploss2

plot(w,Totalprem1000, type="l",main="Total premium according w for n=1000", ylab="Total Premium")

tab1000 <- data.frame(w,Totalprem1000)

pi1000 <- 40*n*0.25 + 50*n*0.5 + 70*n*0.25

w1000 <- tab1000[which(abs(Totalprem1000-pi1000)==min(abs(Totalprem1000-pi1000))),1] #Selecting the w so that it correspond to the pure premium
w1000 <- round(w1000,-3) #Can only be in thousand otherwise the probability will be 0

p(bin)(w1000)

w11000 <- (1000*d(binm1)(w1000-1000)*0.03)/d(bin)(w1000)
w21000 <- (2000*d(binm2)(w1000-2000)*0.025)/d(bin)(w1000)
w31000 <- (3000*d(binm3)(w1000-3000)*0.015)/d(bin)(w1000)

SL11000 <- 1.5*sum(f(x1,w11000)*d(bin)(k))
SL21000 <- 1.5*sum(f(x2,w21000)*d(bin)(k))
SL31000 <- 1.5*sum(f(x3,w31000)*d(bin)(k))

P2P_1 <- 40 - SL11000
P2P_2 <- 50 - SL21000
P2P_3 <- 70 - SL31000

CB11000 <- sum(f(w11000,x1)*d(bin)(k))
CB21000 <- sum(f(w21000,x2)*d(bin)(k))
CB31000 <- sum(f(w31000,x3)*d(bin)(k))

#n=10000
n=10000
lmax=n*0.25*1000+n*0.5*2000+n*0.25*3000

k<-seq(1000,lmax,by=1000)

bin1 <- 1000*Binom(n*0.25,0.03)
bin2 <- 2000*Binom(n*0.5,0.025)
bin3 <- 3000*Binom(n*0.25,0.015)

bin <- bin1 + bin2 + bin3
plot(k,p(bin)(k),type="s",main="Distribution function for n=10000",
     xlab="Cost",ylab="Cumulative probability", col = "blue")

#Removing extrem values
eps <- 10^(-6)/2
low <- q(bin)(eps)
up <- q(bin)(1-eps)

k <- seq(low,up,1000)

plot(k,d(bin)(k),type="l",main="Density function for n=10000",
     xlab="Cost",ylab="Probability")

stoploss <- c()

for (i in 1:(up/1000-low/1000+1)) {
  stoploss[i] <- sum(f(k,(i-1)*1000+low)*d(bin)(k))
}

plot(k,stoploss,type="b",main="Stop-Loss for n=10000",
     xlab="K",ylab="Pure Premium")

bin1m <- 1000*Binom(floor((n-1)*0.25),0.03)
binm1 <- bin1m + bin2 + bin3

x1 <- c()

for (i in 1:((up-low)/1000+1)) {
  x1[i] <- (1000*d(binm1)((i-1)*1000+low-1000)*0.03)/d(bin)((i-1)*1000+low)
}

bin2m <- 2000*Binom(floor((n-1)*0.5),0.025)
binm2 <- bin1 + bin2m + bin3

x2 <- c()

for (i in 1:((up-low)/1000+1)) {
  x2[i] <- (2000*d(binm2)((i-1)*1000+low-2000)*0.025)/d(bin)((i-1)*1000+low)
}

bin3m <- 3000*Binom(floor((n-1)*0.25),0.015)
binm3 <- bin1 + bin2 + bin3m

x3 <- c()

for (i in 1:((up-low)/1000+1)) {
  x3[i] <- (3000*d(binm3)((i-1)*1000+low-3000)*0.015)/d(bin)((i-1)*1000+low)
}

plot(k,x1,type="l",main="Conditional mean risk sharing for n=10000",
     xlab="Total Cost",ylab="Contribution", col = "black",ylim=c(0,80))
lines(k, x2, col = "blue", lty = 1)
lines(k,x3, col = "red", lty = 1)
lines(k,d(bin)(k)*1000,type="h",col="grey")
legend("topleft", legend = c("p=1000", "p=2000", "p=3000"), col = c("black", "blue", "red"), lty=1)

w10000 <- q(bin)(0.8)

w110000 <- (1000*d(binm1)(w10000-1000)*0.03)/d(bin)(w10000)
w210000 <- (2000*d(binm2)(w10000-2000)*0.025)/d(bin)(w10000)
w310000 <- (3000*d(binm3)(w10000-3000)*0.015)/d(bin)(w10000)

SL110000 <- 1.5*sum(f(x1,w110000)*d(bin)(k))
SL210000 <- 1.5*sum(f(x2,w210000)*d(bin)(k))
SL310000 <- 1.5*sum(f(x3,w310000)*d(bin)(k))

purepremium10000 <- c(w110000+SL110000,w210000+SL210000,w310000+SL310000)

CB110000 <- sum(f(w110000,x1)*d(bin)(k))
CB210000 <- sum(f(w210000,x2)*d(bin)(k))
CB310000 <- sum(f(w310000,x3)*d(bin)(k))

#Second method

w <- seq(low,up,10)

stoploss2 <- c()

for (i in 1:(up/10-low/10+1)) {
  stoploss2[i] <- sum(f(w,(i-1)*10+low)*d(bin)(w))
}

Totalprem10000 <- w + 1.5*stoploss2

plot(w,Totalprem10000, type="l",main="Total premium according w for n=10000", ylab="Total Premium")

tab10000 <- data.frame(w,Totalprem10000)

pi10000 <- 40*n*0.25 + 50*n*0.5 + 70*n*0.25

w10000 <- tab10000[which(abs(Totalprem10000-pi10000)==min(abs(Totalprem10000-pi10000))),1] #Selecting the w so that it correspond to the pure premium
w10000 <- round(w10000,-3) #Can only be in thousand otherwise the probability will be 0

p(bin)(w10000)

w110000 <- (1000*d(binm1)(w10000-1000)*0.03)/d(bin)(w10000)
w210000 <- (2000*d(binm2)(w10000-2000)*0.025)/d(bin)(w10000)
w310000 <- (3000*d(binm3)(w10000-3000)*0.015)/d(bin)(w10000)

SL110000 <- 1.5*sum(f(x1,w110000)*d(bin)(k))
SL210000 <- 1.5*sum(f(x2,w210000)*d(bin)(k))
SL310000 <- 1.5*sum(f(x3,w310000)*d(bin)(k))

P2P_1 <- 40 - SL110000
P2P_2 <- 50 - SL210000
P2P_3 <- 70 - SL310000

CB110000 <- sum(f(w110000,x1)*d(bin)(k))
CB210000 <- sum(f(w210000,x2)*d(bin)(k))
CB310000 <- sum(f(w310000,x3)*d(bin)(k))

#Variance E[(Xi|Sn)] en fonction de n
#x1
var1 <- c()
for (n in seq(1000,30000,by=1000)) {
  bin1 <- 1000*Binom(n*0.25,0.03)
  bin2 <- 2000*Binom(n*0.5,0.025)
  bin3 <- 3000*Binom(n*0.25,0.015)
  
  bin <- bin1 + bin2 + bin3
  #Removing extreme values
  low <- q(bin)(eps)
  up <- q(bin)(1-eps)
  k <- seq(low,up,1000)
  bin1m <- 1000*Binom(round((n-1)*0.25,0),0.03)
  binm1 <- bin1m + bin2 + bin3
  esp1 <- c()
  for (i in 1:((up-low)/1000+1)) {
    esp1[i] <- (1000*d(binm1)((i-1)*1000+low-1000)*0.03)/d(bin)((i-1)*1000+low)
  }
  var1[(n/1000)] <- sum(esp1^2 * d(bin)(k)) - (sum(esp1*d(bin)(k)))^2
}

#x2
var2 <- c()
for (n in seq(1000,30000,by=1000)) {
  bin1 <- 1000*Binom(n*0.25,0.03)
  bin2 <- 2000*Binom(n*0.5,0.025)
  bin3 <- 3000*Binom(n*0.25,0.015)
  
  bin <- bin1 + bin2 + bin3
  #Removing extreme values
  low <- q(bin)(eps)
  up <- q(bin)(1-eps)
  k <- seq(low,up,1000)
  bin2m <- 2000*Binom(round((n-1)*0.5,0),0.025)
  binm2 <- bin1 + bin2m + bin3
  esp1 <- c()
  for (i in 1:((up-low)/1000+1)) {
    esp1[i] <- (2000*d(binm2)((i-1)*1000+low-2000)*0.025)/d(bin)((i-1)*1000+low)
  }
  var2[(n/1000)] <- sum(esp1^2 * d(bin)(k)) - (sum(esp1*d(bin)(k)))^2
}

#x3
var3 <- c()
for (n in seq(1000,30000,by=1000)) {
  bin1 <- 1000*Binom(n*0.25,0.03)
  bin2 <- 2000*Binom(n*0.5,0.025)
  bin3 <- 3000*Binom(n*0.25,0.015)
  
  bin <- bin1 + bin2 + bin3
  #Removing extreme values
  low <- q(bin)(eps)
  up <- q(bin)(1-eps)
  k <- seq(low,up,1000)
  bin3m <- 3000*Binom(round((n-1)*0.25,0),0.015)
  binm3 <- bin1 + bin2 + bin3m
  esp1 <- c()
  for (i in 1:((up-low)/1000+1)) {
    esp1[i] <- (3000*d(binm3)((i-1)*1000+low-3000)*0.015)/d(bin)((i-1)*1000+low)
  }
  var3[(n/1000)] <- sum(esp1^2 * d(bin)(k)) - (sum(esp1*d(bin)(k)))^2
}

k <- seq(1000,30000,by=1000)

plot(k,var1,type="l",main="var(E[(Xi|Sn)])",
     xlab="Number of insured",ylab="var(E[(Xi|Sn)])", col = "black",ylim=c(0,200))
lines(k, var2, col = "blue", lty = 1)
lines(k,var3, col = "red", lty = 1)
legend("topright", legend = c("p=1000", "p=2000", "p=3000"), col = c("black", "blue", "red"), lty=1)

plot(k,sqrt(var1),type="l",main="Standard deviation of E[(Xi|Sn))",
     xlab="Number of insured",ylab="Standard deviation", col = "black",ylim=c(0,40))
lines(k, sqrt(var2), col = "blue", lty = 1)
lines(k,sqrt(var3), col = "red", lty = 1)
legend("topright", legend = c("p=1000", "p=2000", "p=3000"), col = c("black", "blue", "red"), lty=1)

#Quantiles plots 
p25_1 <- c()
p50_1 <- c()
p75_1 <- c()
p25_2 <- c()
p50_2 <- c()
p75_2 <- c()
p25_3 <- c()
p50_3 <- c()
p75_3 <- c()
for (n in seq(1000,30000,by=1000)) {
  bin1 <- 1000*Binom(n*0.25,0.03)
  bin2 <- 2000*Binom(n*0.5,0.025)
  bin3 <- 3000*Binom(n*0.25,0.015)
  bin <- bin1 + bin2 + bin3
  q25 <- q(bin)(0.25)
  q50 <- q(bin)(0.5)
  q75 <- q(bin)(0.75)
  bin1m <- 1000*Binom(floor((n-1)*0.25),0.03)
  binm1 <- bin1m + bin2 + bin3
  p25_1 <- c(p25_1,(1000*d(binm1)(q25-1000)*0.03)/d(bin)(q25))
  p50_1 <- c(p50_1,(1000*d(binm1)(q50-1000)*0.03)/d(bin)(q50))
  p75_1 <- c(p75_1,(1000*d(binm1)(q75-1000)*0.03)/d(bin)(q75))
  bin2m <- 2000*Binom(floor((n-1)*0.5),0.025)
  binm2 <- bin1 + bin2m + bin3
  p25_2 <- c(p25_2,(2000*d(binm2)(q25-2000)*0.025)/d(bin)(q25))
  p50_2 <- c(p50_2,(2000*d(binm2)(q50-2000)*0.025)/d(bin)(q50))
  p75_2 <- c(p75_2,(2000*d(binm2)(q75-2000)*0.025)/d(bin)(q75))
  bin3m <- 3000*Binom(floor((n-1)*0.25),0.015)
  binm3 <- bin1 + bin2 + bin3m
  p25_3 <- c(p25_3,(3000*d(binm3)(q25-3000)*0.015)/d(bin)(q25))
  p50_3 <- c(p50_3,(3000*d(binm3)(q50-3000)*0.015)/d(bin)(q50))
  p75_3 <- c(p75_3,(3000*d(binm3)(q75-3000)*0.015)/d(bin)(q75))
}

n <- seq(1000, 30000, by = 1000)

quantiles_data_1 <- data.frame(n, p25_1, p50_1, p75_1)

ggplot(quantiles_data_1, aes(x = n)) +
  geom_ribbon(aes(ymin = p25_1, ymax = p75_1), fill = "blue", alpha = 0.3) +
  geom_line(aes(y = p50_1), color = "blue") +
  labs(title = "Evolution of quantiles as a function of the number of policyholders for p=1000",
       x = "Number of insured",
       y = "Participation") +
  theme_minimal()


quantiles_data_2 <- data.frame(n, p25_2, p50_2, p75_2)

ggplot(quantiles_data_2, aes(x = n)) +
  geom_ribbon(aes(ymin = p25_2, ymax = p75_2), fill = "blue", alpha = 0.3) +
  geom_line(aes(y = p50_2), color = "blue") +
  labs(title = "Evolution of quantiles as a function of the number of policyholders for p=2000",
       x = "Number of insured",
       y = "Participation") +
  theme_minimal()

quantiles_data_3 <- data.frame(n, p25_3, p50_3, p75_3)

ggplot(quantiles_data_3, aes(x = n)) +
  geom_ribbon(aes(ymin = p25_3, ymax = p75_3), fill = "blue", alpha = 0.3) +
  geom_line(aes(y = p50_3), color = "blue") +
  labs(title = "Evolution of quantiles as a function of the number of policyholders for p=3000",
       x = "Number of insured",
       y = "Participation") +
  theme_minimal()


#Illustration

#Uniform distribution
#10000 insured
#Regular bikes

n=10000 #Number of insured
lmax=n/10*sum(seq(500,5000,500)) #Maximum cost

k<-seq(500,lmax,by=500) #All cost values

bin1 <- 500*Binom(n/10,0.0306)
bin2 <- 1000*Binom(n/10,0.0315)
bin3 <- 1500*Binom(n/10,0.024)
bin4 <- 2000*Binom(n/10,0.0218)
bin5 <- 2500*Binom(n/10,0.0204)
bin6 <- 3000*Binom(n/10,0.0195)
bin7 <- 3500*Binom(n/10,0.0189)
bin8 <- 4000*Binom(n/10,0.018)
bin9 <- 4500*Binom(n/10,0.018)
bin10 <- 5000*Binom(n/10,0.018)

bin <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10 #Portfolio distribution

plot(k,p(bin)(k),type="s",main="Distribution function for n=10000 (uniform)",
     xlab="Cost",ylab="Cumulative probability", col = "blue")

#Removing extrem values
eps <- 10^(-6)/2
low <- q(bin)(eps)
up <- q(bin)(1-eps)

k <- seq(low,up,500)

#Premium calculation

bin1m <- 500*Binom(floor((n-1)/10),0.0306)
binm1 <- bin1m + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin2m <- 1000*Binom(floor((n-1)/10),0.0315)
binm2 <- bin1 + bin2m + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin3m <- 1500*Binom(floor((n-1)/10),0.024)
binm3 <- bin1 + bin2 + bin3m + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin4m <- 2000*Binom(floor((n-1)/10),0.0218)
binm4 <- bin1 + bin2 + bin3 + bin4m + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin5m <- 2500*Binom(floor((n-1)/10),0.0204)
binm5 <- bin1 + bin2 + bin3 + bin4 + bin5m + bin6 + bin7 + bin8 + bin9 + bin10

bin6m <- 3000*Binom(floor((n-1)/10),0.0195)
binm6 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6m + bin7 + bin8 + bin9 + bin10

bin7m <- 3500*Binom(floor((n-1)/10),0.0189)
binm7 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7m + bin8 + bin9 + bin10

bin8m <- 4000*Binom(floor((n-1)/10),0.018)
binm8 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8m + bin9 + bin10

bin9m <- 4500*Binom(floor((n-1)/10),0.018)
binm9 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9m + bin10

bin10m <- 5000*Binom(floor((n-1)/10),0.018)
binm10 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10m

x1 <- c()

for (i in 1:((up-low)/500+1)) {
  x1[i] <- (500*d(binm1)((i-1)*500+low-500)*0.0306)/d(bin)((i-1)*500+low)
}

x2 <- c()

for (i in 1:((up-low)/500+1)) {
  x2[i] <- (1000*d(binm2)((i-1)*500+low-1000)*0.0315)/d(bin)((i-1)*500+low)
}

x3 <- c()

for (i in 1:((up-low)/500+1)) {
  x3[i] <- (1500*d(binm3)((i-1)*500+low-1500)*0.024)/d(bin)((i-1)*500+low)
}

x4 <- c()

for (i in 1:((up-low)/500+1)) {
  x4[i] <- (2000*d(binm4)((i-1)*500+low-2000)*0.0218)/d(bin)((i-1)*500+low)
}

x5 <- c()

for (i in 1:((up-low)/500+1)) {
  x5[i] <- (2500*d(binm5)((i-1)*500+low-2500)*0.0204)/d(bin)((i-1)*500+low)
}

x6 <- c()

for (i in 1:((up-low)/500+1)) {
  x6[i] <- (3000*d(binm6)((i-1)*500+low-3000)*0.0195)/d(bin)((i-1)*500+low)
}

x7 <- c()

for (i in 1:((up-low)/500+1)) {
  x7[i] <- (3500*d(binm7)((i-1)*500+low-3500)*0.0189)/d(bin)((i-1)*500+low)
}

x8 <- c()

for (i in 1:((up-low)/500+1)) {
  x8[i] <- (4000*d(binm8)((i-1)*500+low-4000)*0.018)/d(bin)((i-1)*500+low)
}

x9 <- c()

for (i in 1:((up-low)/500+1)) {
  x9[i] <- (4500*d(binm9)((i-1)*500+low-4500)*0.018)/d(bin)((i-1)*500+low)
}

x10 <- c()

for (i in 1:((up-low)/500+1)) {
  x10[i] <- (5000*d(binm10)((i-1)*500+low-5000)*0.018)/d(bin)((i-1)*500+low)
}

#Premium calculation
w <- seq(low,up,10)

stoploss2 <- c()

for (i in 1:(up/10-low/10+1)) {
  stoploss2[i] <- sum(f(w,(i-1)*10+low)*d(bin)(w))
}

Totalprem10000 <- w + 1.5*stoploss2
options(scipen = 10)
plot(w,Totalprem10000, type="l",main="Total premium according the retention level for n=10000 (uniform)", ylab="Total Premium",xlab="Retention level")

tab10000 <- data.frame(w,Totalprem10000)

w10000_1 <- tab10000[Totalprem10000==min(Totalprem10000),1] #Selecting the w that minimize the premium value

w10000_2 <- q(bin)(0.7)

w10000_3 <- q(bin)(0.95)

purepremium10000_4 <- 1.2*c(500*0.0306,1000*0.0315 , 1500*0.024 , 2000*0.0218 , 2500*0.0204 , 
                         3000*0.0195 , 3500*0.0189 , 4000*0.018 , 4500*0.018 , 5000*0.018)

pi10000 <- 1000*sum(purepremium10000_4) #Total premium

w10000_4 <- tab10000[which(abs(Totalprem10000-pi10000)==min(abs(Totalprem10000-pi10000))),1]#Selecting the w so that the premium is equal to the pure premium

w10000_4 <- round(w10000_4,-3)

w110000_1 <- (500*d(binm1)(w10000_1-500)*0.0306)/d(bin)(w10000_1)
w210000_1 <- (1000*d(binm2)(w10000_1-1000)*0.0315)/d(bin)(w10000_1)
w310000_1 <- (1500*d(binm3)(w10000_1-1500)*0.024)/d(bin)(w10000_1)
w410000_1 <- (2000*d(binm4)(w10000_1-2000)*0.0218)/d(bin)(w10000_1)
w510000_1 <- (2500*d(binm5)(w10000_1-2500)*0.0204)/d(bin)(w10000_1)
w610000_1 <- (3000*d(binm6)(w10000_1-3000)*0.0195)/d(bin)(w10000_1)
w710000_1 <- (3500*d(binm7)(w10000_1-3500)*0.0189)/d(bin)(w10000_1)
w810000_1 <- (4000*d(binm8)(w10000_1-4000)*0.018)/d(bin)(w10000_1)
w910000_1 <- (4500*d(binm9)(w10000_1-4500)*0.018)/d(bin)(w10000_1)
w1010000_1 <- (5000*d(binm10)(w10000_1-5000)*0.018)/d(bin)(w10000_1)

SL110000_1 <- 1.5*sum(f(x1,w110000_1)*d(bin)(k))
SL210000_1 <- 1.5*sum(f(x2,w210000_1)*d(bin)(k))
SL310000_1 <- 1.5*sum(f(x3,w310000_1)*d(bin)(k))
SL410000_1 <- 1.5*sum(f(x4,w410000_1)*d(bin)(k))
SL510000_1 <- 1.5*sum(f(x5,w510000_1)*d(bin)(k))
SL610000_1 <- 1.5*sum(f(x6,w610000_1)*d(bin)(k))
SL710000_1 <- 1.5*sum(f(x7,w710000_1)*d(bin)(k))
SL810000_1 <- 1.5*sum(f(x8,w810000_1)*d(bin)(k))
SL910000_1 <- 1.5*sum(f(x9,w910000_1)*d(bin)(k))
SL1010000_1 <- 1.5*sum(f(x10,w1010000_1)*d(bin)(k))

purepremium10000_1 <- c(w110000_1+SL110000_1,w210000_1+SL210000_1,w310000_1+SL310000_1,w410000_1+SL410000_1,
                      w510000_1+SL510000_1,w610000_1+SL610000_1,w710000_1+SL710000_1,w810000_1+SL810000_1,
                      w910000_1+SL910000_1,w1010000_1+SL1010000_1)

#Expected cash-back

CB110000_1 <- sum(f(w110000_1,x1)*d(bin)(k))
CB210000_1 <- sum(f(w210000_1,x2)*d(bin)(k))
CB310000_1 <- sum(f(w310000_1,x3)*d(bin)(k))
CB410000_1 <- sum(f(w410000_1,x4)*d(bin)(k))
CB510000_1 <- sum(f(w510000_1,x5)*d(bin)(k))
CB610000_1 <- sum(f(w610000_1,x6)*d(bin)(k))
CB710000_1 <- sum(f(w710000_1,x7)*d(bin)(k))
CB810000_1 <- sum(f(w810000_1,x8)*d(bin)(k))
CB910000_1 <- sum(f(w910000_1,x9)*d(bin)(k))
CB1010000_1 <- sum(f(w1010000_1,x10)*d(bin)(k))

w110000_2 <- (500*d(binm1)(w10000_2-500)*0.0306)/d(bin)(w10000_2)
w210000_2 <- (1000*d(binm2)(w10000_2-1000)*0.0315)/d(bin)(w10000_2)
w310000_2 <- (1500*d(binm3)(w10000_2-1500)*0.024)/d(bin)(w10000_2)
w410000_2 <- (2000*d(binm4)(w10000_2-2000)*0.0218)/d(bin)(w10000_2)
w510000_2 <- (2500*d(binm5)(w10000_2-2500)*0.0204)/d(bin)(w10000_2)
w610000_2 <- (3000*d(binm6)(w10000_2-3000)*0.0195)/d(bin)(w10000_2)
w710000_2 <- (3500*d(binm7)(w10000_2-3500)*0.0189)/d(bin)(w10000_2)
w810000_2 <- (4000*d(binm8)(w10000_2-4000)*0.018)/d(bin)(w10000_2)
w910000_2 <- (4500*d(binm9)(w10000_2-4500)*0.018)/d(bin)(w10000_2)
w1010000_2 <- (5000*d(binm10)(w10000_2-5000)*0.018)/d(bin)(w10000_2)

SL110000_2 <- 1.5*sum(f(x1,w110000_2)*d(bin)(k))
SL210000_2 <- 1.5*sum(f(x2,w210000_2)*d(bin)(k))
SL310000_2 <- 1.5*sum(f(x3,w310000_2)*d(bin)(k))
SL410000_2 <- 1.5*sum(f(x4,w410000_2)*d(bin)(k))
SL510000_2 <- 1.5*sum(f(x5,w510000_2)*d(bin)(k))
SL610000_2 <- 1.5*sum(f(x6,w610000_2)*d(bin)(k))
SL710000_2 <- 1.5*sum(f(x7,w710000_2)*d(bin)(k))
SL810000_2 <- 1.5*sum(f(x8,w810000_2)*d(bin)(k))
SL910000_2 <- 1.5*sum(f(x9,w910000_2)*d(bin)(k))
SL1010000_2 <- 1.5*sum(f(x10,w1010000_2)*d(bin)(k))

purepremium10000_2 <- c(w110000_2+SL110000_2,w210000_2+SL210000_2,w310000_2+SL310000_2,w410000_2+SL410000_2,w510000_2+SL510000_2,
                      w610000_2+SL610000_2,w710000_2+SL710000_2,w810000_2+SL810000_2,w910000_2+SL910000_2,w1010000_2+SL1010000_2)

#Expected cash-back

CB110000_2 <- sum(f(w110000_2,x1)*d(bin)(k))
CB210000_2 <- sum(f(w210000_2,x2)*d(bin)(k))
CB310000_2 <- sum(f(w310000_2,x3)*d(bin)(k))
CB410000_2 <- sum(f(w410000_2,x4)*d(bin)(k))
CB510000_2 <- sum(f(w510000_2,x5)*d(bin)(k))
CB610000_2 <- sum(f(w610000_2,x6)*d(bin)(k))
CB710000_2 <- sum(f(w710000_2,x7)*d(bin)(k))
CB810000_2 <- sum(f(w810000_2,x8)*d(bin)(k))
CB910000_2 <- sum(f(w910000_2,x9)*d(bin)(k))
CB1010000_2 <- sum(f(w1010000_2,x10)*d(bin)(k))

w110000_3 <- (500*d(binm1)(w10000_3-500)*0.0306)/d(bin)(w10000_3)
w210000_3 <- (1000*d(binm2)(w10000_3-1000)*0.0315)/d(bin)(w10000_3)
w310000_3 <- (1500*d(binm3)(w10000_3-1500)*0.024)/d(bin)(w10000_3)
w410000_3 <- (2000*d(binm4)(w10000_3-2000)*0.0218)/d(bin)(w10000_3)
w510000_3 <- (2500*d(binm5)(w10000_3-2500)*0.0204)/d(bin)(w10000_3)
w610000_3 <- (3000*d(binm6)(w10000_3-3000)*0.0195)/d(bin)(w10000_3)
w710000_3 <- (3500*d(binm7)(w10000_3-3500)*0.0189)/d(bin)(w10000_3)
w810000_3 <- (4000*d(binm8)(w10000_3-4000)*0.018)/d(bin)(w10000_3)
w910000_3 <- (4500*d(binm9)(w10000_3-4500)*0.018)/d(bin)(w10000_3)
w1010000_3 <- (5000*d(binm10)(w10000_3-5000)*0.018)/d(bin)(w10000_3)

SL110000_3 <- 1.5*sum(f(x1,w110000_3)*d(bin)(k))
SL210000_3 <- 1.5*sum(f(x2,w210000_3)*d(bin)(k))
SL310000_3 <- 1.5*sum(f(x3,w310000_3)*d(bin)(k))
SL410000_3 <- 1.5*sum(f(x4,w410000_3)*d(bin)(k))
SL510000_3 <- 1.5*sum(f(x5,w510000_3)*d(bin)(k))
SL610000_3 <- 1.5*sum(f(x6,w610000_3)*d(bin)(k))
SL710000_3 <- 1.5*sum(f(x7,w710000_3)*d(bin)(k))
SL810000_3 <- 1.5*sum(f(x8,w810000_3)*d(bin)(k))
SL910000_3 <- 1.5*sum(f(x9,w910000_3)*d(bin)(k))
SL1010000_3 <- 1.5*sum(f(x10,w1010000_3)*d(bin)(k))

purepremium10000_3 <- c(w110000_3+SL110000_3,w210000_3+SL210000_3,w310000_3+SL310000_3,w410000_3+SL410000_3,w510000_3+SL510000_3,
                        w610000_3+SL610000_3,w710000_3+SL710000_3,w810000_3+SL810000_3,w910000_3+SL910000_3,w1010000_3+SL1010000_3)

#Expected cash-back

CB110000_3 <- sum(f(w110000_3,x1)*d(bin)(k))
CB210000_3 <- sum(f(w210000_3,x2)*d(bin)(k))
CB310000_3 <- sum(f(w310000_3,x3)*d(bin)(k))
CB410000_3 <- sum(f(w410000_3,x4)*d(bin)(k))
CB510000_3 <- sum(f(w510000_3,x5)*d(bin)(k))
CB610000_3 <- sum(f(w610000_3,x6)*d(bin)(k))
CB710000_3 <- sum(f(w710000_3,x7)*d(bin)(k))
CB810000_3 <- sum(f(w810000_3,x8)*d(bin)(k))
CB910000_3 <- sum(f(w910000_3,x9)*d(bin)(k))
CB1010000_3 <- sum(f(w1010000_3,x10)*d(bin)(k))

w110000_4 <- (500*d(binm1)(w10000_4-500)*0.0306)/d(bin)(w10000_4)
w210000_4 <- (1000*d(binm2)(w10000_4-1000)*0.0315)/d(bin)(w10000_4)
w310000_4 <- (1500*d(binm3)(w10000_4-1500)*0.024)/d(bin)(w10000_4)
w410000_4 <- (2000*d(binm4)(w10000_4-2000)*0.0218)/d(bin)(w10000_4)
w510000_4 <- (2500*d(binm5)(w10000_4-2500)*0.0204)/d(bin)(w10000_4)
w610000_4 <- (3000*d(binm6)(w10000_4-3000)*0.0195)/d(bin)(w10000_4)
w710000_4 <- (3500*d(binm7)(w10000_4-3500)*0.0189)/d(bin)(w10000_4)
w810000_4 <- (4000*d(binm8)(w10000_4-4000)*0.018)/d(bin)(w10000_4)
w910000_4 <- (4500*d(binm9)(w10000_4-4500)*0.018)/d(bin)(w10000_4)
w1010000_4 <- (5000*d(binm10)(w10000_4-5000)*0.018)/d(bin)(w10000_4)

SL110000_4 <- 1.5*sum(f(x1,w110000_4)*d(bin)(k))
SL210000_4 <- 1.5*sum(f(x2,w210000_4)*d(bin)(k))
SL310000_4 <- 1.5*sum(f(x3,w310000_4)*d(bin)(k))
SL410000_4 <- 1.5*sum(f(x4,w410000_4)*d(bin)(k))
SL510000_4 <- 1.5*sum(f(x5,w510000_4)*d(bin)(k))
SL610000_4 <- 1.5*sum(f(x6,w610000_4)*d(bin)(k))
SL710000_4 <- 1.5*sum(f(x7,w710000_4)*d(bin)(k))
SL810000_4 <- 1.5*sum(f(x8,w810000_4)*d(bin)(k))
SL910000_4 <- 1.5*sum(f(x9,w910000_4)*d(bin)(k))
SL1010000_4 <- 1.5*sum(f(x10,w1010000_4)*d(bin)(k))

P2P110000_4 <- purepremium10000_4[1] - SL110000_4
P2P210000_4 <- purepremium10000_4[2] - SL210000_4
P2P310000_4 <- purepremium10000_4[3] - SL310000_4
P2P410000_4 <- purepremium10000_4[4] - SL410000_4
P2P510000_4 <- purepremium10000_4[5] - SL510000_4
P2P610000_4 <- purepremium10000_4[6] - SL610000_4
P2P710000_4 <- purepremium10000_4[7] - SL710000_4
P2P810000_4 <- purepremium10000_4[8] - SL810000_4
P2P910000_4 <- purepremium10000_4[9] - SL910000_4
P2P1010000_4 <- purepremium10000_4[10] - SL1010000_4

#Expected cash-back

CB110000_4 <- sum(f(w110000_4,x1)*d(bin)(k))
CB210000_4 <- sum(f(w210000_4,x2)*d(bin)(k))
CB310000_4 <- sum(f(w310000_4,x3)*d(bin)(k))
CB410000_4 <- sum(f(w410000_4,x4)*d(bin)(k))
CB510000_4 <- sum(f(w510000_4,x5)*d(bin)(k))
CB610000_4 <- sum(f(w610000_4,x6)*d(bin)(k))
CB710000_4 <- sum(f(w710000_4,x7)*d(bin)(k))
CB810000_4 <- sum(f(w810000_4,x8)*d(bin)(k))
CB910000_4 <- sum(f(w910000_4,x9)*d(bin)(k))
CB1010000_4 <- sum(f(w1010000_4,x10)*d(bin)(k))


#Binomial (0.5) distribution
#10000 insured
#Regular bikes

n=10000 #Number of insured
lmax=sum(seq(500,5000,500)*round(dbinom(1:10,10,0.5)*n)) #Maximum cost

k<-seq(500,lmax,by=500) #All cost values

bin1 <- 500*Binom(round(dbinom(1,10,0.5)*n),0.0306)
bin2 <- 1000*Binom(round(dbinom(2,10,0.5)*n),0.0315)
bin3 <- 1500*Binom(round(dbinom(3,10,0.5)*n),0.024)
bin4 <- 2000*Binom(round(dbinom(4,10,0.5)*n),0.0218)
bin5 <- 2500*Binom(round(dbinom(5,10,0.5)*n),0.0204)
bin6 <- 3000*Binom(round(dbinom(6,10,0.5)*n),0.0195)
bin7 <- 3500*Binom(round(dbinom(7,10,0.5)*n),0.0189)
bin8 <- 4000*Binom(round(dbinom(8,10,0.5)*n),0.018)
bin9 <- 4500*Binom(round(dbinom(9,10,0.5)*n),0.018)
bin10 <- 5000*Binom(round(dbinom(10,10,0.5)*n),0.018)

bin <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10 #Portfolio distribution

plot(k,p(bin)(k),type="s",main="Distribution function for n=10000 (binomial)",
     xlab="Cost",ylab="Cumulative probability", col = "blue")

#Removing extrem values
eps <- 10^(-6)/2
low <- q(bin)(eps)
up <- q(bin)(1-eps)

k <- seq(low,up,500)

#Premium calculation

bin1m <- 500*Binom(round(dbinom(1,10,0.5)*n)-1,0.0306)
binm1 <- bin1m + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin2m <- 1000*Binom(round(dbinom(2,10,0.5)*n)-1,0.0315)
binm2 <- bin1 + bin2m + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin3m <- 1500*Binom(round(dbinom(3,10,0.5)*n)-1,0.024)
binm3 <- bin1 + bin2 + bin3m + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin4m <- 2000*Binom(round(dbinom(4,10,0.5)*n)-1,0.0218)
binm4 <- bin1 + bin2 + bin3 + bin4m + bin5 + bin6 + bin7 + bin8 + bin9 + bin10

bin5m <- 2500*Binom(round(dbinom(5,10,0.5)*n)-1,0.0204)
binm5 <- bin1 + bin2 + bin3 + bin4 + bin5m + bin6 + bin7 + bin8 + bin9 + bin10

bin6m <- 3000*Binom(round(dbinom(6,10,0.5)*n)-1,0.0195)
binm6 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6m + bin7 + bin8 + bin9 + bin10

bin7m <- 3500*Binom(round(dbinom(7,10,0.5)*n)-1,0.0189)
binm7 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7m + bin8 + bin9 + bin10

bin8m <- 4000*Binom(round(dbinom(8,10,0.5)*n)-1,0.018)
binm8 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8m + bin9 + bin10

bin9m <- 4500*Binom(round(dbinom(9,10,0.5)*n)-1,0.018)
binm9 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9m + bin10

bin10m <- 5000*Binom(round(dbinom(10,10,0.5)*n)-1,0.018)
binm10 <- bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9 + bin10m

x1 <- c()

for (i in 1:((up-low)/500+1)) {
  x1[i] <- (500*d(binm1)((i-1)*500+low-500)*0.0306)/d(bin)((i-1)*500+low)
}

x2 <- c()

for (i in 1:((up-low)/500+1)) {
  x2[i] <- (1000*d(binm2)((i-1)*500+low-1000)*0.0315)/d(bin)((i-1)*500+low)
}

x3 <- c()

for (i in 1:((up-low)/500+1)) {
  x3[i] <- (1500*d(binm3)((i-1)*500+low-1500)*0.024)/d(bin)((i-1)*500+low)
}

x4 <- c()

for (i in 1:((up-low)/500+1)) {
  x4[i] <- (2000*d(binm4)((i-1)*500+low-2000)*0.0218)/d(bin)((i-1)*500+low)
}

x5 <- c()

for (i in 1:((up-low)/500+1)) {
  x5[i] <- (2500*d(binm5)((i-1)*500+low-2500)*0.0204)/d(bin)((i-1)*500+low)
}

x6 <- c()

for (i in 1:((up-low)/500+1)) {
  x6[i] <- (3000*d(binm6)((i-1)*500+low-3000)*0.0195)/d(bin)((i-1)*500+low)
}

x7 <- c()

for (i in 1:((up-low)/500+1)) {
  x7[i] <- (3500*d(binm7)((i-1)*500+low-3500)*0.0189)/d(bin)((i-1)*500+low)
}

x8 <- c()

for (i in 1:((up-low)/500+1)) {
  x8[i] <- (4000*d(binm8)((i-1)*500+low-4000)*0.018)/d(bin)((i-1)*500+low)
}

x9 <- c()

for (i in 1:((up-low)/500+1)) {
  x9[i] <- (4500*d(binm9)((i-1)*500+low-4500)*0.018)/d(bin)((i-1)*500+low)
}

x10 <- c()

for (i in 1:((up-low)/500+1)) {
  x10[i] <- (5000*d(binm10)((i-1)*500+low-5000)*0.018)/d(bin)((i-1)*500+low)
}

#Premium calculation
w <- seq(low,up,10)

stoploss2 <- c()

for (i in 1:(up/10-low/10+1)) {
  stoploss2[i] <- sum(f(w,(i-1)*10+low)*d(bin)(w))
}

Totalprem10000 <- w + 1.5*stoploss2

plot(w,Totalprem10000, type="l",main="Total premium according the retention level for n=10000 (binomial)", ylab="Total Premium",xlab="Retention level")

tab10000 <- data.frame(w,Totalprem10000)

w10000_1 <- tab10000[Totalprem10000==min(Totalprem10000),1] #Selecting the w that minimize the premium value

w10000_2 <- q(bin)(0.7)

w10000_3 <- q(bin)(0.95)

purepremium10000_4 <- 1.2*c(500*0.0306,1000*0.0315 , 1500*0.024 , 2000*0.0218 , 2500*0.0204 , 
                            3000*0.0195 , 3500*0.0189 , 4000*0.018 , 4500*0.018 , 5000*0.018)

pi10000 <- sum(1.2*c(98*500*0.0306,439*1000*0.0315 , 1172*1500*0.024 , 2051*2000*0.0218 , 2461*2500*0.0204 , 
                 2051*3000*0.0195 , 1172*3500*0.0189 , 439*4000*0.018 , 98*4500*0.018 , 10*5000*0.018)) #Total premium

w10000_4 <- tab10000[which(abs(Totalprem10000-pi10000)==min(abs(Totalprem10000-pi10000))),1]#Selecting the w so that the premium is equal to the pure premium

w10000_4 <- 612500

w110000_1 <- (500*d(binm1)(w10000_1-500)*0.0306)/d(bin)(w10000_1)
w210000_1 <- (1000*d(binm2)(w10000_1-1000)*0.0315)/d(bin)(w10000_1)
w310000_1 <- (1500*d(binm3)(w10000_1-1500)*0.024)/d(bin)(w10000_1)
w410000_1 <- (2000*d(binm4)(w10000_1-2000)*0.0218)/d(bin)(w10000_1)
w510000_1 <- (2500*d(binm5)(w10000_1-2500)*0.0204)/d(bin)(w10000_1)
w610000_1 <- (3000*d(binm6)(w10000_1-3000)*0.0195)/d(bin)(w10000_1)
w710000_1 <- (3500*d(binm7)(w10000_1-3500)*0.0189)/d(bin)(w10000_1)
w810000_1 <- (4000*d(binm8)(w10000_1-4000)*0.018)/d(bin)(w10000_1)
w910000_1 <- (4500*d(binm9)(w10000_1-4500)*0.018)/d(bin)(w10000_1)
w1010000_1 <- (5000*d(binm10)(w10000_1-5000)*0.018)/d(bin)(w10000_1)

SL110000_1 <- 1.5*sum(f(x1,w110000_1)*d(bin)(k))
SL210000_1 <- 1.5*sum(f(x2,w210000_1)*d(bin)(k))
SL310000_1 <- 1.5*sum(f(x3,w310000_1)*d(bin)(k))
SL410000_1 <- 1.5*sum(f(x4,w410000_1)*d(bin)(k))
SL510000_1 <- 1.5*sum(f(x5,w510000_1)*d(bin)(k))
SL610000_1 <- 1.5*sum(f(x6,w610000_1)*d(bin)(k))
SL710000_1 <- 1.5*sum(f(x7,w710000_1)*d(bin)(k))
SL810000_1 <- 1.5*sum(f(x8,w810000_1)*d(bin)(k))
SL910000_1 <- 1.5*sum(f(x9,w910000_1)*d(bin)(k))
SL1010000_1 <- 1.5*sum(f(x10,w1010000_1)*d(bin)(k))

purepremium10000_1 <- c(w110000_1+SL110000_1,w210000_1+SL210000_1,w310000_1+SL310000_1,w410000_1+SL410000_1,
                        w510000_1+SL510000_1,w610000_1+SL610000_1,w710000_1+SL710000_1,w810000_1+SL810000_1,
                        w910000_1+SL910000_1,w1010000_1+SL1010000_1)

#Expected cash-back

CB110000_1 <- sum(f(w110000_1,x1)*d(bin)(k))
CB210000_1 <- sum(f(w210000_1,x2)*d(bin)(k))
CB310000_1 <- sum(f(w310000_1,x3)*d(bin)(k))
CB410000_1 <- sum(f(w410000_1,x4)*d(bin)(k))
CB510000_1 <- sum(f(w510000_1,x5)*d(bin)(k))
CB610000_1 <- sum(f(w610000_1,x6)*d(bin)(k))
CB710000_1 <- sum(f(w710000_1,x7)*d(bin)(k))
CB810000_1 <- sum(f(w810000_1,x8)*d(bin)(k))
CB910000_1 <- sum(f(w910000_1,x9)*d(bin)(k))
CB1010000_1 <- sum(f(w1010000_1,x10)*d(bin)(k))


w110000_2 <- (500*d(binm1)(w10000_2-500)*0.0306)/d(bin)(w10000_2)
w210000_2 <- (1000*d(binm2)(w10000_2-1000)*0.0315)/d(bin)(w10000_2)
w310000_2 <- (1500*d(binm3)(w10000_2-1500)*0.024)/d(bin)(w10000_2)
w410000_2 <- (2000*d(binm4)(w10000_2-2000)*0.0218)/d(bin)(w10000_2)
w510000_2 <- (2500*d(binm5)(w10000_2-2500)*0.0204)/d(bin)(w10000_2)
w610000_2 <- (3000*d(binm6)(w10000_2-3000)*0.0195)/d(bin)(w10000_2)
w710000_2 <- (3500*d(binm7)(w10000_2-3500)*0.0189)/d(bin)(w10000_2)
w810000_2 <- (4000*d(binm8)(w10000_2-4000)*0.018)/d(bin)(w10000_2)
w910000_2 <- (4500*d(binm9)(w10000_2-4500)*0.018)/d(bin)(w10000_2)
w1010000_2 <- (5000*d(binm10)(w10000_2-5000)*0.018)/d(bin)(w10000_2)

SL110000_2 <- 1.5*sum(f(x1,w110000_2)*d(bin)(k))
SL210000_2 <- 1.5*sum(f(x2,w210000_2)*d(bin)(k))
SL310000_2 <- 1.5*sum(f(x3,w310000_2)*d(bin)(k))
SL410000_2 <- 1.5*sum(f(x4,w410000_2)*d(bin)(k))
SL510000_2 <- 1.5*sum(f(x5,w510000_2)*d(bin)(k))
SL610000_2 <- 1.5*sum(f(x6,w610000_2)*d(bin)(k))
SL710000_2 <- 1.5*sum(f(x7,w710000_2)*d(bin)(k))
SL810000_2 <- 1.5*sum(f(x8,w810000_2)*d(bin)(k))
SL910000_2 <- 1.5*sum(f(x9,w910000_2)*d(bin)(k))
SL1010000_2 <- 1.5*sum(f(x10,w1010000_2)*d(bin)(k))

purepremium10000_2 <- c(w110000_2+SL110000_2,w210000_2+SL210000_2,w310000_2+SL310000_2,w410000_2+SL410000_2,w510000_2+SL510000_2,
                        w610000_2+SL610000_2,w710000_2+SL710000_2,w810000_2+SL810000_2,w910000_2+SL910000_2,w1010000_2+SL1010000_2)

#Expected cash-back

CB110000_2 <- sum(f(w110000_2,x1)*d(bin)(k))
CB210000_2 <- sum(f(w210000_2,x2)*d(bin)(k))
CB310000_2 <- sum(f(w310000_2,x3)*d(bin)(k))
CB410000_2 <- sum(f(w410000_2,x4)*d(bin)(k))
CB510000_2 <- sum(f(w510000_2,x5)*d(bin)(k))
CB610000_2 <- sum(f(w610000_2,x6)*d(bin)(k))
CB710000_2 <- sum(f(w710000_2,x7)*d(bin)(k))
CB810000_2 <- sum(f(w810000_2,x8)*d(bin)(k))
CB910000_2 <- sum(f(w910000_2,x9)*d(bin)(k))
CB1010000_2 <- sum(f(w1010000_2,x10)*d(bin)(k))

w110000_3 <- (500*d(binm1)(w10000_3-500)*0.0306)/d(bin)(w10000_3)
w210000_3 <- (1000*d(binm2)(w10000_3-1000)*0.0315)/d(bin)(w10000_3)
w310000_3 <- (1500*d(binm3)(w10000_3-1500)*0.024)/d(bin)(w10000_3)
w410000_3 <- (2000*d(binm4)(w10000_3-2000)*0.0218)/d(bin)(w10000_3)
w510000_3 <- (2500*d(binm5)(w10000_3-2500)*0.0204)/d(bin)(w10000_3)
w610000_3 <- (3000*d(binm6)(w10000_3-3000)*0.0195)/d(bin)(w10000_3)
w710000_3 <- (3500*d(binm7)(w10000_3-3500)*0.0189)/d(bin)(w10000_3)
w810000_3 <- (4000*d(binm8)(w10000_3-4000)*0.018)/d(bin)(w10000_3)
w910000_3 <- (4500*d(binm9)(w10000_3-4500)*0.018)/d(bin)(w10000_3)
w1010000_3 <- (5000*d(binm10)(w10000_3-5000)*0.018)/d(bin)(w10000_3)

SL110000_3 <- 1.5*sum(f(x1,w110000_3)*d(bin)(k))
SL210000_3 <- 1.5*sum(f(x2,w210000_3)*d(bin)(k))
SL310000_3 <- 1.5*sum(f(x3,w310000_3)*d(bin)(k))
SL410000_3 <- 1.5*sum(f(x4,w410000_3)*d(bin)(k))
SL510000_3 <- 1.5*sum(f(x5,w510000_3)*d(bin)(k))
SL610000_3 <- 1.5*sum(f(x6,w610000_3)*d(bin)(k))
SL710000_3 <- 1.5*sum(f(x7,w710000_3)*d(bin)(k))
SL810000_3 <- 1.5*sum(f(x8,w810000_3)*d(bin)(k))
SL910000_3 <- 1.5*sum(f(x9,w910000_3)*d(bin)(k))
SL1010000_3 <- 1.5*sum(f(x10,w1010000_3)*d(bin)(k))

purepremium10000_3 <- c(w110000_3+SL110000_3,w210000_3+SL210000_3,w310000_3+SL310000_3,w410000_3+SL410000_3,w510000_3+SL510000_3,
                        w610000_3+SL610000_3,w710000_3+SL710000_3,w810000_3+SL810000_3,w910000_3+SL910000_3,w1010000_3+SL1010000_3)


#Expected cash-back

CB110000_3 <- sum(f(w110000_3,x1)*d(bin)(k))
CB210000_3 <- sum(f(w210000_3,x2)*d(bin)(k))
CB310000_3 <- sum(f(w310000_3,x3)*d(bin)(k))
CB410000_3 <- sum(f(w410000_3,x4)*d(bin)(k))
CB510000_3 <- sum(f(w510000_3,x5)*d(bin)(k))
CB610000_3 <- sum(f(w610000_3,x6)*d(bin)(k))
CB710000_3 <- sum(f(w710000_3,x7)*d(bin)(k))
CB810000_3 <- sum(f(w810000_3,x8)*d(bin)(k))
CB910000_3 <- sum(f(w910000_3,x9)*d(bin)(k))
CB1010000_3 <- sum(f(w1010000_3,x10)*d(bin)(k))


w110000_4 <- (500*d(binm1)(w10000_4-500)*0.0306)/d(bin)(w10000_4)
w210000_4 <- (1000*d(binm2)(w10000_4-1000)*0.0315)/d(bin)(w10000_4)
w310000_4 <- (1500*d(binm3)(w10000_4-1500)*0.024)/d(bin)(w10000_4)
w410000_4 <- (2000*d(binm4)(w10000_4-2000)*0.0218)/d(bin)(w10000_4)
w510000_4 <- (2500*d(binm5)(w10000_4-2500)*0.0204)/d(bin)(w10000_4)
w610000_4 <- (3000*d(binm6)(w10000_4-3000)*0.0195)/d(bin)(w10000_4)
w710000_4 <- (3500*d(binm7)(w10000_4-3500)*0.0189)/d(bin)(w10000_4)
w810000_4 <- (4000*d(binm8)(w10000_4-4000)*0.018)/d(bin)(w10000_4)
w910000_4 <- (4500*d(binm9)(w10000_4-4500)*0.018)/d(bin)(w10000_4)
w1010000_4 <- (5000*d(binm10)(w10000_4-5000)*0.018)/d(bin)(w10000_4)

SL110000_4 <- 1.5*sum(f(x1,w110000_4)*d(bin)(k))
SL210000_4 <- 1.5*sum(f(x2,w210000_4)*d(bin)(k))
SL310000_4 <- 1.5*sum(f(x3,w310000_4)*d(bin)(k))
SL410000_4 <- 1.5*sum(f(x4,w410000_4)*d(bin)(k))
SL510000_4 <- 1.5*sum(f(x5,w510000_4)*d(bin)(k))
SL610000_4 <- 1.5*sum(f(x6,w610000_4)*d(bin)(k))
SL710000_4 <- 1.5*sum(f(x7,w710000_4)*d(bin)(k))
SL810000_4 <- 1.5*sum(f(x8,w810000_4)*d(bin)(k))
SL910000_4 <- 1.5*sum(f(x9,w910000_4)*d(bin)(k))
SL1010000_4 <- 1.5*sum(f(x10,w1010000_4)*d(bin)(k))

P2P110000_4 <- purepremium10000_4[1] - SL110000_4
P2P210000_4 <- purepremium10000_4[2] - SL210000_4
P2P310000_4 <- purepremium10000_4[3] - SL310000_4
P2P410000_4 <- purepremium10000_4[4] - SL410000_4
P2P510000_4 <- purepremium10000_4[5] - SL510000_4
P2P610000_4 <- purepremium10000_4[6] - SL610000_4
P2P710000_4 <- purepremium10000_4[7] - SL710000_4
P2P810000_4 <- purepremium10000_4[8] - SL810000_4
P2P910000_4 <- purepremium10000_4[9] - SL910000_4
P2P1010000_4 <- purepremium10000_4[10] - SL1010000_4

#Expected cash-back

CB110000_4 <- sum(f(w110000_4,x1)*d(bin)(k))
CB210000_4 <- sum(f(w210000_4,x2)*d(bin)(k))
CB310000_4 <- sum(f(w310000_4,x3)*d(bin)(k))
CB410000_4 <- sum(f(w410000_4,x4)*d(bin)(k))
CB510000_4 <- sum(f(w510000_4,x5)*d(bin)(k))
CB610000_4 <- sum(f(w610000_4,x6)*d(bin)(k))
CB710000_4 <- sum(f(w710000_4,x7)*d(bin)(k))
CB810000_4 <- sum(f(w810000_4,x8)*d(bin)(k))
CB910000_4 <- sum(f(w910000_4,x9)*d(bin)(k))
CB1010000_4 <- sum(f(w1010000_4,x10)*d(bin)(k))