## ----eval10-------------------------------------------------------------------
d <- readRDS("pred_mort1.rds")
d <- na.omit(d)


## ----eval11-------------------------------------------------------------------
trigger <- 0.20
d$loss <- pmax(0, d$mortality_rate - trigger)
d$pay  <- pmax(0, d$predicted_mortality - trigger)
head(d, n=2)


## ----eval12-------------------------------------------------------------------
plot(d$loss, d$pay, xlab="Loss", ylab="Pay")
abline(0,1, col="red")


## ----ev1----------------------------------------------------------------------
plot(pay~NDVI, data=d, col="red")


## ----util---------------------------------------------------------------------
library(agro)
agro::ce_income


## ----eval20-------------------------------------------------------------------
markup <- 0.25
premium <- mean(d$pay) * (1 + markup)
premium

base <- 1-d$mortality_rate
insurance_nomarkup <- base + d$pay
insurance <- insurance_nomarkup - premium

base <- base * 100
insurance <- insurance * 100
insurance_nomarkup <- insurance_nomarkup * 100


## ----eval21-------------------------------------------------------------------
rho <- 2
ce_base <- ce_income(base, rho)
ce_ins25  <- ce_income(insurance, rho)
ce_ins_nomarkup  <- ce_income(na.omit(insurance_nomarkup), rho)

ce_base 
ce_ins25
ce_ins_nomarkup

ce_ins25 / ce_base
ce_ins_nomarkup/ce_base
head(d, n=2)



## ----eval22-------------------------------------------------------------------
rhos <- seq(0, 3, .1)
ce_base <- ce_ins <- ce_ins_nomarkup <- rep(NA, length(rhos))
for(i in 1:length(rhos)){
  ce_base[i] <- ce_income(base, rhos[i])
  ce_ins[i]  <- ce_income(insurance, rhos[i])
  ce_ins_nomarkup[i] <- ce_income(insurance_nomarkup, rhos[i])
}

inc <- seq(0.1, 1, 0.1)
plot(rhos, ce_base, type="l", col="blue", ylab= "Pay (%)", xlab="CRRA", cex=2, ylim=c(70, 90))

lines(rhos, ce_ins, col="red")
lines(rhos, ce_ins_nomarkup, col="green")

legend("topright", c("No insurance", "25% Marked-up Insurance", "No mark up"), lty=1, col=c("blue", "red", "green"), title = "IBLI contract type", bty = "n")



## ----eval23-------------------------------------------------------------------
ce_base <- ce_ins <- mqs <- rep(NA, length(rhos))

for(i in 1:length(rhos)){
  ce_base[i] <- ce_income(base, rhos[i])
  ce_ins[i]  <- ce_income(insurance, rhos[i])
  mqs[i] <- ce_ins[i] - ce_base[i]
}

inc <- seq(0.1, 1, 0.1)
plot(rhos, mqs*100, type="l", col="blue", ylab= "MQS (insurance benefit in %)", xlab="CRRA", cex=2)
abline(h=0, col="red")



## ----eval100, message=FALSE---------------------------------------------------
library(msir)
dd <- na.omit(d)



## ----eval200------------------------------------------------------------------
m <- loess.sd(dd$NDVI, dd$mortality_rate)
fitsd <- cbind(fit=m$model$fitted, sd=m$sd)



## ----eval3--------------------------------------------------------------------
ns <- 1000
sample <- apply(fitsd, 1, function(i) {
  pmin(1, pmax(0, rnorm(ns, i[1], i[2])))
})
sample <- t(sample)


## ----eval4--------------------------------------------------------------------
markup <- 0.2
premium <- mean(dd$pay) * (1 + markup)


## ----eval6--------------------------------------------------------------------
out_nomarkup <- out_ins <- out_base <- rep(NA, ns)
rho = 2

for (i in 1:ns) {
  base <- 1- sample[,i]
  insurance_nomarkup <- base + dd$pay
  insurance <- insurance_nomarkup - premium
  out_base[i] <- ce_income(100*base, rho)
  out_ins[i]  <- ce_income(100*insurance, rho)
  out_nomarkup[i]  <- ce_income(100*insurance_nomarkup, rho)
}

mean(out_base)
mean(out_ins)
mean(out_nomarkup)



## ----eval5--------------------------------------------------------------------
benefit_ins <- out_ins - out_base
benefit_nomu <- out_nomarkup - out_base
b <- cbind(benefit_ins, benefit_nomu)
boxplot(b, ylim=c(-1,10))