/* Multivariate Farlie-Gumbel-Morgenstern Copula  */
new; cls;
theta=0.9;               /* -1<=theta<=1 */
dim=3;                   /* dim=2,3,4,... */
n=1000;
u=MCfgm(theta,dim,n);

library pgraph;
graphset;
_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 MCfgm(theta,dim,n);
   local u,q,a,b,c,um,up,i,j,k;
   if theta<-1 or theta>1;
      errorlog "ERROR: Parameter theta must be -1<=theta<=1.";
      retp(".");
   endif;
   u=rndu(n,dim);
/* generation of u[.,k] for k=2,3,... */
   k=2;
   do while k<=dim;
      q=u[.,k];
      c=1;
      i=1;
      do while i<=k-2;
         j=i+1;
         do while j<=k-1;
            c=c+theta*(1-2*u[.,i]).*(1-2*u[.,j]);
            j=j+1;
         endo;
         i=i+1;
      endo;
      a=0;
      b=c;
      i=1;
      do while i<=k-1;
         a=a-theta*(1-2*u[.,i]);
         b=b+theta*(1-2*u[.,i]);
         i=i+1;
      endo;
      um=-b./(2*a)-sqrt((b./(2*a))^2+(c./a).*q); up=-b./(2*a)+sqrt((b./(2*a))^2+(c./a).*q);
      u[.,k]=((um.>=0).and (um.<=1)).*um+((up.>=0).and (up.<=1)).*up;
      k=k+1;
   endo;
   retp(u);
endp;