/* Joe's BB5 Copula                                                                                      */
/* C(u,v)=exp(-((-ln(u))^k+(-ln(v))^k-((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta))^(1/k))       */
/* Solve  q=exp(-((-ln(u))^k+(-ln(v))^k-((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta))^(1/k))     */
/*        .*(-1/k).*((-ln(u))^k+(-ln(v))^k-((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta))^(1/k-1) */
/*        .*( k*(-ln(u))^(k-1).*(-1/u)+1/theta*((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta-1)    */
/*        .*(-k*theta).*(-ln(u))^(-k*theta-1).*(-1/u) )  (=Cu) in terms of v numerically.  Very Slow.    */
new; cls;
theta=1.5;
k=2;
n=1000;
U=Cbb5(theta,k,n);

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



proc Cbb5(theta,k,n);
   local u,q,v;
   if theta<=0;
      errorlog "ERROR: Parameter theta must be positive.";
      retp(".");
   endif;
   if k<1;
      errorlog "ERROR: Parameter k must be k>=1.";
      retp(".");
   endif;
   u=rndu(n,1);
   q=rndu(n,1);
   v=calcv(u,q,theta,k);
   retp(u~v);
endp;

proc calcv(u,q,theta,k);
   local vstar,vstarstar,i,v,x;
   vstar=zeros(rows(q),1);
   i=1;
   do while i<=rows(q);
      v=seqa(0.0001,0.0001,10000);
      x=abs(q[i]-Cu(u[i],v,theta,k));
      vstar[i]=v[minindc(x)];
      v=seqa(vstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-Cu(u[i],v,theta,k));
      vstarstar=v[minindc(x)];
      if vstarstar==vstar[i]-0.00005;
         v=seqa(0.00000001,0.00000001,10000);
         x=abs(q[i]-Cu(u[i],v,theta,k));
         vstar[i]=v[minindc(x)];
      else;
         vstar[i]=vstarstar;
      endif;
      if vstarstar>1;
         vstar[i]=1;
      endif;
      i=i+1;
   endo;
   retp(vstar);
endp;

fn Cu(u,v,theta,k)=exp(-((-ln(u))^k+(-ln(v))^k-((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta))^(1/k))
                   .*(-1/k).*((-ln(u))^k+(-ln(v))^k-((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta))^(1/k-1)
                   .*( k*(-ln(u))^(k-1).*(-1/u)+1/theta*((-ln(u))^(-k*theta)+(-ln(v))^(-k*theta))^(-1/theta-1)
                       .*(-k*theta).*(-ln(u))^(-k*theta-1).*(-1/u) );