上級マクロ経済学

第12回 ARモデル

Author

荻巣嘉高

1 AR(1)モデル

AR(1)モデル

自己回帰モデル（ARモデル）は次のように定められる。 y_t = \phi_0 + \phi_1 y_{t-1} + u_t ここで、u_tはiidで、 E(u_t) = 0, E(u_t^2) = \sigma^2 このu_tを、ホワイトノイズと呼んだ。

t期の値y_tに対してt-1期のy_{t-1}をラグ変数と呼んだりする。
t期の値を1期のラグ変数で説明するモデルは、1次の自己回帰モデル（AR(1)）と呼ばれる。

定常性

変数y_tの定常性は、次の3つの性質を満たすこと。

平均（E(y_t)）が時間を通じて一定
分散（Var(y_t)）が時間を通じて一定
自己共分散（Cov(y_t,y_{t-j})）が時間を通じて一定：

y_t = \phi_0 + \phi_1 y_{t-1} + u_t

AR(1)モデルの定常性は、|\phi_1| < 1によって保証される。

AR(1)のダイナミクス

y_t = \phi_0 + \phi_1 y_{t-1} + u_t を。\phi_0=0のもとで、いくつかの\phi_1の値で比較してみよう。

Code

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1234) # ランダムシードを固定
u = np.random.normal(size=100) # 標準正規分布でホワイトノイズを発生させる

plt.rcParams["font.size"] = 16
fig = plt.figure() # 図を設定

phi0 = 0 # phi0を指定
phi1_list = [-0.5,0.5,0.9,1.1] # phi1のリストを作成

for j in range(4):
    phi1 = phi1_list[j] # phi1を指定
    res = [] # 結果の保存用リストを作成
    y0 = 0 # yの初期値を0に指定
    for i in range(100):
        y1 = phi0 + phi1*y0 + u[i] # 計算
        res.append(y1)
        y0 = y1 # y0をアップデート
    ax = fig.add_subplot(2,2,j+1)
    ax.plot(np.arange(100),res, label=f"phi1={phi1}")
    ax.legend()

AR(1)の推定

推定は、通常の最小二乗法（OLS）で行えば良い。

OLSで一致推計量を得られる条件: E(u_t | y_{t-1} ) = 0

t期の誤差項u_tとy_{t-1}が独立ならOK。
その条件は満たされる？

展開

AR(1)モデルでは、どの連続する2期をとってもパラメータは変化しないので、 \begin{aligned} y_{t-1} &= \phi_0 + \phi_1 y_{t-2} + u_{t-1} \\ &= \phi_0 + u_{t-1} + \phi_1 y_{t-2} \\ &= \phi_0 + \phi_1 (\phi_0 + \phi_1 y_{t-3} + u_{t-2}) + u_{t-1} \\ &= \underbrace{\phi_0 + \phi_1 \phi_0}_{定数} + \underbrace{\phi_1 u_{t-2} + u_{t-1}}_{ホワイトノイズ} + \underbrace{\phi_1^2 y_{t-3}}_{過去の変数} \\ & \hspace{24pt} \vdots \end{aligned}

AR(1)の無限加重和表現

これを続けていくと、次のように表現される。 y_{t-1} = \underbrace{\frac{\phi_0}{1-\phi_1}}_{定数} + \underbrace{\sum_{j=1}^{\infty} \phi_1^{j-1} u_{t-j}}_{過去のホワイトノイズの加重和}

y_{t-1}は、u_{t-1}より前の期のホワイトノイズの加重和に依存している。
u_tは、y_{t-1}を構築するu_{t-1}, u_{t-2}, \cdotsのいずれとも独立。
- u_{t}はiidだったため、独立性が仮定されている。

　AR(1)の実装

pythonでAR(1)を実装するためには、statsmodels.tsa.apiモジュールの中のAutoRegを用いる。

lags=1で、ラグの数を1に指定。
y_tをGDPギャップとしてAR(1)の推計。

import pandas as pd
import statsmodels.tsa.api as ts

df = pd.read_excel("gdp_gap.xlsx", sheet_name="Data")

df["time"] = pd.to_datetime(df["年"].astype(str) + "-" + df["月"].astype(str))

model = ts.AutoReg(df["GDP_gap"], lags=1)
res = model.fit(cov_type="HC0")
print(res.summary())

                            AutoReg Model Results                             
==============================================================================
Dep. Variable:                GDP_gap   No. Observations:                  178
Model:                     AutoReg(1)   Log Likelihood                -262.961
Method:               Conditional MLE   S.D. of innovations              1.069
Date:                Sat, 13 Sep 2025   AIC                            531.922
Time:                        13:05:53   BIC                            541.450
Sample:                             1   HQIC                           535.786
                                  178                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.1607      0.083     -1.939      0.053      -0.323       0.002
GDP_gap.L1     0.8140      0.056     14.515      0.000       0.704       0.924
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.2285           +0.0000j            1.2285            0.0000
-----------------------------------------------------------------------------

constが\phi_0（定数項）の推計値\hat{\phi}_0を表す。
GDP_gap.L1がGDPギャップの1次ラグy_{t-1}の係数\phi_1の推定値\hat{\phi}_1を表す。

モデルの結果の表示

AR(1)でGDPギャップの推計をすると、以下のように推計できた。 y_t = \underset{(0.083)}{-0.161} + \underset{(0.056)}{0.814} y_{t-1} + u_t

u_tを除く \hat{y}_t = \underset{(0.083)}{-0.161} + \underset{(0.056)}{0.814} y_{t-1} 部分が、モデルの予測値として解釈可能。

y_tのフィットの確認

モデルのフィッティングを確認してみる。

phi0 = res.params[0] # const
phi1 = res.params[1] # L1

## modelの値を計算
ar1 = phi0 + phi1*(df["GDP_gap"].iloc[:-1].values) # 1期ラグを用いるので、最後の値は使えない。

## プロットする
plt.plot(df["time"].iloc[1:],ar1, label="model")
plt.plot(df["time"], df["GDP_gap"], label="data", ls=":", color="gray")
plt.legend()
plt.show()

残差のプロット

u = res.resid # 残差の系列

## プロットする
plt.plot(df["time"].iloc[1:], u, label=r"$\hat{u}$")
plt.plot(df["time"].iloc[1:], np.zeros(len(u)), color="gray", alpha=0.3)
plt.legend()
plt.show()

残差はホワイトノイズか？

誤差項u_tはホワイトノイズを仮定されていた。
誤差項の推定値である残差\hat{u}_tもホワイトノイズであることが望まれる。
コレログラムで確認。
- 系列相関が出現

Code

rho = ts.acf(u, nlags = 20) # 20期までの自己相関を計算
SE = (1/len(u)) # 標準誤差の計算
## 区間の生成
top = np.ones(20)*SE*1.96
bottom = -np.ones(20)*SE*1.96

plt.plot(np.arange(1,21),rho[1:], label="correlation")
plt.fill_between(np.arange(1,21), y1=top, y2=bottom, color="gray", alpha=0.3, label=r"$\pm1.96\times$SE")
plt.legend()
plt.show()

誤差の系列相関とバイアス

系列相関が存在する時、ARモデルの推定係数\hat{\phi}_1の分散にはバイアスが発生することが知られている。
\widehat{SE}(\hat{\phi}_1) = \sqrt{Var(\hat{\phi}_1)}なので、標準誤差にもバイアスが生じる。
- \hat{\phi}_1の検定をする場合に問題になる。
誤差項に系列相関が存在する場合には、系列相関に頑健な標準誤差である、HAC標準誤差を利用する。
- Heteroskedasticity- and Autocorrelation-Consistentの略。
  - Heteroskedasticity：不均一分散
  - Autocorrelation：系列相関
- ちなみにpythonでは、これも簡単に計算できる。

HAC標準誤差

HAC標準誤差（Newey-Westの標準誤差）: \widehat{SE}_{HAC} (\hat{\phi}_1) = \frac{1}{\hat{\sigma}_{y}^2} \sqrt{\frac{1}{T} \left[ \hat{\sigma}^2 + 2 \sum_{j = 1}^{m-1} \left( 1 - \frac{j}{m} \right) \hat{\gamma}_j \right]}

ここで、

\hat{\sigma}_y^2はyの分散の推定値
\hat{\sigma}^2と\hat{\gamma}_jは次の\hat{v}_tの分散と自己共分散の推定値 \hat{v}_t = (y_t - \bar{y})\hat{u}_t
mは適切に選択されたラグ数

HAC標準誤差のラグ選択

mは大きすぎないことが（統計的な性質として）重要とされる。しかし少なすぎるとそもそも系列相関に頑健にならない。
mは経験的に、以下のように選ばれることが多い。
- 四半期データならばm=4
- 月次データならばm=12
- あるいは、

m = \left\lfloor 4 \left(\frac{T}{100} \right)^{\frac{1}{3}} \right\rfloor

統計的な裏付けに基けば、\lim_{T \to \infty} m/T = 0かつ\lim_{T \to \infty} m = \inftyとなって欲しい。
保守的に推定・検定を行なっていくならば、時系列データの解析ではHAC標準誤差を最初から使うのが望ましい。

HAC標準誤差の実装

fit("HAC", cov_kwds={"maxlags":m})を指定。
GDPギャップは四半期データなので、m=4で実装。

model = ts.AutoReg(df["GDP_gap"], lags=1)
res_hc0 = model.fit(cov_type="HC0") # 比較用：不均一分散に頑健なSEを用いた演算
res_hac = model.fit(cov_type="HAC", cov_kwds={"maxlags":4}) # HACを使ってSEを計算（ラグm=4）

## modelの値を計算
ar1 = phi0 + phi1*(df["GDP_gap"].iloc[:-1].values) # 1期ラグを用いるので、最後の値は使えない。

SE_hc0 = res_hc0.bse
SE_hac = res_hac.bse
print("HAC標準誤差")
print(SE_hac)
print("------------")
print(f"不均一分散に頑健な標準誤差")
print(SE_hc0)

HAC標準誤差
const         0.105690
GDP_gap.L1    0.056558
dtype: float64
------------
不均一分散に頑健な標準誤差
const         0.082885
GDP_gap.L1    0.056083
dtype: float64

表示の比較

不均一分散に頑健な標準誤差

\hat{y}_t = \underset{(0.083)}{-0.161} + \underset{(0.056)}{0.814} y_{t-1}

HAC標準誤差 \hat{y}_t = \underset{(0.011)}{-0.161} + \underset{(0.057)}{0.814} y_{t-1}
推定係数の値自体は変わっていないことに注意。
統計検定する際に使う\hat{SE}の値が変わるため、t統計量が変わる。

2 AR(p)モデル

AR(p)モデル

ARモデルは、ラグの次数を増やすこともできる。

ラグ次数をp次まで取ったものを、AR(p)モデルと呼ぶ。

y_t = \phi_0 + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + u_t

ただし、u_tはホワイトノイズで、y_tは定常であるとする。

AR(p)モデルの推計

ラグの最大次数pはどうやって決める？

情報量基準と呼ばれる基準を用いて決めることが一般的。
- 赤池情報量基準
- ベイズ情報量基準

直観的には: 情報量 = (a)フィッティング誤差 + (b)ラグペナルティで表し、情報量を最小にするようなラグ次数を選択する。

フィッティング誤算の評価

ラグ次数をpに設定するとする。
データの流列(y_1, y_2, \cdots, y_T)
モデルの予測値の流列(\hat{y}_{1+p}, \hat{y}_{2+p}, \cdots, \hat{y}_T)

残差: \hat{u}_{t} = y_t - \hat{y}_t

残差平方和（Residal sum of squeres; RSS）を次のように計算して、全体のフィッティングの誤差の大きさを計算 RSS(p) = \sum_{t = 1+p}^{T} \hat{u}_{t}^2

赤池情報量基準（AIC）

AIC(p) = \log \left( \frac{RSS(p)}{T-p} \right) + (p+1)\frac{2}{T-p}

分母のT-pはTになっているケースや、他に自由度の調整を行なったケースもある。

ベイズ情報量基準（BIC）

BIC(p) = \log \left( \frac{RSS(p)}{T-p} \right) + (p+1)\frac{\log(T-p)}{T-p}

分母のT-pはTになっているケースや、他に自由度の調整を行なったケースもある。
統計ソフトによって定義されている算出法に若干のばらつきは存在。
- pythonのtsaモジュールで出すものは、少し定義が違うので、自分で計算する必要がある。

AIC？BIC？

AICとBICで選択されるラグ次数は異なりうる。
どちらを使っても基本的にはよい。
- 目安はサンプルサイズが大きい場合（1000以上とか）はBICの方が望ましいというイメージ。
- サンプルが小さければAICでOK。
重要なのは、「恣意的にラグを選んだわけではない」と主張できること。

ラグ選択

GDPギャップデータでは：
- AICなら3
- BICなら1

Code

def aic(resid,T,p):
    """
    赤池情報量の算出
    -------------
    input:
        resid (ndarray): モデルの残差の列
        T (int): サンプルサイズ
        p (int): ラグの数
    output : 
        (int) AICの値
    """
    return np.log(np.sum(resid**2)/(T-p)) + (p+1)*2/(T-p)

def bic(resid,T,p):
    """
    ベイズ情報量の算出
    -------------
    input:
        resid (ndarray): モデルの残差の列
        T (int): サンプルサイズ
        p (int): ラグの数
    output : 
        (int) BICの値
    """
    return np.log(np.sum(resid**2)/(T-p)) + (p+1)*np.log(T-p)/(T-p)

AIC = []
BIC = []

T = len(df)
maxlag = 20
for p in range(1,maxlag):
    res = ts.AutoReg(df["GDP_gap"],lags=p).fit(cov_type="HC0")
    AIC.append(aic(res.resid.values,T,p))
    BIC.append(bic(res.resid.values,T,p))
    
plt.plot(np.arange(1,maxlag), AIC, label="AIC", marker="o")
plt.plot(np.arange(1,maxlag), BIC, label="BIC", marker="*")
plt.legend()
plt.show()

AR(3)とAR(1)の比較

Code

## AR(1)
model1 = ts.AutoReg(df["GDP_gap"], lags=1)
res1 = model1.fit(cov_type="HAC", cov_kwds={"maxlags":4})
phi0_1 = res1.params[0] # AR(1)のphi0
phi1_1 = res1.params[1] # AR(1)のphi1
y1 = phi0_1 + phi1_1 * df["GDP_gap"].values[:-1] # AR(1)の理論値

## AR(3)
model3 = ts.AutoReg(df["GDP_gap"], lags=3)
res3 = model3.fit(cov_type="HAC", cov_kwds={"maxlags":4})
phi0_3 = res3.params[0] # AR(3)のphi0
phi1_3 = res3.params[1] # AR(3)のphi1
phi2_3 = res3.params[2] # AR(3)のphi2
phi3_3 = res3.params[3] # AR(3)のphi3
y3 = phi0_3 + phi1_3 * df["GDP_gap"].values[2:-1] + phi2_3 * df["GDP_gap"].values[1:-2] + phi3_3 * df["GDP_gap"].values[:-3] #  AR(1)の理論値

## プロットする
plt.plot(df["time"].iloc[1:],y1, label="AR(1)", alpha=0.8)
plt.plot(df["time"].iloc[3:],y3, label="AR(3)", alpha=0.8)
plt.plot(df["time"],df["GDP_gap"], ls = ":", color="gray", alpha=0.5, label="data")
plt.legend()
plt.show()

まとめ

ARモデルを概説した。
- 推計は一般的なOLSでOK。
HAC標準誤差を紹介した。
- 誤差項に系列相関があると考えられる場合には、これを用いた方がよい。
ARモデルはラグを増やすこともできる。
- 最適なラグ数は情報量基準で選択。
- 客観性を担保することが重要。