new; cls;
let data[50,6]=
 117.58990   90.32080  223.27970  105.08481   94.15478  149.78601 
 117.41000   90.33525  223.32556  104.30881   94.11623  149.47467 
 116.94000   90.49987  223.05136  104.20602   94.01833  149.22713 
 117.05990   90.69801  223.37370  103.31854   93.99382  149.28649 
 117.09990   90.65874  223.45003  103.28091   94.03349  149.59512 
 118.26000   91.21394  224.32739  104.70120   94.18605  150.41489 
 117.69990   90.95848  224.26539  104.25146   94.17499  150.29100 
 117.83000   90.78802  224.44258  104.31126   94.18865  150.29217 
 117.91990   90.57428  224.59024  104.43707   94.87481  151.23227 
 117.58990   90.23849  224.73782  103.92391   94.75415  151.00895 
 118.15000   90.04211  224.56770  103.82250   94.89198  151.23200 
 117.58990   89.97979  222.93869  103.38483   94.57126  150.63266 
 117.95990   90.26292  222.80266  103.57354   94.68606  151.28357 
 118.13990   90.57786  223.10720  103.46812   94.70892  151.18363 
 117.66000   90.45701  223.01276  102.63433   94.57439  150.91072 
 118.05000   90.92211  223.85822  102.89375   94.87262  151.22205 
 117.83990   90.87813  223.78975  102.87202   94.97856  151.34178 
 116.66000   90.48150  223.38057  102.31538   95.15498  150.96971 
 115.86000   90.19701  223.78359  102.10628   95.83919  151.66074 
 116.06990   90.37202  224.86222  102.50808   96.21179  152.41139 
 116.05990   90.97936  226.61856  102.69879   96.55566  153.19907 
 116.37990   91.14874  226.41709  102.31200   96.19764  153.00466 
 115.77990   91.25772  227.51908  101.45452   96.58789  153.29259 
 115.41990   91.02013  228.45061  100.82982   96.72329  153.90089 
 115.33990   91.00318  228.28073  101.13100   96.56723  153.62121 
 114.81990   90.29437  226.62004  100.59567   96.37393  152.96307 
 115.18000   90.33567  226.50147  100.39223   96.36075  153.08574 
 115.19000   90.87339  226.24468  100.26985   96.40944  153.07599 
 116.40000   91.42056  227.38740  101.32312   96.61355  153.63636 
 116.97000   91.72787  228.87520  101.92576   97.24809  154.80980 
 116.77990   91.80068  230.13815  101.37144   97.36527  155.07203 
 117.51000   92.36286  231.17742  101.60830   97.45397  155.28947 
 117.80990   92.22159  231.07234  101.80600   97.02677  154.84933 
 118.12000   92.29897  230.55843  102.05633   96.74830  154.50096 
 118.11000   92.11399  230.05466  102.08297   96.65303  154.68867 
 118.09990   92.53127  232.44422  102.50838   97.24959  155.75015 
 118.38990   93.01894  232.62432  103.07322   97.25614  155.99053 
 118.33990   92.80215  232.02904  102.49428   97.27900  155.86548 
 118.75000   93.05250  232.67875  102.61839   97.46389  156.06125 
 118.76990   93.10372  234.55868  102.63559   97.44023  155.92112 
 119.12000   93.04463  232.77239  102.65426   97.43170  155.97573 
 118.76990   93.15123  232.38519  102.26442   96.88384  155.80236 
 118.88990   93.91113  233.28576  102.49129   97.33904  156.28077 
 119.01000   93.88699  233.10489  102.16328   97.67728  157.03369 
 119.01990   93.31160  231.23186  101.12141   96.56004  155.71374 
 118.58990   92.36967  228.86665  101.14277   95.95428  154.17873 
 119.40000   93.20364  231.54048  101.51335   96.21273  155.17224 
 119.65000   92.89626  231.10397  101.78647   95.97337  154.74335 
 120.44990   94.07137  234.37142  102.39726   96.44479  155.27197 
 120.27990   94.19119  235.58021  102.88247   96.40903  155.34149 
;
label={"USDJPY","AUDJPY","GBPJPY","CADJPY","CHFJPY","EURJPY"};
data=ln(data[2:rows(data),.])-ln(data[1:rows(data)-1,.]);
lag=2;
d=0;    /* fractional difference */
pp=0.10;
call DAGs_VARrho(data,lag,d,pp,label);



proc DAGs_VARrho(data,lag,d,pp,label);
   local nn,k,j,ddata,corr,xlag,y,x1,b,e,maxdim,mindim,dim,N,low,up,x,p0,p1,p,m0,m1,m,i,l,m3;
   nn=rows(data);
   k=cols(data);
   if rows(label)<k;
      errorlog "ERROR: Set the same number of names as label vector.";
      retp(".");
   endif;
   if rows(d)==1;
      d=d*ones(k,1);
   endif;
   if d==zeros(k,1);   
      ddata=data;
   else;
      ddata={};
      j=1;
      do while j<=k;
         ddata=ddata~fracdiff(data[.,j],d[j]);
         j=j+1;
      endo;
      ddata=ddata[2:nn,.];
   endif;
/* VAR */
   nn=rows(ddata); 
   xlag=ddata[lag:nn-1,.];
   i=2; 
   do while i<=lag;
      xlag=xlag~ddata[lag+1-i:nn-i,.];
      i=i+1;
   endo;
   y=ddata[lag+1:nn,.];
   x1=ones(nn-lag,1)~xlag;
   b=inv(x1'x1)*x1'y;
   e=y-x1*b;
   corr=eye(k);
   j=1;
   do while j<=k;
      i=j+1;
      do while i<=k;
         corr[i,j]=spearman1(e[.,j],e[.,i]);
         corr[j,i]=corr[i,j];
         i=i+1;
      endo;
      j=j+1;
   endo;
/* DAGs */
   format /rd 10,4;
   print "Unconditional correlation:";
   print $label';;
   p0=zerocorrrho(corr,nn);
   m0=p0.<pp;
   print;
   print "Conditional correlations:";
   maxdim=cols(corr); mindim=3;
   dim=mindim;
   do while dim<=maxdim;
      N=maxdim-(dim-1)*1;
      low=1+seqa(0,1,dim);
      up=low+N-1;
      x=low;
      print $label[x]';;
      p1=mparcorrrho(corr,nn,x);
      print;
      p=zeros(maxdim,maxdim);
      p[x,x]=p1;
      m1=p.<pp;
      m=m0.*m1;
      do while not x==up;
         x=recDindex(x,up,1);
         print $label[x]';;
         p1=mparcorrrho(corr,nn,x);
         print;
         p=zeros(maxdim,maxdim);
         p[x,x]=p1;
         m1=p.<pp;
         m=m.*m1;
      endo;
      dim=dim+1;
   endo;
   k=maxdim;
/* directing threesome */
   format /mb1 /ros 16,8;
   print;
   print "DAGs from Spearman's correlation matrix of residuals of VAR:";
   print $label';
   print "d=";
   print/rz d';
   print/lz "n=" nn "as in VAR";; print/lz "   lag=" lag;
   print/rz pp*100 "% significance level:";
   print;
   m3=zeros(k,k);
   i=1;
   do while i<=k-1;
      j=i+1;
      do while j<=k;
         l=1;
         do while l<=k;
            if i/=l and j/=l;
               if m[i,l]==1 and m[l,j]==1 and m0[i,j]==0;
                  print $label[i] "            ->" $label[l] "           <-" $label[j];
                  m3[i,l]=1; m3[l,i]=1; m3[l,j]=1; m3[j,l]=1;
               endif;
            endif;
            l=l+1;
         endo;
         j=j+1;
      endo;
      i=i+1;
   endo;
/* the other edges left */
   i=1;
   do while i<=k-1;
      j=i+1;
      do while j<=k;
         if m[i,j]==1 and m3[i,j]==0;
            print $label[i] "            - - -" $label[j];
         endif;
         j=j+1;
      endo;
      i=i+1;
   endo;
   retp(m);
endp;

proc zerocorrrho(corr,n);
   local k,j,i,order,t0,p0;
   k=cols(corr);
   order=0;
   t0=1./zeros(k,k);
   j=1;
   do while j<=k;
      i=1;
      do while i<=k;
         if i/=j;
            t0[i,j]=corr[i,j]*sqrt((n-order-2)/(1-corr[i,j]^2));
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   p0=2*cdftc(abs(t0),n-order-2); 
   print corr;
   print "p-values=" p0;
   retp(p0);
endp;

proc mparcorrrho(corr,n,index);
   local k,j,i,c,parcorr,order,t1,p1,cindex;
   corr=corr[index,index];
   k=rows(index);
   c=inv(corr);
   parcorr=zeros(k,k);
   j=1;
   do while j<=k;
      i=1;
      do while i<=k;
         parcorr[i,j]=c[i,j]/(sqrt(c[i,i])*sqrt(c[j,j]));
         i=i+1;
      endo;
      j=j+1;
   endo;
   parcorr=-parcorr;
   i=1;
   do while i<=k;
      parcorr[i,i]=-parcorr[i,i];
      i=i+1;
   endo;
   order=k-2;
   t1=1./zeros(k,k);
   j=1;
   do while j<=k;
      i=1;
      do while i<=k;
         if i/=j;
            t1[i,j]=parcorr[i,j]*sqrt((n-order-2)/(1-parcorr[i,j]^2));
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   p1=2*cdftc(abs(t1),n-order-2); 
   print parcorr;
   print "p-values=" p1;

   retp(p1);
endp;

proc spearman1(y,x);
   local n,u,v,t1,t2,Tx,Ty,Sx2,Sy2,rho;
   n=rows(x);
   u=ranksame(x);
   v=ranksame(y);
   t1=calct(x);
   t2=calct(y);
   Tx=sumc((t1^3-t1)/12);
   Ty=sumc((t2^3-t2)/12);
   Sx2=(n^3-n)/12-Tx;
   Sy2=(n^3-n)/12-Ty;
   rho=(Sx2+Sy2-sumc((u-v)^2))/(2*sqrt(Sx2*Sy2));
   retp(rho);
endp;

proc ranksame(x);
   local n,y,i,j,k,index;
   n=rows(x);
   y=zeros(n,1);
   i=1;
   do while i<=n;
      index=sumc(x.==minc(x));
      k=0;
      j=0;
      do while j<=index-1;
         k=k+(i+j);
         j=j+1;
      endo;
      k=k/index;
      y=y+k*(x.==minc(x));
      x=x+(x.==minc(x))*1e256;
      i=i+index;
   endo;
   retp(y);
endp;

proc calct(x);
   local n,i,j,t;
   x=sortc(x,1);
   n=rows(x);
   t=1;
   i=1; 
   j=1;
   do while i<=n-1;
      if x[i]==x[i+1];
         t[j]=t[j]+1;
      else;
         t=t|1;
         j=j+1;
      endif;
      i=i+1;
   endo;
   retp(t);
endp;


proc recDindex(x,up,mindist);
   local dim,j;
   dim=rows(x);
   x[dim]=x[dim]+1;
   j=dim;
   do until j==0;
      if x[j]==up[j]+1;
         x[j-1]=x[j-1]+1;
         if x[j-1]/=up[j-1]+1;
            x[j-1:dim]=x[j-1]+seqa(0,mindist,dim+1-(j-1));
         endif;
      endif;
      j=j-1;
   endo;
   retp(x);
endp;

proc fracdiff(y,d);
   local n,t,y1,x;
   n=rows(y);
   t=seqa(1,1,n-1);
   y1=1|cumprodc((t-1-d)./t);
   x=conv(y,y1,0,n);
   retp(x);
endp;