
#parameter setting for 4 time points within a time block. Only marginal parameters are changed, other parts are identical to scenario 5. 
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
library(VineCopula)
library(nlme)
library(forecast)
#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)



#different copula forms in different time blocks

RVinef1<-c(rep(0,17),5,6,5,rep(0,17),5,4,1,rep(0,17),1,3,6,rep(0,17),1,5,4,rep(0,17),1,1,3,rep(0,17),1,1,5,rep(0,17),1,1,1,rep(0,17),1,1,1,rep(0,17),6,3,4,rep(0,17),6,3,4,rep(0,17),1,6,3,rep(0,17),1,6,3,rep(0,17),1,1,6,rep(0,17),1,1,6,rep(0,17),1,1,1,rep(0,17),1,1,1,rep(0,17),5,4,1,rep(0,18),5,4,rep(0,19),5,rep(0,20))
Rvinef1<-matrix(RVinef1,20,20)
RVinepar11<-c(rep(0,17),3,3,3,rep(0,17),2,2,0.5,rep(0,17),0.5,4,3,rep(0,17),0.5,2,2,rep(0,17),0.5,0.5,4,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),4,3,3,rep(0,17),3,2,4,rep(0,17),0.5,4,3,rep(0,17),0.5,3,2,rep(0,17),0.5,0.5,4,rep(0,17),0.5,0.5,3,rep(0,17),0.5,0.5,0.5,rep(0,17),0.5,0.5,0.5,rep(0,17),3,3,0.8,rep(0,18),3,3,rep(0,19),3,rep(0,20))
Rvinepar11<-matrix(RVinepar11,20,20)
RVM111<-RVineMatrix(Rvinem1,Rvinef1,Rvinepar11,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)

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


MAE_CVINE<-rep(0,200)
MAE_CVINE1<-rep(0,200)
MAE_MRM<-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)
  simulate1<-RVineSim(n,RVM111)
  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)*4+1],0,sigma)
      epsilon2[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*4+2],0,sigma)
      epsilon3[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*4+3],0,sigma)
      epsilon4[(j-1)*n+i]=qnorm(simulate1[i,(j-1)*4+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 16 time points to predict the future 4 time points of the 500 subjects
  ytrain<-c(y1[1:((k-1)*n)],y2[1:((k-1)*n)],y3[1:((k-1)*n)],y4[1:((k-1)*n)])
  y1train<-y1[1:((k-1)*n)]
  y2train<-y2[1:((k-1)*n)]
  y3train<-y3[1:((k-1)*n)]
  y4train<-y4[1:((k-1)*n)]
  z0train<-c(rep(1,n),rep(2,n),rep(3,n),rep(4,n))
  
  yvalid<-c(y1[(n*(k-1)+1):(n*k)],y2[(n*(k-1)+1):(n*k)],y3[(n*(k-1)+1):(n*k)],y4[(n*(k-1)+1):(n*k)])
  y1valid<-y1[(n*(k-1)+1):(n*k)]
  y2valid<-y2[(n*(k-1)+1):(n*k)]
  y3valid<-y3[(n*(k-1)+1):(n*k)]
  y4valid<-y4[(n*(k-1)+1):(n*k)]
  
  xtrain<-c(x1[1:((k-1)*n)],x2[1:((k-1)*n)],x3[1:((k-1)*n)],x4[1:((k-1)*n)])
  x1train<-x1[1:((k-1)*n)]
  x2train<-x2[1:((k-1)*n)]
  x3train<-x3[1:((k-1)*n)]
  x4train<-x4[1:((k-1)*n)]
  
  xvalid<-c(x1[(n*(k-1)+1):(n*k)],x2[(n*(k-1)+1):(n*k)],x3[(n*(k-1)+1):(n*k)],x4[(n*(k-1)+1):(n*k)])
  x1valid<-x1[(n*(k-1)+1):(n*k)]
  x2valid<-x2[(n*(k-1)+1):(n*k)]
  x3valid<-x3[(n*(k-1)+1):(n*k)]
  x4valid<-x4[(n*(k-1)+1):(n*k)]
  z0valid<-rep(5,n)
  zvalid<-rep(5,(n*(k-1)))
  
  #fit CVM model
  #Marginal Regression
  fit1<-lm(y1train~x1train+z0train)
  
  fit2<-lm(y2train~x2train+z0train)
  
  fit3<-lm(y3train~x3train+z0train)
  
  fit4<-lm(y4train~x4train+z0train)
  
  
  
  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((y1train[r]-param[7]-param[8]*x1train[r]-param[9]*z0train[r]),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((y2train[r]-param[11]-param[12]*x2train[r]-param[13]*z0train[r]),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((y3train[r]-param[15]-param[16]*x3train[r]-param[17]*z0train[r]),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((y4train[r]-param[19]-param[20]*x4train[r]-param[21]*z0train[r]),mean=0,sd=param[22],log=TRUE)
    }
    
    u1<-pnorm((y1train-param[7]-param[8]*x1train-param[9]*z0train),mean=0,sd=param[10])
    u2<-pnorm((y2train-param[11]-param[12]*x2train-param[13]*z0train),mean=0,sd=param[14])
    u3<-pnorm((y3train-param[15]-param[16]*x3train-param[17]*z0train),mean=0,sd=param[18])
    u4<-pnorm((y4train-param[19]-param[20]*x4train-param[21]*z0train),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)
  
  
  #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((y1train-fit1$coefficients[1]-fit1$coefficients[2]*x1train-fit1$coefficients[3]*z0train),mean=0,sd=sd(fit1$residuals))
    u2<-pnorm((y2train-fit2$coefficients[1]-fit2$coefficients[2]*x2train-fit2$coefficients[3]*z0train),mean=0,sd=sd(fit2$residuals))
    u3<-pnorm((y3train-fit3$coefficients[1]-fit3$coefficients[2]*x3train-fit3$coefficients[3]*z0train),mean=0,sd=sd(fit3$residuals))
    u4<-pnorm((y4train-fit4$coefficients[1]-fit4$coefficients[2]*x4train-fit4$coefficients[3]*z0train),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])
  
  #fit a regression model without seperation of time points
  ktrain<-c(rep(1,((k-1)*n)),rep(2,((k-1)*n)),rep(3,((k-1)*n)),rep(4,((k-1)*n)))
  ztrain<-rep(c(rep(1,n),rep(2,n),rep(3,n),rep(4,n)),(k-1))
  fit<-lm(ytrain~xtrain+ztrain+ktrain)
  kvalid<-c(rep(1,n),rep(2,n),rep(3,n),rep(4,n))
  
  #fit a ARMA model
  jtrain=NULL
  
  for(j in 1:(k-1)){
    for(i in 1:l){
      jtrain<-c(jtrain,rep((j-1)*l+i,n))
    }
  }
  
  jvalid=NULL
  
  for(i in 1:l){
    jvalid<-c(jvalid,rep((k-1)*l+i,n))
  }
  
  
  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*l)
  prediction2<-rep(0,n*l)
  prediction3<-rep(0,n*l)
  prediction4<-rep(0,n*l)
  #Sequential Estimation Selection Prediction for MRM and CVine
  prediction7<-rep(0,n*l)
  prediction8<-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
  
  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]
  prediction8[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)
  
  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)
  
  #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)
  
  #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)
  
  resi1<-y1valid-prediction2[1:n]
  resi2<-y2valid-prediction2[(n+1):(2*n)]
  resi3<-y3valid-prediction2[(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(j in 1:n){
    a<-seq(-30,30,0.01)
    
    d<-0
    db<-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)
    }
    
    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)
    }
    
    f1<-function(k){
      dnorm(k,0,out1$par[14])
    }
    
    f1b<-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
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01
    }
    
    prediction1[(n+j)]<- sum(d)+prediction2[(n+j)]
    prediction8[(n+j)]<- db+prediction7[(n+j)]
  }
  
  
  #for time point 3
  
  c1<-rep(0,n)
  c1b<-rep(0,n)
  
  for(j in 1:n){
    a<-seq(-30,30,0.01)
    
    d<-0
    db<-0
    
    c1[j]<-VineCopula::RVinePDF(c(pnorm(resi1[j],0,out1$par[10]),pnorm(resi2[j],0,out1$par[14])),RVM4, verbose = TRUE)
    c1b[j]<-VineCopula::RVinePDF(c(pnorm(resi1b[j],0,sd(fit1$residuals)),pnorm(resi2b[j],0,sd(fit2$residuals))),RVM10, 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)
    }
    
    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)
    }
    
    f1<-function(k){
      dnorm(k,0,out1$par[18])
    }
    
    f1b<-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]
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01/c1b[j]
    }
    
    prediction1[2*n+j]<- sum(d)+prediction2[2*n+j]
    prediction8[2*n+j]<- db+prediction7[2*n+j]
  }
  
  
  
  #for time point 4
  
  c1<-rep(0,n)
  c1b<-rep(0,n)
  
  for(j in 1:n){
    a<-seq(-30,30,0.01)
    
    d<-0
    db<-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)
    
    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)
    
    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)
    }
    
    
    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)
    }
    
    
    
    f1<-function(k){
      dnorm(k,0,out1$par[22])
    }
    
    
    f1b<-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]
      
      db<-db+a[i]*c2b(a[i])*f1b(a[i])*0.01/c1b[j]
      
    }
    
    prediction1[3*n+j]<- d+prediction2[3*n+j]
    prediction8[3*n+j]<- db+prediction7[3*n+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 ARMA 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-prediction8))/length(yvalid)
  MAE_MRM1[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


percentage1
percentage2
percentage3