#in this part, only the true vine family and corresponding parameters are changed, other parts are identical to the full likelihood simulatneous. 

library(VineCopula)
library(Matrix)
library(matrixcalc)
library(matlib)
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<-750
k<-4
l<-4

repeattime<-500

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

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


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)]
  }
  
  
  #this is the likelihood function, where the copula family matrix is Rvinef, the misspecified matrix
  optimfun<-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<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(k in 1:n){
      logp1<-logp1+dnorm((y1[k]-param[25]-param[26]*x1[k]-param[27]),mean=0,sd=param[28],log=TRUE)+dnorm((y1[n+k]-param[25]-param[26]*x1[n+k]-param[27]*2),mean=0,sd=param[28],log=TRUE)+dnorm((y1[2*n+k]-param[25]-param[26]*x1[2*n+k]-param[27]*3),mean=0,sd=param[28],log=TRUE)+dnorm((y1[3*n+k]-param[25]-param[26]*x1[3*n+k]-param[27]*4),mean=0,sd=param[28],log=TRUE)
      logp2<-logp2+dnorm((y2[k]-param[29]-param[30]*x2[k]-param[31]),mean=0,sd=param[32],log=TRUE)+dnorm((y2[n+k]-param[29]-param[30]*x2[n+k]-param[31]*2),mean=0,sd=param[32],log=TRUE)+dnorm((y2[2*n+k]-param[29]-param[30]*x2[2*n+k]-param[31]*3),mean=0,sd=param[32],log=TRUE)+dnorm((y2[3*n+k]-param[29]-param[30]*x2[3*n+k]-param[31]*4),mean=0,sd=param[32],log=TRUE)
      logp3<-logp3+dnorm((y3[k]-param[33]-param[34]*x3[k]-param[35]),mean=0,sd=param[36],log=TRUE)+dnorm((y3[n+k]-param[33]-param[34]*x3[n+k]-param[35]*2),mean=0,sd=param[36],log=TRUE)+dnorm((y3[2*n+k]-param[33]-param[34]*x3[2*n+k]-param[35]*3),mean=0,sd=param[36],log=TRUE)+dnorm((y3[3*n+k]-param[33]-param[34]*x3[3*n+k]-param[35]*4),mean=0,sd=param[36],log=TRUE)
      logp4<-logp4+dnorm((y4[k]-param[37]-param[38]*x4[k]-param[39]),mean=0,sd=param[40],log=TRUE)+dnorm((y4[n+k]-param[37]-param[38]*x4[n+k]-param[39]*2),mean=0,sd=param[40],log=TRUE)+dnorm((y4[2*n+k]-param[37]-param[38]*x4[2*n+k]-param[39]*3),mean=0,sd=param[40],log=TRUE)+dnorm((y4[3*n+k]-param[37]-param[38]*x4[3*n+k]-param[39]*4),mean=0,sd=param[40],log=TRUE)
    }
    
    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)
    -(logp1+logp2+logp3+logp4+RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik)
  }
  
  
  
  #start parameters
  startpar<-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,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),lower=lo,upper=up)
  
  
  
  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)
    
    logp1<-0
    logp2<-0
    logp3<-0
    logp4<-0
    for(k in 1:n){
      logp1<-logp1+dnorm((y1[k]-param[25]-param[26]*x1[k]-param[27]),mean=0,sd=param[28],log=TRUE)+dnorm((y1[n+k]-param[25]-param[26]*x1[n+k]-param[27]*2),mean=0,sd=param[28],log=TRUE)+dnorm((y1[2*n+k]-param[25]-param[26]*x1[2*n+k]-param[27]*3),mean=0,sd=param[28],log=TRUE)+dnorm((y1[3*n+k]-param[25]-param[26]*x1[3*n+k]-param[27]*4),mean=0,sd=param[28],log=TRUE)
      logp2<-logp2+dnorm((y2[k]-param[29]-param[30]*x2[k]-param[31]),mean=0,sd=param[32],log=TRUE)+dnorm((y2[n+k]-param[29]-param[30]*x2[n+k]-param[31]*2),mean=0,sd=param[32],log=TRUE)+dnorm((y2[2*n+k]-param[29]-param[30]*x2[2*n+k]-param[31]*3),mean=0,sd=param[32],log=TRUE)+dnorm((y2[3*n+k]-param[29]-param[30]*x2[3*n+k]-param[31]*4),mean=0,sd=param[32],log=TRUE)
      logp3<-logp3+dnorm((y3[k]-param[33]-param[34]*x3[k]-param[35]),mean=0,sd=param[36],log=TRUE)+dnorm((y3[n+k]-param[33]-param[34]*x3[n+k]-param[35]*2),mean=0,sd=param[36],log=TRUE)+dnorm((y3[2*n+k]-param[33]-param[34]*x3[2*n+k]-param[35]*3),mean=0,sd=param[36],log=TRUE)+dnorm((y3[3*n+k]-param[33]-param[34]*x3[3*n+k]-param[35]*4),mean=0,sd=param[36],log=TRUE)
      logp4<-logp4+dnorm((y4[k]-param[37]-param[38]*x4[k]-param[39]),mean=0,sd=param[40],log=TRUE)+dnorm((y4[n+k]-param[37]-param[38]*x4[n+k]-param[39]*2),mean=0,sd=param[40],log=TRUE)+dnorm((y4[2*n+k]-param[37]-param[38]*x4[2*n+k]-param[39]*3),mean=0,sd=param[40],log=TRUE)+dnorm((y4[3*n+k]-param[37]-param[38]*x4[3*n+k]-param[39]*4),mean=0,sd=param[40],log=TRUE)
    }
    
    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)
    logp1+logp2+logp3+logp4+RVineLogLik(simulate2,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
  }
  out<-fdHess(out1$par,optimfun1)
  H<-out$Hessian/n
  
  if(any(is.nan(sqrt(diag(inv(-H))))==TRUE)){
    next
  }else{
    marginalbeta01[i]<-out1$par[25]
    marginalbeta02[i]<-out1$par[29]
    marginalbeta03[i]<-out1$par[33]
    marginalbeta04[i]<-out1$par[37]
    marginalbeta11[i]<-out1$par[26]
    marginalbeta12[i]<-out1$par[30]
    marginalbeta13[i]<-out1$par[34]
    marginalbeta14[i]<-out1$par[38]
    marginalbeta21[i]<-out1$par[27]
    marginalbeta22[i]<-out1$par[31]
    marginalbeta23[i]<-out1$par[35]
    marginalbeta24[i]<-out1$par[39]
    sigma1[i]<-out1$par[28]
    sigma2[i]<-out1$par[32]
    sigma3[i]<-out1$par[36]
    sigma4[i]<-out1$par[40]
    par12[i]<-out1$par[10]
    par13[i]<-out1$par[6]
    par14[i]<-out1$par[3]
    par23[i]<-out1$par[5]
    par24[i]<-out1$par[2]
    par34[i]<-out1$par[1]
    mlefullpar[4,3]<-mlefullpar[4,3]+out1$par[10]
    mlefullpar[4,2]<-mlefullpar[4,2]+out1$par[6]
    mlefullpar[4,1]<-mlefullpar[4,1]+out1$par[3]
    mlefullpar[3,2]<-mlefullpar[3,2]+out1$par[5]
    mlefullpar[3,1]<-mlefullpar[3,1]+out1$par[2]
    mlefullpar[2,1]<-mlefullpar[2,1]+out1$par[1]
    
    
    
    Fullstd<-matrix(0,4,4)
    Fullstd[4,3]<-(sqrt(diag(inv(-H)))/sqrt(n))[10]
    Fullstd[4,2]<-(sqrt(diag(inv(-H)))/sqrt(n))[6]
    Fullstd[4,1]<-(sqrt(diag(inv(-H)))/sqrt(n))[3]
    Fullstd[3,2]<-(sqrt(diag(inv(-H)))/sqrt(n))[5]
    Fullstd[3,1]<-(sqrt(diag(inv(-H)))/sqrt(n))[2]
    Fullstd[2,1]<-(sqrt(diag(inv(-H)))/sqrt(n))[1]
    sdmarginalbeta01[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[25]
    sdmarginalbeta02[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[29]
    sdmarginalbeta03[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[33]
    sdmarginalbeta04[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[37]
    sdmarginalbeta11[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[26]
    sdmarginalbeta12[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[30]
    sdmarginalbeta13[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[34]
    sdmarginalbeta14[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[38]
    sdmarginalbeta21[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[27]
    sdmarginalbeta22[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[31]
    sdmarginalbeta23[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[35]
    sdmarginalbeta24[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[39]
    sdsigma1[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[28]
    sdsigma2[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[32]
    sdsigma3[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[36]
    sdsigma4[i]<-(sqrt(diag(inv(-H)))/sqrt(n))[40]
    asefullstd<-asefullstd+Fullstd
  }
  
  
}

mean(marginalbeta01)
mean(marginalbeta02)
mean(marginalbeta03)
mean(marginalbeta04)
mean(marginalbeta11)
mean(marginalbeta12)
mean(marginalbeta13)
mean(marginalbeta14)
mean(marginalbeta21)
mean(marginalbeta22)
mean(marginalbeta23)
mean(marginalbeta24)
mean(sigma1)
mean(sigma2)
mean(sigma3)
mean(sigma4)
mean(sdmarginalbeta01)
mean(sdmarginalbeta02)
mean(sdmarginalbeta03)
mean(sdmarginalbeta04)
mean(sdmarginalbeta11)
mean(sdmarginalbeta12)
mean(sdmarginalbeta13)
mean(sdmarginalbeta14)
mean(sdmarginalbeta21)
mean(sdmarginalbeta22)
mean(sdmarginalbeta23)
mean(sdmarginalbeta24)
mean(sdsigma1)
mean(sdsigma2)
mean(sdsigma3)
mean(sdsigma4)
mean(par12)
mean(par13)
mean(par14)
mean(par23)
mean(par24)
mean(par34)
asefullstd/repeattime

sd(marginalbeta01)
sd(marginalbeta02)
sd(marginalbeta03)
sd(marginalbeta04)
sd(marginalbeta11)
sd(marginalbeta12)
sd(marginalbeta13)
sd(marginalbeta14)
sd(marginalbeta21)
sd(marginalbeta22)
sd(marginalbeta23)
sd(marginalbeta24)
sd(sigma1)
sd(sigma2)
sd(sigma3)
sd(sigma4)

sd(par12)
sd(par13)
sd(par14)
sd(par23)
sd(par24)
sd(par34)


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

cisigma4up<-sigma4+1.96*sd(sigma4)
cisigma4lo<-sigma4-1.96*sd(sigma4)
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
