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<-matrix(RVinem,16,16)

#Rvinerobust and RVineparrobust are the true family matrix and parameters matrix, where the copula functions connecting time blocks are all Gumbel(3)
RVinerobust<-c(rep(0,13),3,4,1,rep(0,13),3,4,1,rep(0,13),4,3,4,rep(0,13),4,3,4,rep(0,13),4,4,5,rep(0,13),4,4,5,rep(0,13),4,4,4,rep(0,13),4,4,4,rep(0,13),3,4,1,rep(0,13),3,4,1,rep(0,13),4,3,4,rep(0,13),4,3,4,rep(0,13),4,4,5,rep(0,14),4,5,rep(0,15),4,rep(0,16))
Rvinerobust<-matrix(RVinerobust,16,16)
RVineparrobust<-c(rep(0,13),2,3,0.5,rep(0,13),2,3,0.5,rep(0,13),3,3,2,rep(0,13),3,3,2,rep(0,13),3,3,2,rep(0,13),3,3,2,rep(0,13),3,3,3,rep(0,13),3,3,3,rep(0,13),2,3,0.5,rep(0,13),2,3,0.5,rep(0,13),3,3,2,rep(0,13),3,3,2,rep(0,13),3,3,2,rep(0,14),3,2,rep(0,15),3,rep(0,16))
Rvineparrobust<-matrix(RVineparrobust,16,16)

RVinepar2<-rep(0,256)
Rvinepar2<-matrix(RVinepar2,16,16)

RVMrobust<-RVineMatrix(Rvinem,Rvinerobust,Rvineparrobust,Rvinepar2,check.pars=TRUE)

#Rvinef is the misspecified family matrix used to fit the model, where the copula functions connecting time blocks are specified as Gaussian
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<-matrix(RVinef,16,16)

n<-500
k<-4
l<-4

repeattime<-500

#lower
lo<-c(-Inf,-Inf,-Inf,0)
lo1<-c(0,1,-1,-1,0,1,rep(-1,3),-Inf,rep(-1,14))

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


marginalbeta01<-rep(0,repeattime)
marginalbeta02<-rep(0,repeattime)
marginalbeta03<-rep(0,repeattime)
marginalbeta04<-rep(0,repeattime)
marginalbeta11<-rep(0,repeattime)
marginalbeta12<-rep(0,repeattime)
marginalbeta13<-rep(0,repeattime)
marginalbeta14<-rep(0,repeattime)
marginalbeta21<-rep(0,repeattime)
marginalbeta22<-rep(0,repeattime)
marginalbeta23<-rep(0,repeattime)
marginalbeta24<-rep(0,repeattime)

sigma1<-rep(0,repeattime)
sigma2<-rep(0,repeattime)
sigma3<-rep(0,repeattime)
sigma4<-rep(0,repeattime)
sdmarginalbeta01<-rep(0,repeattime)
sdmarginalbeta02<-rep(0,repeattime)
sdmarginalbeta03<-rep(0,repeattime)
sdmarginalbeta04<-rep(0,repeattime)
sdmarginalbeta11<-rep(0,repeattime)
sdmarginalbeta12<-rep(0,repeattime)
sdmarginalbeta13<-rep(0,repeattime)
sdmarginalbeta14<-rep(0,repeattime)
sdmarginalbeta21<-rep(0,repeattime)
sdmarginalbeta22<-rep(0,repeattime)
sdmarginalbeta23<-rep(0,repeattime)
sdmarginalbeta24<-rep(0,repeattime)

sdsigma1<-rep(0,repeattime)
sdsigma2<-rep(0,repeattime)
sdsigma3<-rep(0,repeattime)
sdsigma4<-rep(0,repeattime)

par12<-rep(0,repeattime)
par13<-rep(0,repeattime)
par14<-rep(0,repeattime)
par23<-rep(0,repeattime)
par24<-rep(0,repeattime)
par34<-rep(0,repeattime)
mlefullpar<-matrix(0,4,4)
esefullstd<-matrix(0,4,4)
asefullstd<-matrix(0,4,4)

for(i in 1:repeattime){
  set.seed(i)
  #generate the data from the true structure RVMrobust
  simulate1<-RVineSim(n,RVMrobust)
  
  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)
  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)
  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, here we use the misspecified copula family matrix Rvinef)
  
  optimfun5<-function(param){
    RVinepar3<-c(rep(0,13),param[1:3],rep(0,13),param[1],param[2],param[3],rep(0,13),param[4:6],rep(0,13),param[7],param[5],param[6],rep(0,13),param[8:10],rep(0,13),param[11],param[12],param[10],rep(0,13),param[13:15],rep(0,13),param[16:18],rep(0,13),param[1],param[2],param[3],rep(0,13),param[1],param[2],param[3],rep(0,13),param[19],param[5],param[6],rep(0,13),param[20],param[5],param[6],rep(0,13),param[21],param[22],param[10],rep(0,14),param[23],param[10],rep(0,15),param[24],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,0.5,3,2,0.5,0.5,0.5,2,0.5,0.5,rep(0.5,6),0.5,0.5,0.5,0.5,0.5,0.5)
  out5<-nlminb(objective=optimfun5,start=startpar5,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo1,upper=up1)
  
  
  
  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)
  H1<-out$Hessian
  
  
  optimfun7<-function(param){
    RVinepar3<-c(rep(0,13),param[1:3],rep(0,13),param[1],param[2],param[3],rep(0,13),param[4:6],rep(0,13),param[7],param[5],param[6],rep(0,13),param[8:10],rep(0,13),param[11],param[12],param[10],rep(0,13),param[13:15],rep(0,13),param[16:18],rep(0,13),param[1],param[2],param[3],rep(0,13),param[1],param[2],param[3],rep(0,13),param[19],param[5],param[6],rep(0,13),param[20],param[5],param[6],rep(0,13),param[21],param[22],param[10],rep(0,14),param[23],param[10],rep(0,15),param[24],rep(0,16))
    Rvinepar3<-matrix(RVinepar3,16,16)
    RVM2<-RVineMatrix(Rvinem,Rvinef,Rvinepar3,Rvinepar2,check.pars=FALSE)
    
    
    u1<-pnorm((y1[1:n]-param[25]-param[26]*x1[1:n]-param[27]),mean=0,sd=param[28])
    u2<-pnorm((y2[1:n]-param[29]-param[30]*x2[1:n]-param[31]),mean=0,sd=param[32])
    u3<-pnorm((y3[1:n]-param[33]-param[34]*x3[1:n]-param[35]),mean=0,sd=param[36])
    u4<-pnorm((y4[1:n]-param[37]-param[38]*x4[1:n]-param[39]),mean=0,sd=param[40])
    u5<-pnorm((y1[(n+1):(2*n)]-param[25]-param[26]*x1[(n+1):(2*n)]-param[27]*2),mean=0,sd=param[28])
    u6<-pnorm((y2[(n+1):(2*n)]-param[29]-param[30]*x2[(n+1):(2*n)]-param[31]*2),mean=0,sd=param[32])
    u7<-pnorm((y3[(n+1):(2*n)]-param[33]-param[34]*x3[(n+1):(2*n)]-param[35]*2),mean=0,sd=param[36])
    u8<-pnorm((y4[(n+1):(2*n)]-param[37]-param[38]*x4[(n+1):(2*n)]-param[39]*2),mean=0,sd=param[40])
    u9<-pnorm((y1[(2*n+1):(3*n)]-param[25]-param[26]*x1[(2*n+1):(3*n)]-param[27]*3),mean=0,sd=param[28])
    u10<-pnorm((y2[(2*n+1):(3*n)]-param[29]-param[30]*x2[(2*n+1):(3*n)]-param[31]*3),mean=0,sd=param[32])
    u11<-pnorm((y3[(2*n+1):(3*n)]-param[33]-param[34]*x3[(2*n+1):(3*n)]-param[35]*3),mean=0,sd=param[36])
    u12<-pnorm((y4[(2*n+1):(3*n)]-param[37]-param[38]*x4[(2*n+1):(3*n)]-param[39]*3),mean=0,sd=param[40])
    u13<-pnorm((y1[(3*n+1):(4*n)]-param[25]-param[26]*x1[(3*n+1):(4*n)]-param[27]*4),mean=0,sd=param[28])
    u14<-pnorm((y2[(3*n+1):(4*n)]-param[29]-param[30]*x2[(3*n+1):(4*n)]-param[31]*4),mean=0,sd=param[32])
    u15<-pnorm((y3[(3*n+1):(4*n)]-param[33]-param[34]*x3[(3*n+1):(4*n)]-param[35]*4),mean=0,sd=param[36])
    u16<-pnorm((y4[(3*n+1):(4*n)]-param[37]-param[38]*x4[(3*n+1):(4*n)]-param[39]*4),mean=0,sd=param[40])
    
    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)
  H2<-out$Hessian
  
  
  
  
  #calculate omega matrix
  Q<-H2[1:24,25:40]%*%inv(H1)
  
  J<-matrix(0,16,16)
  omega<-matrix(0,24,24)
  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,13),param[4:6],rep(0,13),param[7],param[5],param[6],rep(0,13),param[8:10],rep(0,13),param[11],param[12],param[10],rep(0,13),param[13:15],rep(0,13),param[16:18],rep(0,13),param[1],param[2],param[3],rep(0,13),param[1],param[2],param[3],rep(0,13),param[19],param[5],param[6],rep(0,13),param[20],param[5],param[6],rep(0,13),param[21],param[22],param[10],rep(0,14),param[23],param[10],rep(0,15),param[24],rep(0,16))
      Rvinepar3<-matrix(RVinepar3,16,16)
      RVM2<-RVineMatrix(Rvinem,Rvinef,Rvinepar3,Rvinepar2,check.pars=FALSE)
      
      
      logp1<-dnorm((y1[r]-param[25]-param[26]*x1[r]-param[27]),mean=0,sd=param[28],log=TRUE)+dnorm((y1[n+r]-param[25]-param[26]*x1[n+r]-param[27]*2),mean=0,sd=param[28],log=TRUE)+dnorm((y1[2*n+r]-param[25]-param[26]*x1[2*n+r]-param[27]*3),mean=0,sd=param[28],log=TRUE)+dnorm((y1[3*n+r]-param[25]-param[26]*x1[3*n+r]-param[27]*4),mean=0,sd=param[28],log=TRUE)
      logp2<-dnorm((y2[r]-param[29]-param[30]*x2[r]-param[31]),mean=0,sd=param[32],log=TRUE)+dnorm((y2[n+r]-param[29]-param[30]*x2[n+r]-param[31]*2),mean=0,sd=param[32],log=TRUE)+dnorm((y2[2*n+r]-param[29]-param[30]*x2[2*n+r]-param[31]*3),mean=0,sd=param[32],log=TRUE)+dnorm((y2[3*n+r]-param[29]-param[30]*x2[3*n+r]-param[31]*4),mean=0,sd=param[32],log=TRUE)
      logp3<-dnorm((y3[r]-param[33]-param[34]*x3[r]-param[35]),mean=0,sd=param[36],log=TRUE)+dnorm((y3[n+r]-param[33]-param[34]*x3[n+r]-param[35]*2),mean=0,sd=param[36],log=TRUE)+dnorm((y3[2*n+r]-param[33]-param[34]*x3[2*n+r]-param[35]*3),mean=0,sd=param[36],log=TRUE)+dnorm((y3[3*n+r]-param[33]-param[34]*x3[3*n+r]-param[35]*4),mean=0,sd=param[36],log=TRUE)
      logp4<-dnorm((y4[r]-param[37]-param[38]*x4[r]-param[39]),mean=0,sd=param[40],log=TRUE)+dnorm((y4[n+r]-param[37]-param[38]*x4[n+r]-param[39]*2),mean=0,sd=param[40],log=TRUE)+dnorm((y4[2*n+r]-param[37]-param[38]*x4[2*n+r]-param[39]*3),mean=0,sd=param[40],log=TRUE)+dnorm((y4[3*n+r]-param[37]-param[38]*x4[3*n+r]-param[39]*4),mean=0,sd=param[40],log=TRUE)
      
      
      u1<-pnorm((y1[r]-param[25]-param[26]*x1[r]-param[27]),mean=0,sd=param[28])
      u2<-pnorm((y2[r]-param[29]-param[30]*x2[r]-param[31]),mean=0,sd=param[32])
      u3<-pnorm((y3[r]-param[33]-param[34]*x3[r]-param[35]),mean=0,sd=param[36])
      u4<-pnorm((y4[r]-param[37]-param[38]*x4[r]-param[39]),mean=0,sd=param[40])
      u5<-pnorm((y1[n+r]-param[25]-param[26]*x1[n+r]-param[27]*2),mean=0,sd=param[28])
      u6<-pnorm((y2[n+r]-param[29]-param[30]*x2[n+r]-param[31]*2),mean=0,sd=param[32])
      u7<-pnorm((y3[n+r]-param[33]-param[34]*x3[n+r]-param[35]*2),mean=0,sd=param[36])
      u8<-pnorm((y4[n+r]-param[37]-param[38]*x4[n+r]-param[39]*2),mean=0,sd=param[40])
      u9<-pnorm((y1[2*n+r]-param[25]-param[26]*x1[2*n+r]-param[27]*3),mean=0,sd=param[28])
      u10<-pnorm((y2[2*n+r]-param[29]-param[30]*x2[2*n+r]-param[31]*3),mean=0,sd=param[32])
      u11<-pnorm((y3[2*n+r]-param[33]-param[34]*x3[2*n+r]-param[35]*3),mean=0,sd=param[36])
      u12<-pnorm((y4[2*n+r]-param[37]-param[38]*x4[2*n+r]-param[39]*3),mean=0,sd=param[40])
      u13<-pnorm((y1[3*n+r]-param[25]-param[26]*x1[3*n+r]-param[27]*4),mean=0,sd=param[28])
      u14<-pnorm((y2[3*n+r]-param[29]-param[30]*x2[3*n+r]-param[31]*4),mean=0,sd=param[32])
      u15<-pnorm((y3[3*n+r]-param[33]-param[34]*x3[3*n+r]-param[35]*4),mean=0,sd=param[36])
      u16<-pnorm((y4[3*n+r]-param[37]-param[38]*x4[3*n+r]-param[39]*4),mean=0,sd=param[40])
      
      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:24]-Q%*%out6
    omega<-omega+ome%*%t(ome)
  }
  
  omega<-omega/n
  H3<-H2/n
  robustase1<-inv(H3[1:24,1:24])%*%omega%*%inv(H3[1:24,1:24])
  
  
  if(any(is.nan(sqrt(diag(robustase1)))==TRUE)){
    next
  }else{
    
    par12[i]<-out5$par[10]
    par13[i]<-out5$par[6]
    par14[i]<-out5$par[3]
    par23[i]<-out5$par[5]
    par24[i]<-out5$par[2]
    par34[i]<-out5$par[1]
    mlefullpar[4,3]<-mlefullpar[4,3]+out5$par[10]
    mlefullpar[4,2]<-mlefullpar[4,2]+out5$par[6]
    mlefullpar[4,1]<-mlefullpar[4,1]+out5$par[3]
    mlefullpar[3,2]<-mlefullpar[3,2]+out5$par[5]
    mlefullpar[3,1]<-mlefullpar[3,1]+out5$par[2]
    mlefullpar[2,1]<-mlefullpar[2,1]+out5$par[1]
    
    
    
    
    Fullstd<-matrix(0,4,4)
    Fullstd[4,3]<-(sqrt(diag(robustase1))/sqrt(n))[10]
    Fullstd[4,2]<-(sqrt(diag(robustase1))/sqrt(n))[6]
    Fullstd[4,1]<-(sqrt(diag(robustase1))/sqrt(n))[3]
    Fullstd[3,2]<-(sqrt(diag(robustase1))/sqrt(n))[5]
    Fullstd[3,1]<-(sqrt(diag(robustase1))/sqrt(n))[2]
    Fullstd[2,1]<-(sqrt(diag(robustase1))/sqrt(n))[1]
    
    asefullstd<-asefullstd+Fullstd
  }
  
  
}