CRPS_calculation_script.R

#####Debugging###
#1. add: in_sample<-in_sample[order(as.numeric(as.character(in_sample$origin))),]
#2. add: out_sample<-out_sample[order(as.numeric(as.character(out_sample$origin))),]
#3. change gamlss to ###GAMLSS
#gamlss2_GA_In<-gamlss(formula=(value+tau_Ga)~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,sigma.formula=~cs(as.numeric(as.character(dev))),family=GA(mu.link="log", sigma.link ="log"),trace=FALSE)
#gamlss2_LN_In<-gamlss(formula=(value+tau_LN)~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,sigma.formula=~cs(as.numeric(as.character(dev))),family=LOGNO(mu.link="identity",sigma.link="log"),trace=FALSE)
#4. Change: 
#in_sample_numeric<-in_sample
#in_sample_numeric$origin=as.numeric(as.character(in_sample$origin))
#in_sample_numeric$dev=as.numeric(as.character(in_sample$dev))
#in_sample_numeric$Calendar=as.numeric(as.character(in_sample$Calendar))

#5. Change: out_sample_numeric<-out_sample
#out_sample_numeric$origin=as.numeric(as.character(out_sample$origin))
#out_sample_numeric$dev=as.numeric(as.character(out_sample$dev))
#out_sample_numeric$Calendar=as.numeric(as.character(out_sample$Calendar))


ntri<-100
SpLN_crps<-matrix(NA,nrow=780,ncol=ntri)
EqEns_crps<-matrix(NA,nrow=780,ncol=ntri)

SLP_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP8_crps<-matrix(NA,nrow=780,ncol=ntri)

ADLP1_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP2_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP3_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP4_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP5_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP6_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP7_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP9_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP10_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP11_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP12_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP13_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP14_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP15_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP16_crps<-matrix(NA,nrow=780,ncol=ntri)
ADLP17_crps<-matrix(NA,nrow=780,ncol=ntri)

for (D in 1:ntri){
  
  set.seed(20200130+D)
  # Package-wise global parameters
  set_parameters(ref_claim = 200000, time_unit = 1/4)
  ref_claim = return_parameters()[1]
  time_unit = return_parameters()[2]
  
  # Claim occurrence
  # Example 1.1: Claim occurrence: Constant exposure and frequency
  
  # Input parameters
  years = 10
  I = years / time_unit
  E = c(rep(12000, I)) # effective annual exposure rates
  lambda = c(rep(0.03, I))
  
  # Number of claims occurring for each period i
  n_vector = claim_frequency(I = I, E = E, freq = lambda)
  
  # Occurrence time of each claim r, for each period i
  occurrence_times = claim_occurrence(frequency_vector = n_vector)
  
  ##############Claim size: Default power normal
  # S^0.2 ~ N(9.5, 3), left truncated at 30
  
  S_df <- function(s) {
    # truncate and rescale
    if (s < 1) {
      return(0)
    } else {
      p_trun <- pnorm(s^0.2, 9.5, 3) - pnorm(1^0.2, 9.5, 3)
      p_rescaled <- p_trun/(1 - pnorm(1^0.2, 9.5, 3))
      return(p_rescaled)
    }
  }
  
  claim_sizes <- claim_size(frequency_vector = n_vector, simfun = S_df, type = "p", range = c(0, 1e24))
  
  ###Notification Delay
  notidel_param <- function(claim_size, occurrence_period){
    mean <- 2.47; cv <- 1.54
    shape <- get_Weibull_parameters(mean, cv)[1, ]
    scale <- get_Weibull_parameters(mean, cv)[2, ]
    c(shape = shape, scale = scale)
  }
  
  ## output
  notidel <- claim_notification(n_vector, claim_size(n_vector),
                                paramfun = notidel_param)
  #By changing claim notification to transformed Gamma, the first development periods contain no zero value claims 
  
  
  # Claim closure: Default Weibull
  # Claim closure: Default Weibull
  setldel_param <- function(claim_size, occurrence_period) {
    target_mean <- 11.74
    
    # specify the target Weibull coefficient of variation
    target_cv <- 0.61
    
    c(shape = get_Weibull_parameters(target_mean, target_cv)[1, ],
      scale = get_Weibull_parameters(target_mean, target_cv)[2, ])
  }
  
  ## output
  # simulate the settlement delays from the Weibull with parameters above
  setldel <- claim_closure(n_vector, claim_sizes, rfun = rweibull, paramfun = setldel_param)
  
  # Claim partial payment: Default mixture distribution
  
  no_payments <- claim_payment_no(n_vector, claim_sizes)
  
  # Claim Partial payment (Without Inflation): Default Distribution
  
  payment_sizes <- claim_payment_size(n_vector, claim_sizes,no_payments)
  
  # Claim payment time: Default Distribution
  
  ## input
  r_pmtdel <- function(n, claim_size, setldel, setldel_mean) {
    result <- c(rep(NA, n))
    
    # First simulate the unnormalised values of d, sampled from a Weibull distribution
    if (n >= 4) {
      # 1) Simulate the last payment delay
      unnorm_d_mean <- (1 / 4) / time_unit
      unnorm_d_cv <- 0.20
      parameters <- get_Weibull_parameters(target_mean = unnorm_d_mean, target_cv = unnorm_d_cv)
      result[n] <- stats::rweibull(1, shape = parameters[1], scale = parameters[2])
      
      # 2) Simulate all the other payment delays
      for (i in 1:(n - 1)) {
        unnorm_d_mean <- setldel_mean / n
        unnorm_d_cv <- 0.35
        parameters <- get_Weibull_parameters(target_mean = unnorm_d_mean, target_cv = unnorm_d_cv)
        result[i] <- stats::rweibull(1, shape = parameters[1], scale = parameters[2])
      }
      
    } else {
      for (i in 1:n) {
        unnorm_d_mean <- setldel_mean / n
        unnorm_d_cv <- 0.35
        parameters <- get_Weibull_parameters(target_mean = unnorm_d_mean, target_cv = unnorm_d_cv)
        result[i] <- stats::rweibull(1, shape = parameters[1], scale = parameters[2])
      }
    }
    
    # Normalise d such that sum(inter-partial delays) = settlement delay
    # To make sure that the pmtdels add up exactly to setldel, we treat the last one separately
    result[1:n-1] <- (setldel/sum(result)) * result[1:n-1]
    result[n] <- setldel - sum(result[1:n-1])
    
    return(result)
  }
  
  param_pmtdel <- function(claim_size, setldel, occurrence_period) {
    # mean settlement delay
    if (claim_size < (0.10 * ref_claim) & occurrence_period >= 21) {
      a <- min(0.85, 0.65 + 0.02 * (occurrence_period - 21))
    } else {
      a <- max(0.85, 1 - 0.0075 * occurrence_period)
    }
    mean_quarter <- a * min(25, max(1, 6 + 4*log(claim_size/(0.10 * ref_claim))))
    target_mean <- mean_quarter / 4 / time_unit
    
    c(claim_size = claim_size,
      setldel = setldel,
      setldel_mean = target_mean)
  }
  
  ## output
  payment_delays <- claim_payment_delay(
    n_vector, claim_sizes, no_payments, setldel,
    rfun = r_pmtdel, paramfun = param_pmtdel,
    occurrence_period = rep(1:I, times = n_vector))
  
  # payment times on a continuous time scale
  payment_times <- claim_payment_time(n_vector, occurrence_times, notidel, payment_delays)
  # payment times in periods
  payment_periods <- claim_payment_time(n_vector, occurrence_times, notidel, payment_delays,discrete = TRUE)
  
  
  # Claim inflation
  # Base inflation: a vector of quarterly rates
  # In this data set we set base inflation to be at 2% p.a. constant for both past and future
  
  demo_rate <- (1 + 0.02)^(1/4) - 1
  base_inflation_past <- rep(demo_rate, times = 40)
  base_inflation_future <- rep(demo_rate, times = 40)
  base_inflation_vector <- c(base_inflation_past, base_inflation_future)
  
  # Remove the Superimposed inflation
  
  SI_occurrence <- function(occurrence_time, claim_size){1}
  # 2) With respect to payment "time" (continuous scale)
  # -> compounding by user-defined time unit
  SI_payment <-  function(payment_time, claim_size){1}
  
  payment_inflated <- claim_payment_inflation(
    n_vector,
    payment_sizes,
    payment_times,
    occurrence_times,
    claim_sizes,
    base_inflation_vector,
    SI_occurrence,
    SI_payment
  )
  agg_sqrt_withTail<-claim_output(n_vector, payment_times, payment_inflated, incremental = TRUE,future=TRUE,adjust=FALSE)
  agg_sqrt<-agg_sqrt_withTail[,-ncol(agg_sqrt_withTail)]
  colnames(agg_sqrt)<-c(1:40)
  rownames(agg_sqrt)<-c(1:40)
  #Generate a Claim Triangle
  agg_tri_withInflation <- claim_output(n_vector, payment_times, payment_inflated,incremental = TRUE)
  colnames(agg_tri_withInflation)<-c(1:40)
  rownames(agg_tri_withInflation)<-c(1:40)
  full_dat<-as.data.frame(as.triangle(agg_sqrt))
  full_dat$value=full_dat$value/10000
  full_dat$Calendar=as.numeric(full_dat$origin)+as.numeric(full_dat$dev)
  full_dat$origin=as.factor(full_dat$origin)
  full_dat$dev=as.factor(full_dat$dev)
  in_sample<-full_dat[full_dat$Calendar<=41,]
  in_sample<-in_sample[order(as.numeric(as.character(in_sample$origin))),]
  out_sample<-full_dat[full_dat$Calendar>41,]
  out_sample<-out_sample[order(as.numeric(as.character(out_sample$origin))),]
  
  z_u<-2*round(max(out_sample$value),0)
  z_l<-0
  z<-z_l:z_u
  N<-length(z_u:z_l)
  I<-function(y,z) ifelse(y<=z,1,0)
  
  
  tri_withInflation<-as.triangle(agg_tri_withInflation)
  
  # convert the aggregate data generated by SynthETIC into triangle format: value, development year and accident year
  dt_withInflation<-as.data.frame(tri_withInflation)
  # Convert the triangle into dataframe format
  dt_withInflation_past<-as.data.frame(tri_withInflation,na.rm=TRUE)
  # only contains the past observations (the upper triangle) as the lower triangle are coded as N.A
  dt_withInflation_past$Calendar=as.character(as.numeric(dt_withInflation_past$origin)+as.numeric(dt_withInflation_past$dev))
  #add calendar periods in the data frame
  
  #Convert the square into data frame
  dt_sqrt<-as.data.frame(as.triangle(agg_sqrt))
  
  #Rescale the Claims Value in the unit of 0000's
  dt_withInflation_past$value=dt_withInflation_past$value/10000
  
  # Construction of Training Set:
  train_1<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==1,]
  train_2<-dt_withInflation_past[as.numeric(dt_withInflation_past$Calendar)<=34&as.numeric(dt_withInflation_past$origin)<=32&as.numeric(dt_withInflation_past$origin)>=2,]
  train_3<-dt_withInflation_past[dt_withInflation_past$dev==1&as.numeric(dt_withInflation_past$origin)>=33&as.numeric(dt_withInflation_past$origin)<=40,]
  train_4<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)>=33&as.numeric(dt_withInflation_past$origin)<=34&as.numeric(dt_withInflation_past$dev)==2,]
  train_5<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==40&as.numeric(dt_withInflation_past$dev)==1,]
  train<-rbind(train_1,train_2,train_3,train_4,train_5)
  # Construction of Validation set:
  valid_1<-dt_withInflation_past[as.numeric(dt_withInflation_past$Calendar)>=35&as.numeric(dt_withInflation_past$Calendar)<=38&as.numeric(dt_withInflation_past$origin)>=2&as.numeric(dt_withInflation_past$origin)<=32,]
  valid_2<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==33&as.numeric(dt_withInflation_past$dev)>=3&as.numeric(dt_withInflation_past$dev)<=5,]
  valid_3<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==34&as.numeric(dt_withInflation_past$dev)>=3&as.numeric(dt_withInflation_past$dev)<=4,]
  valid_4<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)>=35&as.numeric(dt_withInflation_past$origin)<=37&as.numeric(dt_withInflation_past$dev)>=2&as.numeric(dt_withInflation_past$dev)<=3,]
  valid_5<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)>=38&as.numeric(dt_withInflation_past$origin)<=39&as.numeric(dt_withInflation_past$dev)==2,]
  valid_6<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==40&as.numeric(dt_withInflation_past$dev)==1,]
  valid<-rbind(valid_1,valid_2,valid_3,valid_4,valid_5,valid_6)
  # Construction of Test Data
  test_1<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)>=2&as.numeric(dt_withInflation_past$origin)<=35&as.numeric(dt_withInflation_past$Calendar)>=39&as.numeric(dt_withInflation_past$Calendar)<=41,]
  test_2<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==36&as.numeric(dt_withInflation_past$dev)>=4&as.numeric(dt_withInflation_past$dev)<=5,]
  test_3<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==37&as.numeric(dt_withInflation_past$dev)==4,]
  test_4<-dt_withInflation_past[as.numeric(dt_withInflation_past$origin)==38&as.numeric(dt_withInflation_past$dev)==3,]
  test<-rbind(test_1,test_2,test_3,test_4)
  test$origin=as.factor(test$origin)
  test$dev=as.factor(test$dev)
  test<-test[order(as.numeric(as.character(test$origin))),]
  train$origin=as.factor(train$origin)
  train$dev=as.factor(train$dev)
  train$Calendar=as.numeric(train$Calendar)
  #######Generate Predictive Density for Validation Sets
  test<-test[order(as.numeric(as.character(test$origin))),]
  aug_valid<-rbind(valid,test)
  aug_valid<-aug_valid[order(as.numeric(as.character(aug_valid$origin))),]
  aug_valid$origin<-as.factor(aug_valid$origin)
  aug_valid$dev<-as.factor(aug_valid$dev)
  aug_valid$Calendar<-as.numeric(aug_valid$Calendar)
  aug_valid_numeric<-aug_valid
  aug_valid_numeric$origin=as.numeric(as.character(aug_valid$origin))
  aug_valid_numeric$dev=as.numeric(as.character(aug_valid$dev))
  aug_valid_numeric$Calendar=as.numeric(as.character(aug_valid$Calendar))
  ###Fitting of all the models in the training Sets
  ####Basic Models 
  #ODP
  ODP_GLM_train<-glm(formula=value~factor(origin)+factor(dev),family=quasipoisson(link="log"),data=train)
  #ZAGA
  gamma_1_train<-gamlss(formula=value~factor(origin)+factor(dev),nu.formula=~as.numeric(as.character(dev)),data=train,family=ZAGA(mu.link="log",sigma.link = "log", nu.link = "logit"))
  #ZALN
  LN_1_train<-gamlssZadj(y=value,mu.formula = ~factor(origin)+factor(dev),xi0.formula=~as.numeric(as.character(dev)),data=train,family=LOGNO(mu.link="identity",sigma.link="log"))
  tau_LN<-5
  tau_Ga<-5
  #Gamma 
  Ga_optimTau_train<-gamlss(formula=(value+tau_Ga)~factor(origin)+factor(dev),data=train,family=GA(mu.link="log", sigma.link ="log"))
  #Log-Normal
  LN_optimTau_train<-gamlss(formula=(value+tau_LN)~factor(origin)+factor(dev),data=train,family=LOGNO(mu.link="identity",sigma.link="log"))
  dens_ODPGLM2<-cal_dens_ODP(ODP_GLM_train,aug_valid)
  dens_GAGLM2<-cal_dens_GA(Ga_optimTau_train,tau_Ga,aug_valid)
  dens_LNGLM2<-cal_dens_LN(LN_optimTau_train,tau_LN,aug_valid)
  dens_ZAGA2<-cal_dens_ZAGA(gamma_1_train,aug_valid)
  dens_ZALN2<-cal_dens_ZALN(LN_1_train,aug_valid)
  
  train_new<-train
  train_new$origin=as.numeric(as.character(train$origin))
  train_new$dev=as.numeric(as.character(train$dev))
  
  ####With Hoerl Curve
  #Under ODP Assumption
  train_Ho<-train
  train_Ho$dev=as.numeric(as.character(train$dev))
  glm_ODP_Ho_tr1<-glm(formula=value~factor(origin)+log(dev)+dev,family=quasipoisson(link="log"),data=train_Ho)
  #R retun warning for rank deficient fit for Hoerl Curve ODP
  #Under Gamma Assumption
  glm_Ga_Ho_tr1<-gamlss(formula=(value+tau_Ga)~factor(origin)+log(dev)+dev,data=train_Ho,family=GA(mu.link="log", sigma.link ="log"))
  #Under Log-Normal Assumption
  glm_LN_Ho_tr1<-gamlss(formula=(value+tau_LN)~factor(origin)+log(dev)+dev,data=train_Ho,family=LOGNO(mu.link="identity",sigma.link="log"))
  #Calculate the density for Hoerl Curve
  #Density for Hoerl curve
  aug_valid_Ho<-aug_valid
  aug_valid_Ho$dev<-as.numeric(as.character(aug_valid$dev))
  dens_ODPHo2<-cal_dens_ODP(glm_ODP_Ho_tr1,aug_valid_Ho)
  dens_GaHo2<-cal_dens_GA(glm_Ga_Ho_tr1,tau_Ga,aug_valid_Ho)
  dens_LNHo2<-cal_dens_LN(glm_LN_Ho_tr1,tau_LN,aug_valid_Ho)
  
  ####With Calendar Periods
  #Under ODP Assumption
  glm_ODP_Cal_tr1<-glm(formula=value~factor(dev)+Calendar,family=quasipoisson(link="log"),data=train)
  #Under Gamma Assumption
  glm_Ga_Cal_tr1<-gamlss(formula=(value+tau_Ga)~factor(dev)+Calendar,data=train,family=GA(mu.link="log", sigma.link ="log"))
  #Under Log-Normal Assumption
  glm_LN_Cal_tr1<-gamlss(formula=(value+tau_LN)~factor(dev)+Calendar,data=train,family=LOGNO(mu.link="identity",sigma.link="log"))
 
  
  
 
  
  ####PPCI
  all_claims <- claims(
    frequency_vector = n_vector,
    occurrence_list = occurrence_times,
    claim_size_list = claim_sizes,
    notification_list = notidel,
    settlement_list = setldel,
    no_payments_list = no_payments,
    payment_size_list = payment_sizes,
    payment_delay_list = payment_delays,
    payment_time_list = payment_times,
    payment_inflated_list = payment_inflated
  )
  transaction_dataset <- generate_transaction_dataset(all_claims,adjust=TRUE)
  transaction_dataset$ReportingPeriod=ceiling(transaction_dataset$notidel)
  transaction_dataset$Calendar=transaction_dataset$ReportingPeriod+transaction_dataset$occurrence_period
  matri_NotiCount<-matrix(NA,nrow=40,ncol=40)
  
  for (i in 1:40){
    for (j in 1:40){
      if (i+j<=41){matri_NotiCount[i,j]=nrow(transaction_dataset[transaction_dataset$occurrence_period==i&transaction_dataset$ReportingPeriod==j&transaction_dataset$Calendar<=41,])}
      else{matri_NotiCount[i,j]='NA'}
    }
  }
  
  #Using the Chain-Ladder package to generate a triangular object:
  NotiCount_tri<-as.triangle(matri_NotiCount)
  #Convert it to a dataframe format
  NotiCount_dt<-as.data.frame(NotiCount_tri,na.rm=TRUE)
  NotiCount_dt$Calendar=as.numeric(NotiCount_dt$origin)+as.numeric(NotiCount_dt$dev)
  #Construct a full claim count square
  matri_NotiCount_sqrt<-matrix(NA,nrow=40,ncol=40)
  
  for (i in 1:40){
    for (j in 1:40){
      matri_NotiCount_sqrt[i,j]=nrow(transaction_dataset[transaction_dataset$occurrence_period==i&transaction_dataset$ReportingPeriod==j,])
    }
  }
  NotiCount_full_dt<-as.data.frame(as.triangle(matri_NotiCount_sqrt))
  NotiCount_full_dt$Calendar=as.numeric(NotiCount_full_dt$origin)+as.numeric(NotiCount_full_dt$dev)
  insample_NC<-NotiCount_full_dt[NotiCount_full_dt$Calendar<=41,]
  test_NC<-NotiCount_full_dt[NotiCount_full_dt$Calendar>41,]
  # Construction of Training Set for Notification Counts(to avoid data leakage):
  train_1_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)==1,]
  train_2_N<-NotiCount_dt[as.numeric(NotiCount_dt$Calendar)<=34&as.numeric(NotiCount_dt$origin)<=32&as.numeric(NotiCount_dt$origin)>=2,]
  train_3_N<-NotiCount_dt[NotiCount_dt$dev==1&as.numeric(NotiCount_dt$origin)>=33&as.numeric(NotiCount_dt$origin)<=40,]
  train_4_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)>=33&as.numeric(NotiCount_dt$origin)<=34&as.numeric(NotiCount_dt$dev)==2,]
  train_5_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)==40&as.numeric(NotiCount_dt$dev)==1,]
  train_NC<-rbind(train_1_N,train_2_N,train_3_N,train_4_N,train_5_N)
  #Construct of Validation set(identical to Paid Loss to avoid data leakage)
  valid_1_N<-NotiCount_dt[as.numeric(NotiCount_dt$Calendar)>=35&as.numeric(NotiCount_dt$Calendar)<=38&as.numeric(NotiCount_dt$origin)>=2&as.numeric(NotiCount_dt$origin)<=32,]
  valid_2_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)==33&as.numeric(NotiCount_dt$dev)>=3&as.numeric(NotiCount_dt$dev)<=5,]
  valid_3_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)==34&as.numeric(NotiCount_dt$dev)>=3&as.numeric(NotiCount_dt$dev)<=4,]
  valid_4_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)>=35&as.numeric(NotiCount_dt$origin)<=37&as.numeric(NotiCount_dt$dev)>=2&as.numeric(NotiCount_dt$dev)<=3,]
  valid_5_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)>=38&as.numeric(NotiCount_dt$origin)<=39&as.numeric(NotiCount_dt$dev)==2,]
  valid_6_N<-NotiCount_dt[as.numeric(NotiCount_dt$origin)==40&as.numeric(NotiCount_dt$dev)==1,]
  valid_N<-rbind(valid_1_N,valid_2_N,valid_3_N,valid_4_N,valid_5_N,valid_6_N)
  ##Fit a ODP/Chain-Ladder model to predict total claim notification count at each accident period 
  fit_nc<-glm(value~factor(origin)+factor(dev),data=train_NC,family=quasipoisson(link="log"))
  ##Obtain Estimated notification claims count for each accident period: Nk
  N<-c()
  N[1]<-sum(train_NC[train_NC$origin==1,]$value)
  for (i in 2:40){
    N[i]<-sum(train_NC[train_NC$origin==i,]$value)+sum(round(predict(fit_nc,newdata=test_NC[test_NC$origin==i,],type="response"),0))+sum(round(predict(fit_nc,newdata=valid_N[valid_N$origin==i,],type="response"),0))
  }
  ##Fit a second model for paid loss
  
  ##Payment per Notified Claim
  ###Create a function that counts the number of entries in each accident period
  count_by_origin<-function(data,start,end){
    cbo<-c()
    for (i in start:end){
      cbo[i]<-nrow(data[data$origin==i,]) 
    }
    return(cbo)
  }
  ###Sort the data by accident period
  new_train<-train[order(as.numeric(train$origin)),]
  ###Repete each element in N_rep by the number of entries in each accident period
  N_rep<-rep(N,times=count_by_origin(train,1,40))
  new_train$PPCI=new_train$value/N_rep
  
  ##Fitting a ODP model on PPCI
  fit_ODP_ppci<-glm(PPCI~factor(dev),family=quasipoisson(link="log"),data=new_train)
  
  #####PPCF
  transaction_dataset$SettlementPeriod=ceiling(transaction_dataset$setldel)
  matri_SettlementCount<-matrix(NA,nrow=40,ncol=40)
  
  for (i in 1:40){
    for (j in 1:40){
      if (i+j<=41){matri_SettlementCount[i,j]=nrow(transaction_dataset[transaction_dataset$occurrence_period==i&transaction_dataset$SettlementPeriod==j&transaction_dataset$Calendar<=41,])}
      else{matri_SettlementCount[i,j]='NA'}
    }
  }
  
  sqrt_SettlementCount<-matrix(NA,nrow=40,ncol=40)
  
  for (i in 1:40){
    for (j in 1:40){
      sqrt_SettlementCount[i,j]=nrow(transaction_dataset[transaction_dataset$occurrence_period==i&transaction_dataset$SettlementPeriod==j,])
    }
  }
  #Using the Chain-Ladder package to generate a triangular object:
  SettlCount_tri<-as.triangle(matri_SettlementCount)
  SettlCount_sqrt<-as.triangle(sqrt_SettlementCount)
  #Convert it into a dataframe
  SettlCount_dt<-as.data.frame(SettlCount_tri,na.rm=TRUE)
  SettlCount_dt$Calendar=as.numeric(SettlCount_dt$origin)+as.numeric(SettlCount_dt$dev)
  SettlCount_in<-SettlCount_dt
  SettlCount_full<-as.data.frame(SettlCount_sqrt)
  SettlCount_full$Calendar=as.numeric(SettlCount_full$origin)+as.numeric(SettlCount_full$dev)
  SettlCount_out<-SettlCount_full[SettlCount_full$Calendar>41,]
  # Construction of Training Set for Finalized Counts(to avoid data leakage):
  train_1_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)==1,]
  train_2_F<-SettlCount_dt[as.numeric(SettlCount_dt$Calendar)<=34&as.numeric(SettlCount_dt$origin)<=32&as.numeric(SettlCount_dt$origin)>=2,]
  train_3_F<-SettlCount_dt[SettlCount_dt$dev==1&as.numeric(SettlCount_dt$origin)>=33&as.numeric(SettlCount_dt$origin)<=40,]
  train_4_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)>=33&as.numeric(SettlCount_dt$origin)<=34&as.numeric(SettlCount_dt$dev)==2,]
  train_5_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)==40&as.numeric(SettlCount_dt$dev)==1,]
  train_F<-rbind(train_1_F,train_2_F,train_3_F,train_4_F,train_5_F)
  train_F<-train_F[order(as.numeric(train_F$origin)),]
  
  
  #Construct of Validation set(identical to Paid Loss to avoid data leakage)
  valid_1_F<-SettlCount_dt[as.numeric(SettlCount_dt$Calendar)>=35&as.numeric(SettlCount_dt$Calendar)<=38&as.numeric(SettlCount_dt$origin)>=2&as.numeric(SettlCount_dt$origin)<=32,]
  valid_2_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)==33&as.numeric(SettlCount_dt$dev)>=3&as.numeric(SettlCount_dt$dev)<=5,]
  valid_3_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)==34&as.numeric(SettlCount_dt$dev)>=3&as.numeric(SettlCount_dt$dev)<=4,]
  valid_4_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)>=35&as.numeric(SettlCount_dt$origin)<=37&as.numeric(SettlCount_dt$dev)>=2&as.numeric(SettlCount_dt$dev)<=3,]
  valid_5_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)>=38&as.numeric(SettlCount_dt$origin)<=39&as.numeric(SettlCount_dt$dev)==2,]
  valid_6_F<-SettlCount_dt[as.numeric(SettlCount_dt$origin)==40&as.numeric(SettlCount_dt$dev)==1,]
  valid_F<-rbind(valid_1_F,valid_2_F,valid_3_F,valid_4_F,valid_5_F,valid_6_F)
  valid_F<-valid_F[order(as.numeric(valid_F$origin)),]
  #Create a triangle for Cummulative Claims Count
  cum_Count_tri<-incr2cum(NotiCount_tri)
  cum_Count_dt<-as.data.frame(cum_Count_tri,na.rm=TRUE)
  cum_Count_dt<-cum_Count_dt[order(as.numeric(cum_Count_dt$origin)),]
  cum_Count_dt$Calendar=as.numeric(as.character(cum_Count_dt$origin))+as.numeric(cum_Count_dt$dev)
  #Create a triangle for Cummulative finalised claims count
  cum_F_tri<-incr2cum(SettlCount_tri)
  cum_F_dt<-as.data.frame(cum_F_tri,na.rm=TRUE)
  cum_F_dt<-cum_F_dt[order(as.numeric(cum_F_dt$origin)),]
  cum_F_dt$Calendar=as.numeric(as.character(cum_F_dt$origin))+as.numeric(cum_F_dt$dev)
  # Construction of Training Set for Cummulative Notified Claim Counts(to avoid data leakage):
  train_1_cumN<-cum_Count_dt[as.numeric(as.character(cum_Count_dt$origin))==1,]
  train_2_cumN<-cum_Count_dt[as.numeric(cum_Count_dt$Calendar)<=34&as.numeric(as.character(cum_Count_dt$origin))<=32&as.numeric(as.character(cum_Count_dt$origin))>=2,]
  train_3_cumN<-cum_Count_dt[cum_Count_dt$dev==1&as.numeric(as.character(cum_Count_dt$origin))>=33&as.numeric(as.character(cum_Count_dt$origin))<=40,]
  train_4_cumN<-cum_Count_dt[as.numeric(as.character(cum_Count_dt$origin))>=33&as.numeric(as.character(cum_Count_dt$origin))<=34&as.numeric(cum_Count_dt$dev)==2,]
  train_5_cumN<-cum_Count_dt[as.numeric(as.character(cum_Count_dt$origin))==40&as.numeric(cum_Count_dt$dev)==1,]
  train_cumN<-rbind(train_1_cumN,train_2_cumN,train_3_cumN,train_4_cumN,train_5_cumN)
  train_cumN<-train_cumN[order(as.numeric(train_cumN$origin)),]
  # Construction of Training Set for Cummulative Finalized Counts(to avoid data leakage):
  train_1_cumF<-cum_F_dt[as.numeric(cum_F_dt$origin)==1,]
  train_2_cumF<-cum_F_dt[as.numeric(cum_F_dt$Calendar)<=34&as.numeric(cum_F_dt$origin)<=32&as.numeric(cum_F_dt$origin)>=2,]
  train_3_cumF<-cum_F_dt[cum_F_dt$dev==1&as.numeric(cum_F_dt$origin)>=33&as.numeric(cum_F_dt$origin)<=40,]
  train_4_cumF<-cum_F_dt[as.numeric(cum_F_dt$origin)>=33&as.numeric(cum_F_dt$origin)<=34&as.numeric(cum_F_dt$dev)==2,]
  train_5_cumF<-cum_F_dt[as.numeric(cum_F_dt$origin)==40&as.numeric(cum_F_dt$dev)==1,]
  train_cumF<-rbind(train_1_cumF,train_2_cumF,train_3_cumF,train_4_cumF,train_5_cumF)
  train_cumF<-train_cumF[order(as.numeric(train_cumF$origin)),]
  cum_U_dt<-data.frame(origin=cum_F_dt$origin,dev=cum_F_dt$dev,value=cum_Count_dt$value-cum_F_dt$value)
  #Generate a full data frame containing finalized claim count, unfinalized claim count and notification claim count
  NotiCount_dt<-NotiCount_dt[order(as.numeric(NotiCount_dt$origin)),]
  SettlCount_dt<-SettlCount_dt[order(as.numeric(SettlCount_dt$origin)),]
  UplusN<-list()
  
  
  #UplusN[1]<-NA
  #UplusN[2:nrow(cum_U_dt)]<-cum_U_dt$value[1:(nrow(cum_U_dt)-1)]+NotiCount_dt$value[2:nrow(NotiCount_dt)]
  
  for (i in 1:40){
    UplusN[[i]]<-cum_U_dt[cum_U_dt$origin==i,]$value[1:(nrow(cum_U_dt[cum_U_dt$origin==i,])-1)]+NotiCount_dt[NotiCount_dt$origin==i,]$value[2:nrow(NotiCount_dt[NotiCount_dt$origin==i,])]
  }
  
  UplusN<-as.vector(unlist(UplusN))
  FCount<-SettlCount_dt[SettlCount_dt$dev!=1,]$value
  FUN_dat<-data.frame(origin=SettlCount_dt[SettlCount_dt$dev!=1,]$origin,dev=SettlCount_dt[SettlCount_dt$dev!=1,]$dev,UplusN=UplusN[1:(length(UplusN)-2)],FCount=FCount)
  FUN_dat$Calendar=as.numeric(as.character(FUN_dat$origin))+as.numeric(FUN_dat$dev)
  #full_FUN<-data.frame(origin=NotiCount_dt$origin,dev=NotiCount_dt$dev,U=cum_U_dt$value,N=NotiCount_dt$value,UplusN=UplusN,F=SettlCount_dt$value)
  #View(full_FUN)
  addrow<-data.frame(origin=as.factor(40),dev=1,UplusN=UplusN[length(UplusN)],FCount=SettlCount_dt[SettlCount_dt$origin==40,]$value,Calendar=41)
  FUN_dat<-rbind(FUN_dat,addrow)
  train_1_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==1,]
  train_2_FC<-FUN_dat[as.numeric(FUN_dat$Calendar)<=34&as.numeric(as.character(FUN_dat$origin))<=32&as.numeric(as.character(FUN_dat$origin))>=2,]
  train_3_FC<-FUN_dat[FUN_dat$dev==1&as.numeric(as.character(FUN_dat$origin))>=33&as.numeric(as.character(FUN_dat$origin))<=40,]
  train_4_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))>=33&as.numeric(as.character(FUN_dat$origin))<=34&as.numeric(FUN_dat$dev)==2,]
  train_5_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==40&as.numeric(FUN_dat$dev)==1,]
  train_FC<-rbind(train_1_FC,train_2_FC,train_3_FC,train_4_FC,train_5_FC)
  train_FC<-train_FC[order(as.numeric(as.character(train_FC$origin))),]
  
  valid_1_FC<-FUN_dat[as.numeric(FUN_dat$Calendar)>=35&as.numeric(FUN_dat$Calendar)<=38&as.numeric(as.character(FUN_dat$origin))>=2&as.numeric(as.character(FUN_dat$origin))<=32,]
  valid_2_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==33&as.numeric(FUN_dat$dev)>=3&as.numeric(FUN_dat$dev)<=5,]
  valid_3_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==34&as.numeric(FUN_dat$dev)>=3&as.numeric(FUN_dat$dev)<=4,]
  valid_4_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))>=35&as.numeric(as.character(FUN_dat$origin))<=37&as.numeric(FUN_dat$dev)>=2&as.numeric(FUN_dat$dev)<=3,]
  valid_5_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))>=38&as.numeric(as.character(FUN_dat$origin))<=39&as.numeric(FUN_dat$dev)==2,]
  valid_6_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==40&as.numeric(FUN_dat$dev)==1,]
  valid_FC<-rbind(valid_1_FC,valid_2_FC,valid_3_FC,valid_4_FC,valid_5_FC,valid_6_FC)
  valid_FC<-valid_FC[order(as.numeric(as.character(valid_FC$origin))),]
  
  # Construction of Test Data
  test_1_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))>=2&as.numeric(as.character(FUN_dat$origin))<=35&as.numeric(as.character(FUN_dat$Calendar))>=39&as.numeric(as.character(FUN_dat$Calendar))<=41,]
  test_2_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==36&as.numeric(as.character(FUN_dat$dev))>=4&as.numeric(as.character(FUN_dat$dev))<=5,]
  test_3_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==37&as.numeric(as.character(FUN_dat$dev))==4,]
  test_4_FC<-FUN_dat[as.numeric(as.character(FUN_dat$origin))==38&as.numeric(as.character(FUN_dat$dev))==3,]
  test_FC<-rbind(test_1_FC,test_2_FC,test_3_FC,test_4_FC)
  test_FC<-test_FC[order(as.numeric(as.character(test_FC$origin))),]
  aug_valid_FC<-rbind(valid_FC,test_FC)
  aug_valid_FC<-aug_valid_FC[order(as.numeric(as.character(aug_valid_FC$origin))),]
  #Alternatively, fit a ODP on finalization claims that depends on development periods
  odp_FC<-glm(FCount~factor(dev),data=train_FC,family=quasipoisson(link="log"))
  new_train$PPCF=new_train$value/train_F$value
  #Create a column for operation time
  new_train$OT=train_cumF$value/N_rep
  #Remove the cells with zero finalized claim count but with positive payments; 
  new_train_PPCF<-na.omit(new_train[new_train$PPCF!=Inf,])
  ODP_PPCF_train<-glm(PPCF~OT,family=quasipoisson(link="log"),data=new_train_PPCF)
  #Density for PPCI and PPCF
  dens_PPCI2<-cal_dens_PPCI(fit_ODP_ppci,N,aug_valid,2,40)
  aug_valid_FC$dev=as.factor(aug_valid_FC$dev)
  dens_PPCF2<-cal_PPCF_dens(odp_FC,ODP_PPCF_train,train_cumF,aug_valid_FC,aug_valid,N,2,40) 
  
  
  
  mu_PPCI2<-cal_mu_PPCI(fit_ODP_ppci,N,aug_valid,2,40)
  mu_PPCF2<-cal_PPCF_mu(odp_FC,ODP_PPCF_train,train_cumF,aug_valid_FC,aug_valid,N,2,40)
  
  ###Fitting of all the models in the in-sample 
  #####Models with Basic Structures
  #Fitting of ODP GLM
  ODP_GLM_in<-glm(formula=value~factor(origin)+factor(dev),family=quasipoisson(link="log"),data=in_sample)
  #Fitting of ZAGA
  
  ZAGA_in<-gamlss(formula=value~factor(origin)+factor(dev),nu.formula=~as.numeric(as.character(dev)),data=in_sample,family=ZAGA(mu.link="log",sigma.link = "log", nu.link = "logit"),trace=FALSE)
  #Fitting of ZALN
  LN_in<-gamlssZadj(y=value,mu.formula =~factor(origin)+factor(dev),xi0.formula=~as.numeric(as.character(dev)),data=in_sample,family=LOGNO(mu.link="identity",sigma.link="log"))
  #Fitting of opt_LN:
  tau_LN<-5
  LN_optimTau<-gamlss(formula=(value+tau_LN)~factor(origin)+factor(dev),data=in_sample,family=LOGNO(mu.link="identity",sigma.link="log"),trace=FALSE)
  tau_Ga<-5
  Ga_optimTau<-gamlss(formula=(value+tau_Ga)~factor(origin)+factor(dev),data=in_sample,family=GA(mu.link="log", sigma.link ="log"),trace=FALSE)
  
  ###Hoerl Curve
  #Under ODP Assumption
  glm_ODP_Ho_In<-glm(formula=value~factor(origin)+log(as.numeric(as.character(dev)))+as.numeric(as.character(dev)),family=quasipoisson(link="log"),data=in_sample,trace=FALSE)
  
  #Density for Hoerl curve
  out_sample_Ho<-out_sample
  out_sample_Ho$dev<-as.numeric(as.character(out_sample$dev))
  
  
  #Under Gamma Assumption
  in_sample_Ho<-in_sample
  in_sample_Ho$dev<-as.numeric(as.character(in_sample$dev))
  glm_Ga_Ho_In<-gamlss(formula=(value+tau_Ga)~factor(origin)+log(dev)+dev,data=in_sample_Ho,family=GA(mu.link="log", sigma.link ="log"),trace=FALSE)
  
  #Under Log-Normal Assumption
  glm_LN_Ho_In<-gamlss(formula=(value+tau_LN)~factor(origin)+log(dev)+dev,data=in_sample_Ho,family=LOGNO(mu.link="identity",sigma.link="log"),trace=FALSE)
  
  ###GLM with Calendar Periods Effects
  
  #Under ODP Assumption
  
  glm_ODP_Cal_In<-glm(formula=value~factor(dev)+Calendar,family=quasipoisson(link="log"),data=in_sample)
  
  #Under Gamma Assumption
  
  glm_Ga_Cal_In<-gamlss(formula=(value+tau_Ga)~factor(dev)+Calendar,data=in_sample,family=GA(mu.link="log", sigma.link ="log"),trace=FALSE)
  #predict(glm_Ga_Cal_In,what="mu",newdata=out_sample,type="response")
  #Under Log-Normal Assumption
  glm_LN_Cal_In<-gamlss(formula=(value+tau_LN)~factor(dev)+Calendar,data=in_sample,family=LOGNO(mu.link="identity",sigma.link="log"),trace=FALSE)
  in_sample_numeric<-in_sample
  in_sample_numeric$origin=as.numeric(as.character(in_sample$origin))
  in_sample_numeric$dev=as.numeric(as.character(in_sample$dev))
  in_sample_numeric$Calendar=as.numeric(as.character(in_sample$Calendar))
  
  ###Smoothing Spline
  sp_Normal_In<-gamlss(formula=value~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,family=NO(),trace=FALSE)
  sp_Gamma_In<-gamlss(formula=(value+tau_Ga)~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,family=GA(mu.link="log", sigma.link ="log"),trace=FALSE)
  sp_LN_In<-gamlss(formula=(value+tau_LN)~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,family=LOGNO(mu.link="identity",sigma.link="log"),trace=FALSE)
  
  ###GAMLSS
  gamlss2_GA_In<-gamlss(formula=(value+tau_Ga)~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,sigma.formula=~cs(as.numeric(as.character(dev))),family=GA(mu.link="log", sigma.link ="log"),trace=FALSE)
  gamlss2_LN_In<-gamlss(formula=(value+tau_LN)~scs(as.numeric(as.character(origin)))+scs(as.numeric(as.character(dev))),data=in_sample_numeric,sigma.formula=~cs(as.numeric(as.character(dev))),family=LOGNO(mu.link="identity",sigma.link="log"),trace=FALSE)
  
  

  
  out_sample_numeric<-out_sample
  out_sample_numeric$origin=as.numeric(as.character(out_sample$origin))
  out_sample_numeric$dev=as.numeric(as.character(out_sample$dev))
  out_sample_numeric$Calendar=as.numeric(as.character(out_sample$Calendar))
  
  
  #Fit a PPCI model
  
  ##Fit a ODP/Chain-Ladder model to predict total claim notification count at each accident period 
  fit_nc_In<-glm(value~factor(origin)+factor(dev),data=insample_NC,family=quasipoisson(link="log"))
  ##Obtain Estimated notification claims count for each accident period: Nk
  N_in<-c()
  N_in[1]<-sum(insample_NC[insample_NC$origin==1,]$value)
  for (i in 2:40){
    N_in[i]<-sum(insample_NC[insample_NC$origin==i,]$value)+sum(round(predict(fit_nc_In,newdata=test_NC[test_NC$origin==i,],type="response"),0))
  }
  
  ##Fit a second model for paid loss
  
  ##Payment per Notified Claim
  ###Create a function that counts the number of entries in each accident period
  count_by_origin<-function(data,start,end){
    cbo<-c()
    for (i in start:end){
      cbo[i]<-nrow(data[data$origin==i,]) 
    }
    return(cbo)
  }
  ###Sort the data by accident period
  in_sample<-in_sample[order(as.numeric(as.character(in_sample$origin))),]
  ###Repete each element in N_rep by the number of entries in each accident period
  N_rep_In<-rep(N_in,times=count_by_origin(in_sample,1,40))
  in_sample$PPCI=in_sample$value/N_rep_In
  
  ##Fitting a ODP model on PPCI
  fit_ODP_ppci_In<-glm(PPCI~factor(dev),family=quasipoisson(link="log"),data=in_sample)
  
  ###For PPCF model
  transaction_dataset$SettlementPeriod=ceiling(transaction_dataset$setldel)
  transaction_dataset$ReportingPeriod=ceiling(transaction_dataset$notidel)
  transaction_dataset$Calendar=transaction_dataset$ReportingPeriod+transaction_dataset$occurrence_period
  matri_SettlementCount<-matrix(NA,nrow=40,ncol=40)
  transaction_dataset$Calendar=as.numeric(transaction_dataset$occurrence_period)+as.numeric(transaction_dataset$SettlementPeriod)
  for (i in 1:40){
    for (j in 1:40){
      if (i+j<=41){matri_SettlementCount[i,j]=nrow(transaction_dataset[transaction_dataset$occurrence_period==i&transaction_dataset$SettlementPeriod==j&transaction_dataset$Calendar<=41,])}
      else{matri_SettlementCount[i,j]='NA'}
    }
  }
  
  sqrt_SettlementCount<-matrix(NA,nrow=40,ncol=40)
  
  for (i in 1:40){
    for (j in 1:40){
      sqrt_SettlementCount[i,j]=nrow(transaction_dataset[transaction_dataset$occurrence_period==i&transaction_dataset$SettlementPeriod==j,])
    }
  }
  
  #Using the Chain-Ladder package to generate a triangular object:
  SettlCount_tri<-as.triangle(matri_SettlementCount)
  SettlCount_sqrt<-as.triangle(sqrt_SettlementCount)
  
  #Convert it into a dataframe
  SettlCount_dt<-as.data.frame(SettlCount_tri,na.rm=TRUE)
  SettlCount_dt$Calendar=as.numeric(SettlCount_dt$origin)+as.numeric(SettlCount_dt$dev)
  SettlCount_in<-SettlCount_dt
  
  SettlCount_full<-as.data.frame(SettlCount_sqrt)
  SettlCount_full$Calendar=as.numeric(SettlCount_full$origin)+as.numeric(SettlCount_full$dev)
  SettlCount_out<-SettlCount_full[SettlCount_full$Calendar>41,]
  
  
  #Create a triangle for Full Cummulative Claims Count
  cum_Count_tri_full<-incr2cum(as.triangle(matri_NotiCount_sqrt))
  cum_Count_dt_full<-as.data.frame(cum_Count_tri_full,na.rm=TRUE)
  cum_Count_dt_full<-cum_Count_dt_full[order(as.numeric(cum_Count_dt_full$origin)),]
  cum_Count_dt_full$Calendar=as.numeric(as.character(cum_Count_dt_full$origin))+as.numeric(cum_Count_dt_full$dev)
  
  #Create a triangle for In-Sample Cummulative finalised claims count
  cum_F_tri<-incr2cum(SettlCount_tri)
  cum_F_dt<-as.data.frame(cum_F_tri,na.rm=TRUE)
  cum_F_dt<-cum_F_dt[order(as.numeric(cum_F_dt$origin)),]
  cum_F_dt$Calendar=as.numeric(as.character(cum_F_dt$origin))+as.numeric(cum_F_dt$dev)
  
  #Create a triangle for Out-Sample Cummulative finalised claims count
  cum_F_tri_full<-incr2cum(SettlCount_sqrt)
  cum_F_dt_full<-as.data.frame(cum_F_tri_full)
  cum_F_dt_full$Calendar=as.numeric(as.character(cum_F_dt_full$origin))+as.numeric(as.character(cum_F_dt_full$dev))
  cum_F_dt_full<-cum_F_dt_full[order(as.numeric(as.character(cum_F_dt_full$origin))),]
  cum_F_dt_out<-cum_F_dt_full[cum_F_dt_full$Calendar>41,]
  
  #Create a triangle for Unfinalized claim count triangle (full)
  cum_U_dt_full<-data.frame(origin=cum_F_dt_full$origin,dev=cum_F_dt_full$dev,Calendar=cum_F_dt_full$Calendar,value=cum_Count_dt_full$value-cum_F_dt_full$value)
  cum_U_dt_full<-cum_U_dt_full[order(as.numeric(as.character(cum_U_dt_full$origin))),]
  cum_U_dt_In<-cum_U_dt_full[cum_U_dt_full$Calendar<=41,]
  cum_U_dt_Out<-cum_U_dt_full[cum_U_dt_full$Calendar>41,]
  
  #Generate a full data frame containing finalized claim count, unfinalized claim count and notification claim count
  
  NotiCount_full_dt<-NotiCount_full_dt[order(as.numeric(NotiCount_full_dt$origin)),]
  UplusN_full<-list()
  SettlCount_full<-SettlCount_full[order(as.numeric(as.character(SettlCount_full$origin))),]
  
  
  for (i in 1:40){
    UplusN_full[[i]]<-cum_U_dt_full[cum_U_dt_full$origin==i,]$value[1:(nrow(cum_U_dt_full[cum_U_dt_full$origin==i,])-1)]+NotiCount_full_dt[NotiCount_full_dt$origin==i,]$value[2:nrow(NotiCount_full_dt[NotiCount_full_dt$origin==i,])]
  }
  
  UplusN_full<-as.vector(unlist(UplusN_full))
  FCount_full<-SettlCount_full[SettlCount_full$dev!=1,]$value
  
  FUN_dat_full<-data.frame(origin=SettlCount_full[SettlCount_full$dev!=1,]$origin,dev=SettlCount_full[SettlCount_full$dev!=1,]$dev,UplusN=UplusN_full,FCount=FCount_full)
  
  FUN_dat_full$Calendar=as.numeric(as.character(FUN_dat_full$origin))+as.numeric(FUN_dat_full$dev)
  FUN_dat_full<-FUN_dat_full[order(as.numeric(as.character(FUN_dat_full$origin))),]
  addrow<-data.frame(origin=as.factor(40),dev=1,UplusN=UplusN[length(UplusN)],FCount=SettlCount_dt[SettlCount_dt$origin==40,]$value,Calendar=41)
  FUN_dat_full<-rbind(FUN_dat_full,addrow)
  
  FUN_dat_In<-FUN_dat_full[FUN_dat_full$Calendar<=41,]
  FUN_dat_Out<-FUN_dat_full[FUN_dat_full$Calendar>41,]
  
  in_sample$PPCF=in_sample$value/SettlCount_in$value
  odp_FC_In<-glm(FCount~factor(dev),data=FUN_dat_In,family=quasipoisson(link="log"))
  #Create a column for operation time
  in_sample$OT=cum_F_dt$value/N_rep_In
  #Remove the cells with zero finalized claim count but with positive payments; 
  new_in_PPCF<-na.omit(in_sample[in_sample$PPCF!=Inf,])
  ODP_PPCF_In<-glm(PPCF~OT,family=quasipoisson(link="log"),data=new_in_PPCF)
  
  #CRPS for BMV
  SpLN_crps[,D]<-cal_crps_LN(z,sp_LN_In,out_sample_numeric)
  #CRPS for EW
  EqEns_crps[,D]<-cal_CRPS_ensemble(w=rep(1/18,18),z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  #CRPS for SLP
  SLP_crps[,D]<-cal_CRPS_ensemble(w=as.vector(unlist(model_weights_simul_par0[[D]])),z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  
  #CRPS for ADLP8
  ADLP8_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par7_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par7_new[[D]][2])),index_subset1=ind_subset1_par7_new,index_subset2=ind_subset2_par7_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  
  #CRPS for others ADLP
  ADLP1_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par1_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par1_new[[D]][2])),index_subset1=ind_subset1_par1_new,index_subset2=ind_subset2_par1_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP2_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par2_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par2_new[[D]][2])),index_subset1=ind_subset1_par2_new,index_subset2=ind_subset2_par2_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP3_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par3_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par3_new[[D]][2])),index_subset1=ind_subset1_par3_new,index_subset2=ind_subset2_par3_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP4_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par4_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par4_new[[D]][2])),index_subset1=ind_subset1_par4_new,index_subset2=ind_subset2_par4_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP5_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par5_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par5_new[[D]][2])),index_subset1=ind_subset1_par5_new,index_subset2=ind_subset2_par5_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP6_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par6_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par6_new[[D]][2])),index_subset1=ind_subset1_par6_new,index_subset2=ind_subset2_par6_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP7_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par8_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par8_new[[D]][2])),index_subset1=ind_subset1_par8_new,index_subset2=ind_subset2_par8_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP9_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par9_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par9_new[[D]][2])),index_subset1=ind_subset1_par9_new,index_subset2=ind_subset2_par9_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  
  ADLP10_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par10_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par10_new[[D]][2])),index_subset1=ind_subset1_par10_new,index_subset2=ind_subset2_par10_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP11_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par11_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par11_new[[D]][2])),index_subset1=ind_subset1_par11_new,index_subset2=ind_subset2_par11_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP12_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par12_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par12_new[[D]][2])),index_subset1=ind_subset1_par12_new,index_subset2=ind_subset2_par12_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP13_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par13_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par13_new[[D]][2])),index_subset1=ind_subset1_par13_new,index_subset2=ind_subset2_par13_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP14_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par14_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par14_new[[D]][2])),index_subset1=ind_subset1_par14_new,index_subset2=ind_subset2_par14_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP15_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par15_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par15_new[[D]][2])),index_subset1=ind_subset1_par15_new,index_subset2=ind_subset2_par15_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  
  ADLP16_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par16_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par16_new[[D]][2])),index_subset1=ind_subset1_par16_new,index_subset2=ind_subset2_par16_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  ADLP17_crps[,D]<-cal_CRPS_ADLPensemble(w1=as.vector(unlist(model_weights_simul_par17_new[[D]][1])),w2=as.vector(unlist(model_weights_simul_par17_new[[D]][2])),index_subset1=ind_subset1_par17_new,index_subset2=ind_subset2_par17_new,tau_Ga=tau_Ga,tau_LN=tau_LN,z=z,newdata=out_sample,ODPGLM=ODP_GLM_in,GAGLM=Ga_optimTau,LNGLM=LN_optimTau,ZAGA=ZAGA_in,ZALN=LN_in,ODPHo=glm_ODP_Ho_In,GaHo=glm_Ga_Ho_In,LNHo=glm_LN_Ho_In,ODPCal=glm_ODP_Cal_In,GaCAL=glm_Ga_Cal_In,LNCAL=glm_LN_Cal_In,NoSp=sp_Normal_In,GaSp=sp_Gamma_In,LNSp=sp_LN_In,GaGAMLSS=gamlss2_GA_In,LNGAMLSS=gamlss2_LN_In,PPCI=fit_ODP_ppci_In,NO_PPCI=N_in,index_start=2,index_end=40,odp_FC=odp_FC_In,model_subPayments=ODP_PPCF_In,train_cumF=cum_F_dt,newdataFC=FUN_dat_Out,newdata_Pay=out_sample,NO_PPCF=N_in)
  
  
  
  
}  

save(SpLN_crps,file="SpLN_crps")
save(EqEns_crps,file="SpLN_crps")

save(SLP_crps,file="SpLN_crps")
save(ADLP8_crps,file="SpLN_crps")

save(ADLP1_crps,file="SpLN_crps")
save(ADLP2_crps,file="SpLN_crps")
save(ADLP3_crps,file="SpLN_crps")
save(ADLP4_crps,file="SpLN_crps")
save(ADLP5_crps,file="SpLN_crps")
save(ADLP6_crps,file="SpLN_crps")
save(ADLP7_crps,file="SpLN_crps")
save(ADLP9_crps,file="SpLN_crps")
save(ADLP10_crps,file="SpLN_crps")
save(ADLP11_crps,file="SpLN_crps")
save(ADLP12_crps,file="SpLN_crps")
save(ADLP13_crps,file="SpLN_crps")
save(ADLP14_crps,file="SpLN_crps")
save(ADLP15_crps,file="SpLN_crps")
save(ADLP16_crps,file="SpLN_crps")
save(ADLP17_crps,file="SpLN_crps")


mean(apply(SpLN_crps,MARGIN=2,FUN=mean))
mean(apply(EqEns_crps,MARGIN=2,FUN=mean))
mean(apply(ADLP7_crps,MARGIN=2,FUN=mean))
mean(apply(ADLP5_crps,MARGIN=2,FUN=mean))
mean(apply(ADLP13_crps,MARGIN=2,FUN=mean))