/* Multivariate Fang-Fang-von Rosen Copula in Archimedean Family */
/* phi(u)=ln((1-theta*(1-u^(1/k)))./u^(1/k))                     */
/*                 (Latin Hypercube version) Wait for a while... */
new; cls;
theta=0.99;
k=0.5;
dim=3;
nn=10;              /* n=nn^dim */
U=MCffrh(theta,k,dim,nn);

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 MCffrh(theta,k,dim,nn);
   local s,q,t,U,i,a;
   if theta<-1 or theta>=1;
      errorlog "ERROR: Parameter theta must be -1<=theta<1.";
      retp(".");
   endif;
   if k<=0;
      errorlog "ERROR: Parameter k must be positive.";
      retp(".");
   endif;
   s=LHss(nn,dim);
   q=s[.,dim];
   U=zeros(nn^dim,dim);
   i=1;
      t=calct2(q,theta,k);
      a=phi_1(s[.,i].*phi(t,theta,k),theta,k);
      U[.,dim-i+1]=phi_1((1-s[.,i]).*phi(t,theta,k),theta,k);
   i=2;
   do while i<=dim-1;
      t=calct2(a,theta,k);
      a=phi_1(s[.,i].*phi(t,theta,k),theta,k);
      U[.,dim-i+1]=phi_1((1-s[.,i]).*phi(t,theta,k),theta,k);
      i=i+1;
   endo;
   U[.,1]=a;
   retp(U);
endp;

proc calct2(q,theta,k);
   local tstar,tstarstar,i,t,x;
   tstar=zeros(rows(q),1);
   i=1;
   do while i<=rows(q);
      t=seqa(0.0001,0.0001,10000);
      x=abs(q[i]-(t-phi(t,theta,k)./phi1(t,theta,k)));
      tstar[i]=t[minindc(x)];
      t=seqa(tstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-(t-phi(t,theta,k)./phi1(t,theta,k)));
      tstarstar=t[minindc(x)];
      if tstarstar==tstar[i]-0.00005;
         t=seqa(0.00000001,0.00000001,10000);
         x=abs(q[i]-(t-phi(t,theta,k)./phi1(t,theta,k)));
         tstar[i]=t[minindc(x)];
      else;
         tstar[i]=tstarstar;
      endif;
      if tstarstar>1;
         tstar[i]=1;
      endif;
      i=i+1;
   endo;
   retp(tstar);
endp;


fn phi(u,theta,k)=ln((1-theta*(1-u^(1/k)))./u^(1/k));

fn phi1(u,theta,k)=1/k*(theta*u^(1/k-1)./(1-theta*(1-u^(1/k)))-1/u);

fn phi_1(u,theta,k)=((1-theta)/(exp(u)-theta))^k;


/*
**  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;