/* Archimedean Copula ( phi(u)=(u^(-1)-1)^theta ) */
new; cls;
theta=1.5;
n=1000;
U=Carch(theta,n);

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




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

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

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