RとPythonでの自己相関関数（ACF）の信頼区間の違い

はじめに
自己相関関数（ACF）と偏自己相関関数（PACF）
Rの信頼区間
Pythonの信頼区間
どちらを用いれば良いのか
参考文献

はじめに

自己相関関数（ACF）のコレログラムをRとPythonで作成すると、自己相関係数や偏自己相関係数がゼロである、という帰無仮説の95%信頼区間も合わせて出力されますが、その信頼区間の算出方法の違いについてのメモです。結論から言うと、Rはラグ $h$ の値に関わらず信頼区間が一定（直線）になり、Pythonではラグ $h$ が大きくなるにつれて信頼区間が大きく（曲線状に広がる）なります。という計算方法の違いです。これは自己相関係数の分散がラグ $h$ に依存するかどうかという前提条件の違いです。

以下のstackoverflowのページにRとPythonでの自己相関係数の分散の違いについて簡潔に書かれています
stackoverflow.com

自己相関関数（ACF）と偏自己相関関数（PACF）

自己相関関数（auto corelation function; ACF）は、自己相関係数をラグの関数としてみたもの。 $h$ 次の（ラグが $h$ の）自己相関は、 $y_t$ と $y_{t-h}$ の間の自己相関になる。偏自己相関関数（partial auto corelation function; PACF）は、同じく偏自己相関係数をラグの関数としてみたもの。 $h$ 次の偏自己相関は $y_t$ を定数と $y_{t-1}, \cdots, y_{t-h}$ に線形射影した時の射影係数で定義される、すなわち $y_t$ と $y_{t-h}$ から $y_{t-1}, \cdots, y_{t-h+1}$ の影響を取り除いたものとして定義される。自己相関関数（ACF）と偏自己相関関数（PACF）をグラフで描いたものをコレログラムと呼ぶ。

Rの信頼区間

それではAirPassengersのデータを使ってコレログラムを書いてみます。

# ライブラリのインストールと読み込み
install.packages("forecast")
library(forecast)

# データの読み込み（例：AirPassengersデータ）
data <- AirPassengers

#acf
mx=20
acf(AirPassengers, lag.max=mx, xaxt="n", xlab="Lag (months)")
axis(1, at=0:mx/12, labels=0:mx)

Rのacf 関数では95%信頼区間は $\pm \frac{1.96}{\sqrt{n}}$ で計算される。つまり信頼区間はサンプルサイズ $n$ のみに依存する。

Pythonの信頼区間

# ライブラリのインストールと読み込み
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# データの読み込み（例：AirPassengersデータ）
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
data = pd.read_csv(url, header=0, parse_dates=[0], index_col=0)

# acf
plot_acf(data,lags=20)
plt.title("ACF")
plt.xlabel("Lag")
plt.show()

Pythonのstatsmodelsライブラリのplot_acf()だと、ラグ $h$ に対する95%信頼区間は $\pm \frac{1.96}{\sqrt{n-h}}$ で計算される。つまりラグが増加するにつれて標準誤差が大きくなることを前提としている。plot_acf()は、デフォルトで

bartlett_confint=True

が設定されており、Bartlett's formulaを使用して信頼区間が計算されているらしい。冒頭のstackoverflowの回答によると、 $h$ 次の自己相関の標準誤差は $1, \cdots, h-1$ 次の自己相関の関数として表される。なのでこのオプションをFalseに設定すると、Rと同じ出力結果が得られる。

どちらを用いれば良いのか

そもそもACFのコレログラムを出力する目的の一つに、AR過程、MA過程、ARMA過程の次数を選択するというものがある。次数を多く取りすぎると推定するモデルが複雑になってしまうため、やはりラグが大きくなるにつれて信頼区間の幅が大きくなる（＝過去の値との関係が徐々になくなる）前提で、次数の大きさにペナルティをかける考え方が良いのではないかと思う。

次数選択と自己相関・偏自己相関の関係

自己相関係数	偏自己相関係数	選択モデル
2次以降0	ゆっくり減衰	MA(1)
3次以降0	ゆっくり減衰	MA(2)
ゆっくり減衰	2次以降0	AR(1)
ゆっくり減衰	3次以降0	AR(2)
ゆっくり減衰	ゆっくり減衰	ARMA(1,1)

また、pythonのstatsmodelの公式ドキュメントのbartlett_confint の箇所には以下のようにある。

（略）
The formula used for standard error depends upon the situation. If the autocorrelations are being used to test for randomness of residuals as part of the ARIMA routine, the standard errors are determined assuming the residuals are white noise. The approximate formula for any lag is that standard error of each r_k = 1/sqrt(N). See section 9.4 of [2] for more details on the 1/sqrt(N) result. For more elementary discussion, see section 5.3.2 in [3]. For the ACF of raw data, the standard error at a lag k is found as if the right model was an MA(k-1). This allows the possible interpretation that if all autocorrelations past a certain lag are within the limits, the model might be an MA of order defined by the last significant autocorrelation. In this case, a moving average model is assumed for the data and the standard errors for the confidence intervals should be generated using Bartlett's formula.
（略）
標準誤差に使用される公式は、状況によって異なる。自己相関がARIMAルーチンの一部として残差のランダム性を検定するために使用されている場合、標準誤差は残差がホワイトノイズであると仮定して決定されます。任意のラグの近似式は、各r_kの標準誤差 = 1/sqrt(N) です。1/sqrt(N)の結果の詳細については[2]の9.4節を参照。より初歩的な議論については、[3]の5.3.2節を参照のこと。生データのACFでは、ラグkでの標準誤差は、正しいモデルがMA(k-1)であるかのように求められます。これは、あるラグを過ぎたすべての自己相関が限界内にある場合、モデルは最後の有意な自己相関で定義される次数のMAであるかもしれないという解釈を可能にします。この場合、データについて移動平均モデルが仮定され、信頼区間の標準誤差がBartlettの公式を用いて生成されなければならない

www.statsmodels.org

これを踏まえると、一般的な時系列データだと誤差項がホワイトノイズであるとは限らない、なのでBartlettの公式を用いて、 $h$ 次の自己相関の標準誤差には $1, \cdots, h-1$ 次の自己相関が誤差項を経由して影響すべきだと仮定した方が自然では、と解釈した。ちなみに偏自己相関係数では他のラグの影響を受けないため、信頼区間は一定となる。