library(VineCopula)
library(Matrix)
library(matrixcalc)
library(matlib)
library(maxLik)
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,-1.0,1,rep(-1,3),-Inf,rep(-1,14))

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


#set up the matrix that record the simulation results
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)

#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)
  #this is the likelihood function for marginal parameters with negative sign (so that the nlminb function can minimize it)
  optimfun<-function(param){
  
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-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)
      logp2<-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<-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<-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)
  }
  
  
  
  #start parameters
  startpar<-c(beta01,beta11,beta21,sigma,beta02,beta12,beta22,sigma,beta03,beta13,beta23,sigma,beta04,beta14,beta24,sigma)
  out1<-nlminb(objective=optimfun,start=startpar,control=list(trace=10,eval.max=5000,iter.max=5000,rel.tol=1e-12),lower=lo,upper=up)
  
  
  #calculate the Hessian matrix and J matrix of the marginal parameters
  H<-matrix(0,16,16)
  J<-matrix(0,16,16)
  for(k in 1:n){
    
  #this is the likelihood function without the negative sign
  optimfun2<-function(param){
    
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    
      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<-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<-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<-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(out1$par,optimfun2)
  H<-H+out$Hessian
  J<-J+out$gradient%*%t(out$gradient)
  }
  H<-H/n
  J<-J/n
  robustase<-inv(H)%*%J%*%inv(H)
  
  #record the results for marginal parameters
  marginalbeta01[i]<-out1$par[1]
  marginalbeta02[i]<-out1$par[5]
  marginalbeta03[i]<-out1$par[9]
  marginalbeta04[i]<-out1$par[13]
  marginalbeta11[i]<-out1$par[2]
  marginalbeta12[i]<-out1$par[6]
  marginalbeta13[i]<-out1$par[10]
  marginalbeta14[i]<-out1$par[14]
  marginalbeta21[i]<-out1$par[3]
  marginalbeta22[i]<-out1$par[7]
  marginalbeta23[i]<-out1$par[11]
  marginalbeta24[i]<-out1$par[15]
  sigma1[i]<-out1$par[4]
  sigma2[i]<-out1$par[8]
  sigma3[i]<-out1$par[12]
  sigma4[i]<-out1$par[16]
  

  sdmarginalbeta01[i]<-(sqrt(diag(robustase))/sqrt(n))[1]
  sdmarginalbeta02[i]<-(sqrt(diag(robustase))/sqrt(n))[5]
  sdmarginalbeta03[i]<-(sqrt(diag(robustase))/sqrt(n))[9]
  sdmarginalbeta04[i]<-(sqrt(diag(robustase))/sqrt(n))[13]
  sdmarginalbeta11[i]<-(sqrt(diag(robustase))/sqrt(n))[2]
  sdmarginalbeta12[i]<-(sqrt(diag(robustase))/sqrt(n))[6]
  sdmarginalbeta13[i]<-(sqrt(diag(robustase))/sqrt(n))[10]
  sdmarginalbeta14[i]<-(sqrt(diag(robustase))/sqrt(n))[14]
  sdmarginalbeta21[i]<-(sqrt(diag(robustase))/sqrt(n))[3]
  sdmarginalbeta22[i]<-(sqrt(diag(robustase))/sqrt(n))[7]
  sdmarginalbeta23[i]<-(sqrt(diag(robustase))/sqrt(n))[11]
  sdmarginalbeta24[i]<-(sqrt(diag(robustase))/sqrt(n))[15]
  sdsigma1[i]<-(sqrt(diag(robustase))/sqrt(n))[4]
  sdsigma2[i]<-(sqrt(diag(robustase))/sqrt(n))[8]
  sdsigma3[i]<-(sqrt(diag(robustase))/sqrt(n))[12]
  sdsigma4[i]<-(sqrt(diag(robustase))/sqrt(n))[16]
  
  
  
  #estimate the vine parameters(second stage estimation)
  
  #this is the likelihood function for the whole vine structure with the negative sign
  optimfun1<-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]-out1$par[5]-out1$par[6]*x2[1:n]-out1$par[7]),mean=0,sd=out1$par[8])
    u3<-pnorm((y3[1:n]-out1$par[9]-out1$par[10]*x3[1:n]-out1$par[11]),mean=0,sd=out1$par[12])
    u4<-pnorm((y4[1:n]-out1$par[13]-out1$par[14]*x4[1:n]-out1$par[15]),mean=0,sd=out1$par[16])
    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)]-out1$par[5]-out1$par[6]*x2[(n+1):(2*n)]-out1$par[7]*2),mean=0,sd=out1$par[8])
    u7<-pnorm((y3[(n+1):(2*n)]-out1$par[9]-out1$par[10]*x3[(n+1):(2*n)]-out1$par[11]*2),mean=0,sd=out1$par[12])
    u8<-pnorm((y4[(n+1):(2*n)]-out1$par[13]-out1$par[14]*x4[(n+1):(2*n)]-out1$par[15]*2),mean=0,sd=out1$par[16])
    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)]-out1$par[5]-out1$par[6]*x2[(2*n+1):(3*n)]-out1$par[7]*3),mean=0,sd=out1$par[8])
    u11<-pnorm((y3[(2*n+1):(3*n)]-out1$par[9]-out1$par[10]*x3[(2*n+1):(3*n)]-out1$par[11]*3),mean=0,sd=out1$par[12])
    u12<-pnorm((y4[(2*n+1):(3*n)]-out1$par[13]-out1$par[14]*x4[(2*n+1):(3*n)]-out1$par[15]*3),mean=0,out1$par[16])
    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)]-out1$par[5]-out1$par[6]*x2[(3*n+1):(4*n)]-out1$par[7]*4),mean=0,sd=out1$par[8])
    u15<-pnorm((y3[(3*n+1):(4*n)]-out1$par[9]-out1$par[10]*x3[(3*n+1):(4*n)]-out1$par[11]*4),mean=0,sd=out1$par[12])
    u16<-pnorm((y4[(3*n+1):(4*n)]-out1$par[13]-out1$par[14]*x4[(3*n+1):(4*n)]-out1$par[15]*4),mean=0,sd=out1$par[16])
    
    
    
    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
  }
  
  #minimize the likelihood function with negative sign
  startpar1<-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)
  out2<-nlminb(objective=optimfun1,start=startpar1,control=list(trace=10,eval.max=5000,iter.max=5000,rel.tol=1e-12),lower=lo1,upper=up1)
  
  
  #optimfun3 is the likelihood function used to calculate the Hessian matrix for the dependence parameters. Since in the second stage, we use the full likelihood(the whole vine structure), we don't need to calculate the J matrix, and just use the inverse of the hessian matrix
  optimfun3<-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]-out1$par[5]-out1$par[6]*x2[1:n]-out1$par[7]),mean=0,sd=out1$par[8])
    u3<-pnorm((y3[1:n]-out1$par[9]-out1$par[10]*x3[1:n]-out1$par[11]),mean=0,sd=out1$par[12])
    u4<-pnorm((y4[1:n]-out1$par[13]-out1$par[14]*x4[1:n]-out1$par[15]),mean=0,sd=out1$par[16])
    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)]-out1$par[5]-out1$par[6]*x2[(n+1):(2*n)]-out1$par[7]*2),mean=0,sd=out1$par[8])
    u7<-pnorm((y3[(n+1):(2*n)]-out1$par[9]-out1$par[10]*x3[(n+1):(2*n)]-out1$par[11]*2),mean=0,sd=out1$par[12])
    u8<-pnorm((y4[(n+1):(2*n)]-out1$par[13]-out1$par[14]*x4[(n+1):(2*n)]-out1$par[15]*2),mean=0,sd=out1$par[16])
    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)]-out1$par[5]-out1$par[6]*x2[(2*n+1):(3*n)]-out1$par[7]*3),mean=0,sd=out1$par[8])
    u11<-pnorm((y3[(2*n+1):(3*n)]-out1$par[9]-out1$par[10]*x3[(2*n+1):(3*n)]-out1$par[11]*3),mean=0,sd=out1$par[12])
    u12<-pnorm((y4[(2*n+1):(3*n)]-out1$par[13]-out1$par[14]*x4[(2*n+1):(3*n)]-out1$par[15]*3),mean=0,out1$par[16])
    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)]-out1$par[5]-out1$par[6]*x2[(3*n+1):(4*n)]-out1$par[7]*4),mean=0,sd=out1$par[8])
    u15<-pnorm((y3[(3*n+1):(4*n)]-out1$par[9]-out1$par[10]*x3[(3*n+1):(4*n)]-out1$par[11]*4),mean=0,sd=out1$par[12])
    u16<-pnorm((y4[(3*n+1):(4*n)]-out1$par[13]-out1$par[14]*x4[(3*n+1):(4*n)]-out1$par[15]*4),mean=0,sd=out1$par[16])
    
    
    
    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
  }
  
  out3<-fdHess(out2$par,optimfun3)
  H1<-out3$Hessian/n
  
  
  #record the results for the dependence parameters
  par12[i]<-out2$par[10]
  par13[i]<-out2$par[6]
  par14[i]<-out2$par[3]
  par23[i]<-out2$par[5]
  par24[i]<-out2$par[2]
  par34[i]<-out2$par[1]
  mlefullpar[4,3]<-mlefullpar[4,3]+out2$par[10]
  mlefullpar[4,2]<-mlefullpar[4,2]+out2$par[6]
  mlefullpar[4,1]<-mlefullpar[4,1]+out2$par[3]
  mlefullpar[3,2]<-mlefullpar[3,2]+out2$par[5]
  mlefullpar[3,1]<-mlefullpar[3,1]+out2$par[2]
  mlefullpar[2,1]<-mlefullpar[2,1]+out2$par[1]
  
  
  
  
  Fullstd<-matrix(0,4,4)
  Fullstd[4,3]<-(sqrt(diag(inv(-H1)))/sqrt(n))[10]
  Fullstd[4,2]<-(sqrt(diag(inv(-H1)))/sqrt(n))[6]
  Fullstd[4,1]<-(sqrt(diag(inv(-H1)))/sqrt(n))[3]
  Fullstd[3,2]<-(sqrt(diag(inv(-H1)))/sqrt(n))[5]
  Fullstd[3,1]<-(sqrt(diag(inv(-H1)))/sqrt(n))[2]
  Fullstd[2,1]<-(sqrt(diag(inv(-H1)))/sqrt(n))[1]
  
  asefullstd<-asefullstd+Fullstd
}
