/* Empirical Copula */
/* Notice that the plots of pdf and cdf have no points between the adjacent two points. */
/* We just connect the points just for convenience in the surface plots.                */
new; cls;
n=100;
mu={5 4};
Sig={1   0.9,
     0.9   1};
data=mu+rndn(n,2)*chol(Sig);
x=data[.,1]; y=data[.,2];
call cdfCemp(x,y);
call pdfCemp(x,y);

nn=1000;
U=Cemp(x,y);
print U;
library pgraph;
graphset;
pqgwin auto;
_plctrl=-1;
xy(U[.,1],U[.,2]);


proc cdfCemp(x,y);
   local n,xi,yj,m,i,j;
   n=rows(x);
   xi=sortc(x,1); yj=sortc(y,1);
   m=zeros(n,n);
   j=1;
   do while j<=n;
      i=1;
      do while i<=n;
         m[i,j]=1/n*sumc(x.<=xi[i] .and y.<=yj[j]);
         i=i+1;
      endo;
      j=j+1;
   endo;
   library pgraph;
   graphset;
   pqgwin auto;
   title("cdf");
   surface(seqa(1,1,n)'/n,seqa(1,1,n)/n,m');
   retp(m);
endp;

proc pdfCemp(x,y);
   local n,xi,yj,k,m;
   n=rows(x);
   xi=rankindx(x,1); yj=rankindx(y,1);
   m=zeros(n,n);
   k=1;
   do while k<=n;
      m[xi[k],yj[k]]=1/n;
      k=k+1;
   endo;
   library pgraph;
   graphset;
   pqgwin auto;
   title("pdf");
   surface(seqa(1,1,n)'/n,seqa(1,1,n)/n,m');
   retp(m);
endp;

proc Cemp(x,y);
   local n,xi,yj,index,U;
   n=rows(x);
   xi=rankindx(x,1); yj=rankindx(y,1);
   index=xi~yj;
   U=index[rankindx(rndu(n,1),1),.]/n;
   retp(U);
endp;