Assignment

##Install necessary packages
install.packages("markovchain")
##Set up packages and data
library(markovchain)
library(readr)
data<-read.csv("claimsdata.csv")

============================================================================================================================================

##Task 1
##Consider that the MLE of one step transition probability is nothing more than the number of policy holders who moved from level i to j over the total number of policy holders at level i during 2018.
##Number of policy holders in state i in 2018
N2<-length(which(data$ncdlevel18==-2))
N1<-length(which(data$ncdlevel18==-1))
N0<-length(which(data$ncdlevel18==0))
n1<-length(which(data$ncdlevel18==1))
n2<-length(which(data$ncdlevel18==2))

##Number of policy holdies transitioning from state i to j
##Note that since states can only shift to the left by 1 or 2, or 1 to the right means that there can only be 2 or 3 transitions from a specified state
N2_N2<-length(which(data$ncdlevel18==-2&data$ncdlevel19==-2))
N2_N1<-length(which(data$ncdlevel18==-2&data$ncdlevel19==-1))

N1_N2<-length(which(data$ncdlevel18==-1&data$ncdlevel19==-2))
N1_N0<-length(which(data$ncdlevel18==-1&data$ncdlevel19==0))

N0_N2<-length(which(data$ncdlevel18==0&data$ncdlevel19==-2))
N0_N1<-length(which(data$ncdlevel18==0&data$ncdlevel19==-1))
N0_n1<-length(which(data$ncdlevel18==0&data$ncdlevel19==1))

n1_N1<-length(which(data$ncdlevel18==1&data$ncdlevel19==-1))
n1_N0<-length(which(data$ncdlevel18==1&data$ncdlevel19==0))
n1_n2<-length(which(data$ncdlevel18==1&data$ncdlevel19==2))

n2_N0<-length(which(data$ncdlevel18==2&data$ncdlevel19==0))
n2_n1<-length(which(data$ncdlevel18==2&data$ncdlevel19==1))
n2_n2<-length(which(data$ncdlevel18==2&data$ncdlevel19==2))

##Solving for probability
pN2_N2<-N2_N2/N2
pN2_N1<-N2_N1/N2

pN1_N2<-N1_N2/N1
pN1_N0<-N1_N0/N1

pN0_N2<-N0_N2/N0
pN0_N1<-N0_N1/N0
pN0_n1<-N0_n1/N0

pn1_N1<-n1_N1/n1
pn1_N0<-n1_N0/n1
pn1_n2<-n1_n2/n1

pn2_N0<-n2_N0/n2
pn2_n1<-n2_n1/n2
pn2_n2<-n2_n2/n2

##Transitional Matrix
DiscountLevel <- new("markovchain", states = c("-2", "-1", "0", "1", "2"), 
transitionMatrix = matrix(data = c(pN2_N2, pN2_N1, 0, 0, 0, 
pN1_N2, 0, pN1_N0, 0, 0, 
pN0_N2, pN0_N1, 0, pN0_n1, 0, 
0, pn1_N1, pn1_N0, 0, pn1_n2, 
0, 0, pn2_N0, pn2_n1, pn2_n2), 
byrow = TRUE, nrow = 5), name = "Discount Level")

============================================================================================================================================

##Task 2 
##Set up empty vectors
ncdlevel20<-c()
N_2<-c()
N_1<-c()
N_0<-c()
n_1<-c()
n_2<-c()
##Create loop to simulate 2020 states 1000 times
for(n in 1:1000){
    for(t in 1:10000){
        ncdlevel20[t]=rmarkovchain(1, object=DiscountLevel,t0=data$ncdlevel19[t])[1]
    }
    ncdlevel20<-as.integer(ncdlevel20)
    N_2[n]<-length(which(ncdlevel20==-2))
    N_1[n]<-length(which(ncdlevel20==-1))
    N_0[n]<-length(which(ncdlevel20==0))
    n_1[n]<-length(which(ncdlevel20==1))
    n_2[n]<-length(which(ncdlevel20==2))
}
##Create a string of values with Total Expected Premium for each of the 1000 simulations
Simulations<-300*(N_2*1.2+N_1*1.1+N_0+n_1*0.9+n_2*0.8)

##Total Expected Premium
mean(Simulations)
>[1] 2697954
##Solve for Expected Value of a single premium
mean(Simulations)/10000
>[1] 269.7954
##General Statistics
summary(Simulations)
##SD
sd(Simulations)


===========================================================================================================================================

##Task 3
##Total premiums in 2018
Premium18 <- 300*(N2*1.2+N1*1.1+N0+n1*0.9+n2*0.8)
>[1] 3004680
##Amount of claims
sum(data$numberofclaims)
>[1] 1452
##Profit
300*(N2*1.2+N1*1.1+N0+n1*0.9+n2*0.8) - (sum(data$numberofclaims))*2000
>[1] 100680

##Total Premiums in 2019
N__2<-length(which(data$ncdlevel19==-2))
N__1<-length(which(data$ncdlevel19==-1))
N__0<-length(which(data$ncdlevel19==0))
n__1<-length(which(data$ncdlevel19==1))
n__2<-length(which(data$ncdlevel19==2))
Premium19 <- 300*(N__2*1.2+N__1*1.1+N__0+n__1*0.9+n__2*0.8)
>[1] 2829870

##Total Premiums in 2021
ncdlevel21<-c()
N_21<-c()
N_11<-c()
N_01<-c()
n_11<-c()
n_21<-c()
##Create loop to simulate 2020 states 1000 times
for(n in 1:1000){
    for(t in 1:10000){
        ncdlevel21[t]=rmarkovchain(1, object=DiscountLevel,t0=ncdlevel20[t])
    }
    ncdlevel21<-as.integer(ncdlevel21)
    N_21[n]<-length(which(ncdlevel21==-2))
    N_11[n]<-length(which(ncdlevel21==-1))
    N_01[n]<-length(which(ncdlevel21==0))
    n_11[n]<-length(which(ncdlevel21==1))
    n_21[n]<-length(which(ncdlevel21==2))
}
##Create a string of values with Total Expected Premium for each of the 1000 simulations
Simulations21<-300*(N_21*1.2+N_11*1.1+N_01+n_11*0.9+n_21*0.8)
mean(Simulations)
>[1] 2600021

##Total Premiums in 2022
ncdlevel22<-c()
N_22<-c()
N_12<-c()
N_02<-c()
n_12<-c()
n_22<-c()
##Create loop to simulate 2020 states 1000 times
for(n in 1:1000){
    for(t in 1:10000){
        ncdlevel22[t]=rmarkovchain(1, object=DiscountLevel,t0=ncdlevel21[t])
    }
    ncdlevel22<-as.integer(ncdlevel22)
    N_22[n]<-length(which(ncdlevel22==-2))
    N_12[n]<-length(which(ncdlevel22==-1))
    N_02[n]<-length(which(ncdlevel22==0))
    n_12[n]<-length(which(ncdlevel22==1))
    n_22[n]<-length(which(ncdlevel22==2))
}
##Create a string of values with Total Expected Premium for each of the 1000 simulations
Simulations22<-300*(N_22*1.2+N_12*1.1+N_02+n_12*0.9+n_22*0.8)

##Plotting Premiums
Premiums <- c(Premium18, Premium19, mean(Simulations), mean(Simulations21), mean(Simulations22))
EPremium<-data.frame(c(2018, 2019, 2020, 2021, 2022), Premiums)
names(EPremium) <- c("Year", "Expected Premium")
plot(EPremium,type = "o",col = "red", main = "Total Expected Premium")

##Producing Premiums
Profit<-data.frame(c(2018, 2019, 2020, 2021, 2022), Premiums, c(2000*sum(data$numberofclaims), 2000*sum(data$numberofclaims), 2000*sum(data$numberofclaims), 2000*sum(data$numberofclaims), 2000*sum(data$numberofclaims)))
names(Profit) <- c("Year", "Expected Premium", "Expenses")
Profit$Profit <- (Profit$`Expected Premium` - Profit$Expenses)
Profit<-floor(Profit)

#Plotting States
States<-data.frame(c(-2, -1, 0, 1, 2), c(N2, N1, N0, n1, n2),
c(length(which(data$ncdlevel19==-2)), length(which(data$ncdlevel19==-1)), length(which(data$ncdlevel19==0)), length(which(data$ncdlevel19==1)), length(which(data$ncdlevel19==2))),
c(mean(N_2), mean(N_1), mean(N_0), mean(n_1), mean(n_2)),
c(mean(N_21), mean(N_11), mean(N_01), mean(n_11), mean(n_21)),
c(mean(N_22), mean(N_12), mean(N_02), mean(n_12), mean(n_22)))
names(States) <- c("State", "2018", "2019", "2020", "2021", "2022")
States<-t(States)
States<-floor(States)
States<-States[-c(1), ]
States<-as.data.frame(States)
names(States) <- c("-2", "-1", "0", "1", "2")
State123<-c(-2, -1, 0, 1, 2)
matplot(States, type="b", pch=15:19, xlab="Year",ylab="No. of Policyholders",main="No. of Policyholders in a state")
legend("bottomleft", inset=0.005, legend=State123, col=c(1:5),ncol=2, pch=15:19, bg= ("white"), horiz=F, cex = 0.42)

steadyStates(DiscountLevel)