/* Multivariate Nelsen's no.13 Copula in Archimedean Family ( phi(u)=(1-ln(u))^theta-1 ) */
/*                                                                   Wait for a while... */
new; cls;
theta=3; /* theta>0 */
dim=3;
n=1000;
U=MCno13(theta,dim,n);

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

proc calct2(q,theta);
   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)./phi1(t,theta)));
      tstar[i]=t[minindc(x)];
      t=seqa(tstar[i]-0.00005,0.00000001,10000);
      x=abs(q[i]-(t-phi(t,theta)./phi1(t,theta)));
      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)./phi1(t,theta)));
         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)=(1-ln(u))^theta-1;

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

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