/* Gumbel-Barnett Copula (by acceptance-rejection method)           */
/* cdf: C(u,v)=u.*v.*exp(-theta*(-ln(u)).*(-ln(v)))                 */
/* Notice: This simulation is still rough.(Latin Hypercube version) */
new; cls;
theta=1;
nn=14;                  /* max nn=14 in Light version of GAUSS */
X=Cgbh(theta,nn);

library pgraph;
graphset;
pqgwin auto;
u=seqa(0.01,0.01,99); v=seqa(0.01,0.01,99);
title("target density");
xlabel("u");
ylabel("v");
contour(u',v,pdf(u,v',theta));
surface(u',v,pdf(u,v',theta));

title("Gumbel-Barnett Copula");
_plctrl=-1;
xy(X[.,1],X[.,2]);



proc Cgbh(theta,nn);
   local u,v,c,h,i,us,vs,x,z,m;
   if theta<0 or theta>1;
      errorlog "ERROR: Parameter theta must be 0<=theta<=1.";
      retp(".");
   endif;
   u=seqa(0.01,0.01,100);
   v=seqa(0.01,0.01,100);
   c=pdf(u',v,theta);
   h=maxc(maxc(c));
   m=LHss(nn,3);
   u=miss(0,0); v=miss(0,0);
   i=1;
   do while i<=nn^3;
      us=m[i,1];
      vs=m[i,2];
      x=m[i,3];
      z=pdf(us,vs,theta);
      if x<=(z/h);
         u=u|us;
         v=v|vs;
      endif;
      i=i+1;
   endo;
   u=u[2:rows(u)];
   v=v[2:rows(v)];
   print/lz "# of pair=" rows(u);
   retp(u~v);
endp;

proc pdf(u,v,theta);
   retp( (1-theta-theta*(ln(u)+ln(v))+theta^2*ln(u).*ln(v)).*exp(-theta*ln(u).*ln(v)) );
endp;


/*
**  LHss.txt - Latin Hypercube U(0,1).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Latin Hypercube numbers in a very easy way. new algorithm.
**
**  Format:     x=LHss(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     x        matrix ,   (nn^dim) x (dim) of resulting U(0,1) numbers
**
**  Notice:     It takes lots of memory to run. Light version of GAUSS will not work in most cases.
*/
proc LHss(nn,dim);
   local x,b,n,z,y,indexvec,d,U,i,index;
/* convert x(base 10) to y(base b) */
   x=seqa(1,1,nn^dim-1); 
   b=nn;
   n=maxc(log(x)/log(b))+1;
   z=reshape(b,n,rows(x));
   y=rev(recserrc(x,z))';
   indexvec=y+1;                     /* shift all elements by 1    */
   indexvec=ones(1,dim)|indexvec;    /* insert 1's at the 1st row  */
/* divide between 0 and 1 into nn segments and sample from each segment */
   d=1/nn;
   U=(d*indexvec-d*(indexvec-1)).*rndu(nn^dim,dim)+d*(indexvec-1);
   do while NOT(U>0);                                                      /* check if U contains 0's */
      U=(d*indexvec-d*(indexvec-1)).*rndu(nn^dim,dim)+d*(indexvec-1);
   endo;
/* randomize the order as a row */
   x=zeros(nn^dim,dim);
   i=1;
   do while i<=nn^dim-1;
      index=ceil((nn^dim-(i-1))*rndu(1,1));
      x[i,.]=U[index,.];
      if index==1;
         U=U[2:rows(U),.];
      elseif index==rows(U);
         U=U[1:rows(U)-1,.];
      else;
         U=U[1:(index-1) (index+1):rows(U),.];
      endif;
      i=i+1;
   endo;
   x[nn^dim,.]=U;
   retp(x);
endp;