/* Galambos Copula (by acceptance-rejection method)                        */
/* cdf: C(u,v)=u.*v.*exp(((-ln(u))^(-theta)+(-ln(v))^(-theta))^(-1/theta)) */
/* Notice: This simulation is rough. For large theta, very slow.           */
new; cls;
theta=1.5;
n=1000;
X=Cgalambos(theta,n);

library pgraph;
graphset;
pqgwin auto;
u=seqa(0.01,0.01,99); v=seqa(0.01,0.01,99);
title("target density");
xlabel("u");
ylabel("v");
contour(u',v,pdf(u,v',theta));
surface(u',v,pdf(u,v',theta));

title("Galambos Copula");
_plctrl=-1;
xy(X[.,1],X[.,2]);



proc Cgalambos(theta,n);
   local u,v,c,h,i,us,vs,x,z;
   if theta<0;
      errorlog "ERROR: Parameter theta must be nonnegative.";
      retp(".");
   endif;
   u=seqa(0.01,0.01,100);
   v=seqa(0.01,0.01,100);
   c=pdf(u',v,theta);
   h=maxc(maxc(c));
   u=zeros(n,1); v=zeros(n,1);
   i=1;
   do while i<=n;
   label:
      us=rndu(1,1);
      vs=rndu(1,1);
      x=rndu(1,1);
      z=pdf(us,vs,theta);
      if x<=(z/h);
         u[i]=us;
         v[i]=vs;
      else;
         goto label;
      endif;
      i=i+1;
   endo;
   retp(u~v);
endp;


proc pdf(u,v,theta);
  local p;
  p=1-((-ln(u))^(-theta)+(-ln(v))^(-theta))^(-1/theta-1).*((-ln(u))^(-theta-1)+(-ln(v))^(-theta-1)) 
     +((-ln(u))^(-theta)+(-ln(v))^(-theta))^(-1/theta-2).*((-ln(u)).*(-ln(v)))^(-1-theta) 
     .*(1+theta+((-ln(u))^(-theta)+(-ln(v))^(-theta))^(-1/theta));
  retp( u.*v.*exp(((-ln(u))^(-theta)+(-ln(v))^(-theta))^(-1/theta))./(u.*v).*p );
endp;