
#marginal parameters. In this setting, I only change the periodic terms sine and cosine in the marginal model into unperiodic terms, other parts are identical to scenario 3. 
beta0<-2.5
beta1<-3.5
beta2<-4.5


#generate random numbers from a specific R-Vine Structure
library(VineCopula)
library(nlme)
library(forecast)


#time block vine structure
RVineblock<-c(4,3,2,1,0,3,2,1,0,0,2,1,0,0,0,1)
Rvineblock<-matrix(RVineblock,4,4)

#Time block family and parameters
RVineblockf<-c(0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0)
Rvineblockf<-matrix(RVineblockf,4,4)

RVineblockpar<-c(0,0.2,0.3,0.5,0,0,0.3,0.5,0,0,0,0.5,0,0,0,0)
Rvineblockpar<-matrix(RVineblockpar,4,4)

#lower
lo<-c(rep(-1,6),rep(c(-Inf,-Inf,-Inf,0),4))

#upper
up<-c(rep(1,6),rep(c(Inf,Inf,Inf,Inf),4))


#specify sample size and time block
n<-500
k<-5
l<-4


MAE_CVINE<-rep(0,200)
MAE_CVINE1<-rep(0,200)
MAE_CVINE2<-rep(0,200)
MAE_CVINE3<-rep(0,200)
MAE_MRM<-rep(0,200)
MAE_MRM1<-rep(0,200)
MAE_MRM2<-rep(0,200)
MAE_LRM<-rep(0,200)
MAE_AR<-rep(0,200)
percentage1<-rep(0,200)
percentage2<-rep(0,200)
percentage3<-rep(0,200)
#the simulation is repeated 200 times
for(m in 1:200){
  set.seed(m)
  ar_simM<-matrix(0,k*l,n)
  for(i in 1:n){
    ar_simM[,i]<-arima.sim(list(ar=c(0.5)),k*l)
  }
  ar_sim<-as.vector(ar_simM)
  
  #generate x from uniform(0,100)
  x<-runif(n*k*l,0,100)
  
  #label each data point with time j
  j<-rep(seq(1,(k*l),1),n)
  z<-rep(c(rep(1,4),rep(2,4),rep(3,4),rep(4,4),rep(5,4)),n)
  t<-rep(seq(1,5,1),(n*k))
  
  #create response variables(remove periodic terms and add unperiodic terms j)
  
  y<-beta0+beta1*x+beta2*j+ar_sim
  
  #seperate the data into training set and validation set 
  ytrain<-NULL
  xtrain<-NULL
  jtrain<-NULL
  ztrain<-NULL
  ttrain<-NULL
  
  yvalid<-NULL
  xvalid<-NULL
  jvalid<-NULL
  zvalid<-NULL
  tvalid<-NULL
  for(i in 1:n){
    ytrain=c(ytrain,y[((i-1)*k*l+1):((i-1)*k*l+(k-1)*l)])
    xtrain=c(xtrain,x[((i-1)*k*l+1):((i-1)*k*l+(k-1)*l)])
    jtrain=c(jtrain,j[((i-1)*k*l+1):((i-1)*k*l+(k-1)*l)])
    ztrain=c(ztrain,z[((i-1)*k*l+1):((i-1)*k*l+(k-1)*l)])
    ttrain=c(ttrain,t[((i-1)*k*l+1):((i-1)*k*l+(k-1)*l)])
    
    yvalid=c(yvalid,y[((i-1)*k*l+(k-1)*l+1):(i*k*l)])
    xvalid=c(xvalid,x[((i-1)*k*l+(k-1)*l+1):(i*k*l)])
    jvalid=c(jvalid,j[((i-1)*k*l+(k-1)*l+1):(i*k*l)])
    zvalid=c(zvalid,z[((i-1)*k*l+(k-1)*l+1):(i*k*l)])
    tvalid=c(tvalid,t[((i-1)*k*l+(k-1)*l+1):(i*k*l)])
  }
  
  
  
  
  ytrainM<-matrix(ytrain,l,(length(ytrain)/l))
  ytrain1<-ytrainM[1,]
  ytrain2<-ytrainM[2,]
  ytrain3<-ytrainM[3,]
  ytrain4<-ytrainM[4,]
  
  yvalidM<-matrix(yvalid,l,(length(yvalid)/l))
  yvalid1<-yvalidM[1,]
  yvalid2<-yvalidM[2,]
  yvalid3<-yvalidM[3,]
  yvalid4<-yvalidM[4,]
  yvalidadjusted<-c(yvalid1,yvalid2,yvalid3,yvalid4)
  
  xtrainM<-matrix(xtrain,l,length(xtrain)/l)
  xtrain1<-xtrainM[1,]
  xtrain2<-xtrainM[2,]
  xtrain3<-xtrainM[3,]
  xtrain4<-xtrainM[4,]
  
  xvalidM<-matrix(xvalid,l,(length(xvalid)/l))
  xvalid1<-xvalidM[1,]
  xvalid2<-xvalidM[2,]
  xvalid3<-xvalidM[3,]
  xvalid4<-xvalidM[4,]
  
  ztrainM<-matrix(ztrain,l,length(ztrain)/l)
  ztrain1<-ztrainM[1,]
  ztrain2<-ztrainM[2,]
  ztrain3<-ztrainM[3,]
  ztrain4<-ztrainM[4,]
  
  zvalidM<-matrix(zvalid,l,(length(zvalid)/l))
  zvalid1<-zvalidM[1,]
  zvalid2<-zvalidM[2,]
  zvalid3<-zvalidM[3,]
  zvalid4<-zvalidM[4,]
  
  #fit AR1 model
  sin<-sin(jtrain/l*2*pi)
  cos<-cos(jtrain/l*2*pi)
  fit0<-gls(ytrain~xtrain+jtrain,correlation = corARMA(form=~1|jtrain,p=1,q=0))
  
  #fit LRM
  fit<-lm(ytrain~xtrain+ztrain+ttrain)
  
  #fit CVM and MRM
  fit1<-lm(ytrain1~xtrain1+ztrain1)
  fit2<-lm(ytrain2~xtrain2+ztrain2)
  fit3<-lm(ytrain3~xtrain3+ztrain3)
  fit4<-lm(ytrain4~xtrain4+ztrain4)
  
  u1<-pnorm(fit1$residuals,0,(var(fit1$residuals))^0.5)
  u2<-pnorm(fit2$residuals,0,(var(fit2$residuals))^0.5)
  u3<-pnorm(fit3$residuals,0,(var(fit3$residuals))^0.5)
  u4<-pnorm(fit4$residuals,0,(var(fit4$residuals))^0.5)
  uM<-cbind(u1,u2,u3,u4)
  
  e<-seq(1,4,1)
  f<-rep(0,6)
  g<-rep(0,6)
  h<-rep(0,6)
  CVINEMATRIX<-C2RVine(order=e,family=f,par=g,par2=h)$Matrix
  RVM2<-RVineCopSelect(uM,family=c(1,3,4,5,6,13,14,26,36),Matrix=CVINEMATRIX,selectioncrit = "BIC",indeptest = TRUE)
  
  #corresponding upper and lower bound
  lo1<-c(rep(0,6),rep(c(-Inf,-Inf,-Inf,0),4))
  up1<-c(rep(0,6),rep(c(Inf,Inf,Inf,Inf),4))
  num<-1
  for(s in 1:3){
    for(q in (s+1):4){
      if(RVM2$family[q,s]==1){
        lo1[num]<--1
        up1[num]<-1
      }else if (RVM2$family[q,s] %in% c(3, 13)) {
        # clayton
        lo1[num] <- 1e-04
        up1[num] <- 100
      }else if (RVM2$family[q,s] %in% c(4, 14)) {
        # gumbel
        lo1[num] <- 1.0001
        up1[num] <- 100
      }else if (RVM2$family[q,s] %in% c(5)) {
        # gumbel
        lo1[num] <- -100
        up1[num] <- 100
      }else if (RVM2$family[q,s] %in% c(6, 16)) {
        # joe
        lo1[num] <- 1.0001
        up1[num] <- 50
      }else if (RVM2$family[q,s] %in% c(26, 36)) {
        # joe
        lo1[num] <- -50
        up1[num] <- -1.0001
      }
      num=num+1
    }
  }
  
  #MLE
  
  
  #start parameters (with selection)
  startpar1<-c(RVM2$par[2,1],RVM2$par[3,1],RVM2$par[4,1],RVM2$par[3,2],RVM2$par[4,2],RVM2$par[4,3],fit1$coefficients[1:3],sd(fit1$residuals),fit2$coefficients[1:3],sd(fit2$residuals),fit3$coefficients[1:3],sd(fit3$residuals),fit4$coefficients[1:3],sd(fit4$residuals))
  
  
  optimfun1<-function(param=startpar1){
    RVinepar3<-c(0,param[1:3],0,0,param[4:5],rep(0,3),param[6],rep(0,4))
    Rvinepar3<-matrix(RVinepar3,4,4)
    RVM3<-VineCopula::RVineMatrix(Rvineblock,RVM2$family,Rvinepar3,Rvinepar2,check.pars=FALSE)
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(r in 1:((k-1)*n)){
      logp1<-logp1+dnorm((ytrain1[r]-param[7]-param[8]*xtrain1[r]-param[9]*ztrain1[r]),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((ytrain2[r]-param[11]-param[12]*xtrain2[r]-param[13]*ztrain2[r]),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((ytrain3[r]-param[15]-param[16]*xtrain3[r]-param[17]*ztrain3[r]),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((ytrain4[r]-param[19]-param[20]*xtrain4[r]-param[21]*ztrain4[r]),mean=0,sd=param[22],log=TRUE)
    }
    
    u1<-pnorm((ytrain1-param[7]-param[8]*xtrain1-param[9]*ztrain1),mean=0,sd=param[10])
    u2<-pnorm((ytrain2-param[11]-param[12]*xtrain2-param[13]*ztrain2),mean=0,sd=param[14])
    u3<-pnorm((ytrain3-param[15]-param[16]*xtrain3-param[17]*ztrain3),mean=0,sd=param[18])
    u4<-pnorm((ytrain4-param[19]-param[20]*xtrain4-param[21]*ztrain4),mean=0,sd=param[22])
    
    simulate2<-cbind(u1,u2,u3,u4)
    
    
    -(logp1+logp2+logp3+logp4+VineCopula::RVineLogLik(simulate2,RVM3,check.pars = FALSE,calculate.V = FALSE)$loglik)
  }
  
  
  out1<-nlminb(objective=optimfun1,start=startpar1,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo1,upper=up1)
  
  #without selection(force to the right copula)
  startpar2<-c(Rvineblockpar[2,1],Rvineblockpar[3,1],Rvineblockpar[4,1],Rvineblockpar[3,2],Rvineblockpar[4,2],Rvineblockpar[4,3],fit1$coefficients[1:3],sd(fit1$residuals),fit2$coefficients[1:3],sd(fit2$residuals),fit3$coefficients[1:3],sd(fit3$residuals),fit4$coefficients[1:3],sd(fit4$residuals))
  
  
  optimfun2<-function(param=startpar1){
    RVinepar3<-c(0,param[1:3],0,0,param[4:5],rep(0,3),param[6],rep(0,4))
    Rvinepar3<-matrix(RVinepar3,4,4)
    RVM3<-VineCopula::RVineMatrix(Rvineblock,Rvineblockf,Rvinepar3,Rvinepar2,check.pars=FALSE)
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(r in 1:((k-1)*n)){
      logp1<-logp1+dnorm((ytrain1[r]-param[7]-param[8]*xtrain1[r]-param[9]*ztrain1[r]),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((ytrain2[r]-param[11]-param[12]*xtrain2[r]-param[13]*ztrain2[r]),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((ytrain3[r]-param[15]-param[16]*xtrain3[r]-param[17]*ztrain3[r]),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((ytrain4[r]-param[19]-param[20]*xtrain4[r]-param[21]*ztrain4[r]),mean=0,sd=param[22],log=TRUE)
    }
    
    u1<-pnorm((ytrain1-param[7]-param[8]*xtrain1-param[9]*ztrain1),mean=0,sd=param[10])
    u2<-pnorm((ytrain2-param[11]-param[12]*xtrain2-param[13]*ztrain2),mean=0,sd=param[14])
    u3<-pnorm((ytrain3-param[15]-param[16]*xtrain3-param[17]*ztrain3),mean=0,sd=param[18])
    u4<-pnorm((ytrain4-param[19]-param[20]*xtrain4-param[21]*ztrain4),mean=0,sd=param[22])
    
    simulate2<-cbind(u1,u2,u3,u4)
    
    
    -(logp1+logp2+logp3+logp4+VineCopula::RVineLogLik(simulate2,RVM3,check.pars = FALSE,calculate.V = FALSE)$loglik)
  }
  
  
  out2<-nlminb(objective=optimfun2,start=startpar2,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo,upper=up)
  
  #sequential(with selection)
  startpar3<-c(RVM2$par[2,1],RVM2$par[3,1],RVM2$par[4,1],RVM2$par[3,2],RVM2$par[4,2],RVM2$par[4,3])
  
  
  optimfun3<-function(param=startpar1){
    RVinepar3<-c(0,param[1:3],0,0,param[4:5],rep(0,3),param[6],rep(0,4))
    Rvinepar3<-matrix(RVinepar3,4,4)
    RVM3<-VineCopula::RVineMatrix(Rvineblock,RVM2$family,Rvinepar3,Rvinepar2,check.pars=FALSE)
    
    
    u1<-pnorm((ytrain1-fit1$coefficients[1]-fit1$coefficients[2]*xtrain1-fit1$coefficients[3]*ztrain1),mean=0,sd=sd(fit1$residuals))
    u2<-pnorm((ytrain2-fit2$coefficients[1]-fit2$coefficients[2]*xtrain2-fit2$coefficients[3]*ztrain2),mean=0,sd=sd(fit2$residuals))
    u3<-pnorm((ytrain3-fit3$coefficients[1]-fit3$coefficients[2]*xtrain3-fit3$coefficients[3]*ztrain3),mean=0,sd=sd(fit3$residuals))
    u4<-pnorm((ytrain4-fit4$coefficients[1]-fit4$coefficients[2]*xtrain4-fit4$coefficients[3]*ztrain4),mean=0,sd=sd(fit4$residuals))
    
    simulate2<-cbind(u1,u2,u3,u4)
    
    
    -(VineCopula::RVineLogLik(simulate2,RVM3,check.pars = FALSE,calculate.V = FALSE)$loglik)
  }
  
  
  out3<-nlminb(objective=optimfun3,start=startpar3,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo1[1:6],upper=up1[1:6])
  
  #sequential(without selection)
  startpar4<-c(Rvineblockpar[2,1],Rvineblockpar[3,1],Rvineblockpar[4,1],Rvineblockpar[3,2],Rvineblockpar[4,2],Rvineblockpar[4,3])
  
  
  optimfun4<-function(param=startpar1){
    RVinepar3<-c(0,param[1:3],0,0,param[4:5],rep(0,3),param[6],rep(0,4))
    Rvinepar3<-matrix(RVinepar3,4,4)
    RVM3<-VineCopula::RVineMatrix(Rvineblock,Rvineblockf,Rvinepar3,Rvinepar2,check.pars=FALSE)
    
    
    u1<-pnorm((ytrain1-fit1$coefficients[1]-fit1$coefficients[2]*xtrain1-fit1$coefficients[3]*ztrain1),mean=0,sd=sd(fit1$residuals))
    u2<-pnorm((ytrain2-fit2$coefficients[1]-fit2$coefficients[2]*xtrain2-fit2$coefficients[3]*ztrain2),mean=0,sd=sd(fit2$residuals))
    u3<-pnorm((ytrain3-fit3$coefficients[1]-fit3$coefficients[2]*xtrain3-fit3$coefficients[3]*ztrain3),mean=0,sd=sd(fit3$residuals))
    u4<-pnorm((ytrain4-fit4$coefficients[1]-fit4$coefficients[2]*xtrain4-fit4$coefficients[3]*ztrain4),mean=0,sd=sd(fit4$residuals))
    
    simulate2<-cbind(u1,u2,u3,u4)
    
    
    -(VineCopula::RVineLogLik(simulate2,RVM3,check.pars = FALSE,calculate.V = FALSE)$loglik)
  }
  
  
  out4<-nlminb(objective=optimfun4,start=startpar4,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo[1:6],upper=up[1:6])
  
  #prediction 
  prediction1<-rep(0,n*l)
  prediction2<-rep(0,n*l)
  prediction3<-rep(0,n*l)
  prediction4<-rep(0,n*l)
  #MLE NO selection prediction for MRM and CVine
  prediction5<-rep(0,n*l)
  prediction6<-rep(0,n*l)
  
  #Sequential Estimation Selection Prediction for MRM and CVine
  prediction7<-rep(0,n*l)
  prediction8<-rep(0,n*l)
  
  #sequential Estimation no selection prediction for CVine (MRM is same as prediction 7)
  prediction9<-rep(0,n*l)
  
  #CVM and MRM 
  #MRM predictions
  prediction2[1:n]=out1$par[7]+out1$par[8]*x1valid+out1$par[9]*z0valid
  prediction2[(n+1):(2*n)]=out1$par[11]+out1$par[12]*x2valid+out1$par[13]*z0valid
  prediction2[(2*n+1):(3*n)]=out1$par[15]+out1$par[16]*x3valid+out1$par[17]*z0valid
  prediction2[(3*n+1):(4*n)]=out1$par[19]+out1$par[20]*x4valid+out1$par[21]*z0valid
  
  prediction5[1:n]=out2$par[7]+out2$par[8]*x1valid+out2$par[9]*z0valid
  prediction5[(n+1):(2*n)]=out2$par[11]+out2$par[12]*x2valid+out2$par[13]*z0valid
  prediction5[(2*n+1):(3*n)]=out2$par[15]+out2$par[16]*x3valid+out2$par[17]*z0valid
  prediction5[(3*n+1):(4*n)]=out2$par[19]+out2$par[20]*x4valid+out2$par[21]*z0valid
  
  newdata1<-data.frame(x1train=x1valid,z0train=z0valid)
  prediction7[1:n]=predict(fit1,newdata1)
  newdata2<-data.frame(x2train=x2valid,z0train=z0valid)
  prediction7[(n+1):(2*n)]=predict(fit2,newdata2)
  newdata3<-data.frame(x3train=x3valid,z0train=z0valid)
  prediction7[(2*n+1):(3*n)]=predict(fit3,newdata3)
  newdata4<-data.frame(x4train=x4valid,z0train=z0valid)
  prediction7[(3*n+1):(4*n)]=predict(fit4,newdata4)
  
  #CVM predictions
  #for time point 1, it is just the marginal prediction without the C-VINE structure
  prediction1[1:n]=prediction2[1:n]
  prediction6[1:n]=prediction5[1:n]
  prediction8[1:n]=prediction7[1:n]
  prediction9[1:n]=prediction7[1:n]
  
  #sub-parameters
  RVinepar4<-c(0,out1$par[1:3],0,0,out1$par[4:5],rep(0,3),out1$par[6],rep(0,4))
  Rvinepar4<-matrix(RVinepar4,4,4)
  
  RVinepar5<-c(0,out2$par[1:3],0,0,out2$par[4:5],rep(0,3),out2$par[6],rep(0,4))
  Rvinepar5<-matrix(RVinepar5,4,4)
  
  RVinepar6<-c(0,out3$par[1:3],0,0,out3$par[4:5],rep(0,3),out3$par[6],rep(0,4))
  Rvinepar6<-matrix(RVinepar6,4,4)
  
  RVinepar7<-c(0,out4$par[1:3],0,0,out4$par[4:5],rep(0,3),out4$par[6],rep(0,4))
  Rvinepar7<-matrix(RVinepar7,4,4)
  
  #for time point 2
  #MLE Selection
  RVM4<-VineCopula::RVineMatrix(RVM2$Matrix[3:4,3:4], family = RVM2$family[3:4,3:4], par = Rvinepar4[3:4,3:4], par2 = RVM2$par2[3:4,3:4],check.pars = TRUE)
  RVM5<-VineCopula::RVineMatrix(RVM2$Matrix[2:4,2:4], family = RVM2$family[2:4,2:4], par = Rvinepar4[2:4,2:4], par2 = RVM2$par2[2:4,2:4],check.pars = TRUE)
  RVM6<-VineCopula::RVineMatrix(RVM2$Matrix, family = RVM2$family, par = Rvinepar4,check.pars = TRUE)
  
  #MLE NO selection
  RVM7<-VineCopula::RVineMatrix(RVM2$Matrix[3:4,3:4], family = Rvineblockf[3:4,3:4], par = Rvinepar5[3:4,3:4], par2 = Rvinepar2[3:4,3:4],check.pars = TRUE)
  RVM8<-VineCopula::RVineMatrix(RVM2$Matrix[2:4,2:4], family = Rvineblockf[2:4,2:4], par = Rvinepar5[2:4,2:4], par2 = Rvinepar2[2:4,2:4],check.pars = TRUE)
  RVM9<-VineCopula::RVineMatrix(RVM2$Matrix, family = Rvineblockf, par = Rvinepar5,check.pars = TRUE)
  
  #Sequential with selection
  RVM10<-VineCopula::RVineMatrix(RVM2$Matrix[3:4,3:4], family = RVM2$family[3:4,3:4], par = Rvinepar6[3:4,3:4], par2 = Rvinepar2[3:4,3:4],check.pars = TRUE)
  RVM11<-VineCopula::RVineMatrix(RVM2$Matrix[2:4,2:4], family = RVM2$family[2:4,2:4], par = Rvinepar6[2:4,2:4], par2 = Rvinepar2[2:4,2:4],check.pars = TRUE)
  RVM12<-VineCopula::RVineMatrix(RVM2$Matrix, family = RVM2$family, par = Rvinepar6,check.pars = TRUE)
  
  #Sequential NO selection
  RVM13<-VineCopula::RVineMatrix(RVM2$Matrix[3:4,3:4], family = Rvineblockf[3:4,3:4], par = Rvinepar7[3:4,3:4], par2 = Rvinepar2[3:4,3:4],check.pars = TRUE)
  RVM14<-VineCopula::RVineMatrix(RVM2$Matrix[2:4,2:4], family = Rvineblockf[2:4,2:4], par = Rvinepar7[2:4,2:4], par2 = Rvinepar2[2:4,2:4],check.pars = TRUE)
  RVM15<-VineCopula::RVineMatrix(RVM2$Matrix, family = Rvineblockf, par = Rvinepar7,check.pars = TRUE)
  
  
  resi1<-y1valid-prediction2[1:n]
  resi2<-y2valid-prediction2[(n+1):(2*n)]
  resi3<-y3valid-prediction2[(2*n+1):(3*n)]
  
  resi1a<-y1valid-prediction5[1:n]
  resi2a<-y2valid-prediction5[(n+1):(2*n)]
  resi3a<-y3valid-prediction5[(2*n+1):(3*n)]
  
  resi1b<-y1valid-prediction7[1:n]
  resi2b<-y2valid-prediction7[(n+1):(2*n)]
  resi3b<-y3valid-prediction7[(2*n+1):(3*n)]
  
  for(o in 1:n){
    a<-seq(-30,30,0.01)
    
    d<-0
    da<-0
    db<-0
    dc<-0
    
    c2<-function(k){
      cop<-pnorm(k,0,out1$par[14])
      dat<-c(pnorm(resi1[o],0,out1$par[10]),cop)
      VineCopula::RVinePDF(dat,RVM4,verbose=TRUE)
    }
    
    c2a<-function(k){
      cop<-pnorm(k,0,out2$par[14])
      dat<-c(pnorm(resi1a[o],0,out2$par[10]),cop)
      VineCopula::RVinePDF(dat,RVM7,verbose=TRUE)
    }
    
    c2b<-function(k){
      cop<-pnorm(k,0,sd(fit2$residuals))
      dat<-c(pnorm(resi1b[o],0,sd(fit1$residuals)),cop)
      VineCopula::RVinePDF(dat,RVM10,verbose=TRUE)
    }
    
    c2c<-function(k){
      cop<-pnorm(k,0,sd(fit2$residuals))
      dat<-c(pnorm(resi1b[o],0,sd(fit1$residuals)),cop)
      VineCopula::RVinePDF(dat,RVM13,verbose=TRUE)
    }
    
    f1<-function(k){
      dnorm(k,0,out1$par[14])
    }
    
    f1a<-function(k){
      dnorm(k,0,out2$par[14])
    }
    
    f1b<-function(k){
      dnorm(k,0,sd(fit2$residuals))
    }
    
    f1c<-function(k){
      dnorm(k,0,sd(fit2$residuals))
    }
    
    for(i in 1:6001){
      d<-d+a[i]*c2(a[i])*f1(a[i])*0.01
      
      da<-da+a[i]*c2a(a[i])*f1a(a[i])*0.01
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01
      
      dc<-dc+a[i]*c2c(a[i])*f1c(a[i])*0.01
      
    }
    
    prediction1[(n+o)]<- d+prediction2[(n+o)]
    prediction6[(n+o)]<- da+prediction5[(n+o)]
    prediction8[(n+o)]<- db+prediction7[(n+o)]
    prediction9[(n+o)]<- dc+prediction7[(n+o)]
  }
  
  
  #for time point 3
  
  c1<-rep(0,n)
  c1a<-rep(0,n)
  c1b<-rep(0,n)
  c1c<-rep(0,n)
  
  for(o in 1:n){
    a<-seq(-30,30,0.01)
    
    d<-0
    da<-0
    db<-0
    dc<-0
    
    c1[o]<-VineCopula::RVinePDF(c(pnorm(resi1[o],0,out1$par[10]),pnorm(resi2[o],0,out1$par[14])),RVM4, verbose = TRUE)
    c1a[o]<-VineCopula::RVinePDF(c(pnorm(resi1a[o],0,out2$par[10]),pnorm(resi2a[o],0,out2$par[14])),RVM7, verbose = TRUE)
    c1b[o]<-VineCopula::RVinePDF(c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals))),RVM10, verbose = TRUE)
    c1c[o]<-VineCopula::RVinePDF(c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals))),RVM13, verbose = TRUE)
    
    c2<-function(k){
      cop<-pnorm(k,0,out1$par[18])
      dat<-c(pnorm(resi1[o],0,out1$par[10]),pnorm(resi2[o],0,out1$par[14]),cop)
      VineCopula::RVinePDF(dat,RVM5,verbose=TRUE)
    }
    
    c2a<-function(k){
      cop<-pnorm(k,0,out2$par[18])
      dat<-c(pnorm(resi1a[o],0,out2$par[10]),pnorm(resi2a[o],0,out2$par[14]),cop)
      VineCopula::RVinePDF(dat,RVM8,verbose=TRUE)
    }
    
    c2b<-function(k){
      cop<-pnorm(k,0,sd(fit3$residuals))
      dat<-c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals)),cop)
      VineCopula::RVinePDF(dat,RVM11,verbose=TRUE)
    }
    
    c2c<-function(k){
      cop<-pnorm(k,0,sd(fit3$residuals))
      dat<-c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals)),cop)
      VineCopula::RVinePDF(dat,RVM14,verbose=TRUE)
    }
    
    
    
    f1<-function(k){
      dnorm(k,0,out1$par[18])
    }
    
    f1a<-function(k){
      dnorm(k,0,out2$par[18])
    }
    
    f1b<-function(k){
      dnorm(k,0,sd(fit3$residuals))
    }
    
    f1c<-function(k){
      dnorm(k,0,sd(fit3$residuals))
    }
    
    for(i in 1:6001){
      d<-d+a[i]*c2(a[i])*f1(a[i])*0.01/c1[o]
      
      da<-da+a[i]*c2a(a[i])*f1a(a[i])*0.01/c1a[o]
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01/c1b[o]
      
      dc<-dc+a[i]*c2c(a[i])*f1c(a[i])*0.01/c1c[o]
    }
    
    prediction1[2*n+o]<- d+prediction2[2*n+o]
    prediction6[2*n+o]<- da+prediction5[2*n+o]
    prediction8[2*n+o]<- db+prediction7[2*n+o]
    prediction9[2*n+o]<- dc+prediction7[2*n+o]
  }
  
  
  
  #for time point 4
  
  c1<-rep(0,n)
  c1a<-rep(0,n)
  c1b<-rep(0,n)
  c1c<-rep(0,n)
  
  for(o in 1:n){
    a<-seq(-30,30,0.01)
    b<-rep(0,6001)
    c<-rep(0,6001)
    
    d<-0
    da<-0
    db<-0
    dc<-0
    
    c1[o]<-VineCopula::RVinePDF(c(pnorm(resi1[o],0,out1$par[10]),pnorm(resi2[o],0,out1$par[14]),pnorm(resi3[o],0,out1$par[18])),RVM5, verbose = TRUE)
    
    c1a[o]<-VineCopula::RVinePDF(c(pnorm(resi1a[o],0,out2$par[10]),pnorm(resi2a[o],0,out2$par[14]),pnorm(resi3a[o],0,out2$par[18])),RVM8, verbose = TRUE)
    c1b[o]<-VineCopula::RVinePDF(c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals)),pnorm(resi3b[o],0,sd(fit3$residuals))),RVM11, verbose = TRUE)
    c1c[o]<-VineCopula::RVinePDF(c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals)),pnorm(resi3b[o],0,sd(fit3$residuals))),RVM14, verbose = TRUE)
    
    c2<-function(k){
      cop<-pnorm(k,0,out1$par[22])
      dat<-c(pnorm(resi1[o],0,out1$par[10]),pnorm(resi2[o],0,out1$par[14]),pnorm(resi3[o],0,out1$par[18]),cop)
      VineCopula::RVinePDF(dat,RVM6,verbose=TRUE)
    }
    
    c2a<-function(k){
      cop<-pnorm(k,0,out2$par[22])
      dat<-c(pnorm(resi1a[o],0,out2$par[10]),pnorm(resi2a[o],0,out2$par[14]),pnorm(resi3a[o],0,out2$par[18]),cop)
      VineCopula::RVinePDF(dat,RVM9,verbose=TRUE)
    }
    
    c2b<-function(k){
      cop<-pnorm(k,0,sd(fit3$residuals))
      dat<-c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals)),pnorm(resi3b[o],0,sd(fit3$residuals)),cop)
      VineCopula::RVinePDF(dat,RVM12,verbose=TRUE)
    }
    
    c2c<-function(k){
      cop<-pnorm(k,0,sd(fit3$residuals))
      dat<-c(pnorm(resi1b[o],0,sd(fit1$residuals)),pnorm(resi2b[o],0,sd(fit2$residuals)),pnorm(resi3b[o],0,sd(fit3$residuals)),cop)
      VineCopula::RVinePDF(dat,RVM15,verbose=TRUE)
    }
    
    
    f1<-function(k){
      dnorm(k,0,out1$par[22])
    }
    
    f1a<-function(k){
      dnorm(k,0,out2$par[22])
    }
    
    f1b<-function(k){
      dnorm(k,0,sd(fit4$residuals))
    }
    
    f1c<-function(k){
      dnorm(k,0,sd(fit4$residuals))
    }
    
    for(i in 1:6001){
      d<-d+a[i]*c2(a[i])*f1(a[i])*0.01/c1[o]
      
      da<-da+a[i]*c2a(a[i])*f1a(a[i])*0.01/c1a[o]
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01/c1b[o]
      
      dc<-dc+a[i]*c2c(a[i])*f1c(a[i])*0.01/c1c[o]
    }
    
    prediction1[3*n+o]<- d+prediction2[3*n+o]
    prediction6[3*n+o]<- da+prediction5[3*n+o]
    prediction8[3*n+o]<- db+prediction7[3*n+o]
    prediction9[3*n+o]<- dc+prediction7[3*n+o]
  }
  
  #fit independent regression model with seperating time block and time points
  newdata5<-data.frame(xtrain=xvalid,ztrain=zvalid,ttrain=tvalid)
  prediction3<-predict(fit,newdata5)
  
  #fit AR model and prediction
  newdata6<-data.frame(xtrain=xvalid,jtrain=jvalid)
  prediction4<-predict(fit0,newdata6)
  
  MAE_CVINE[m]<-sum(abs(yvalidadjusted-prediction1))/length(yvalid)
  MAE_MRM[m]<-sum(abs(yvalidadjusted-prediction2))/length(yvalid)
  MAE_LRM[m]<-sum(abs(yvalid-prediction3))/length(yvalid)
  MAE_AR[m]<-sum(abs(yvalid-prediction4))/length(yvalid)
  MAE_CVINE1[m]<-sum(abs(yvalidadjusted-prediction6))/length(yvalid)
  MAE_MRM1[m]<-sum(abs(yvalidadjusted-prediction5))/length(yvalid)
  MAE_CVINE2[m]<-sum(abs(yvalidadjusted-prediction8))/length(yvalid)
  MAE_CVINE3[m]<-sum(abs(yvalidadjusted-prediction9))/length(yvalid)
  MAE_MRM2[m]<-sum(abs(yvalidadjusted-prediction7))/length(yvalid)
  
  count=0
  for(i in 1:length(prediction1)){
    if(abs(yvalidadjusted-prediction1)[i]<=abs(yvalidadjusted-prediction2)[i]){
      count=count+1
    }
  }
  percentage1[m]<-count/length(prediction1)
  
  count=0
  for(i in 1:length(prediction1)){
    if(abs(yvalidadjusted-prediction1)[i]<=abs(yvalid-prediction3)[i]){
      count=count+1
    }
  }
  percentage2[m]<-count/length(prediction1)
  
  count=0
  for(i in 1:length(prediction1)){
    if(abs(yvalidadjusted-prediction1)[i]<=abs(yvalid-prediction4)[i]){
      count=count+1
    }
  }
  percentage3[m]<-count/length(prediction1)
}






MAE_CVINE
MAE_MRM
MAE_LRM
MAE_AR
MAE_CVINE1
MAE_MRM1
MAE_CVINE2
MAE_CVINE3
MAE_MRM2
percentage1
percentage2
percentage3