/* Joe's BB6 Copula in Archimedean Family ( phi(u)=(-ln(1-(1-u)^k))^theta ) */
new; cls;
theta=1.5;
k=2;
n=1000;
U=Cbb6(theta,k,n);

library pgraph;
graphset;
_plctrl=-1;
xy(U[.,1],U[.,2]);



proc Cbb6(theta,k,n);
   local s,q,t,u,v;
   if theta<1;
      errorlog "ERROR: Parameter theta must be greater than or equal to 1.";
      retp(".");
   endif;
   if k<1;
      errorlog "ERROR: Parameter k must be greater than or equal to 1.";
      retp(".");
   endif;
   s=rndu(n,1);
   q=rndu(n,1);
   t=calct2(q,theta,k);
   u=phi_1(s.*phi(t,theta,k),theta,k);
   v=phi_1((1-s).*phi(t,theta,k),theta,k);
   retp(u~v);
endp;

proc calct2(q,theta,k);
   local tstar,tstarstar,i,t,x;
   tstar=zeros(rows(q),1);
   i=1;
   do while i<=rows(q);
      t=seqa(0.0001,0.0001,10000);
      x=abs(q[i]-(t-phi(t,theta,k)./phi1(t,theta,k)));
      tstar[i]=t[minindc(x)];
      t=seqa(tstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-(t-phi(t,theta,k)./phi1(t,theta,k)));
      tstarstar=t[minindc(x)];
      if tstarstar==tstar[i]-0.00005;
         t=seqa(0.00000001,0.00000001,10000);
         x=abs(q[i]-(t-phi(t,theta,k)./phi1(t,theta,k)));
         tstar[i]=t[minindc(x)];
      else;
         tstar[i]=tstarstar;
      endif;
      if tstarstar>1;
         tstar[i]=1;
      endif;
      i=i+1;
   endo;
   retp(tstar);
endp;


fn phi(u,theta,k)=(-ln(1-(1-u)^k))^theta;

fn phi1(u,theta,k)=theta*(-ln(1-(1-u)^k))^(theta-1).*((-1)/(1-(1-u)^k)).*((-k)*(1-u)^(k-1))*(-1);

fn phi_1(u,theta,k)=1-(1-exp(-u^(1/theta)))^(1/k);