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


library(foreach)
library(doParallel)
library(parallel)
library(tictoc)
library(DEoptim)

rm(list=ls(all=TRUE))

k1=3
k2=2

ff=function(n, m){
  a1=18
  a2=21
  b1=0.2
  b2=0.3
  p=6
  u=array(0,c(m,2,n))
  DX=matrix(0,m,n)
  DY=matrix(0,m,n)
  t=seq(0,0.3,0.01)
  
  for(j in 1:n){
    for(i in 1:m){
      library(copula)
      library(VineCopula)
      library(DEoptim)
      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=5
  d2=12
  tX=array(m,n)
  tY=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)
  }
  
 m1=sort(tX)
 m2=sort(tY)
 jst=min(m1[n-k1+1],m2[n-k2+1])
 
L=function(tt) -sum(sapply(1:n,function(j) sum(sapply(1:min(tX[j],tY[j],jst),function(i) log(dCopula(u[i,,j], gumbelCopula(tt))) ))))
  
  
  Lx=function(aa){
    a=aa[1]
    b=aa[2]
    sum(sapply(1:n,function(j) sum(sapply(1:min(tX[j],jst),function(i)  (exp(a*t[i+1])-exp(a*t[i]))*log(b)+DX[i,j]/b-(exp(a*t[i+1])-exp(a*t[i])-1)*log(DX[i,j])+log(gamma(exp(a*t[i+1])-exp(a*t[i]))) ))))
  }
 
  

  Ly=function(cc){
    c=cc[1]
    d=cc[2]
    sum(sapply(1:n,function(j) sum(sapply(1:min(tY[j],jst),function(i)  (exp(c*t[i+1])-exp(c*t[i]))*log(d)+DY[i,j]/d-(exp(c*t[i+1])-exp(c*t[i])-1)*log(DY[i,j])+log(gamma(exp(c*t[i+1])-exp(c*t[i]))) ))))
  }
  
 
  theta=DEoptim(L,lower=4,upper = 8,control=list(itermax=10))$optim$bestmem
  theta1=DEoptim(Lx,lower=c(16,0.025),upper=c(22,0.4),control=list(itermax=10))$optim$bestmem
  theta2=DEoptim(Ly,lower=c(17,0.025),upper=c(22,0.4),control=list(itermax=10))$optim$bestmem
  return(c(theta,theta1,theta2))
}



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=30, m=10---------

A<- matrix(0,1000,10)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(ff(30,10),(ff(30, 10)-c(6,18,0.2,21,0.3))^2)
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))
theta.hat[1]=mean(A[,1])
a1.hat[1]=mean(A[,2])
b1.hat[1]=mean(A[,3])
a2.hat[1]=mean(A[,4])
b2.hat[1]=mean(A[,5])

mse.theta[1]=mean(A[,6])
mse.a1[1]=mean(A[,7])
mse.b1[1]=mean(A[,8])
mse.a2[1]=mean(A[,9])
mse.b2[1]=mean(A[,10])


###----n=30, m=15---------


A<- matrix(0,1000,10)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(ff(30,15),(ff(30,15)-c(6,18,0.2,21,0.3))^2)
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))
theta.hat[2]=mean(A[,1])
a1.hat[2]=mean(A[,2])
b1.hat[2]=mean(A[,3])
a2.hat[2]=mean(A[,4])
b2.hat[2]=mean(A[,5])

mse.theta[2]=mean(A[,6])
mse.a1[2]=mean(A[,7])
mse.b1[2]=mean(A[,8])
mse.a2[2]=mean(A[,9])
mse.b2[2]=mean(A[,10])


###----n=30, m=20---------


A<- matrix(0,1000,10)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(ff(30,20),(ff(30,20)-c(6,18,0.2,21,0.3))^2)
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))
theta.hat[3]=mean(A[,1])
a1.hat[3]=mean(A[,2])
b1.hat[3]=mean(A[,3])
a2.hat[3]=mean(A[,4])
b2.hat[3]=mean(A[,5])

mse.theta[3]=mean(A[,6])
mse.a1[3]=mean(A[,7])
mse.b1[3]=mean(A[,8])
mse.a2[3]=mean(A[,9])
mse.b2[3]=mean(A[,10])

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


A<- matrix(0,1000,10)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(ff(40,10),(ff(40,10)-c(6,18,0.2,21,0.3))^2)
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))
theta.hat[4]=mean(A[,1])
a1.hat[4]=mean(A[,2])
b1.hat[4]=mean(A[,3])
a2.hat[4]=mean(A[,4])
b2.hat[4]=mean(A[,5])

mse.theta[4]=mean(A[,6])
mse.a1[4]=mean(A[,7])
mse.b1[4]=mean(A[,8])
mse.a2[4]=mean(A[,9])
mse.b2[4]=mean(A[,10])



###----n=40, m=15---------


A<- matrix(0,1000,10)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(ff(40,15),(ff(40,15)-c(6,18,0.2,21,0.3))^2)
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))
theta.hat[5]=mean(A[,1])
a1.hat[5]=mean(A[,2])
b1.hat[5]=mean(A[,3])
a2.hat[5]=mean(A[,4])
b2.hat[5]=mean(A[,5])

mse.theta[5]=mean(A[,6])
mse.a1[5]=mean(A[,7])
mse.b1[5]=mean(A[,8])
mse.a2[5]=mean(A[,9])
mse.b2[5]=mean(A[,10])



###----n=40, m=20---------


A<- matrix(0,1000,10)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:1000,.combine='rbind') %dopar% c(ff(40, 20),(ff(40,20)-c(6,18,0.2,21,0.3))^2)
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))
theta.hat[6]=mean(A[,1])
a1.hat[6]=mean(A[,2])
b1.hat[6]=mean(A[,3])
a2.hat[6]=mean(A[,4])
b2.hat[6]=mean(A[,5])

mse.theta[6]=mean(A[,6])
mse.a1[6]=mean(A[,7])
mse.b1[6]=mean(A[,8])
mse.a2[6]=mean(A[,9])
mse.b2[6]=mean(A[,10])
write(rbind(theta.hat,a1.hat,b1.hat,a2.hat,b2.hat),"C:/Users/admin/Desktop/Result  simulation  gumbel new/ mle2-gumbel.txt")
write(rbind(mse.theta,mse.a1,mse.b1,mse.a2,mse.b2),"C:/Users/admin/Desktop/Result  simulation  gumbel new/ mse2-gumbel.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>