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

#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

# vine structure within a time block
RVineblock<-c(4,3,2,1,0,3,2,1,0,0,2,1,0,0,0,1)
Rvineblock<-matrix(RVineblock,4,4)

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

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


#the upper and lower bounds for the corresponding copula parameters
#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))


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


#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)]
  }
  
  #start parameters
  startpar1<-c(2,3,0.5,3,2,2,beta01,beta11,beta21,sigma,beta02,beta12,beta22,sigma,beta03,beta13,beta23,sigma,beta04,beta14,beta24,sigma)
  
  
  #this is the composite likelihood function, where the copula parameters outside the time blocks are all specified to be 0, so we only have 6 parameters in the RVinepar3
  optimfun1<-function(param=startpar1){
    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,Rvinef1,Rvinepar3,Rvinepar2,check.pars=FALSE)
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(r in 1:n){
      logp1<-logp1+dnorm((y1[r]-param[7]-param[8]*x1[r]-param[9]),mean=0,sd=param[10],log=TRUE)+dnorm((y1[n+r]-param[7]-param[8]*x1[n+r]-param[9]*2),mean=0,sd=param[10],log=TRUE)+dnorm((y1[2*n+r]-param[7]-param[8]*x1[2*n+r]-param[9]*3),mean=0,sd=param[10],log=TRUE)+dnorm((y1[3*n+r]-param[7]-param[8]*x1[3*n+r]-param[9]*4),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((y2[r]-param[11]-param[12]*x2[r]-param[13]),mean=0,sd=param[14],log=TRUE)+dnorm((y2[n+r]-param[11]-param[12]*x2[n+r]-param[13]*2),mean=0,sd=param[14],log=TRUE)+dnorm((y2[2*n+r]-param[11]-param[12]*x2[2*n+r]-param[13]*3),mean=0,sd=param[14],log=TRUE)+dnorm((y2[3*n+r]-param[11]-param[12]*x2[3*n+r]-param[13]*4),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((y3[r]-param[15]-param[16]*x3[r]-param[17]),mean=0,sd=param[18],log=TRUE)+dnorm((y3[n+r]-param[15]-param[16]*x3[n+r]-param[17]*2),mean=0,sd=param[18],log=TRUE)+dnorm((y3[2*n+r]-param[15]-param[16]*x3[2*n+r]-param[17]*3),mean=0,sd=param[18],log=TRUE)+dnorm((y3[3*n+r]-param[15]-param[16]*x3[3*n+r]-param[17]*4),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((y4[r]-param[19]-param[20]*x4[r]-param[21]),mean=0,sd=param[22],log=TRUE)+dnorm((y4[n+r]-param[19]-param[20]*x4[n+r]-param[21]*2),mean=0,sd=param[22],log=TRUE)+dnorm((y4[2*n+r]-param[19]-param[20]*x4[2*n+r]-param[21]*3),mean=0,sd=param[22],log=TRUE)+dnorm((y4[3*n+r]-param[19]-param[20]*x4[3*n+r]-param[21]*4),mean=0,sd=param[22],log=TRUE)
    }
    
    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)
    
    #take negative sign here since the nlminb function minimize the target function
    -(logp1+logp2+logp3+logp4+RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik)
  }
  
  
  out1<-nlminb(objective=optimfun1,start=startpar1,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo,upper=up)
  
  
  #this is the composite likelihood function without the negative sign to calculate the hessian matrix
  optimfun3<-function(param=startpar1){
    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,Rvinef1,Rvinepar3,Rvinepar2,check.pars=FALSE)
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(r in 1:n){
      logp1<-logp1+dnorm((y1[r]-param[7]-param[8]*x1[r]-param[9]),mean=0,sd=param[10],log=TRUE)+dnorm((y1[n+r]-param[7]-param[8]*x1[n+r]-param[9]*2),mean=0,sd=param[10],log=TRUE)+dnorm((y1[2*n+r]-param[7]-param[8]*x1[2*n+r]-param[9]*3),mean=0,sd=param[10],log=TRUE)+dnorm((y1[3*n+r]-param[7]-param[8]*x1[3*n+r]-param[9]*4),mean=0,sd=param[10],log=TRUE)
      logp2<-logp2+dnorm((y2[r]-param[11]-param[12]*x2[r]-param[13]),mean=0,sd=param[14],log=TRUE)+dnorm((y2[n+r]-param[11]-param[12]*x2[n+r]-param[13]*2),mean=0,sd=param[14],log=TRUE)+dnorm((y2[2*n+r]-param[11]-param[12]*x2[2*n+r]-param[13]*3),mean=0,sd=param[14],log=TRUE)+dnorm((y2[3*n+r]-param[11]-param[12]*x2[3*n+r]-param[13]*4),mean=0,sd=param[14],log=TRUE)
      logp3<-logp3+dnorm((y3[r]-param[15]-param[16]*x3[r]-param[17]),mean=0,sd=param[18],log=TRUE)+dnorm((y3[n+r]-param[15]-param[16]*x3[n+r]-param[17]*2),mean=0,sd=param[18],log=TRUE)+dnorm((y3[2*n+r]-param[15]-param[16]*x3[2*n+r]-param[17]*3),mean=0,sd=param[18],log=TRUE)+dnorm((y3[3*n+r]-param[15]-param[16]*x3[3*n+r]-param[17]*4),mean=0,sd=param[18],log=TRUE)
      logp4<-logp4+dnorm((y4[r]-param[19]-param[20]*x4[r]-param[21]),mean=0,sd=param[22],log=TRUE)+dnorm((y4[n+r]-param[19]-param[20]*x4[n+r]-param[21]*2),mean=0,sd=param[22],log=TRUE)+dnorm((y4[2*n+r]-param[19]-param[20]*x4[2*n+r]-param[21]*3),mean=0,sd=param[22],log=TRUE)+dnorm((y4[3*n+r]-param[19]-param[20]*x4[3*n+r]-param[21]*4),mean=0,sd=param[22],log=TRUE)
    }
    
    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)
    
    
    logp1+logp2+logp3+logp4+RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
  }
  H<-fdHess(out1$par,optimfun3)$Hessian
  H<-H/n
  
  
  #here we calculate the variability matrix J for each data point from 1 to n
  J<-matrix(0,22,22)
  for(j in 1:n){
    optimfun2<-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,Rvinef1,Rvinepar3,Rvinepar2,check.pars=FALSE)
      logp1<-dnorm((y1[j]-param[7]-param[8]*x1[j]-param[9]),mean=0,sd=param[10],log=TRUE)+dnorm((y1[n+j]-param[7]-param[8]*x1[n+j]-param[9]*2),mean=0,sd=param[10],log=TRUE)+dnorm((y1[2*n+j]-param[7]-param[8]*x1[2*n+j]-param[9]*3),mean=0,sd=param[10],log=TRUE)+dnorm((y1[3*n+j]-param[7]-param[8]*x1[3*n+j]-param[9]*4),mean=0,sd=param[10],log=TRUE)
      logp2<-dnorm((y2[j]-param[11]-param[12]*x2[j]-param[13]),mean=0,sd=param[14],log=TRUE)+dnorm((y2[n+j]-param[11]-param[12]*x2[n+j]-param[13]*2),mean=0,sd=param[14],log=TRUE)+dnorm((y2[2*n+j]-param[11]-param[12]*x2[2*n+j]-param[13]*3),mean=0,sd=param[14],log=TRUE)+dnorm((y2[3*n+j]-param[11]-param[12]*x2[3*n+j]-param[13]*4),mean=0,sd=param[14],log=TRUE)
      logp3<-dnorm((y3[j]-param[15]-param[16]*x3[j]-param[17]),mean=0,sd=param[18],log=TRUE)+dnorm((y3[n+j]-param[15]-param[16]*x3[n+j]-param[17]*2),mean=0,sd=param[18],log=TRUE)+dnorm((y3[2*n+j]-param[15]-param[16]*x3[2*n+j]-param[17]*3),mean=0,sd=param[18],log=TRUE)+dnorm((y3[3*n+j]-param[15]-param[16]*x3[3*n+j]-param[17]*4),mean=0,sd=param[18],log=TRUE)
      logp4<-dnorm((y4[j]-param[19]-param[20]*x4[j]-param[21]),mean=0,sd=param[22],log=TRUE)+dnorm((y4[n+j]-param[19]-param[20]*x4[n+j]-param[21]*2),mean=0,sd=param[22],log=TRUE)+dnorm((y4[2*n+j]-param[19]-param[20]*x4[2*n+j]-param[21]*3),mean=0,sd=param[22],log=TRUE)+dnorm((y4[3*n+j]-param[19]-param[20]*x4[3*n+j]-param[21]*4),mean=0,sd=param[22],log=TRUE)
      
      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)
      logp1+logp2+logp3+logp4+RVineLogLik(simulate2[j,],RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
    }
    out<-genGradient(x=out1$par,func=optimfun2)
    J=J+out%*%t(out)
    
  }
  J<-J/n
  
  
  #output the results to the record matrix
  cmarginalbeta01[i]<-out1$par[7]
  cmarginalbeta02[i]<-out1$par[11]
  cmarginalbeta03[i]<-out1$par[15]
  cmarginalbeta04[i]<-out1$par[19]
  cmarginalbeta11[i]<-out1$par[8]
  cmarginalbeta12[i]<-out1$par[12]
  cmarginalbeta13[i]<-out1$par[16]
  cmarginalbeta14[i]<-out1$par[20]
  cmarginalbeta21[i]<-out1$par[9]
  cmarginalbeta22[i]<-out1$par[13]
  cmarginalbeta23[i]<-out1$par[17]
  cmarginalbeta24[i]<-out1$par[21]
  csigma1[i]<-out1$par[10]
  csigma2[i]<-out1$par[14]
  csigma3[i]<-out1$par[18]
  csigma4[i]<-out1$par[22]
  cpar12[i]<-out1$par[6]
  cpar13[i]<-out1$par[5]
  cpar14[i]<-out1$par[3]
  cpar23[i]<-out1$par[4]
  cpar24[i]<-out1$par[2]
  cpar34[i]<-out1$par[1]
  mlecomppar[4,3]<-mlecomppar[4,3]+out1$par[6]
  mlecomppar[4,2]<-mlecomppar[4,2]+out1$par[5]
  mlecomppar[4,1]<-mlecomppar[4,1]+out1$par[3]
  mlecomppar[3,2]<-mlecomppar[3,2]+out1$par[4]
  mlecomppar[3,1]<-mlecomppar[3,1]+out1$par[2]
  mlecomppar[2,1]<-mlecomppar[2,1]+out1$par[1]
  
 
  #calculate the Godambe information matrix, which is robustase here
  robustase<-inv(H)%*%J%*%inv(H)
  
  Compstd<-matrix(0,4,4)
  Compstd[4,3]<-(sqrt(diag(robustase))/sqrt(n))[6]
  Compstd[4,2]<-(sqrt(diag(robustase))/sqrt(n))[5]
  Compstd[4,1]<-(sqrt(diag(robustase))/sqrt(n))[3]
  Compstd[3,2]<-(sqrt(diag(robustase))/sqrt(n))[4]
  Compstd[3,1]<-(sqrt(diag(robustase))/sqrt(n))[2]
  Compstd[2,1]<-(sqrt(diag(robustase))/sqrt(n))[1]
  asecompstd<-asecompstd+Compstd
  
  csdmarginalbeta01[i]<-(sqrt(diag(robustase))/sqrt(n))[7]
  csdmarginalbeta02[i]<-(sqrt(diag(robustase))/sqrt(n))[11]
  csdmarginalbeta03[i]<-(sqrt(diag(robustase))/sqrt(n))[15]
  csdmarginalbeta04[i]<-(sqrt(diag(robustase))/sqrt(n))[19]
  csdmarginalbeta11[i]<-(sqrt(diag(robustase))/sqrt(n))[8]
  csdmarginalbeta12[i]<-(sqrt(diag(robustase))/sqrt(n))[12]
  csdmarginalbeta13[i]<-(sqrt(diag(robustase))/sqrt(n))[16]
  csdmarginalbeta14[i]<-(sqrt(diag(robustase))/sqrt(n))[20]
  csdmarginalbeta21[i]<-(sqrt(diag(robustase))/sqrt(n))[9]
  csdmarginalbeta22[i]<-(sqrt(diag(robustase))/sqrt(n))[13]
  csdmarginalbeta23[i]<-(sqrt(diag(robustase))/sqrt(n))[17]
  csdmarginalbeta24[i]<-(sqrt(diag(robustase))/sqrt(n))[21]
  csdsigma1[i]<-(sqrt(diag(robustase))/sqrt(n))[10]
  csdsigma2[i]<-(sqrt(diag(robustase))/sqrt(n))[14]
  csdsigma3[i]<-(sqrt(diag(robustase))/sqrt(n))[18]
  csdsigma4[i]<-(sqrt(diag(robustase))/sqrt(n))[22]
}


#calculate the mean, standard error and corresponding intervals
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<-cpar12+1.96*sd(cpar12)
ci12lower<-cpar12-1.96*sd(cpar12)
count<-0
for(i in 1:repeattime)
  if(ci12upper[i]<2 | ci12lower[i]>2){
    count=count+1
  }
ecpcomp12<-(repeattime-count)/repeattime
ecpcomp12

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

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

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

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

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