Excel VBAでダービン・ワトソン統計量を求める関数を作ってみた

前回、エクセルExcel分析ツールの回帰分析で出力できるグラフの読み方を学んだわけですが、そのうちの残差グラフには宿題がありました。

「残差の分布から回帰モデルの妥当性を検証」ってどうやるの？　です。

まだ、これについての解説が載っている本に巡り合っていないので、「回帰分析　残差」でググったサイトを手当たり次第に読んでみました。

が、残念ながら、わかりやすい説明がされているところはあまり多くありませんでした。

http://www.aoni.waseda.jp/abek/document/regression-3.html

概論としては、ここの説明が一番分かりやすかったです。かいつまんで紹介すると、

分布に規則性が見られる場合　→　回帰式そのものを見直す必要がある
残差のバラツキが不均一な場合　→　加重最小二乗法，Box-Cox変換などで対応する
データの並び方に意味があり，隣り合うデータから生じた残差に相関が存在する場合　→　系列相関がある　→　これが今回のお題になります。

残差の分布に偏りが見られる場合　→　データが正規分布していない　→　変数変換を検討する

分布に規則性が見られるケースについて、事例に即して詳しく解説されています。理解を深めるために、ぜひ読んでおくべきかと。

ググった程度では、即実行できるほどの知識は得られませんでしたね。
正直言って、「分布に規則性が見られる場合」と「データの並び方に意味がある場合」の違いすら分からないのですが、系列相関を疑うケースというのは時系列データに限りますので、その辺で区別するのかな、と。

系列相関は『マンガでわかる統計学　回帰分析編』にも説明があったのですが、私自身が時系列データを扱っていないので、スルーするつもりでいました。が、今回学習の流れで再びお目見えしましたので、じゃあいっちょ勉強のために、ダービン・ワトソン統計量を求める関数でも作ってみよう、と思い直した次第。

『マンガでわかる統計学　回帰分析編』９７ページによると、ダービン・ワトソン統計量は、(隣り合う残差を引いたもの)²を足したもの ÷ (個々の残差)²を足したもので求められる由。なるほど。残差を求めるには予測値ŷが必要、ŷを求めるには回帰係数と切片が必要ですね。

というわけで、以下のように実装してみました。

'ダービン・ワトソン統計量を求める
Function Durbin_Watson(既知のy, 既知のx)

    Dim c, n, u, Reg, Se, Su, i, j, haty
    
    If 既知のy.Rows.Count = 既知のx.Rows.Count And 既知のy.Columns.Count = 1 Then
        With Application.WorksheetFunction
            c = 既知のx.Columns.Count                   '説明変数の数
            n = 既知のy.Rows.Count                      '観測数
            Reg = .LinEst(既知のy, 既知のx, True, True) '回帰統計量を取得
            Se = .Index(Reg, 5, 2)                      '残差平方和
            ReDim u(1 To n)
            haty = 0                                    '予測値を初期化
            For j = 1 To c
                haty = haty + 既知のx(1, j) * .Index(Reg, 1, c - j + 1)
            Next
            haty = haty + .Index(Reg, 1, c + 1)         '予測値
            u(1) = 既知のy(1) - haty                    '残差
            For i = 2 To n
                haty = 0                                '予測値を初期化
                For j = 1 To c
                    haty = haty + 既知のx(i, j) * .Index(Reg, 1, c - j + 1)
                Next
                haty = haty + .Index(Reg, 1, c + 1)     '予測値
                u(i) = 既知のy(i) - haty                '残差
                '隣り合う残差の差の２乗を加算
                Su = Su + (u(i) - u(i - 1)) ^ 2
            Next
            
            Durbin_Watson = Su / Se
        End With
    Else
        Durbin_Watson = CVErr(xlErrNA)
    End If

End Function

回帰係数と切片をLinest関数を使って取得したのがミソでして、これにより線形単回帰だけでなく、重回帰モデルや非線形回帰モデルにも対応可能です。
注意すべきポイントとしては、Linest関数は説明変数と逆の並びで回帰係数を返すので、
.Index(Reg, 1, c - j + 1) として変数j と逆に［行番号］をデクリメントさせています。

２次多項式回帰モデルを例に使い方を説明すると、