new; cls;
let data[100,5]=
    -1.602      0.684     -0.263     -0.775     -0.289 
     0.490      0.369     -0.530      0.663     -0.979 
    -0.051     -0.527     -0.145     -0.873      1.802 
     1.205      1.464     -0.510      0.452     -1.026 
    -2.806     -0.493      0.555     -0.751      0.931 
    -1.240      0.530      0.178      0.658     -0.743 
     0.216     -0.053      0.220      1.315     -0.525 
    -1.039      0.961     -0.662     -1.329      0.626 
    -0.001      0.698      0.510     -0.001      1.322 
    -1.263     -2.696     -0.814      0.145     -0.257 
    -1.582     -0.056     -0.007     -0.440     -0.266 
     0.140      0.214      0.933     -0.213      0.406 
     0.707      1.299      1.638      1.336     -0.933 
    -1.851      0.383     -1.439      1.070      0.712 
    -0.220      0.116     -0.084      1.687      0.257 
     1.302     -1.466     -1.040     -0.669     -1.051 
     0.708     -0.079      0.716     -0.808      0.589 
    -0.481     -0.006      0.825     -0.167     -0.680 
     0.756     -0.750      2.197      0.135      1.333 
    -1.189      0.236     -0.034      0.519     -1.853 
    -0.098      0.306      0.611     -0.812      0.276 
     0.503      0.350      0.450     -1.392     -1.405 
    -0.034      0.137      1.073     -0.600      0.458 
     0.212     -2.381      2.526      0.439     -0.282 
    -0.010     -0.627     -0.276     -1.924      0.534 
     1.322      0.201     -0.689      0.468     -1.421 
     1.726      0.336      1.028     -1.671     -0.192 
    -0.418     -1.417     -2.673      0.550     -1.035 
    -0.304     -0.182     -0.741     -1.492     -0.232 
    -0.458      0.035      0.232      1.125      0.486 
     0.886     -0.779     -0.443     -1.671     -0.687 
     1.704      0.743      1.792     -0.237     -1.316 
     1.699     -0.182     -0.175      0.151      1.039 
    -1.119      1.103     -0.219      1.444      0.410 
    -1.593     -2.449     -1.246     -0.297      1.495 
    -0.731     -0.992      2.081     -0.811      0.485 
     2.461      1.648     -1.823     -0.365      1.250 
     1.453     -0.755     -0.713     -0.742     -0.560 
    -0.565     -0.544     -0.023     -0.582     -0.861 
    -0.859      0.482     -0.776     -0.506     -0.533 
    -0.377     -0.690     -0.496     -0.153      0.457 
    -1.553      0.335     -0.098     -1.799     -0.019 
    -1.249     -0.606     -0.349      0.282     -1.031 
     1.029      0.064      0.630     -1.654     -0.799 
     0.337      0.080      0.279      0.610      0.272 
     0.848      0.085     -1.448      0.344      0.007 
    -0.226     -1.122      1.101     -0.386      0.210 
     0.098     -0.262     -0.730      0.045     -0.854 
    -0.490      0.163     -0.008      1.114      0.340 
    -0.483      1.215      0.058     -0.081      0.585 
    -0.605      1.681      0.378     -0.109     -0.422 
    -0.525      0.014      0.070      1.335     -0.391 
    -0.154     -0.208      0.671      0.204      1.145 
    -0.389      0.533     -0.815     -0.665      1.056 
    -0.083     -1.154     -0.682      1.231      0.472 
     0.979     -0.110      1.438      1.514     -1.816 
     0.950     -0.076     -1.114     -0.384     -0.363 
    -0.888      1.137     -1.549      1.789     -1.359 
     0.296     -0.138     -0.520     -0.062     -0.359 
     1.270     -0.756     -0.755      0.607     -0.020 
     0.630      0.936      2.540      1.115     -0.684 
     0.483     -0.193      0.430     -0.884     -0.256 
     0.524     -1.192      0.246      0.504     -0.431 
    -0.887     -0.500     -0.756      1.280     -0.174 
    -0.672     -0.819     -0.722      0.362     -1.332 
    -0.627     -0.126     -0.393     -2.411     -0.110 
     0.622     -0.245     -1.898     -0.353      1.062 
     0.728      1.065     -0.345      0.083     -1.960 
    -0.705     -0.662      0.149     -0.225      0.168 
     1.585      1.748      0.130      1.646      0.217 
    -1.223     -0.268      0.074     -0.125      2.589 
    -0.771      0.508     -0.549      0.565      1.048 
     0.346      1.199     -0.407      0.632      1.158 
     1.117     -0.619      1.048      0.586      0.667 
     0.259     -0.834     -1.289      1.047     -0.732 
     0.050     -0.485     -0.374      0.205      0.128 
    -0.145      1.986     -0.148      1.492     -1.131 
    -0.103     -0.849      0.180      1.350     -1.570 
    -0.001     -0.844     -0.742     -0.452     -1.305 
     0.811     -0.285      0.644      0.395     -0.040 
    -0.410     -0.523     -0.852     -0.251     -0.751 
     1.514      1.354     -0.862      0.249      0.568 
    -0.787     -1.879     -1.330      1.143     -1.010 
    -0.767      0.811      0.862      1.037      1.083 
     1.042     -0.147     -1.601     -0.668      1.113 
    -0.208      0.494      0.274      0.882      1.368 
     0.870      0.726      2.219      0.755     -0.001 
    -0.466      0.125     -1.385     -0.161     -0.142 
     0.779     -0.455      0.183     -1.151      0.092 
     0.200      0.462     -0.979     -0.638      0.841 
    -0.831     -1.302     -0.987      0.861     -1.200 
     0.580      0.947      1.998     -0.120     -0.832 
    -0.688     -0.832      0.880      0.507      1.632 
    -0.647     -1.074     -1.172     -0.062      1.673 
    -0.408      1.557      1.177      0.608     -0.346 
    -0.902      0.091     -0.053      1.309     -0.320 
     1.913      0.578      2.655      0.593      0.074 
     1.370      0.208     -0.405     -1.080      1.354 
     1.616     -0.282      0.571      1.259      0.823 
     0.227     -2.435     -1.352     -0.174      0.143
;
d={0.1,0,0.1,0,0.05}; /* d=0 for non time series data */
pp=0.20;
label={"A","B","C","D","E"};
call ToddII(data,d,pp,label);




proc ToddII(data,d,pp,label);
   local nn,k,j,ddata,corr,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;
   corr=eye(k);
   j=1;
   do while j<=k;
      i=j+1;
      do while i<=k;
         corr[i,j]=spearman1(ddata[.,j],ddata[.,i]);
         corr[j,i]=corr[i,j];
         i=i+1;
      endo;
      j=j+1;
   endo;
   nn=rows(ddata);
/* 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 $label';
   print "d=";
   print/rz d';
   print/rz "n=" nn;
   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); 
   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); 
   retp(p1);
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;