

library(VineCopula)
library(forecast)
library(nlme)

#parameter setting for 4 time points within a time block(Marginal Parameters are set to be identical). Other than this, other parts are identical to prediction scenario 1
beta01<-2.5
beta02<-2.5
beta03<-2.5
beta04<-2.5
beta11<-3.5
beta12<-3.5
beta13<-3.5
beta14<-3.5
beta21<-4.5
beta22<-4.5
beta23<-4.5
beta24<-4.5
sigma<-2

#generate random numbers from a specific R-Vine Structure

#20 time points structure
RVinem1<-c(20,4,3,2,1,8,7,6,5,12,11,10,9,16,15,14,13,19,18,17,0,4,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,3,2,1,rep(0,2),19,3,2,1,8,7,6,5,12,11,10,9,16,15,14,13,18,17,rep(0,3),3,18,17,16,15,14,13,12,11,10,9,8,7,6,5,2,1,rep(0,4),18,2,1,8,7,6,5,12,11,10,9,16,15,14,13,17,rep(0,5),2,17,16,15,14,13,12,11,10,9,8,7,6,5,1,rep(0,6),17,1,8,7,6,5,12,11,10,9,16,15,14,13,rep(0,7),1,16,15,14,13,12,11,10,9,8,7,6,5,rep(0,8),16,8,7,6,5,12,11,10,9,15,14,13,rep(0,9),8,15,14,13,12,11,10,9,7,6,5,rep(0,10),15,7,6,5,12,11,10,9,14,13,rep(0,11),7,14,13,12,11,10,9,6,5,rep(0,12),14,6,5,12,11,10,9,13,rep(0,13),6,13,12,11,10,9,5,rep(0,14),13,5,12,11,10,9,rep(0,15),5,12,11,10,9,rep(0,16),12,11,10,9,rep(0,17),11,10,9,rep(0,18),10,9,rep(0,19),9)
Rvinem1<-matrix(RVinem1,20,20)


##Rvine family and parameters
RVinef<-c(rep(0,17),3,4,1,rep(0,17),3,4,1,rep(0,17),1,3,4,rep(0,17),1,3,4,rep(0,17),1,1,5,rep(0,17),1,1,5,rep(0,17),1,1,1,rep(0,17),1,1,1,rep(0,17),3,4,1,rep(0,17),3,4,1,rep(0,17),1,3,4,rep(0,17),1,3,4,rep(0,17),1,1,5,rep(0,17),1,1,5,rep(0,17),1,1,1,rep(0,17),1,1,1,rep(0,17),3,4,1,rep(0,18),3,4,rep(0,19),5,rep(0,20))
Rvinef<-matrix(RVinef,20,20)
RVinepar1<-c(rep(0,17),2,3,0.5,rep(0,17),2,3,0.5,rep(0,17),0.5,3,2,rep(0,17),0.5,3,2,rep(0,17),0.5,0.5,2,rep(0,17),0.5,0.5,2,rep(0,17),0.5,0.5,0.5,rep(0,17),0.5,0.5,0.5,rep(0,17),2,3,0.5,rep(0,17),2,3,0.5,rep(0,17),0.5,3,2,rep(0,17),0.5,3,2,rep(0,17),0.5,0.5,2,rep(0,17),0.5,0.5,2,rep(0,17),0.5,0.5,0.5,rep(0,17),0.5,0.5,0.5,rep(0,17),2,3,0.5,rep(0,18),3,2,rep(0,19),2,rep(0,20))
Rvinepar1<-matrix(RVinepar1,20,20)

RVinepar2<-rep(0,400)
Rvinepar2<-matrix(RVinepar2,20,20)

RVM1<-RVineMatrix(Rvinem1,Rvinef,Rvinepar1,Rvinepar2,check.pars=TRUE)

#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)


#family and parameters within a time block 
RVineblockf<-c(0,3,4,1,0,0,3,4,0,0,0,5,0,0,0,0)
Rvineblockf<-matrix(RVineblockf,4,4)

RVineblockpar<-c(0,2,3,0.5,0,0,3,2,0,0,0,2,0,0,0,0)
Rvineblockpar<-matrix(RVineblockpar,4,4)

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

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

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

#the proportion of training set 
p=0.9


##Case1
#Create response data and covariate data
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)
for(m in 1:200){
  
  set.seed(m)
  simulate1<-RVineSim(n,RVM1)
  epsilon1<-rep(0,n*k)
  epsilon2<-rep(0,n*k)
  epsilon3<-rep(0,n*k)
  epsilon4<-rep(0,n*k)
  
  for (j in 1:k){
    for(i in 1:n){
      epsilon1[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*l+1],0,sigma)
      epsilon2[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*l+2],0,sigma)
      epsilon3[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*l+3],0,sigma)
      epsilon4[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*l+4],0,sigma)
    }
  }
  
  
  #generate covariates x
  x1<-runif(n*k,0,100)
  x2<-runif(n*k,0,100)
  x3<-runif(n*k,0,100)
  x4<-runif(n*k,0,100)
  
  #calculate corresponding response y
  y1<-rep(0,n*k)
  y2<-rep(0,n*k)
  y3<-rep(0,n*k)
  y4<-rep(0,n*k)
  
  z<-c(rep(1,n),rep(2,n),rep(3,n),rep(4,n))
  
  for(j in 1:k){
    y1[((j-1)*n+1):(j*n)]<-beta01+beta11*x1[((j-1)*n+1):(j*n)]+beta21*j+epsilon1[((j-1)*n+1):(j*n)]
    y2[((j-1)*n+1):(j*n)]<-beta02+beta12*x2[((j-1)*n+1):(j*n)]+beta22*j+epsilon2[((j-1)*n+1):(j*n)]
    y3[((j-1)*n+1):(j*n)]<-beta03+beta13*x3[((j-1)*n+1):(j*n)]+beta23*j+epsilon3[((j-1)*n+1):(j*n)]
    y4[((j-1)*n+1):(j*n)]<-beta04+beta14*x4[((j-1)*n+1):(j*n)]+beta24*j+epsilon4[((j-1)*n+1):(j*n)]
  }
  
  
  #use the first 450 data points as training set and the last 50 data points as the validation set
  ytrain<-c(y1[1:(n*p)],y1[(n+1):(n*(1+p))],y1[(2*n+1):(n*(2+p))],y1[(3*n+1):(n*(3+p))],y2[1:(n*p)],y2[(n+1):(n*(1+p))],y2[(2*n+1):(n*(2+p))],y2[(3*n+1):(n*(3+p))],y3[1:(n*p)],y3[(n+1):(n*(1+p))],y3[(2*n+1):(n*(2+p))],y3[(3*n+1):(n*(3+p))],y4[1:(n*p)],y4[(n+1):(n*(1+p))],y4[(2*n+1):(n*(2+p))],y4[(3*n+1):(n*(3+p))])
  xtrain<-c(x1[1:(n*p)],x1[(n+1):(n*(1+p))],x1[(2*n+1):(n*(2+p))],x1[(3*n+1):(n*(3+p))],x2[1:(n*p)],x2[(n+1):(n*(1+p))],x2[(2*n+1):(n*(2+p))],x2[(3*n+1):(n*(3+p))],x3[1:(n*p)],x3[(n+1):(n*(1+p))],x3[(2*n+1):(n*(2+p))],x3[(3*n+1):(n*(3+p))],x4[1:(n*p)],x4[(n+1):(n*(1+p))],x4[(2*n+1):(n*(2+p))],x4[(3*n+1):(n*(3+p))])
  ztrain<-c(rep(1,n*p),rep(2,n*p),rep(3,n*p),rep(4,n*p),rep(1,n*p),rep(2,n*p),rep(3,n*p),rep(4,n*p),rep(1,n*p),rep(2,n*p),rep(3,n*p),rep(4,n*p),rep(1,n*p),rep(2,n*p),rep(3,n*p),rep(4,n*p))
  
  ytrain1<-c(y1[1:(n*p)],y1[(n+1):(n*(1+p))],y1[(2*n+1):(n*(2+p))],y1[(3*n+1):(n*(3+p))])
  ytrain2<-c(y2[1:(n*p)],y2[(n+1):(n*(1+p))],y2[(2*n+1):(n*(2+p))],y2[(3*n+1):(n*(3+p))])
  ytrain3<-c(y3[1:(n*p)],y3[(n+1):(n*(1+p))],y3[(2*n+1):(n*(2+p))],y3[(3*n+1):(n*(3+p))])
  ytrain4<-c(y4[1:(n*p)],y4[(n+1):(n*(1+p))],y4[(2*n+1):(n*(2+p))],y4[(3*n+1):(n*(3+p))])
  xtrain1<-c(x1[1:(n*p)],x1[(n+1):(n*(1+p))],x1[(2*n+1):(n*(2+p))],x1[(3*n+1):(n*(3+p))])
  xtrain2<-c(x2[1:(n*p)],x2[(n+1):(n*(1+p))],x2[(2*n+1):(n*(2+p))],x2[(3*n+1):(n*(3+p))])
  xtrain3<-c(x3[1:(n*p)],x3[(n+1):(n*(1+p))],x3[(2*n+1):(n*(2+p))],x3[(3*n+1):(n*(3+p))])
  xtrain4<-c(x4[1:(n*p)],x4[(n+1):(n*(1+p))],x4[(2*n+1):(n*(2+p))],x4[(3*n+1):(n*(3+p))])
  ztrain0<-c(rep(1,n*p),rep(2,n*p),rep(3,n*p),rep(4,n*p))
  
  yvalid<-c(y1[(n*p+1):(n)],y1[(n*(1+p)+1):(2*n)],y1[(n*(2+p)+1):(n*3)],y1[(n*(3+p)+1):(n*4)],y2[(n*p+1):(n)],y2[(n*(1+p)+1):(2*n)],y2[(n*(2+p)+1):(n*3)],y2[(n*(3+p)+1):(n*4)],y3[(n*p+1):(n)],y3[(n*(1+p)+1):(2*n)],y3[(n*(2+p)+1):(n*3)],y3[(n*(3+p)+1):(n*4)],y4[(n*p+1):(n)],y4[(n*(1+p)+1):(2*n)],y4[(n*(2+p)+1):(n*3)],y4[(n*(3+p)+1):(n*4)])
  xvalid<-c(x1[(n*p+1):(n)],x1[(n*(1+p)+1):(2*n)],x1[(n*(2+p)+1):(n*3)],x1[(n*(3+p)+1):(n*4)],x2[(n*p+1):(n)],x2[(n*(1+p)+1):(2*n)],x2[(n*(2+p)+1):(n*3)],x2[(n*(3+p)+1):(n*4)],x3[(n*p+1):(n)],x3[(n*(1+p)+1):(2*n)],x3[(n*(2+p)+1):(n*3)],x3[(n*(3+p)+1):(n*4)],x4[(n*p+1):(n)],x4[(n*(1+p)+1):(2*n)],x4[(n*(2+p)+1):(n*3)],x4[(n*(3+p)+1):(n*4)])
  zvalid<-rep(c(rep(1,n*(1-p)+1),rep(2,n*(1-p)+1),rep(3,n*(1-p)+1),rep(4,n*(1-p)+1)),4)
  
  yvalid1<-c(y1[(n*p+1):(n)],y1[(n*(1+p)+1):(2*n)],y1[(n*(2+p)+1):(n*3)],y1[(n*(3+p)+1):(n*4)])
  yvalid2<-c(y2[(n*p+1):(n)],y2[(n*(1+p)+1):(2*n)],y2[(n*(2+p)+1):(n*3)],y2[(n*(3+p)+1):(n*4)])
  yvalid3<-c(y3[(n*p+1):(n)],y3[(n*(1+p)+1):(2*n)],y3[(n*(2+p)+1):(n*3)],y3[(n*(3+p)+1):(n*4)])
  yvalid4<-c(y4[(n*p+1):(n)],y4[(n*(1+p)+1):(2*n)],y4[(n*(2+p)+1):(n*3)],y4[(n*(3+p)+1):(n*4)])
  xvalid1<-c(x1[(n*p+1):(n)],x1[(n*(1+p)+1):(2*n)],x1[(n*(2+p)+1):(n*3)],x1[(n*(3+p)+1):(n*4)])
  xvalid2<-c(x2[(n*p+1):(n)],x2[(n*(1+p)+1):(2*n)],x2[(n*(2+p)+1):(n*3)],x2[(n*(3+p)+1):(n*4)])
  xvalid3<-c(x3[(n*p+1):(n)],x3[(n*(1+p)+1):(2*n)],x3[(n*(2+p)+1):(n*3)],x3[(n*(3+p)+1):(n*4)])
  xvalid4<-c(x4[(n*p+1):(n)],x4[(n*(1+p)+1):(2*n)],x4[(n*(2+p)+1):(n*3)],x4[(n*(3+p)+1):(n*4)])
  zvalid0<-c(rep(1,n*(1-p)+1),rep(2,n*(1-p)+1),rep(3,n*(1-p)+1),rep(4,n*(1-p)+1))
  
  #fit CVM model
  #Marginal Regression
  fit1<-lm(ytrain1~xtrain1+ztrain0)
  
  fit2<-lm(ytrain2~xtrain2+ztrain0)
  
  fit3<-lm(ytrain3~xtrain3+ztrain0)
  
  fit4<-lm(ytrain4~xtrain4+ztrain0)
  
  
  
  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<-VineCopula::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,check.pars=FALSE)
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(r in 1:(n*p*k)){
      logp1<-logp1+dnorm((ytrain1[r]-param[7]-param[8]*xtrain1[r]-param[9]*ztrain0[r]),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((ytrain2[r]-param[11]-param[12]*xtrain2[r]-param[13]*ztrain0[r]),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((ytrain3[r]-param[15]-param[16]*xtrain3[r]-param[17]*ztrain0[r]),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((ytrain4[r]-param[19]-param[20]*xtrain4[r]-param[21]*ztrain0[r]),mean=0,sd=param[22],log=TRUE)
    }
    
    u1<-pnorm((ytrain1-param[7]-param[8]*xtrain1-param[9]*ztrain0),mean=0,sd=param[10])
    u2<-pnorm((ytrain2-param[11]-param[12]*xtrain2-param[13]*ztrain0),mean=0,sd=param[14])
    u3<-pnorm((ytrain3-param[15]-param[16]*xtrain3-param[17]*ztrain0),mean=0,sd=param[18])
    u4<-pnorm((ytrain4-param[19]-param[20]*xtrain4-param[21]*ztrain0),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,check.pars=FALSE)
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(r in 1:(n*p*k)){
      logp1<-logp1+dnorm((ytrain1[r]-param[7]-param[8]*xtrain1[r]-param[9]*ztrain0[r]),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((ytrain2[r]-param[11]-param[12]*xtrain2[r]-param[13]*ztrain0[r]),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((ytrain3[r]-param[15]-param[16]*xtrain3[r]-param[17]*ztrain0[r]),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((ytrain4[r]-param[19]-param[20]*xtrain4[r]-param[21]*ztrain0[r]),mean=0,sd=param[22],log=TRUE)
    }
    
    u1<-pnorm((ytrain1-param[7]-param[8]*xtrain1-param[9]*ztrain0),mean=0,sd=param[10])
    u2<-pnorm((ytrain2-param[11]-param[12]*xtrain2-param[13]*ztrain0),mean=0,sd=param[14])
    u3<-pnorm((ytrain3-param[15]-param[16]*xtrain3-param[17]*ztrain0),mean=0,sd=param[18])
    u4<-pnorm((ytrain4-param[19]-param[20]*xtrain4-param[21]*ztrain0),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,check.pars=FALSE)
    
    
    u1<-pnorm((ytrain1-fit1$coefficients[1]-fit1$coefficients[2]*xtrain1-fit1$coefficients[3]*ztrain0),mean=0,sd=sd(fit1$residuals))
    u2<-pnorm((ytrain2-fit2$coefficients[1]-fit2$coefficients[2]*xtrain2-fit2$coefficients[3]*ztrain0),mean=0,sd=sd(fit2$residuals))
    u3<-pnorm((ytrain3-fit3$coefficients[1]-fit3$coefficients[2]*xtrain3-fit3$coefficients[3]*ztrain0),mean=0,sd=sd(fit3$residuals))
    u4<-pnorm((ytrain4-fit4$coefficients[1]-fit4$coefficients[2]*xtrain4-fit4$coefficients[3]*ztrain0),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,check.pars=FALSE)
    
    
    u1<-pnorm((ytrain1-fit1$coefficients[1]-fit1$coefficients[2]*xtrain1-fit1$coefficients[3]*ztrain0),mean=0,sd=sd(fit1$residuals))
    u2<-pnorm((ytrain2-fit2$coefficients[1]-fit2$coefficients[2]*xtrain2-fit2$coefficients[3]*ztrain0),mean=0,sd=sd(fit2$residuals))
    u3<-pnorm((ytrain3-fit3$coefficients[1]-fit3$coefficients[2]*xtrain3-fit3$coefficients[3]*ztrain0),mean=0,sd=sd(fit3$residuals))
    u4<-pnorm((ytrain4-fit4$coefficients[1]-fit4$coefficients[2]*xtrain4-fit4$coefficients[3]*ztrain0),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])
  
  
  
  #fit a regression model without seperation of time points
  ktrain<-c(rep(1,n*p*k),rep(2,n*p*k),rep(3,n*p*k),rep(4,n*p*k))
  fit<-lm(ytrain~xtrain+ztrain+ktrain)
  kvalid<-c(rep(1,n*(1-p)*k+1),rep(2,n*(1-p)*k+1),rep(3,n*(1-p)*k+1),rep(4,n*(1-p)*k+1))
  
  #fit a ARMA model
  
  jtrain=NULL
  
  for(i in 1:k){
    for(j in 1:l){
      jtrain<-c(jtrain,rep((j-1)*l+i,n*p))
    }
  }
  
  jvalid=NULL
  for(i in 1:k){
    for(j in 1:l){
      jvalid<-c(jvalid,rep((j-1)*l+i,n*(1-p)+1))
    }
  }
  
  sin<-sin(jtrain/4*2*pi)
  cos<-cos(jtrain/4*2*pi)
  fit0<-gls(ytrain~xtrain+ceiling(jtrain/4),correlation = corARMA(form=~1|jtrain,p=2,q=0))
  
  #prediction 
  prediction1<-rep(0,n*(1-p)*k*l+1)
  prediction2<-rep(0,n*(1-p)*k*l+1)
  prediction3<-rep(0,n*(1-p)*k*l+1)
  prediction4<-rep(0,n*(1-p)*k*l+1)
  #MLE NO selection prediction for MRM and CVine
  prediction5<-rep(0,n*(1-p)*k*l+1)
  prediction6<-rep(0,n*(1-p)*k*l+1)
  
  #Sequential Estimation Selection Prediction for MRM and CVine
  prediction7<-rep(0,n*(1-p)*k*l+1)
  prediction8<-rep(0,n*(1-p)*k*l+1)
  
  #sequential Estimation no selection prediction for CVine (MRM is same as prediction 7)
  prediction9<-rep(0,n*(1-p)*k*l+1)
  
  
  #CVM and MRM 
  #MRM predictions
  n1<-n*(1-p)*k
  n2<-2*n*(1-p)*k
  n3<-3*n*(1-p)*k
  n4<-4*n*(1-p)*k
  
  prediction2[1:n1]=out1$par[7]+out1$par[8]*xvalid1+out1$par[9]*zvalid0
  prediction2[(n1+2):(n2+1)]=out1$par[11]+out1$par[12]*xvalid2+out1$par[13]*zvalid0
  prediction2[(n2+2):(n3+1)]=out1$par[15]+out1$par[16]*xvalid3+out1$par[17]*zvalid0
  prediction2[(n3+2):(n4+1)]=out1$par[19]+out1$par[20]*xvalid4+out1$par[21]*zvalid0
  
  prediction5[1:n1]=out2$par[7]+out2$par[8]*xvalid1+out2$par[9]*zvalid0
  prediction5[(n1+2):(n2+1)]=out2$par[11]+out2$par[12]*xvalid2+out2$par[13]*zvalid0
  prediction5[(n2+2):(n3+1)]=out2$par[15]+out2$par[16]*xvalid3+out2$par[17]*zvalid0
  prediction5[(n3+2):(n4+1)]=out2$par[19]+out2$par[20]*xvalid4+out2$par[21]*zvalid0
  
  newdata1<-data.frame(xtrain1=xvalid1,ztrain0=zvalid0)
  prediction7[1:n1]=predict(fit1,newdata1)
  newdata2<-data.frame(xtrain2=xvalid2,ztrain0=zvalid0)
  prediction7[(n1+2):(n2+1)]=predict(fit2,newdata2)
  newdata3<-data.frame(xtrain3=xvalid3,ztrain0=zvalid0)
  prediction7[(n2+2):(n3+1)]=predict(fit3,newdata3)
  newdata4<-data.frame(xtrain4=xvalid4,ztrain0=zvalid0)
  prediction7[(n3+2):(n4+1)]=predict(fit4,newdata4)
  
  #CVM predictions
  #for time point 1, it is just the marginal prediction without the C-VINE structure
  prediction1[1:(n*(1-p)*k)]=prediction2[1:(n*(1-p)*k)]
  prediction6[1:(n*(1-p)*k)]=prediction5[1:(n*(1-p)*k)]
  prediction8[1:(n*(1-p)*k)]=prediction7[1:(n*(1-p)*k)]
  prediction9[1:(n*(1-p)*k)]=prediction7[1:(n*(1-p)*k)]
  
  #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, par2 = RVM2$par2,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<-yvalid1-prediction2[1:(n*(1-p)*k)]
  resi2<-yvalid2-prediction2[(n*(1-p)*k+2):(2*n*(1-p)*k+1)]
  resi3<-yvalid3-prediction2[(2*n*(1-p)*k+2):(3*n*(1-p)*k+1)]
  
  resi1a<-yvalid1-prediction5[1:(n*(1-p)*k)]
  resi2a<-yvalid2-prediction5[(n*(1-p)*k+2):(2*n*(1-p)*k+1)]
  resi3a<-yvalid3-prediction5[(2*n*(1-p)*k+2):(3*n*(1-p)*k+1)]
  
  resi1b<-yvalid1-prediction7[1:(n*(1-p)*k)]
  resi2b<-yvalid2-prediction7[(n*(1-p)*k+2):(2*n*(1-p)*k+1)]
  resi3b<-yvalid3-prediction7[(2*n*(1-p)*k+2):(3*n*(1-p)*k+1)]
  
  for(j in 1:(n*(1-p)*k)){
    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[j],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[j],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[j],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[j],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*(1-p)*k+1+j)]<- d+prediction2[(n*(1-p)*k+1+j)]
    prediction6[(n*(1-p)*k+1+j)]<- da+prediction5[(n*(1-p)*k+1+j)]
    prediction8[(n*(1-p)*k+1+j)]<- db+prediction7[(n*(1-p)*k+1+j)]
    prediction9[(n*(1-p)*k+1+j)]<- dc+prediction7[(n*(1-p)*k+1+j)]
  }
  
  
  #for time point 3
  
  c1<-rep(0,(n*(1-p)*k+1))
  c1a<-rep(0,(n*(1-p)*k+1))
  c1b<-rep(0,(n*(1-p)*k+1))
  c1c<-rep(0,(n*(1-p)*k+1))
  
  for(j in 1:(n*(1-p)*k)){
    a<-seq(-30,30,0.01)
    
    d<-0
    da<-0
    db<-0
    dc<-0
    
    c1[j]<-VineCopula::RVinePDF(c(pnorm(resi1[j],0,out1$par[10]),pnorm(resi2[j],0,out1$par[14])),RVM4, verbose = TRUE)
    c1a[j]<-VineCopula::RVinePDF(c(pnorm(resi1a[j],0,out2$par[10]),pnorm(resi2a[j],0,out2$par[14])),RVM7, verbose = TRUE)
    c1b[j]<-VineCopula::RVinePDF(c(pnorm(resi1b[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals))),RVM10, verbose = TRUE)
    c1c[j]<-VineCopula::RVinePDF(c(pnorm(resi1b[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals))),RVM13, verbose = TRUE)
    
    c2<-function(k){
      cop<-pnorm(k,0,out1$par[18])
      dat<-c(pnorm(resi1[j],0,out1$par[10]),pnorm(resi2[j],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[j],0,out2$par[10]),pnorm(resi2a[j],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[j],0,sd(fit1$residuals)),pnorm(resi2b[j],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[j],0,sd(fit1$residuals)),pnorm(resi2b[j],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[j]
      
      da<-da+a[i]*c2a(a[i])*f1a(a[i])*0.01/c1a[j]
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01/c1b[j]
      
      dc<-dc+a[i]*c2c(a[i])*f1c(a[i])*0.01/c1c[j]
    }
    
    prediction1[2*n*(1-p)*k+1+j]<- d+prediction2[2*n*(1-p)*k+1+j]
    prediction6[2*n*(1-p)*k+1+j]<- da+prediction5[2*n*(1-p)*k+1+j]
    prediction8[2*n*(1-p)*k+1+j]<- db+prediction7[2*n*(1-p)*k+1+j]
    prediction9[2*n*(1-p)*k+1+j]<- dc+prediction7[2*n*(1-p)*k+1+j]
  }
  
  
  
  #for time point 4
  
  c1<-rep(0,(n*(1-p)*k+1))
  c1a<-rep(0,(n*(1-p)*k+1))
  c1b<-rep(0,(n*(1-p)*k+1))
  c1c<-rep(0,(n*(1-p)*k+1))
  
  for(j in 1:(n*(1-p)*k)){
    a<-seq(-30,30,0.01)
    b<-rep(0,6001)
    c<-rep(0,6001)
    
    d<-0
    da<-0
    db<-0
    dc<-0
    
    c1[j]<-VineCopula::RVinePDF(c(pnorm(resi1[j],0,out1$par[10]),pnorm(resi2[j],0,out1$par[14]),pnorm(resi3[j],0,out1$par[18])),RVM5, verbose = TRUE)
    
    c1a[j]<-VineCopula::RVinePDF(c(pnorm(resi1a[j],0,out2$par[10]),pnorm(resi2a[j],0,out2$par[14]),pnorm(resi3a[j],0,out2$par[18])),RVM8, verbose = TRUE)
    c1b[j]<-VineCopula::RVinePDF(c(pnorm(resi1b[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals)),pnorm(resi3b[j],0,sd(fit3$residuals))),RVM11, verbose = TRUE)
    c1c[j]<-VineCopula::RVinePDF(c(pnorm(resi1b[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals)),pnorm(resi3b[j],0,sd(fit3$residuals))),RVM14, verbose = TRUE)
    
    c2<-function(k){
      cop<-pnorm(k,0,out1$par[22])
      dat<-c(pnorm(resi1[j],0,out1$par[10]),pnorm(resi2[j],0,out1$par[14]),pnorm(resi3[j],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[j],0,out2$par[10]),pnorm(resi2a[j],0,out2$par[14]),pnorm(resi3a[j],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[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals)),pnorm(resi3b[j],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[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals)),pnorm(resi3b[j],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[j]
      
      da<-da+a[i]*c2a(a[i])*f1a(a[i])*0.01/c1a[j]
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01/c1b[j]
      
      dc<-dc+a[i]*c2c(a[i])*f1c(a[i])*0.01/c1c[j]
    }
    
    prediction1[3*n*(1-p)*k+1+j]<- d+prediction2[3*n*(1-p)*k+1+j]
    prediction6[3*n*(1-p)*k+1+j]<- da+prediction5[3*n*(1-p)*k+1+j]
    prediction8[3*n*(1-p)*k+1+j]<- db+prediction7[3*n*(1-p)*k+1+j]
    prediction9[3*n*(1-p)*k+1+j]<- dc+prediction7[3*n*(1-p)*k+1+j]
  }
  
  #fit independent regression model with seperating time block and time points
  newdata5<-data.frame(xtrain=xvalid,ztrain=zvalid,ktrain=kvalid)
  prediction3<-predict(fit,newdata5)
  
  #fit AR3 model and prediction
  newdata6<-data.frame(xtrain=xvalid,jtrain=jvalid)
  prediction4<-predict(fit0,newdata6)
  
  MAE_CVINE[m]<-sum(abs(yvalid-prediction1))/length(yvalid)
  MAE_MRM[m]<-sum(abs(yvalid-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(yvalid-prediction6))/length(yvalid)
  MAE_MRM1[m]<-sum(abs(yvalid-prediction5))/length(yvalid)
  MAE_CVINE2[m]<-sum(abs(yvalid-prediction8))/length(yvalid)
  MAE_CVINE3[m]<-sum(abs(yvalid-prediction9))/length(yvalid)
  MAE_MRM2[m]<-sum(abs(yvalid-prediction7))/length(yvalid)
  
  count=0
  for(i in 1:length(prediction1)){
    if(abs(yvalid-prediction1)[i]<abs(yvalid-prediction2)[i]){
      count=count+1
    }
  }
  percentage1[m]<-count/length(prediction1)
  
  count=0
  for(i in 1:length(prediction1)){
    if(abs(yvalid-prediction1)[i]<abs(yvalid-prediction3)[i]){
      count=count+1
    }
  }
  percentage2[m]<-count/length(prediction1)
  
  count=0
  for(i in 1:length(prediction1)){
    if(abs(yvalid-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