#Compare ChainLadder and GAMLSS standard model
#Perform the analysis on PPA Data

groups2Scan<-unique(privatePassengersAutoRawData$GRCODE)
#function that returns BE for ChainLadder and GAMLSS standard model

doBEs<-function(group) {
	#Prepare data
	dataGr<-subset(privatePassengersAutoRawData, GRCODE==group, select=c("AccidentYear","DevelopmentYear","DevelopmentLag","CumPaidLoss_B"))
	##fix negative/equal incremental payments
	for(i in 2:nrow(dataGr)) {
		if((dataGr$CumPaidLoss_B[i]<=dataGr$CumPaidLoss_B[i-1])&(dataGr$AccidentYear[i]==dataGr$AccidentYear[i-1])) dataGr$CumPaidLoss_B[i]=dataGr$CumPaidLoss_B[i-1]+1
	}

	#get reserve and ultimate payment
	ultimatePayments<-with(subset(dataGr, DevelopmentLag==10), sum(CumPaidLoss_B))
	lastPayments<-with(subset(dataGr, AccidentYear+DevelopmentLag==1998), sum(CumPaidLoss_B))
	reserve<-ultimatePayments-lastPayments #(true reserve)

	#prepare data
	##collect
	df4Triangles<-subset(dataGr, AccidentYear+DevelopmentLag<=1998)
	names(df4Triangles)=c("ay","dy","dev","cumpaid")
	triGr<-as.triangle(df4Triangles,origin="ay",dev="dev",value="cumpaid")

	#make the triangles
	require(ChainLadder)
	incrementalTriangle<-as.data.frame(cum2incr(triGr))
	incrementalTriangle$ay<-factor(incrementalTriangle$ay,ordered=TRUE)
	incrementalTriangle$dev<-factor(incrementalTriangle$dev,ordered=TRUE)
	triangleDf<-subset(incrementalTriangle, !is.na(value))
	predictionDf<-subset(incrementalTriangle, is.na(value)) #take only the unfilled
	predictionDf<-predictionDf[,c("ay","dev")];rownames(predictionDf)<-NULL

	#estimate the models
	##chain ladder
	chainladderReserve<-summary(MackChainLadder(triGr, alpha=2))$Totals[4,1] #ALPHA = 2 Murphy BLUE
	#gamlss
	gamlssModel<-gamlss(formula=value~ay+dev, sigma.formula = ~ay, data=triangleDf, family="GA", control=con)
	gamlssReserve<-sum(predict(gamlssModel, newdata=predictionDf,type="response"))
	out<-data.frame(group=group,real=reserve, CL=chainladderReserve, GAMLSS=gamlssReserve)
	invisible(out)
}

#perform the BE for all PAP triangles

if(exists("resultsAll")) rm(resultsAll)
for(i in groups2Scan) 
	{ 
		ciao<-tryCatch(doBEs(i),error = function(e){return( out=data.frame(group=i,real=NA, CL=NA, GAMLSS=NA))})
		if(exists("resultsAll")) resultsAll<-rbind(resultsAll,ciao) else resultsAll<-ciao
}

#check the RMSE
resultsAll<-subset(resultsAll, is.finite(real))
cor(resultsAll[,2:4],use="complete")
library(caret)
with(resultsAll, RMSE(pred=CL, obs=real))
with(resultsAll, RMSE(pred=CL, obs=GAMLSS))