/* Multivariate t Copula */
new; cls;
rhos=0.9;
nu=2;              /* degree of freedom */
dim=3;
n=1000;
U=MCt(rhos,nu,dim,n);

library pgraph;
graphset;
pqgwin auto;
_plctrl=-1;
xlabel("u1"); ylabel("u2");
xy(U[.,1],U[.,2]);
xlabel("u2"); ylabel("u3");
xy(U[.,2],U[.,3]);
xlabel("u1"); ylabel("u3");
xy(U[.,1],U[.,3]);



proc MCt(rhos,nu,dim,n);
   local Sig,A,Z,X2,Y,X,U;
   if rows(rhos)==1;
      Sig=rhos*ones(dim,dim); Sig=diagrv(Sig,ones(dim,1));
   else;
      Sig=rhos;
   endif;
   A=chol(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;