install.packages("copula")
install.packages("VineCopula")
install.packages("DEoptim")


library(copula)
library(VineCopula)
library(DEoptim)
#=============================== Data ==========================================

M=c(0.9,0.95,1,1.05,1.12,1.19,1.27,1.35,1.48,1.64,0,0,
0.9,0.94,0.98,1.03,1.08,1.14,1.21,1.28,1.37,1.47,1.6,0,
0.90,0.94,0.98,1.03,1.08,1.13,1.19,1.26,1.35,1.46,1.58,1.77,
0.90,0.94,0.98,1.03,1.07,1.12,1.19,1.25,1.34,1.43,1.55,1.73,
0.90,0.94,0.98,1.03,1.07,1.12,1.19,1.24,1.34,1.43,1.55,1.71,
0.90,0.94,0.98,1.03,1.07,1.12,1.18,1.23,1.33,1.41,1.51,1.68,
0.90,0.94,0.98,1.02,1.07,1.11,1.17,1.23,1.32,1.41,1.52,1.66,
0.90,0.93,0.97,1.00,1.06,1.11,1.17,1.23,1.30,1.39,1.49,1.62,
0.90,0.92,0.97,1.01,1.05,1.09,1.15,1.21,1.28,1.36,1.44,1.55,
0.90,0.92,0.96,1.00,1.04,1.08,1.13,1.19,1.26,1.34,1.42,1.52,
0.90,0.93,0.96,1.00,1.04,1.08,1.13,1.18,1.24,1.31,1.39,1.49,
0.90,0.93,0.97,1.00,1.03,1.07,1.10,1.16,1.22,1.29,1.37,1.48,
0.90,0.92,0.97,0.99,1.03,1.06,1.10,1.14,1.20,1.26,1.31,1.4,
0.90,0.93,0.96,1.00,1.03,1.07,1.12,1.16,1.20,1.26,1.30,1.37,
0.90,0.92,0.96,0.99,1.03,1.06,1.10,1.16,1.21,1.27,1.33,1.4,
0.90,0.92,0.95,0.97,1.00,1.03,1.07,1.11,1.16,1.22,1.26,1.33,
0.90,0.93,0.96,0.97,1.00,1.05,1.08,1.11,1.16,1.20,1.24,1.32,
0.90,0.92,0.94,0.97,1.01,1.04,1.07,1.09,1.14,1.19,1.23,1.28,
0.90,0.92,0.94,0.97,0.99,1.02,1.05,1.08,1.12,1.16,1.20,1.25,
0.90,0.92,0.94,0.97,0.99,1.02,1.05,1.08,1.12,1.16,1.19,1.24)
N=matrix(M,20,12,byrow=T)
dx1=N[1,c(-1,-11,-12)]-N[1,c(-10,-11,-12)]
dx2=N[2,c(-1,-12)]-N[2,c(-11,-12)]
dx3=N[3,-1]-N[3,-12]
dx4=N[4,-1]-N[4,-12]
dx5=N[5,-1]-N[5,-12]
dx6=N[6,-1]-N[6,-12]
dx7=N[7,-1]-N[7,-12]
dx8=N[8,-1]-N[8,-12]
dx9=N[9,-1]-N[9,-12]
dx10=N[10,-1]-N[10,-12]

dy1=N[11,-1]-N[11,-12]
dy2=N[12,-1]-N[12,-12]
dy3=N[13,-1]-N[13,-12]
dy4=N[14,-1]-N[14,-12]
dy5=N[15,-1]-N[15,-12]
dy6=N[16,-1]-N[16,-12]
dy7=N[17,-1]-N[17,-12]
dy8=N[18,-1]-N[18,-12]
dy9=N[19,-1]-N[19,-12]
dy10=N[20,-1]-N[20,-12]


x=matrix(0,10,11)
x[1,]=c(dx1,0,0);x[2,]=c(dx2,0);x[3,]=dx3;x[4,]=dx4;x[5,]=dx5;x[6,]=dx6;x[7,]=dx7;x[8,]=dx8;x[9,]=dx9;x[10,]=dx10
y=matrix(0,10,11)
y[1,]=dy1;y[2,]=dy2;y[3,]=dy3;y[4,]=dy4;y[5,]=dy5;y[6,]=dy6;y[7,]=dy7;y[8,]=dy8;y[9,]=dy9;y[10,]=dy10

#=================== Likelihood function =======================================

Lx=function(aa){
  a=aa[1]
  b=aa[2]
  t=c(seq(0.00,0.11,0.01))
  sum(sapply(1:9,function(i) (exp(a*t[i+1])-exp(a*t[i]))*log(b)+x[1,i]/b-(exp(a*t[i+1])-exp(a*t[i])-1)*log(x[1,i])+log(gamma(exp(a*t[i+1])-exp(a*t[i])))) )
 +sum(sapply(1:10,function(i) (exp(a*t[i+1])-exp(a*t[i]))*log(b)+x[2,i]/b-(exp(a*t[i+1])-exp(a*t[i])-1)*log(x[2,i])+log(gamma(exp(a*t[i+1])-exp(a*t[i])))) )
 +sum(sapply(3:10, function(h) sum(sapply(1:11,function(i) (exp(a*t[i+1])-exp(a*t[i]))*log(b)+x[h,i]/b-(exp(a*t[i+1])-exp(a*t[i])-1)*log(x[h,i])+log(gamma(exp(a*t[i+1])-exp(a*t[i])))))) )
}


Ly=function(cc){
  c=cc[1]
  d=cc[2]
  t=c(seq(0.00,0.11,0.01))
  sum(sapply(1:10, function(h) sum(sapply(1:11,function(i) (exp(c*t[i+1])-exp(c*t[i]))*log(d)+y[h,i]/d-(exp(c*t[i+1])-exp(c*t[i])-1)*log(y[h,i])+log(gamma(exp(c*t[i+1])-exp(c*t[i])))))) )
}


DEoptim(Lx,lower=c(25,0.025),upper=c(37,0.03),control=list(itermax=100))$optim$bestmem
DEoptim(Ly,lower=c(25,0.025),upper=c(37,0.03),control=list(itermax=100))$optim$bestmem


#======================   estimation theta =====================================


n=10
pdx=matrix(0,n,11)
pdy=matrix(0,n,11)

t=c(seq(0,0.11,0.01))
for(h in 1:n){
  for(i in 1:11){
    pdx[h,i]=pgamma(x[h,i],shape=exp(t[i+1]*30.1)-exp(t[i]*30.1),scale=0.028)
    pdy[h,i]=pgamma(y[h,i],shape=exp(t[i+1]*27.06)-exp(t[i]*27.06),scale=0.025)
  }
}

LFg=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), fgmCopula(p)))
  )    ) )
+sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), fgmCopula(p)))
  )    ) )
 + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), fgmCopula(p))))     ) )    ))
}


Lg=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), gumbelCopula(p)))
  )    ) )
+sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), gumbelCopula(p)))
  )    ) )
 + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), gumbelCopula(p))))     ) )    ))
}


Lf=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), frankCopula(p)))
  )    ) )
+sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), frankCopula(p)))
  )    ) )
 + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), frankCopula(p))))     ) )    ))
}



Lj=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), joeCopula(p)))
  )    ) )
+sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), joeCopula(p)))
  )    ) )
 + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i) 
   (  -log (dCopula(c(pdx[h,i],pdy[h,i]), joeCopula(p))))     ) )    ))
}




Lnorm=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), normalCopula(p)))
  )    ) )
  +sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), normalCopula(p)))
  )    ) )
  + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i)  
    (  -log (dCopula(c(pdx[h,i],pdy[h,i]), normalCopula(p))))     ) )    ))
}



Lt=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), tCopula(p)))
  )    ) )
  +sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), tCopula(p)))
  )    ) )
  + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i)  
    (  -log (dCopula(c(pdx[h,i],pdy[h,i]), tCopula(p))))     ) )    ))
}


Lc=function(p){
  sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), claytonCopula(p)))
  )    ) )
+sum(sapply(1:10,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), claytonCopula(p)))
  )    ) )
 + sum(sapply(3:10, function(h)      sum(sapply(1:11,function(i)  
   (  -log (dCopula(c(pdx[h,i],pdy[h,i]), claytonCopula(p))))     ) )    ))
}




DEoptim(LFg,lower=0.8,upper = 1, control=list(itermax=10))$optim$bestmem
DEoptim(Lg,lower=1,upper = 9, control=list(itermax=10))$optim$bestmem
DEoptim(Lf,lower=5,upper = 20,control=list(itermax=20))$optim$bestmem
DEoptim(Lc,lower=1.5,upper = 10,control=list(itermax=20))$optim$bestmem
DEoptim(Lj,lower=3,upper = 15,control=list(itermax=20))$optim$bestmem
DEoptim(Lnorm,lower=0.002,upper = 1,control=list(itermax=20))$optim$bestmem
DEoptim(Lt,lower=0.002,upper = 1,control=list(itermax=20))$optim$bestmem

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