カバン検定を自分で実装してみる

前の記事で、statsmodels に実装されている acorr_ljungbox 関数を使って、カバン検定を行ってみました。
statsmodelsでかばん検定 (自己相関の検定)

今回は学習のため、numpyで実装しました。
実装にあったって数式から説明します。
検定の帰無仮説や対立仮説については前の記事を参照してください。

$\{y_t\}_{0}^{T-1}$を時系列データとします。 (沖本本の定義とインデックスが一つずれているので注意。)

標本平均$\bar{y}$, 標本自己共分散$\hat{\gamma}_k$, 標本自己相関係数$\hat{\rho}_k$を次のように定義します。
\begin{eqnarray}
\bar{y} & = & \frac{1}{T}\sum_{t=0}^{T-1}y_t\\
\hat{\gamma}_k & = & \frac{1}{T}\sum_{t=k}^{T-1}(y_t-\bar{y})(y_{t-k}-\bar{y}) , \hspace{1em} k = 0, 1, 2, \cdots\\
\hat{\rho}_k & = & \frac{\hat{\gamma}_k}{\hat{\gamma}_0},\hspace{1em} k = 1, 2, 3, \cdots
\end{eqnarray}

この時、Ljung と Box が考案した次の統計量$Q(m)$は漸近的に自由度$m$のカイ２乗分布に従います。
$$Q(m) = T(T+2)\sum_{k=1}^{m}\frac{\hat{\rho}_k^2}{T-k}$$

今回も7点周期のデータを準備します。
m = 7 で検定を行いますので、帰無仮説が棄却されれば、
LAG1〜7のいずれかの自己相関係数が0でないことが主張されます。
(p値は5%を基準に判断します。)

import numpy as np
import pandas as pd
from scipy.stats import chi2
# 答え合わせ用
from statsmodels.stats.diagnostic import acorr_ljungbox

# 7点ごとに周期性のあるデータを準備
series = pd.Series([1, 1, 1, 1, 1, 1, 5]*10)
# 乱数追加
series += np.random.randn(70)

# データの個数
T = len(series)
# 検定対象のm
m = 7
# 標本平均
y_ = series.mean()
# 標本自己共分散の列
gamma_list = np.array([
                np.sum([
                    (series[t]-y_)*(series[t-k]-y_) for t in range(k, T)
                ])/T for k in range(m+1)
            ])
# 標本自己相関係数の列
ro_list = gamma_list/gamma_list[0]

# Q(m)の計算
Q = T*(T+2)*np.sum([ro_list[k]**2 / (T-k) for k in range(1, m+1)])
print("lb値:", Q)
print("p値:", chi2.sf(Q, m))

# ライブラリを使った計算結果
lbvalues, pvalues = acorr_ljungbox(series, lags=7)
print(lbvalues[6])
print(pvalues[6])

# 出力結果
lb値: 43.78077348272421
p値: 2.3564368405873674e-07
43.780773482724186
2.3564368405873896e-07

p値が非常に小さな値であることがわかりました。
また、ライブラリを使った結果と四捨五入による誤差がありますが、それ以外は同じ値なのでミスも無さそうです。

この計算の注意点としては、
標本自己相関係数を使うところでしょうか。
pandasのautocorrは自己相関係数を出力してくれますが、
これが標本自己相関係数とは微妙に定義が異なり、autocorrを使うと少し値がずれてしまいます。