#======================  FGM Copula                ============================
#====================== cost in  finite time case  ============================ 
install.packages("foreach")
install.packages("doParallel")
install.packages("parallel")
install.packages("tictoc")
install.packages("copula")
install.packages("DEoptim")

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



#====== a1:=Alpha_1, a1:=Alpha_2,  b1:=Beta1_1,   b2:=Beta1_2    
#====== S1:=omega_1, S2:=omega_2, p:=lambda   ,   t0:=t
#====== m1:=M_1    , m2:=M_2  
#====== x1 and x2 are the initial states of components in each element. 


CCost=function(T){
  
library(copula)
x1=0
x2=0
t0=4
p=2

a1=2
b1=1
a2=3
b2=1
S1=1.2
S2=1.2
aaa=0.2
m1=aaa*S1
m2=0.8
c1=5
c2=5
cr=0.05
cp=0.01
cu=10


Fy1=function(y1,t){
  pgamma(y1,(a1*t),1/b1)
}

fy1=function(y1,t){
  dgamma(y1,(a1*t),1/b1)
}


Fy2=function(y2,t){
  pgamma(y2,(a2*t),1/b2)
}

fy2=function(y2,t){
  dgamma(y2,(a2*t),1/b2)
}


Fc=function(y1,y2,t){
  pCopula(c(Fy1(y1,t),Fy2(y2,t)),claytonCopula(p))
}



fc=function(y1,y2,t){
  if (Fy1(y1,t)<1*10^(-4) | Fy2(y2,t)<1*10^(-4)) 
    0  else 
  dCopula(c(Fy1(y1,t),Fy2(y2,t)),claytonCopula(p))*fy1(y1,t)*fy2(y2,t) }


Fbc=function(y1,y2,t){
  Fc(y1,y2,t)-Fy2(y2,t)+1-Fy1(y1,t)
}

#========================  Cost function finite time   =========================

I0=function(y1,t){
  if (y1==0) 1
  else 0
}

Iy1=function(y1,t){
  if (y1<=m1 & x1<=y1) 1
  else 0
}

Q11=function(y1,t){
  (1-Fy1(m1-x1,t))*I0(y1,t)+fy1(y1-x1,t)*Iy1(y1,t) 
}





IIy2=function(y2,t){
  if (y2<=m2 & x2<=y2) 1
  else 0
}


Q22=function(y2,t){
  (1-Fy2(m2-x2,t))*I0(y2,t)+fy2(y2-x2,t)*IIy2(y2,t) 
}



#===   Density function  ===========================


J1=function(y1,y2,t){
  if (y1<=m1 & x1<=y1 & y2<=m2 & x2<=y2) 1
  else 0
}


Q1=function(y1,y2,t){
  J1(y1,y2,t)*fc(y1-x1,y2-x2,t)
}


J2=function(y1,y2,t){
  if (y2<=m2 & x2<=y2 & y1==0) 1
  else 0
}


D=function(y2,t){
  dd=function(u) {fc(u-x1,y2-x2,t) }
  dd=Vectorize(dd)
  integrate(dd,m1,30)$value
}


Q2=function(y1,y2,t){
  J2(y1,y2,t)*D(y2,t)
}

J3=function(y1,y2,t){
  if (y1<=m1 & x1<=y1 & y2==0) 1
  else 0
}


Dq=function(y1,t){
  dd=function(u) {fc(y1-x1,u-x2,t) }
  dd=Vectorize(dd)
  integrate(dd,m2,30)$value
}


Q3=function(y1,y2,t){
  J3(y1,y2,t)*Dq(y1,t)
}



J4=function(y1,y2,t){
  if (y1==0 & y2==0) 1
  else 0
}

Q4=function(y1,y2,t){
  Fbc(m1-x1,m2-x2,t)*J4(y1,y2,t)  
}

Q=function(y1,y2,t){
  Q1(y1,y2,t)+Q2(y1,y2,t)+Q3(y1,y2,t)+Q4(y1,y2,t)
}


#===============================   n>=2   ======================================


f1=function(T) {
  a=c()
  a[1]=1-Fy1(m1-x1,T) 
  if (floor(t0/T)>=2){
  for(i in 2:floor(t0/T))
    a[i]=integrate(Vectorize(function(y1) a[i-1]*Q11(y1,T)),0.0001,m1)$v}
  return(a)
}


f2=function(T) {
  b=c()
  b[1]=1-Fy2(m2-x2,T) 
  if (floor(t0/T)>=2){
  for(i in 2:floor(t0/T))
    b[i]=integrate(Vectorize(function(y2) b[i-1]*Q22(y2,T)),0.0001,m2)$v}
  return(b)
}



f12=function(T) {
  a=c()
  a[1]=Fbc(m1-x1,m2-x2,T)
  if (floor(t0/T)>=2){
  for(i in 2:floor(t0/T))
    a[i]=integrate(Vectorize(function(y2) 
      integrate(Vectorize(function(y1) a[i-1]*Q(y1,y2,T)),0.0001,m1)$v),0.0001,m2)$v}
  return(a)
}




ff1=function(T){
  n=floor(t0/T)
  return(sum(f1(T)[1:n]))
}


ff2=function(T){
  n=floor(t0/T)
  return(sum(f2(T)[1:n]))
}




ff12=function(T){
  n=floor(t0/T)
  return(sum(f12(T)[1:n]))
}


#=======================     Cost function     ==============================


UT1=function(T,x11,x12){
  LL=function(u) {Fbc(S1-x11,S2-x12,u)}
  LL=Vectorize(LL)
  integrate(LL,0.0001,T)$value
}

UT2=function(T,x21,x22){
  LL=function(u) {Fbc(S1-x21,S2-x22,u)}
  LL=Vectorize(LL)
  integrate(LL,0.0001,T)$value
}





UU0=function(y11,T){
  LL=function(y12) {UT1(t0-T,y11,y12)*Q(y11,y12,T)}
  LL=Vectorize(LL)
  (integrate(LL,0.0001,m2)$value)
}


UU1=function(T){
  LL=function(y11) {UU0(y11,T)}
  LL=Vectorize(LL)
  integrate(LL,0.0001,m1)$value
}



UU00=function(y21,T){
  LL=function(y22) {UT2(t0-T,y21,y22)*Q(y21,y22,T)}
  LL=Vectorize(LL)
  (integrate(LL,0.0001,m2)$value)
}


UU2=function(T){
  LL=function(y21) {UU00(y21,T)}
  LL=Vectorize(LL)
  integrate(LL,0.0001,m1)$value
}


U=function(T){
  UT1(T,0,0)*UT2(T,0,0)+(UU1(T)*UU2(T))
}


Cost=function(T){
  (cu*U(T)+2*((c1+cr)*ff1(T)+(c2+cr)*ff2(T)-cr*ff12(T))+cp*floor(t0/T))/t0
}

Cost=Vectorize(Cost)
return(Cost(T))
}


#============================================================================

TT=seq(0.01,4,0.02)


CCC=c()
registerDoParallel(makeCluster(detectCores()-1))
tic()
CCC=foreach(i=1:length(TT),.combine='rbind') %dopar% CCost(TT[i])
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()-1))



plot(TT,CCC,xlab="T",ylab="Cost function for Clayton copula",
     type="l",lwd=1,lty=1,col=1)
grid( ny = 6, col = "lightgray", lty = 3,
      lwd = 0.0005, equilogs = TRUE)


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