/* t Copula */
new; cls;
rho=0.9;
nu=2;              /* degree of freedom */
n=1000;
U=Ct(rho,nu,n);

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



proc Ct(rho,nu,n);
   local Sig,A,dim,Z,X2,Y,X,U;
   Sig=1~rho|
       rho~1;
   A=chol(Sig);
   dim=cols(Sig);
   Z=rndn(n,dim); 
   Z=(Z-meanc(Z)')*inv(chol(vcx(Z))); /* This moment match is optional. You could erase. */
   X2=rndchisq(nu,n);
   Y=Z*A;
   X=Y./(sqrt(X2)/sqrt(nu));
   U=1-cdftc(X,nu);
   retp(U);
endp;

proc rndchisq(m,nn);
   local x;
   x=rndu(nn,1);
   x=2*gammaii(m/2*ones(nn,1),x);
   retp(x);
endp;