/* Multivariate Gaussian Copula */
new; cls;
rhos=0.9;
dim=3;               
n=1000;
U=MCgaus(rhos,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 MCgaus(rhos,dim,n);
   local Sig,A,Z,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. */
   X=Z*A;
   U=cdfn(X);
   retp(U);
endp;