install.packages("rhmc")
install.packages("copula")
install.packages("VineCopula")
install.packages("foreach")
install.packages("doParallel")
install.packages("parallel")
install.packages("tictoc")

library(rhmc)
library(foreach)
library(doParallel)
library(parallel)
library(tictoc)
library(copula)
library(VineCopula)


k=4

ff=function(n,m){
a1=6
a2=4
b1=0.8
b2=2
p=2
  library(rhmc)
  u=array(0,c(m,2,n))
  DX=matrix(0,m,n)
  DY=matrix(0,m,n)
  t=seq(0,0.45,0.01)
  
  for(j in 1:n){
    for(i in 1:m){
      library(copula)
      library(VineCopula)
      u[i,,j]=rCopula(1,gumbelCopula(p))
      DX[i,j]=qgamma(u[i,1,j] ,shape=exp(t[i+1]*a1)-exp(t[i]*a1), scale = b1)
      DY[i,j]=qgamma(u[i,2,j] ,shape=exp(t[i+1]*a2)-exp(t[i]*a2), scale = b2)
    }
  }
  
  X=matrix(0,m,n)
  Y=matrix(0,m,n)
  for(j in 1:n){
    for(i in 1:m){
      X[i,j]=sum(DX[1:i,j])
      Y[i,j]=sum(DY[1:i,j])
    }
  }
  
  
  
  d1=2
  d2=2.5
  tX=array(m,n)
  tY=array(m,n)
  tsc=array(m,n)
  for(j in 1:n){
    if(X[m,j]<d1) tX[j]=m else tX[j]=which.max(X[,j]>d1)
    if(Y[m,j]<d2) tY[j]=m else tY[j]=which.max(Y[,j]>d2)
    tsc[j]=min(tX[j],tY[j])
  }
  
  
  S=sort(tsc)
  jsts=S[n-k+1] 
  
  Lx=function(aa){
    a1=aa[1]
    b1=aa[2]
    sum(sapply(1:n,function(j) sum(sapply(1:min(tsc[j],jsts),function(i) 
      (exp(a1*t[i+1])-exp(a1*t[i]))*log(b1)+DX[i,j]/b1-(exp(a1*t[i+1])-exp(a1*t[i])-1)*log(DX[i,j])+log(gamma(exp(a1*t[i+1])-exp(a1*t[i]))) ))))
  }
  
  
  Ly=function(bb){
    a2=bb[1]
    b2=bb[2]
    sum(sapply(1:n,function(j)
      sum(sapply(1:min(tsc[j],jsts),function(i)
        (exp(a2*t[i+1])-exp(a2*t[i]))*log(b2)
        +DY[i,j]/b2-(exp(a2*t[i+1])-exp(a2*t[i])-1)*log(DY[i,j])
        +log(gamma(exp(a2*t[i+1])-exp(a2*t[i]))) ))))
  }
  
  
  
  
  
  init=c(6,0.7)
  hat.Lx=hmc(Lx,init,1000,20,eps=0.001,1)
  a1hat=mean(hat.Lx$chain[1,])
  b1hat=mean(hat.Lx$chain[2,])
  
  init=c(4,1.9)
  hat.Ly=hmc(Ly,init,1000,20,eps=0.001,1)
  a2hat=mean(hat.Ly$chain[1,])
  b2hat=mean(hat.Ly$chain[2,])
  
  
  
  pdx=matrix(0,m,n)
  pdy=matrix(0,m,n)
  
  t=c(seq(0,0.45,0.01))
  for(i in 1:n){
    for(j in 1:m){
      pdx[j,i]=pgamma(DX[j,i],shape=exp(t[i+1]*a1hat)-exp(t[i]*a1hat),scale=b1hat)
      pdy[j,i]=pgamma(DY[j,i],shape=exp(t[i+1]*a2hat)-exp(t[i]*a2hat),scale=b2hat)
    }
  }
  
  
  
  Lhg=function(tt){
    sum(sapply(1:n,function(j) 
      sum(sapply(1:min(tX[j],tY[j],jsts),function(i) 
        ((-log(pdx[i,j]))^tt+(-log(pdy[i,j]))^tt)^(1/tt)
        +log(pdx[i,j])+log(pdy[i,j])
        -(-2+2/tt)*log((-log(pdx[i,j]))^tt+(-log(pdy[i,j]))^tt)
        -(tt-1)*log(log(pdx[i,j])*log(pdy[i,j]))
        -log(1+(tt-1)*((-log(pdx[i,j]))^tt+(-log(pdy[i,j]))^tt)^(-1/tt))
      ))))
  }
  
  init=2
  hat.L=hmc(Lhg,init,1000,20,eps=0.001,1)
  tthat=mean(hat.L$chain[1,])
  return(c(a1hat,b1hat,a2hat,b2hat,tthat))
}




theta.hat=c()
a1.hat=c()
b1.hat=c()
a2.hat=c()
b2.hat=c()
mse.theta=c()
mse.a1=c()
mse.b1=c()
mse.a2=c()
mse.b2=c()

#============================  n=10, m=15  ====================================

A<- matrix(0,1000,5)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(tryCatch(ff(10, 15),error=function(e)return(NA)))
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))

B<-A[!rowSums(!is.finite(A)),]

a1.hat[1]=mean(B[,1])
b1.hat[1]=mean(B[,2])
a2.hat[1]=mean(B[,3])
b2.hat[1]=mean(B[,4])
theta.hat[1]=mean(B[,5])


mse.a1[1]=mean((B[,1]-6)^2)
mse.b1[1]=mean((B[,2]-0.8)^2)
mse.a2[1]=mean((B[,3]-4)^2)
mse.b2[1]=mean((B[,4]-2)^2)
mse.theta[1]=mean((B[,5]-2)^2)

#============================  n=10, m=30  ====================================

A<- matrix(0,1000,5)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(tryCatch(ff(10, 30),error=function(e)return(NA)))
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))

B<-A[!rowSums(!is.finite(A)),]

a1.hat[2]=mean(B[,1])
b1.hat[2]=mean(B[,2])
a2.hat[2]=mean(B[,3])
b2.hat[2]=mean(B[,4])
theta.hat[2]=mean(B[,5])


mse.a1[2]=mean((B[,1]-6)^2)
mse.b1[2]=mean((B[,2]-0.8)^2)
mse.a2[2]=mean((B[,3]-4)^2)
mse.b2[2]=mean((B[,4]-2)^2)
mse.theta[2]=mean((B[,5]-2)^2)





#============================  n=10, m=40  ====================================

A<- matrix(0,1000,5)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(tryCatch(ff(10, 40),error=function(e)return(NA)))
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))

B<-A[!rowSums(!is.finite(A)),]

a1.hat[3]=mean(B[,1])
b1.hat[3]=mean(B[,2])
a2.hat[3]=mean(B[,3])
b2.hat[3]=mean(B[,4])
theta.hat[3]=mean(B[,5])


mse.a1[3]=mean((B[,1]-6)^2)
mse.b1[3]=mean((B[,2]-0.8)^2)
mse.a2[3]=mean((B[,3]-4)^2)
mse.b2[3]=mean((B[,4]-2)^2)
mse.theta[3]=mean((B[,5]-2)^2)


write(rbind(a1.hat,b1.hat,a2.hat,b2.hat,theta.hat),"C:/Users/Dodman/Desktop/New folder/ mle-frank.txt")
write(rbind(theta.hat,a.hat,b.hat,c.hat,d.hat),"C:/Users/Dodman/Desktop/result simulation/ mse-frank.txt")



<pre>An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler-&gt;renderException()
#3 [internal function]: yii\base\ErrorHandler-&gt;handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application-&gt;end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files-&gt;getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController-&gt;actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction-&gt;runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller-&gt;runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module-&gt;runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application-&gt;handleRequest()
#10 /var/www/html/prof-homepages/frontend/web/index.php(18): yii\base\Application-&gt;run()
#11 {main}</pre>