/* Tawn Copula(Asymmetric Logistic) in Extreme Value Family                */
/* C(u,v)=exp((1-a)*ln(u)+(1-b)*ln(v)-((-a*ln(u))^c+(-b*ln(v))^c)^(1/c))   */
/* Solve  q=exp((1-a)*ln(u)+(1-b)*ln(v)-((-a*ln(u))^c+(-b*ln(v))^c)^(1/c)) */
/*          .*( (1-a)/u-1/c*((-a*ln(u))^c+(-b*ln(v))^c)^(1/c-1)            */
/*               .*c.*(-a*ln(u))^(c-1).*(-a/u) )                  (=Cu)    */
/* in terms of v numerically. This is very slow. (Latin Hypercube version) */
new; cls;
a=0.9;
b=0.9;
c=2;
nn=30;             /* n=nn^2 */
U=Ctawnh(a,b,c,nn);

library pgraph;
graphset;
_plctrl=-1;
xy(U[.,1],U[.,2]);



proc Ctawnh(a,b,c,nn);
   local dim,UU,u,q,v;
   if a<0 or b<0 or a>1 or b>1;
      errorlog "ERROR: Parameter a and b must be 0<=a<=1 and 0<=b<=1";
      retp(".");
   endif;
   if a==0 and b==0;
      errorlog "ERROR: Parameter a and b must not equal to 0 at the same time.";
      retp(".");
   endif;
   if c<1;
      errorlog "ERROR: Parameter c must be c>=1";
      retp(".");
   endif;
   dim=2;
   UU=LHss(nn,dim);
   u=UU[.,1];
   q=UU[.,2];
   v=calcv(u,q,a,b,c);
   retp(u~v);
endp;

proc calcv(u,q,a,b,c);
   local vstar,vstarstar,i,v,x;
   vstar=zeros(rows(q),1);
   i=1;
   do while i<=rows(q);
      v=seqa(0.0001,0.0001,10000);
      x=abs(q[i]-Cu(u[i],v,a,b,c));
      vstar[i]=v[minindc(x)];
      v=seqa(vstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-Cu(u[i],v,a,b,c));
      vstarstar=v[minindc(x)];
      if vstarstar==vstar[i]-0.00005;
         v=seqa(0.00000001,0.00000001,10000);
         x=abs(q[i]-Cu(u[i],v,a,b,c));
         vstar[i]=v[minindc(x)];
      else;
         vstar[i]=vstarstar;
      endif;
      if vstarstar>1;
         vstar[i]=1;
      endif;
      i=i+1;
   endo;
   retp(vstar);
endp;

fn Cu(u,v,a,b,c)=exp((1-a)*ln(u)+(1-b)*ln(v)-((-a*ln(u))^c+(-b*ln(v))^c)^(1/c))
                 .*( (1-a)/u-1/c*((-a*ln(u))^c+(-b*ln(v))^c)^(1/c-1)
                     .*c.*(-a*ln(u))^(c-1).*(-a/u) );




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