library(VineCopula)
library(Matrix)
library(matrixcalc)
library(matlib)
library(SpatioTemporal)
library(nlme)

#parameter setting for 4 time points within a time block
beta01<-1
beta02<-2
beta03<-3
beta04<-4
beta11<-2
beta12<-3
beta13<-4
beta14<-5
beta21<-3
beta22<-4
beta23<-5
beta24<-6
sigma<-2

#16 time points structure
RVinem<-c(16,4,3,2,1,8,7,6,5,12,11,10,9,15,14,13,0,4,15,14,13,12,11,10,9,8,7,6,5,3,2,1,rep(0,2),3,15,14,13,12,11,10,9,8,7,6,5,2,1,rep(0,3),15,2,1,8,7,6,5,12,11,10,9,14,13,rep(0,4),2,14,13,12,11,10,9,8,7,6,5,1,rep(0,5),14,1,8,7,6,5,12,11,10,9,13,rep(0,6),1,13,12,11,10,9,8,7,6,5,rep(0,7),13,8,7,6,5,12,11,10,9,rep(0,8),8,12,11,10,9,7,6,5,rep(0,9),12,7,6,5,11,10,9,rep(0,10),7,11,10,9,6,5,rep(0,11),11,6,5,10,9,rep(0,12),6,10,9,5,rep(0,13),10,5,9,rep(0,14),5,9,rep(0,15),9)

#Rvinem is the matrix that specify the RVine structure for a=4
Rvinem<-matrix(RVinem,16,16)


RVinef<-c(rep(0,13),3,4,1,rep(0,13),3,4,1,rep(0,13),1,3,4,rep(0,13),1,3,4,rep(0,13),1,1,5,rep(0,13),1,1,5,rep(0,13),1,1,1,rep(0,13),1,1,1,rep(0,13),3,4,1,rep(0,13),3,4,1,rep(0,13),1,3,4,rep(0,13),1,3,4,rep(0,13),1,1,5,rep(0,14),1,5,rep(0,15),1,rep(0,16))
#Rvinef is the matrix that specify the copula forms for each (conditional) bivariate copula
Rvinef<-matrix(RVinef,16,16)

RVinepar1<-c(rep(0,13),2,3,0.5,rep(0,13),2,3,0.5,rep(0,13),0.5,3,2,rep(0,13),0.5,3,2,rep(0,13),0.5,0.5,2,rep(0,13),0.5,0.5,2,rep(0,13),0.5,0.5,0.5,rep(0,13),0.5,0.5,0.5,rep(0,13),2,3,0.5,rep(0,13),2,3,0.5,rep(0,13),0.5,3,2,rep(0,13),0.5,3,2,rep(0,13),0.5,0.5,2,rep(0,14),0.5,2,rep(0,15),0.5,rep(0,16))
#Rvinepar1 is the matrix that specify the corresponding first parameters for the copula forms specified above
Rvinepar1<-matrix(RVinepar1,16,16)

RVinepar2<-rep(0,256)
#Rvinepar2 is the matrix that specify the corresponding second parameters for the copula forms specified above, which is all 0 since we use the one-parameter copulas here 
Rvinepar2<-matrix(RVinepar2,16,16)

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


#n is sample size, k is the number of time block and l is the number of time points within a time block#n is sample size, k is the number of time block and l is the number of time points within a time block
n<-500
k<-4
l<-4

repeattime<-500

#the upper and lower bounds for the marginal parameters(lo and up) and the copula parameters(lo1 and up1)
#lower
lo<-c(rep(c(-Inf,-Inf,-Inf,0),4))
lo1<-c(0,1,-1,0,1,-Inf)

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


#set up the matrix that record the simulation results
cmarginalbeta01<-rep(0,repeattime)
cmarginalbeta02<-rep(0,repeattime)
cmarginalbeta03<-rep(0,repeattime)
cmarginalbeta04<-rep(0,repeattime)
cmarginalbeta11<-rep(0,repeattime)
cmarginalbeta12<-rep(0,repeattime)
cmarginalbeta13<-rep(0,repeattime)
cmarginalbeta14<-rep(0,repeattime)
cmarginalbeta21<-rep(0,repeattime)
cmarginalbeta22<-rep(0,repeattime)
cmarginalbeta23<-rep(0,repeattime)
cmarginalbeta24<-rep(0,repeattime)
cmu1<-rep(0,repeattime)
cmu2<-rep(0,repeattime)
cmu3<-rep(0,repeattime)
cmu4<-rep(0,repeattime)
csigma1<-rep(0,repeattime)
csigma2<-rep(0,repeattime)
csigma3<-rep(0,repeattime)
csigma4<-rep(0,repeattime)
csdmarginalbeta01<-rep(0,repeattime)
csdmarginalbeta02<-rep(0,repeattime)
csdmarginalbeta03<-rep(0,repeattime)
csdmarginalbeta04<-rep(0,repeattime)
csdmarginalbeta11<-rep(0,repeattime)
csdmarginalbeta12<-rep(0,repeattime)
csdmarginalbeta13<-rep(0,repeattime)
csdmarginalbeta14<-rep(0,repeattime)
csdmarginalbeta21<-rep(0,repeattime)
csdmarginalbeta22<-rep(0,repeattime)
csdmarginalbeta23<-rep(0,repeattime)
csdmarginalbeta24<-rep(0,repeattime)
csdsigma1<-rep(0,repeattime)
csdsigma2<-rep(0,repeattime)
csdsigma3<-rep(0,repeattime)
csdsigma4<-rep(0,repeattime)

cpar12<-rep(0,repeattime)
cpar13<-rep(0,repeattime)
cpar14<-rep(0,repeattime)
cpar23<-rep(0,repeattime)
cpar24<-rep(0,repeattime)
cpar34<-rep(0,repeattime)
mlecomppar<-matrix(0,4,4)
esecompstd<-matrix(0,4,4)
asecompstd<-matrix(0,4,4)

#simulation repeat 500 times
for(i in 1:repeattime){
  set.seed(i)
  #generate copula data that is uniformly distributed
  simulate1<-RVineSim(n,RVM1)
  
  #transform the copula data to the error terms using the normal distribution
  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(h in 1:n){
      epsilon1[(j-1)*n+h]=qnorm(simulate1[h,(j-1)*l+1],0,sigma)
      epsilon2[(j-1)*n+h]=qnorm(simulate1[h,(j-1)*l+2],0,sigma)
      epsilon3[(j-1)*n+h]=qnorm(simulate1[h,(j-1)*l+3],0,sigma)
      epsilon4[(j-1)*n+h]=qnorm(simulate1[h,(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)
  
  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)]
  }
  
  #estimate the marginal parameters first(first stage estimation)
  
  #likelihood function for marginal parameters in time block 1,2,3,4, which are optimfun1,optimfun2, optimfun3 and optimfun4, respectively
  optimfun1<-function(param){
    
    logp1<-0
    for(r in 1:n){
      logp1<-logp1+dnorm((y1[r]-param[1]-param[2]*x1[r]-param[3]),mean=0,sd=param[4],log=TRUE)+dnorm((y1[n+r]-param[1]-param[2]*x1[n+r]-param[3]*2),mean=0,sd=param[4],log=TRUE)+dnorm((y1[2*n+r]-param[1]-param[2]*x1[2*n+r]-param[3]*3),mean=0,sd=param[4],log=TRUE)+dnorm((y1[3*n+r]-param[1]-param[2]*x1[3*n+r]-param[3]*4),mean=0,sd=param[4],log=TRUE)
    }
    -logp1
  }
  
  optimfun2<-function(param){
    
    logp2<-0
    for(r in 1:n){
      logp2<-logp2+dnorm((y2[r]-param[1]-param[2]*x2[r]-param[3]),mean=0,sd=param[4],log=TRUE)+dnorm((y2[n+r]-param[1]-param[2]*x2[n+r]-param[3]*2),mean=0,sd=param[4],log=TRUE)+dnorm((y2[2*n+r]-param[1]-param[2]*x2[2*n+r]-param[3]*3),mean=0,sd=param[4],log=TRUE)+dnorm((y2[3*n+r]-param[1]-param[2]*x2[3*n+r]-param[3]*4),mean=0,sd=param[4],log=TRUE)
    }
    -logp2
  }
  
  optimfun3<-function(param){
    
    logp3<-0
    for(r in 1:n){
      logp3<-logp3+dnorm((y3[r]-param[1]-param[2]*x3[r]-param[3]),mean=0,sd=param[4],log=TRUE)+dnorm((y3[n+r]-param[1]-param[2]*x3[n+r]-param[3]*2),mean=0,sd=param[4],log=TRUE)+dnorm((y3[2*n+r]-param[1]-param[2]*x3[2*n+r]-param[3]*3),mean=0,sd=param[4],log=TRUE)+dnorm((y3[3*n+r]-param[1]-param[2]*x3[3*n+r]-param[3]*4),mean=0,sd=param[4],log=TRUE)
    }
    -logp3
  }
  
  
  optimfun4<-function(param){
    
    logp4<-0
    for(r in 1:n){
      logp4<-logp4+dnorm((y4[r]-param[1]-param[2]*x4[r]-param[3]),mean=0,sd=param[4],log=TRUE)+dnorm((y4[n+r]-param[1]-param[2]*x4[n+r]-param[3]*2),mean=0,sd=param[4],log=TRUE)+dnorm((y4[2*n+r]-param[1]-param[2]*x4[2*n+r]-param[3]*3),mean=0,sd=param[4],log=TRUE)+dnorm((y4[3*n+r]-param[1]-param[2]*x4[3*n+r]-param[3]*4),mean=0,sd=param[4],log=TRUE)
    }
    -logp4
  }
  
  
  #start parameters
  startpar1<-c(beta01,beta11,beta21,sigma)
  startpar2<-c(beta02,beta12,beta22,sigma)
  startpar3<-c(beta03,beta13,beta23,sigma)
  startpar4<-c(beta04,beta14,beta24,sigma)
  
  #minimize the likelihood function with negativve sign to obtain the estimate of marginal parameters
  out1<-nlminb(objective=optimfun1,start=startpar1,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo,upper=up)
  out2<-nlminb(objective=optimfun2,start=startpar2,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo,upper=up)
  out3<-nlminb(objective=optimfun3,start=startpar3,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo,upper=up)
  out4<-nlminb(objective=optimfun4,start=startpar4,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo,upper=up)
  
  
  
  
  
  #estimate the vine parameters(second stage estimation)
  #likelihood function with negative sign (here we use composite likelihood instead of the full likelihood, so there is only 6 parameters in the likelihood function)
  optimfun5<-function(param){
    RVinepar3<-c(rep(0,13),param[1:3],rep(0,13),param[1],param[2],param[3],rep(0,14),param[4:5],rep(0,14),param[4],param[5],rep(0,15),param[6],rep(0,15),param[6],rep(0,16),rep(0,16),rep(0,13),param[1],param[2],param[3],rep(0,13),param[1],param[2],param[3],rep(0,14),param[4],param[5],rep(0,14),param[4],param[5],rep(0,15),param[6],rep(0,15),param[6],rep(0,16),rep(0,16))
    Rvinepar3<-matrix(RVinepar3,16,16)
    RVM2<-RVineMatrix(Rvinem,Rvinef,Rvinepar3,Rvinepar2,check.pars=FALSE)
    
    
    u1<-pnorm((y1[1:n]-out1$par[1]-out1$par[2]*x1[1:n]-out1$par[3]),mean=0,sd=out1$par[4])
    u2<-pnorm((y2[1:n]-out2$par[1]-out2$par[2]*x2[1:n]-out2$par[3]),mean=0,sd=out2$par[4])
    u3<-pnorm((y3[1:n]-out3$par[1]-out3$par[2]*x3[1:n]-out3$par[3]),mean=0,sd=out3$par[4])
    u4<-pnorm((y4[1:n]-out4$par[1]-out4$par[2]*x4[1:n]-out4$par[3]),mean=0,sd=out4$par[4])
    u5<-pnorm((y1[(n+1):(2*n)]-out1$par[1]-out1$par[2]*x1[(n+1):(2*n)]-out1$par[3]*2),mean=0,sd=out1$par[4])
    u6<-pnorm((y2[(n+1):(2*n)]-out2$par[1]-out2$par[2]*x2[(n+1):(2*n)]-out2$par[3]*2),mean=0,sd=out2$par[4])
    u7<-pnorm((y3[(n+1):(2*n)]-out3$par[1]-out3$par[2]*x3[(n+1):(2*n)]-out3$par[3]*2),mean=0,sd=out3$par[4])
    u8<-pnorm((y4[(n+1):(2*n)]-out4$par[1]-out4$par[2]*x4[(n+1):(2*n)]-out4$par[3]*2),mean=0,sd=out4$par[4])
    u9<-pnorm((y1[(2*n+1):(3*n)]-out1$par[1]-out1$par[2]*x1[(2*n+1):(3*n)]-out1$par[3]*3),mean=0,sd=out1$par[4])
    u10<-pnorm((y2[(2*n+1):(3*n)]-out2$par[1]-out2$par[2]*x2[(2*n+1):(3*n)]-out2$par[3]*3),mean=0,sd=out2$par[4])
    u11<-pnorm((y3[(2*n+1):(3*n)]-out3$par[1]-out3$par[2]*x3[(2*n+1):(3*n)]-out3$par[3]*3),mean=0,sd=out3$par[4])
    u12<-pnorm((y4[(2*n+1):(3*n)]-out4$par[1]-out4$par[2]*x4[(2*n+1):(3*n)]-out4$par[3]*3),mean=0,sd=out4$par[4])
    u13<-pnorm((y1[(3*n+1):(4*n)]-out1$par[1]-out1$par[2]*x1[(3*n+1):(4*n)]-out1$par[3]*4),mean=0,sd=out1$par[4])
    u14<-pnorm((y2[(3*n+1):(4*n)]-out2$par[1]-out2$par[2]*x2[(3*n+1):(4*n)]-out2$par[3]*4),mean=0,sd=out2$par[4])
    u15<-pnorm((y3[(3*n+1):(4*n)]-out3$par[1]-out3$par[2]*x3[(3*n+1):(4*n)]-out3$par[3]*4),mean=0,sd=out3$par[4])
    u16<-pnorm((y4[(3*n+1):(4*n)]-out4$par[1]-out4$par[2]*x4[(3*n+1):(4*n)]-out4$par[3]*4),mean=0,sd=out4$par[4])
    
    
    
    simulate2<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16)
    -RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
  }
  startpar5<-c(2,3,0.5,3,2,2)
  out5<-nlminb(objective=optimfun5,start=startpar5,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo1,upper=up1)
  
  
  
  #this is the likelihood function for marginal parameters without the negative sign to calculate the hessian matrix and J matrix
  optimfun6<-function(param){
    
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(k in 1:n){
      logp1<-logp1+dnorm((y1[k]-param[1]-param[2]*x1[k]-param[3]),mean=0,sd=param[4],log=TRUE)+dnorm((y1[n+k]-param[1]-param[2]*x1[n+k]-param[3]*2),mean=0,sd=param[4],log=TRUE)+dnorm((y1[2*n+k]-param[1]-param[2]*x1[2*n+k]-param[3]*3),mean=0,sd=param[4],log=TRUE)+dnorm((y1[3*n+k]-param[1]-param[2]*x1[3*n+k]-param[3]*4),mean=0,sd=param[4],log=TRUE)
      logp2<-logp2+dnorm((y2[k]-param[5]-param[6]*x2[k]-param[7]),mean=0,sd=param[8],log=TRUE)+dnorm((y2[n+k]-param[5]-param[6]*x2[n+k]-param[7]*2),mean=0,sd=param[8],log=TRUE)+dnorm((y2[2*n+k]-param[5]-param[6]*x2[2*n+k]-param[7]*3),mean=0,sd=param[8],log=TRUE)+dnorm((y2[3*n+k]-param[5]-param[6]*x2[3*n+k]-param[7]*4),mean=0,sd=param[8],log=TRUE)
      logp3<-logp3+dnorm((y3[k]-param[9]-param[10]*x3[k]-param[11]),mean=0,sd=param[12],log=TRUE)+dnorm((y3[n+k]-param[9]-param[10]*x3[n+k]-param[11]*2),mean=0,sd=param[12],log=TRUE)+dnorm((y3[2*n+k]-param[9]-param[10]*x3[2*n+k]-param[11]*3),mean=0,sd=param[12],log=TRUE)+dnorm((y3[3*n+k]-param[9]-param[10]*x3[3*n+k]-param[11]*4),mean=0,sd=param[12],log=TRUE)
      logp4<-logp4+dnorm((y4[k]-param[13]-param[14]*x4[k]-param[15]),mean=0,sd=param[16],log=TRUE)+dnorm((y4[n+k]-param[13]-param[14]*x4[n+k]-param[15]*2),mean=0,sd=param[16],log=TRUE)+dnorm((y4[2*n+k]-param[13]-param[14]*x4[2*n+k]-param[15]*3),mean=0,sd=param[16],log=TRUE)+dnorm((y4[3*n+k]-param[13]-param[14]*x4[3*n+k]-param[15]*4),mean=0,sd=param[16],log=TRUE)
    }
    
    logp1+logp2+logp3+logp4
  }
  out<-fdHess(c(out1$par,out2$par,out3$par,out4$par),optimfun6)
  
  #hessian matrix for marginal parameters
  H1<-out$Hessian
  
  
  #this is the likelihood function without the negative sign of the dependence parameters to calculate the hessian matrix and J matrix
  optimfun7<-function(param){
    RVinepar3<-c(rep(0,13),param[1:3],rep(0,13),param[1],param[2],param[3],rep(0,14),param[4:5],rep(0,14),param[4],param[5],rep(0,15),param[6],rep(0,15),param[6],rep(0,16),rep(0,16),rep(0,13),param[1],param[2],param[3],rep(0,13),param[1],param[2],param[3],rep(0,14),param[4],param[5],rep(0,14),param[4],param[5],rep(0,15),param[6],rep(0,15),param[6],rep(0,16),rep(0,16))
    Rvinepar3<-matrix(RVinepar3,16,16)
    RVM2<-RVineMatrix(Rvinem,Rvinef,Rvinepar3,Rvinepar2,check.pars=FALSE)
    
    
    u1<-pnorm((y1[1:n]-param[7]-param[8]*x1[1:n]-param[9]),mean=0,sd=param[10])
    u2<-pnorm((y2[1:n]-param[11]-param[12]*x2[1:n]-param[13]),mean=0,sd=param[14])
    u3<-pnorm((y3[1:n]-param[15]-param[16]*x3[1:n]-param[17]),mean=0,sd=param[18])
    u4<-pnorm((y4[1:n]-param[19]-param[20]*x4[1:n]-param[21]),mean=0,sd=param[22])
    u5<-pnorm((y1[(n+1):(2*n)]-param[7]-param[8]*x1[(n+1):(2*n)]-param[9]*2),mean=0,sd=param[10])
    u6<-pnorm((y2[(n+1):(2*n)]-param[11]-param[12]*x2[(n+1):(2*n)]-param[13]*2),mean=0,sd=param[14])
    u7<-pnorm((y3[(n+1):(2*n)]-param[15]-param[16]*x3[(n+1):(2*n)]-param[17]*2),mean=0,sd=param[18])
    u8<-pnorm((y4[(n+1):(2*n)]-param[19]-param[20]*x4[(n+1):(2*n)]-param[21]*2),mean=0,sd=param[22])
    u9<-pnorm((y1[(2*n+1):(3*n)]-param[7]-param[8]*x1[(2*n+1):(3*n)]-param[9]*3),mean=0,sd=param[10])
    u10<-pnorm((y2[(2*n+1):(3*n)]-param[11]-param[12]*x2[(2*n+1):(3*n)]-param[13]*3),mean=0,sd=param[14])
    u11<-pnorm((y3[(2*n+1):(3*n)]-param[15]-param[16]*x3[(2*n+1):(3*n)]-param[17]*3),mean=0,sd=param[18])
    u12<-pnorm((y4[(2*n+1):(3*n)]-param[19]-param[20]*x4[(2*n+1):(3*n)]-param[21]*3),mean=0,sd=param[22])
    u13<-pnorm((y1[(3*n+1):(4*n)]-param[7]-param[8]*x1[(3*n+1):(4*n)]-param[9]*4),mean=0,sd=param[10])
    u14<-pnorm((y2[(3*n+1):(4*n)]-param[11]-param[12]*x2[(3*n+1):(4*n)]-param[13]*4),mean=0,sd=param[14])
    u15<-pnorm((y3[(3*n+1):(4*n)]-param[15]-param[16]*x3[(3*n+1):(4*n)]-param[17]*4),mean=0,sd=param[18])
    u16<-pnorm((y4[(3*n+1):(4*n)]-param[19]-param[20]*x4[(3*n+1):(4*n)]-param[21]*4),mean=0,sd=param[22])
    
    
    
    simulate2<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16)
    RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
  }
  out<-fdHess(c(out5$par,out1$par,out2$par,out3$par,out4$par),optimfun7)
  #hessian matrix of the dependence parameters
  H2<-out$Hessian
  
  
  #calculate omega matrix(a matrix used to calculate the correct ase of dependence parameters)
  Q<-H2[1:6,7:22]%*%inv(H1)
  
  J<-matrix(0,16,16)
  omega<-matrix(0,6,6)
  for(r in 1:n){
    optimfun8<-function(param){
      
      logp1<-dnorm((y1[r]-param[1]-param[2]*x1[r]-param[3]),mean=0,sd=param[4],log=TRUE)+dnorm((y1[n+r]-param[1]-param[2]*x1[n+r]-param[3]*2),mean=0,sd=param[4],log=TRUE)+dnorm((y1[2*n+r]-param[1]-param[2]*x1[2*n+r]-param[3]*3),mean=0,sd=param[4],log=TRUE)+dnorm((y1[3*n+r]-param[1]-param[2]*x1[3*n+r]-param[3]*4),mean=0,sd=param[4],log=TRUE)
      logp2<-dnorm((y2[r]-param[5]-param[6]*x2[r]-param[7]),mean=0,sd=param[8],log=TRUE)+dnorm((y2[n+r]-param[5]-param[6]*x2[n+r]-param[7]*2),mean=0,sd=param[8],log=TRUE)+dnorm((y2[2*n+r]-param[5]-param[6]*x2[2*n+r]-param[7]*3),mean=0,sd=param[8],log=TRUE)+dnorm((y2[3*n+r]-param[5]-param[6]*x2[3*n+r]-param[7]*4),mean=0,sd=param[8],log=TRUE)
      logp3<-dnorm((y3[r]-param[9]-param[10]*x3[r]-param[11]),mean=0,sd=param[12],log=TRUE)+dnorm((y3[n+r]-param[9]-param[10]*x3[n+r]-param[11]*2),mean=0,sd=param[12],log=TRUE)+dnorm((y3[2*n+r]-param[9]-param[10]*x3[2*n+r]-param[11]*3),mean=0,sd=param[12],log=TRUE)+dnorm((y3[3*n+r]-param[9]-param[10]*x3[3*n+r]-param[11]*4),mean=0,sd=param[12],log=TRUE)
      logp4<-dnorm((y4[r]-param[13]-param[14]*x4[r]-param[15]),mean=0,sd=param[16],log=TRUE)+dnorm((y4[n+r]-param[13]-param[14]*x4[n+r]-param[15]*2),mean=0,sd=param[16],log=TRUE)+dnorm((y4[2*n+r]-param[13]-param[14]*x4[2*n+r]-param[15]*3),mean=0,sd=param[16],log=TRUE)+dnorm((y4[3*n+r]-param[13]-param[14]*x4[3*n+r]-param[15]*4),mean=0,sd=param[16],log=TRUE)
      
      
      logp1+logp2+logp3+logp4
    }
    out6<-genGradient(c(out1$par,out2$par,out3$par,out4$par),optimfun8)
    J<-J+out6%*%t(out6)
    
    optimfun9<-function(param){
      RVinepar3<-c(rep(0,13),param[1:3],rep(0,13),param[1],param[2],param[3],rep(0,14),param[4:5],rep(0,14),param[4],param[5],rep(0,15),param[6],rep(0,15),param[6],rep(0,16),rep(0,16),rep(0,13),param[1],param[2],param[3],rep(0,13),param[1],param[2],param[3],rep(0,14),param[4],param[5],rep(0,14),param[4],param[5],rep(0,15),param[6],rep(0,15),param[6],rep(0,16),rep(0,16))
      Rvinepar3<-matrix(RVinepar3,16,16)
      RVM2<-RVineMatrix(Rvinem,Rvinef,Rvinepar3,Rvinepar2,check.pars=FALSE)
      
      
      
      
      u1<-pnorm((y1[r]-param[7]-param[8]*x1[r]-param[9]),mean=0,sd=param[10])
      u2<-pnorm((y2[r]-param[11]-param[12]*x2[r]-param[13]),mean=0,sd=param[14])
      u3<-pnorm((y3[r]-param[15]-param[16]*x3[r]-param[17]),mean=0,sd=param[18])
      u4<-pnorm((y4[r]-param[19]-param[20]*x4[r]-param[21]),mean=0,sd=param[22])
      u5<-pnorm((y1[n+r]-param[7]-param[8]*x1[n+r]-param[9]*2),mean=0,sd=param[10])
      u6<-pnorm((y2[n+r]-param[11]-param[12]*x2[n+r]-param[13]*2),mean=0,sd=param[14])
      u7<-pnorm((y3[n+r]-param[15]-param[16]*x3[n+r]-param[17]*2),mean=0,sd=param[18])
      u8<-pnorm((y4[n+r]-param[19]-param[20]*x4[n+r]-param[21]*2),mean=0,sd=param[22])
      u9<-pnorm((y1[2*n+r]-param[7]-param[8]*x1[2*n+r]-param[9]*3),mean=0,sd=param[10])
      u10<-pnorm((y2[2*n+r]-param[11]-param[12]*x2[2*n+r]-param[13]*3),mean=0,sd=param[14])
      u11<-pnorm((y3[2*n+r]-param[15]-param[16]*x3[2*n+r]-param[17]*3),mean=0,sd=param[18])
      u12<-pnorm((y4[2*n+r]-param[19]-param[20]*x4[2*n+r]-param[21]*3),mean=0,sd=param[22])
      u13<-pnorm((y1[3*n+r]-param[7]-param[8]*x1[3*n+r]-param[9]*4),mean=0,sd=param[10])
      u14<-pnorm((y2[3*n+r]-param[11]-param[12]*x2[3*n+r]-param[13]*4),mean=0,sd=param[14])
      u15<-pnorm((y3[3*n+r]-param[15]-param[16]*x3[3*n+r]-param[17]*4),mean=0,sd=param[18])
      u16<-pnorm((y4[3*n+r]-param[19]-param[20]*x4[3*n+r]-param[21]*4),mean=0,sd=param[22])
      
      simulate2<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16)
      RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
    }
    out7<-genGradient(c(out5$par,out1$par,out2$par,out3$par,out4$par),optimfun9)
    ome<-out7[1:6]-Q%*%out6
    omega<-omega+ome%*%t(ome)
  }
  J<-J/n
  H1<-H1/n
  
  #robustase is the Godambe information matrix for marginal parameters
  robustase<-inv(H1)%*%J%*%inv(H1)
  
  omega<-omega/n
  H3<-H2/n
  #robustase1 is the Godambe information matrix for dependence parameters
  robustase1<-inv(H3[1:6,1:6])%*%omega%*%inv(H3[1:6,1:6])
  
  
  #record the results
  cmarginalbeta01[i]<-out1$par[1]
  cmarginalbeta02[i]<-out2$par[1]
  cmarginalbeta03[i]<-out3$par[1]
  cmarginalbeta04[i]<-out4$par[1]
  cmarginalbeta11[i]<-out1$par[2]
  cmarginalbeta12[i]<-out2$par[2]
  cmarginalbeta13[i]<-out3$par[2]
  cmarginalbeta14[i]<-out4$par[2]
  cmarginalbeta21[i]<-out1$par[3]
  cmarginalbeta22[i]<-out2$par[3]
  cmarginalbeta23[i]<-out3$par[3]
  cmarginalbeta24[i]<-out4$par[3]
  csigma1[i]<-out1$par[4]
  csigma2[i]<-out2$par[4]
  csigma3[i]<-out3$par[4]
  csigma4[i]<-out4$par[4]
  
  
  csdmarginalbeta01[i]<-(sqrt(diag(robustase))/sqrt(n))[1]
  csdmarginalbeta02[i]<-(sqrt(diag(robustase))/sqrt(n))[5]
  csdmarginalbeta03[i]<-(sqrt(diag(robustase))/sqrt(n))[9]
  csdmarginalbeta04[i]<-(sqrt(diag(robustase))/sqrt(n))[13]
  csdmarginalbeta11[i]<-(sqrt(diag(robustase))/sqrt(n))[2]
  csdmarginalbeta12[i]<-(sqrt(diag(robustase))/sqrt(n))[6]
  csdmarginalbeta13[i]<-(sqrt(diag(robustase))/sqrt(n))[10]
  csdmarginalbeta14[i]<-(sqrt(diag(robustase))/sqrt(n))[14]
  csdmarginalbeta21[i]<-(sqrt(diag(robustase))/sqrt(n))[3]
  csdmarginalbeta22[i]<-(sqrt(diag(robustase))/sqrt(n))[7]
  csdmarginalbeta23[i]<-(sqrt(diag(robustase))/sqrt(n))[11]
  csdmarginalbeta24[i]<-(sqrt(diag(robustase))/sqrt(n))[15]
  csdsigma1[i]<-(sqrt(diag(robustase))/sqrt(n))[4]
  csdsigma2[i]<-(sqrt(diag(robustase))/sqrt(n))[8]
  csdsigma3[i]<-(sqrt(diag(robustase))/sqrt(n))[12]
  csdsigma4[i]<-(sqrt(diag(robustase))/sqrt(n))[16]
  
  
  cpar12[i]<-out5$par[6]
  cpar13[i]<-out5$par[5]
  cpar14[i]<-out5$par[3]
  cpar23[i]<-out5$par[4]
  cpar24[i]<-out5$par[2]
  cpar34[i]<-out5$par[1]
  mlecomppar[4,3]<-mlecomppar[4,3]+out5$par[6]
  mlecomppar[4,2]<-mlecomppar[4,2]+out5$par[5]
  mlecomppar[4,1]<-mlecomppar[4,1]+out5$par[3]
  mlecomppar[3,2]<-mlecomppar[3,2]+out5$par[4]
  mlecomppar[3,1]<-mlecomppar[3,1]+out5$par[2]
  mlecomppar[2,1]<-mlecomppar[2,1]+out5$par[1]
  
  
  
  
  compstd<-matrix(0,4,4)
  compstd[4,3]<-(sqrt(diag(robustase1))/sqrt(n))[6]
  compstd[4,2]<-(sqrt(diag(robustase1))/sqrt(n))[5]
  compstd[4,1]<-(sqrt(diag(robustase1))/sqrt(n))[3]
  compstd[3,2]<-(sqrt(diag(robustase1))/sqrt(n))[4]
  compstd[3,1]<-(sqrt(diag(robustase1))/sqrt(n))[2]
  compstd[2,1]<-(sqrt(diag(robustase1))/sqrt(n))[1]
  
  asecompstd<-asecompstd+compstd
}


#summarize the mean, standard error and corresponding intervals for the marginal and dependence parameters
mean(cmarginalbeta01)
mean(cmarginalbeta02)
mean(cmarginalbeta03)
mean(cmarginalbeta04)
mean(cmarginalbeta11)
mean(cmarginalbeta12)
mean(cmarginalbeta13)
mean(cmarginalbeta14)
mean(cmarginalbeta21)
mean(cmarginalbeta22)
mean(cmarginalbeta23)
mean(cmarginalbeta24)
mean(csigma1)
mean(csigma2)
mean(csigma3)
mean(csigma4)
mean(csdmarginalbeta01)
mean(csdmarginalbeta02)
mean(csdmarginalbeta03)
mean(csdmarginalbeta04)
mean(csdmarginalbeta11)
mean(csdmarginalbeta12)
mean(csdmarginalbeta13)
mean(csdmarginalbeta14)
mean(csdmarginalbeta21)
mean(csdmarginalbeta22)
mean(csdmarginalbeta23)
mean(csdmarginalbeta24)
mean(csdsigma1)
mean(csdsigma2)
mean(csdsigma3)
mean(csdsigma4)
mean(cpar12)
mean(cpar13)
mean(cpar14)
mean(cpar23)
mean(cpar24)
mean(cpar34)
asecompstd/repeattime

sd(cmarginalbeta01)
sd(cmarginalbeta02)
sd(cmarginalbeta03)
sd(cmarginalbeta04)
sd(cmarginalbeta11)
sd(cmarginalbeta12)
sd(cmarginalbeta13)
sd(cmarginalbeta14)
sd(cmarginalbeta21)
sd(cmarginalbeta22)
sd(cmarginalbeta23)
sd(cmarginalbeta24)
sd(csigma1)
sd(csigma2)
sd(csigma3)
sd(csigma4)

sd(cpar12)
sd(cpar13)
sd(cpar14)
sd(cpar23)
sd(cpar24)
sd(cpar34)


cibeta01up<-cmarginalbeta01+1.96*sd(cmarginalbeta01)
cibeta01lo<-cmarginalbeta01-1.96*sd(cmarginalbeta01)
count<-0
for(i in 1:repeattime)
  if(cibeta01up[i]<1 | cibeta01lo[i]>1){
    count=count+1
  }
ecpbeta01<-(repeattime-count)/repeattime
ecpbeta01

cibeta02up<-cmarginalbeta02+1.96*sd(cmarginalbeta02)
cibeta02lo<-cmarginalbeta02-1.96*sd(cmarginalbeta02)
count<-0
for(i in 1:repeattime)
  if(cibeta02up[i]<2 | cibeta02lo[i]>2){
    count=count+1
  }
ecpbeta02<-(repeattime-count)/repeattime
ecpbeta02

cibeta03up<-cmarginalbeta03+1.96*sd(cmarginalbeta03)
cibeta03lo<-cmarginalbeta03-1.96*sd(cmarginalbeta03)
count<-0
for(i in 1:repeattime)
  if(cibeta03up[i]<3 | cibeta03lo[i]>3){
    count=count+1
  }
ecpbeta03<-(repeattime-count)/repeattime
ecpbeta03

cibeta04up<-cmarginalbeta04+1.96*sd(cmarginalbeta04)
cibeta04lo<-cmarginalbeta04-1.96*sd(cmarginalbeta04)
count<-0
for(i in 1:repeattime)
  if(cibeta04up[i]<4 | cibeta04lo[i]>4){
    count=count+1
  }
ecpbeta04<-(repeattime-count)/repeattime
ecpbeta04

cibeta11up<-cmarginalbeta11+1.96*sd(cmarginalbeta11)
cibeta11lo<-cmarginalbeta11-1.96*sd(cmarginalbeta11)
count<-0
for(i in 1:repeattime)
  if(cibeta11up[i]<2 | cibeta11lo[i]>2){
    count=count+1
  }
ecpbeta11<-(repeattime-count)/repeattime
ecpbeta11

cibeta12up<-cmarginalbeta12+1.96*sd(cmarginalbeta12)
cibeta12lo<-cmarginalbeta12-1.96*sd(cmarginalbeta12)
count<-0
for(i in 1:repeattime)
  if(cibeta12up[i]<3 | cibeta12lo[i]>3){
    count=count+1
  }
ecpbeta12<-(repeattime-count)/repeattime
ecpbeta12

cibeta13up<-cmarginalbeta13+1.96*sd(cmarginalbeta13)
cibeta13lo<-cmarginalbeta13-1.96*sd(cmarginalbeta13)
count<-0
for(i in 1:repeattime)
  if(cibeta13up[i]<4 | cibeta13lo[i]>4){
    count=count+1
  }
ecpbeta13<-(repeattime-count)/repeattime
ecpbeta13

cibeta14up<-cmarginalbeta14+1.96*sd(cmarginalbeta14)
cibeta14lo<-cmarginalbeta14-1.96*sd(cmarginalbeta14)
count<-0
for(i in 1:repeattime)
  if(cibeta14up[i]<5 | cibeta14lo[i]>5){
    count=count+1
  }
ecpbeta14<-(repeattime-count)/repeattime
ecpbeta14

cibeta21up<-cmarginalbeta21+1.96*sd(cmarginalbeta21)
cibeta21lo<-cmarginalbeta21-1.96*sd(cmarginalbeta21)
count<-0
for(i in 1:repeattime)
  if(cibeta21up[i]<3 | cibeta21lo[i]>3){
    count=count+1
  }
ecpbeta21<-(repeattime-count)/repeattime
ecpbeta21

cibeta22up<-cmarginalbeta22+1.96*sd(cmarginalbeta22)
cibeta22lo<-cmarginalbeta22-1.96*sd(cmarginalbeta22)
count<-0
for(i in 1:repeattime)
  if(cibeta22up[i]<4 | cibeta22lo[i]>4){
    count=count+1
  }
ecpbeta22<-(repeattime-count)/repeattime
ecpbeta22

cibeta23up<-cmarginalbeta23+1.96*sd(cmarginalbeta23)
cibeta23lo<-cmarginalbeta23-1.96*sd(cmarginalbeta23)
count<-0
for(i in 1:repeattime)
  if(cibeta23up[i]<5 | cibeta23lo[i]>5){
    count=count+1
  }
ecpbeta23<-(repeattime-count)/repeattime
ecpbeta23

cibeta24up<-cmarginalbeta24+1.96*sd(cmarginalbeta24)
cibeta24lo<-cmarginalbeta24-1.96*sd(cmarginalbeta24)
count<-0
for(i in 1:repeattime)
  if(cibeta24up[i]<6 | cibeta24lo[i]>6){
    count=count+1
  }
ecpbeta24<-(repeattime-count)/repeattime
ecpbeta24

cisigma1up<-csigma1+1.96*sd(csigma1)
cisigma1lo<-csigma1-1.96*sd(csigma1)
count<-0
for(i in 1:repeattime)
  if(cisigma1up[i]<2 | cisigma1lo[i]>2){
    count=count+1
  }
ecpsigma1<-(repeattime-count)/repeattime
ecpsigma1

cisigma2up<-csigma2+1.96*sd(csigma2)
cisigma2lo<-csigma2-1.96*sd(csigma2)
count<-0
for(i in 1:repeattime)
  if(cisigma2up[i]<2 | cisigma2lo[i]>2){
    count=count+1
  }
ecpsigma2<-(repeattime-count)/repeattime
ecpsigma2

cisigma3up<-csigma3+1.96*sd(csigma3)
cisigma3lo<-csigma3-1.96*sd(csigma3)
count<-0
for(i in 1:repeattime)
  if(cisigma3up[i]<2 | cisigma3lo[i]>2){
    count=count+1
  }
ecpsigma3<-(repeattime-count)/repeattime
ecpsigma3

cisigma4up<-csigma4+1.96*sd(csigma4)
cisigma4lo<-csigma4-1.96*sd(csigma4)
count<-0
for(i in 1:repeattime)
  if(cisigma4up[i]<2 | cisigma4lo[i]>2){
    count=count+1
  }
ecpsigma4<-(repeattime-count)/repeattime
ecpsigma4


ci12upper<-par12+1.96*sd(par12)
ci12lower<-par12-1.96*sd(par12)
count<-0
for(i in 1:repeattime)
  if(ci12upper[i]<2 | ci12lower[i]>2){
    count=count+1
  }
ecpfull12<-(repeattime-count)/repeattime
ecpfull12

ci13upper<-par13+1.96*sd(par13)
ci13lower<-par13-1.96*sd(par13)
count<-0
for(i in 1:repeattime)
  if(ci13upper[i]<2 | ci13lower[i]>2){
    count=count+1
  }
ecpfull13<-(repeattime-count)/repeattime
ecpfull13

ci14upper<-par14+1.96*sd(par14)
ci14lower<-par14-1.96*sd(par14)
count<-0
for(i in 1:repeattime)
  if(ci14upper[i]<0.5 | ci14lower[i]>0.5){
    count=count+1
  }
ecpfull14<-(repeattime-count)/repeattime
ecpfull14

ci23upper<-par23+1.96*sd(par23)
ci23lower<-par23-1.96*sd(par23)
count<-0
for(i in 1:repeattime)
  if(ci23upper[i]<3 | ci23lower[i]>3){
    count=count+1
  }
ecpfull23<-(repeattime-count)/repeattime
ecpfull23

ci24upper<-par24+1.96*sd(par24)
ci24lower<-par24-1.96*sd(par24)
count<-0
for(i in 1:repeattime)
  if(ci24upper[i]<3 | ci24lower[i]>3){
    count=count+1
  }
ecpfull24<-(repeattime-count)/repeattime
ecpfull24

ci34upper<-par34+1.96*sd(par34)
ci34lower<-par34-1.96*sd(par34)
count<-0
for(i in 1:repeattime)
  if(ci34upper[i]<2 | ci34lower[i]>2){
    count=count+1
  }
ecpfull34<-(repeattime-count)/repeattime
ecpfull34
