SASでWald-Wolfowitz test

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

この本の2-3節に出てくる連の検定(runs test)は、Rではtseriesライブラリのruns.test関数を使えば実行できると書いてあります。
SASではどうやるんだろうかと調べたところ、SAS社のサポートページにコードが書いてあったので試してみました。

33092 - Wald-Wolfowitz (or Runs) test for randomness

こちらのページでは検定の名称がWald-Wolfowitz testと紹介されていたので、このエントリのタイトルもそちらを使いました。
コードはほぼコピペで、データセット名と変数名をマクロのキーワードパラメータとして渡せるようにする部分だけ変更してあります。

使用するデータ

本のサポートページからダウンロードしたデータから、4銘柄の収益率の時系列データ return4 を使います。

time	x5202	x7272	x4927	x4502
1	1.80185	-1.47713	3.31831	-0.55788
2	0.00000	-0.59702	0.03932	-0.13996
3	1.76996	0.89419	-1.74469	-1.98026
4	5.12933	5.34423	1.39029	1.27752
5	3.27898	2.08776	-0.51414	0.14094

以下略

コード

%macro Wald_Wolfowitz_Test(in1=, var1=);

    proc standard data = &in1. out = __m1 mean = 0;
        var &var1.;
    run;

    data __runcount;
        set __m1 nobs = nobs;
        if &var1. eq 0 then delete;
        if &var1. >  0 then n + 1;
        if &var1. <  0 then m + 1;
        retain runs 0 numpos 0 numneg 0;
        previous = lag(&var1.);

        if _n_ eq 1 then do;
            runs = 1;
            prevpos = .;
            currpos = .;
            prevneg = .;
            currneg = .;
        end;
        else do;
            prevpos = (previous > 0);
            currpos = (&var1. > 0 );
            prevneg = (previous < 0);
            currneg = (&var1. < 0);

            if _n_ eq 2 and (currpos and prevpos) then numpos + 1;
            else if _n_ eq 2 and (currpos and prevneg) then numneg + 1;
            else if _n_ eq 2 and (currneg and prevpos) then numpos + 1;
            else if _n_ eq 2 and (currneg and prevneg) then numneg + 1;

            if currpos and prevneg then do;
                runs + 1;
                numpos + 1;
            end;

            if currneg and prevpos then do;
                runs + 1;
                numneg + 1;
            end;
        end;
    run;

    data __runcount;
        set __runcount end = last;
        if last;
    run;

    data __waldwolf;
        label z = 'Wald-Wolfowitz Z'
              pvalue = 'Pr > |Z|';
        set __runcount;
        mu = ((2 * n * m) / (n + m)) + 1;
        sigmasq = ((2 * n * m) * (2 * n * m - (n + m))) / (((n + m) ** 2) * (n + m - 1));
        sigma = sqrt(sigmasq);
        drop sigmasq;

        if N GE 50 then Z = (Runs - mu) / sigma;
        else if Runs - mu LT 0 then Z = (Runs - mu + 0.5) / sigma;
        else Z = (Runs - mu - 0.5) / sigma;

        pvalue = 2 * (1 - probnorm(abs(Z)));
    run;

    title  'Wald-Wolfowitz Test for Randomness';
    title2 'H0: The data are random';
    proc print data = __waldwolf label noobs;
        var z pvalue;
        format pvalue pvalue.;
    run;
    title;
    title2;

%mend Wald_Wolfowitz_Test;

%Wald_Wolfowitz_Test(in1=return4, var1=x5202);
%Wald_Wolfowitz_Test(in1=return4, var1=x7272);
%Wald_Wolfowitz_Test(in1=return4, var1=x4927);
%Wald_Wolfowitz_Test(in1=return4, var1=x4502);

実行結果

Wald-Wolfowitz Z	Pr > |Z|
-1.37058	0.1705

Wald-Wolfowitz Z	Pr > |Z|
-0.17132	0.8640

Wald-Wolfowitz Z	Pr > |Z|
1.65012	0.0989

Wald-Wolfowitz Z	Pr > |Z|
-0.12194	0.9029

本と同じ結果が得られました。