/* Runs Test[of NIST SP800-22] for 01 random sequence */
/* Based on Katsuichi Hirose's paper(2004) in Japanese. */
new; cls;
x=rndu(10000,1);
x=(x.>0.5);          /* x=rev(x); if you need the reverse order of it. */
call runstest(x);





proc runstest(x);
   local n,freq,Vn,i,Z,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   freq=sumc(x)/n;
   if abs(freq-0.5)>=2/sqrt(n);
      retp(".");
   endif;
   Vn=1;
   i=1;
   do while i<=n-1;
      if x[i]/=x[i+1];
         Vn=Vn+1;
      endif;
      i=i+1;
   endo;
   Z=(Vn-2*n*freq*(1-freq))/(2*sqrt(n)*freq*(1-freq));
   pval=2*(1-cdfn(abs(Z)));   /* Or, pval=erfc(abs(Vn-2*n*freq*(1-freq))/(2*sqrt(2*n)*freq*(1-freq))); */
/* results */
   print "Runs Test:";
   print/rz "n   =" n;
   print    "Z   =" Z;
   print    "pval=" pval;
   retp(Z);
endp;