/* Husler-Reiss Copula                                                                                */
/* C(u,v)=exp(ln(u).*cdfn(1/r+1/2*r*ln(ln(u)./ln(v)))+ln(v).*cdfn(1/r+1/2*r*ln(ln(v)./ln(u))))        */
/* Solve  q=exp(ln(u).*cdfn(1/r+1/2*r*ln(ln(u)./ln(v)))+ln(v).*cdfn(1/r+1/2*r*ln(ln(v)./ln(u))))      */
/*          .*( 1/u.*cdfn(1/r+1/2*r*ln(ln(u)./ln(v)))                                                 */
/*          +ln(u).*pdfn(1/r+1/2*r*ln(ln(u)./ln(v)))*(1/2)*r.*ln(v)./ln(u).*(1/u)./ln(v)              */
/*          +ln(v).*pdfn(1/r+1/2*r*ln(ln(v)./ln(u)))*(1/2)*r.*ln(u)./ln(v).*(-1/u).*ln(v)./ln(u)^2 ); */
/* in terms of v numerically. This is very slow.                                                      */
new; cls;
r=10;
n=1000;
U=Chr(r,n);

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



proc Chr(r,n);
   local u,q,v;
   if r<0;
      errorlog "ERROR: Parameter r must be r>=0";
      retp(".");
   endif;
   u=rndu(n,1);
   q=rndu(n,1);
   v=calcv(u,q,r);
   retp(u~v);
endp;

proc calcv(u,q,r);
   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,r));
      vstar[i]=v[minindc(x)];
      v=seqa(vstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-Cu(u[i],v,r));
      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,r));
         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,r)=exp(ln(u).*cdfn(1/r+1/2*r*ln(ln(u)./ln(v)))+ln(v).*cdfn(1/r+1/2*r*ln(ln(v)./ln(u))))
            .*( 1/u.*cdfn(1/r+1/2*r*ln(ln(u)./ln(v)))
            +ln(u).*pdfn(1/r+1/2*r*ln(ln(u)./ln(v)))*(1/2)*r.*ln(v)./ln(u).*(1/u)./ln(v)
            +ln(v).*pdfn(1/r+1/2*r*ln(ln(v)./ln(u)))*(1/2)*r.*ln(u)./ln(v).*(-1/u).*ln(v)./ln(u)^2 );