/* 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.                                                */                 
n=1000;
nu=5;      /* nu=1,2,3,... */
rho=0.9;   /* -1<=rho<=1 */
U=Ctev(nu,rho,n);

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



proc Ctev(nu,rho,n);
   local u,q,v;
   if rho<-1 or rho>1;
      errorlog "ERROR: Parameter rho must be -1<=rho<=1";
      retp(".");
   endif;
   u=rndu(n,1);
   q=rndu(n,1);
   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 );