#WROFIN 2016, code: Alicja Wolny-Dominiak, woali@ue.katowice.pl
library(hglm)
library(tweedie)
library(statmod)
library(StatMatch) ##dummy
set.seed(456)

triangle_upper=read.csv2(file="http://web.ue.katowice.pl/woali/triangle_upper.csv")
triangle_lower=read.csv2(file="http://web.ue.katowice.pl/woali/triangle_lower.csv")
X.upper = model.matrix(~as.factor(dev),data=triangle_upper) 
Z.upper = fact2dummy(as.factor(triangle_upper$origin))
X.lower = as.matrix(cbind(rep(1,length(triangle_lower$dev)),fact2dummy(as.factor(triangle_lower$dev)))) 
Z.lower = fact2dummy(as.factor(triangle_lower$origin))

##Data 
p=1 
claim <- triangle_upper$claim

beta_hglm=hglm(fixed=claim~as.factor(triangle_upper$dev),random=~1|as.factor(triangle_upper$origin), 
family=tweedie(var.power=p,link.power=0), rand.family=Gamma(log))

beta=beta_hglm$fixef
u=beta_hglm$ranef
v=log(beta_hglm$ranef)
phi <- beta_hglm$varFix

Y.u.lower <- exp(X.lower%*%beta+Z.lower%*%v[2:10])
Y.u.upper <- exp(X.upper%*%beta+Z.upper%*%v) #beta_hglm$fv 
#################

##Reserves
Ri_hglm <- tapply(Y.u.lower, triangle_lower$origin, sum)
R_hglm <- sum(Ri_hglm)
R_hglm

#Residual Bootstrap HGLM RMSEP
nsim=500
Ri_hglmB <- Ri_hglm_pB <- R_hglmB <- R_hglm_pB <- NULL
resid=(claim-Y.u.upper)/sqrt(Y.u.upper^p) 

for (i in 1:nsim) {
residB <- sample(resid, nrow(triangle_upper), replace = TRUE)
Y.u.upperB <- abs(residB * sqrt(Y.u.upper^p) + Y.u.upper) 
			
beta_hglmB <- hglm(fixed=Y.u.upperB~as.factor(dev),random=~1|as.factor(origin), family=tweedie(var.power=p,link.power=0), 
rand.family=Gamma(log), data=triangle_upper)

betaB=beta_hglmB$fixef
uB=beta_hglmB$ranef
vB=log(uB)
phiB <- beta_hglmB$varFix
Y.u.lowerB <- exp(X.lower%*%betaB+Z.lower%*%vB[2:10])

uBB = rtweedie(length(u), mu = c(u), phi =  beta_hglmB$varRanef, power = 2)
vBB = log(uBB)
Y.u.lowerBB <- rtweedie(length(Y.u.lowerB), mu = c(exp(X.lower%*%betaB+Z.lower%*%vBB[2:10])), phi = phiB, power = p)

Ri_hglmB <- rbind(Ri_hglmB, tapply(Y.u.lowerB, triangle_lower$origin, sum))
Ri_hglm_pB <- rbind(Ri_hglm_pB, tapply(Y.u.lowerBB, triangle_lower$origin, sum))
R_hglmB=rbind(R_hglmB, sum(Ri_hglmB[i,]))
R_hglm_pB=rbind(R_hglm_pB, sum(Ri_hglm_pB[i,]))
}

RMSEPi_hglm <-  NULL
for (k in 1:9){
RMSEPi_hglm <-rbind(RMSEPi_hglm, sqrt(sum((Ri_hglm_pB[,k] - Ri_hglmB[,k])^2)/nsim))
}

RMSEP_hglm <- sqrt(sum((R_hglm_pB - R_hglmB)^2)/nsim)

#### Parametric bootstrap HGLM RMSE
Ri_hglm1B <- Ri_hglm1BB <- R_hglm1B <- R_hglm1BB <- NULL

for(i in 1:nsim){
u1B = rtweedie(length(u), mu = c(u), phi =  beta_hglm$varRanef, power = 2)
v1B=log(u1B)
Y.u.upper1B = rtweedie(nrow(triangle_upper), mu = c(exp(X.upper%*%beta+Z.upper%*%v1B)), phi = phi, power = p)
Y.u.lower1B = rtweedie(nrow(triangle_lower), mu = c(exp(X.lower%*%beta+Z.lower%*%v1B[2:10])), phi = phi, power = p)

beta_hglm1B <- hglm(fixed=Y.u.upper1B~as.factor(dev),random=~1|as.factor(origin), family=tweedie(var.power=p,link.power=0), 
rand.family=Gamma(log), data=triangle_upper, maxit=10)
 
beta1BB=beta_hglm1B$fixef
u1BB=beta_hglm1B$ranef
v1BB=log(u1BB)
   
Y.u.lower1BB <-exp(X.lower%*%beta1BB+Z.lower%*%v1BB[2:10])

Ri_hglm1B <- rbind(Ri_hglm1B, tapply(Y.u.lower1B, triangle_lower$origin, sum))
Ri_hglm1BB <- rbind(Ri_hglm1BB, tapply(Y.u.lower1BB, triangle_lower$origin, sum))
R_hglm1B=rbind(R_hglm1B, sum(Ri_hglm1B[i,]))
R_hglm1BB=rbind(R_hglm1BB, sum(Ri_hglm1BB[i,]))
}
	
RMSEP_hglm1 = sqrt(sum((R_hglm1B - R_hglm1BB)^2)/nsim)

RMSEPi_hglm1 <- NULL
for (k in 1:9){
RMSEPi_hglm1 <-rbind(RMSEPi_hglm1, sqrt(sum((Ri_hglm1B[,k] - Ri_hglm1BB[,k])^2)/nsim))
}

summary <-cbind(rbind(as.matrix(Ri_hglm), R_hglm), rbind(RMSEPi_hglm, RMSEP_hglm), rbind(RMSEPi_hglm1, RMSEP_hglm1))
colnames(summary) <-c("R", "residual boot RMSEP", "parametric boot RMSEP")