/* t-EV Copula                                                                                  */
/* C(u,v)=phi_1((phi(u)+phi(v)).*A(phi(u)./(phi(u)+phi(v))))                                    */
/*  where phi(t)=ln(t) and phi_1(t)=exp(t) for EV Copula in Archimax Family                     */
/*        A(w)=w.*cdft(((w./(1-w))^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)                   */
/*            +(1-w).*cdft((((1-w)./w)^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)               */
/* Solve  q=exp((ln(u)+ln(v)).*A(ln(v)./(ln(u)+ln(v))))                                         */
/*         .*( (1/u).*A(ln(v)./(ln(u)+ln(v)))                                                   */
/*             +(ln(u)+ln(v)).*A1(ln(v)./(ln(u)+ln(v))).*(-1/u).*ln(v)./(ln(u)+ln(v))^2 ) =(Cu) */
/*  where                                                                                       */
/*        A1(w)=cdft(((w./(1-w))^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)                     */
/*             +w.*pdft(((w./(1-w))^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)                  */
/*            .*sqrt((nu+1)/(1-rho^2)).*(1/nu).*(w./(1-w))^(1/nu-1)./(1-w)^2                    */
/*            -cdft((((1-w)./w)^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)                      */
/*            +(1-w).*pdft((((1-w)./w)^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)               */
/*            .*sqrt((nu+1)/(1-rho^2)).*(1/nu).*((1-w)./w)^(1/nu-1).*(-1/w^2)                   */
/* in terms of v numerically. This is very slow.                     (Latin Hypercube version)  */                 
nn=30;     /* n=nn^2 */
nu=5;      /* nu=1,2,3,... */
rho=0.9;   /* -1<=rho<=1 */
U=Ctevh(nu,rho,nn);

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



proc Ctevh(nu,rho,nn);
   local dim,UU,u,q,v;
   if rho<-1 or rho>1;
      errorlog "ERROR: Parameter rho must be -1<=rho<=1";
      retp(".");
   endif;
   dim=2;
   UU=LHss(nn,dim);
   u=UU[.,1];
   q=UU[.,2];
   v=calcv(u,q);
   retp(u~v);
endp;

proc calcv(u,q);
   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));
      vstar[i]=v[minindc(x)];
      v=seqa(vstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-Cu(u[i],v));
      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));
         vstar[i]=v[minindc(x)];
      else;
         vstar[i]=vstarstar;
      endif;
      if vstarstar>1;
         vstar[i]=1;
      endif;
      i=i+1;
   endo;
   retp(vstar);
endp;

proc cdft(x,nu);
   retp( 1-cdftc(x,nu) );
endp;

proc pdft(x,nu);
   retp(gamma((nu+1)/2)/((pi*nu)^(1/2)*gamma(nu/2)*(1+x^2/nu)^((nu+1)/2)));
endp;


fn A(w)=w.*cdft(((w./(1-w))^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)
        +(1-w).*cdft((((1-w)./w)^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1);
fn A1(w)=cdft(((w./(1-w))^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)
        +w.*pdft(((w./(1-w))^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)
         .*sqrt((nu+1)/(1-rho^2)).*(1/nu).*(w./(1-w))^(1/nu-1)./(1-w)^2
        -cdft((((1-w)./w)^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)
        +(1-w).*pdft((((1-w)./w)^(1/nu)-rho)/sqrt(1-rho^2)*sqrt(nu+1),nu+1)
         .*sqrt((nu+1)/(1-rho^2)).*(1/nu).*((1-w)./w)^(1/nu-1).*(-1/w^2);
fn Cu(u,v)=exp((ln(u)+ln(v)).*A(ln(v)./(ln(u)+ln(v))))
          .*( (1/u).*A(ln(v)./(ln(u)+ln(v)))
             +(ln(u)+ln(v)).*A1(ln(v)./(ln(u)+ln(v))).*(-1/u).*ln(v)./(ln(u)+ln(v))^2 );

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