自己回帰モデル¶
このノートブックでは、
AutoReg モデルを使用した自己回帰モデリングを紹介します。また、ar_select_order がAICなどの情報基準を最小化するモデルの選択をサポートする方法についても説明します。自己回帰モデルは以下のダイナミクスで表されます：
\[y_t = \delta + \phi_1 y_{t-1} + \ldots + \phi_p y_{t-p} + \epsilon_t.\]
AutoReg では次のようなモデルも許可されます：
決定論的項（
trend）n: 決定論的項なしc: 定数（デフォルト）ct: 定数と時間トレンドt: 時間トレンドのみ
季節ダミー（
seasonal）Trueは、時間系列の期間 \(s\) に対して \(s-1\) のダミー変数を含みます（例えば、月次の場合は12）
カスタム決定論的項（
deterministic）DeterministicProcessを受け入れます
外生変数（
exog）モデルに含める外生変数の
DataFrameまたはarray
選択された遅延の省略（
lags）lagsが整数の反復可能なオブジェクトの場合、これらのみがモデルに含まれます。
完全な仕様は次の通りです：
\[y_t = \delta_0 + \delta_1 t + \phi_1 y_{t-1} + \ldots + \phi_p y_{t-p} + \sum_{i=1}^{s-1} \gamma_i d_i + \sum_{j=1}^{m} \kappa_j x_{t,j} + \epsilon_t.\]
ここで：
\(d_i\) は季節ダミーで、\(mod(t, period) = i\) の場合に1となります。モデルに定数（
cがtrendに含まれている場合）が含まれている場合、期間0は除外されます。\(t\) は時間トレンド（\(1, 2, \ldots\)）で、最初の観測に1から始まります。
\(x_{t,j}\) は外生回帰変数です。注：これらはモデルを定義する際に左辺変数と時間的に揃えられます。
\(\epsilon_t\) はホワイトノイズ過程と仮定されます。
この最初のセルは標準パッケージをインポートし、プロットをインラインで表示するように設定します。
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import pandas_datareader as pdr
import seaborn as sns
from statsmodels.tsa.api import acf, graphics, pacf
from statsmodels.tsa.ar_model import AutoReg, ar_select_order
このセルは、プロット スタイルを設定し、matplotlib 用の pandas 日付コンバーターを登録し、デフォルトの図のサイズを設定します。
[2]:
sns.set_style("darkgrid")
pd.plotting.register_matplotlib_converters()
# Default figure size
sns.mpl.rc("figure", figsize=(16, 6))
sns.mpl.rc("font", size=14)
最初の例では、季節調整されていない米国の住宅着工件数の月次成長率を使用しています。季節性は、ピークとトラフの定期的なパターンによって明らかです。AutoRegを使用する際に警告を避けるため、時系列の頻度を「MS」（月初）に設定しています。
[3]:
data = pdr.get_data_fred("HOUSTNSA", "1959-01-01", "2019-06-01")
housing = data.HOUSTNSA.pct_change().dropna()
# Scale by 100 to get percentages
housing = 100 * housing.asfreq("MS")
fig, ax = plt.subplots()
ax = housing.plot(ax=ax)
まず、AR(3) から始めましょう。これはこのデータに適したモデルではありませんが、API の基本的な使用方法を示しています。
[4]:
mod = AutoReg(housing, 3, old_names=False)
res = mod.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: HOUSTNSA No. Observations: 725
Model: AutoReg(3) Log Likelihood -2993.442
Method: Conditional MLE S.D. of innovations 15.289
Date: Thu, 28 Nov 2024 AIC 5996.884
Time: 23:10:29 BIC 6019.794
Sample: 05-01-1959 HQIC 6005.727
- 06-01-2019
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.1228 0.573 1.961 0.050 0.000 2.245
HOUSTNSA.L1 0.1910 0.036 5.235 0.000 0.120 0.263
HOUSTNSA.L2 0.0058 0.037 0.155 0.877 -0.067 0.079
HOUSTNSA.L3 -0.1939 0.036 -5.319 0.000 -0.265 -0.122
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 0.9680 -1.3298j 1.6448 -0.1499
AR.2 0.9680 +1.3298j 1.6448 0.1499
AR.3 -1.9064 -0.0000j 1.9064 -0.5000
-----------------------------------------------------------------------------
AutoReg は OLS と同じ共分散推定器をサポートしています。以下では、White の共分散推定器である cov_type="HC0" を使用しています。パラメータの推定値は同じですが、標準誤差に依存するすべての量が変化します。
[5]:
res = mod.fit(cov_type="HC0")
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: HOUSTNSA No. Observations: 725
Model: AutoReg(3) Log Likelihood -2993.442
Method: Conditional MLE S.D. of innovations 15.289
Date: Thu, 28 Nov 2024 AIC 5996.884
Time: 23:10:29 BIC 6019.794
Sample: 05-01-1959 HQIC 6005.727
- 06-01-2019
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.1228 0.601 1.869 0.062 -0.055 2.300
HOUSTNSA.L1 0.1910 0.035 5.499 0.000 0.123 0.259
HOUSTNSA.L2 0.0058 0.039 0.150 0.881 -0.070 0.081
HOUSTNSA.L3 -0.1939 0.036 -5.448 0.000 -0.264 -0.124
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 0.9680 -1.3298j 1.6448 -0.1499
AR.2 0.9680 +1.3298j 1.6448 0.1499
AR.3 -1.9064 -0.0000j 1.9064 -0.5000
-----------------------------------------------------------------------------
[6]:
sel = ar_select_order(housing, 13, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: HOUSTNSA No. Observations: 725
Model: AutoReg(13) Log Likelihood -2676.157
Method: Conditional MLE S.D. of innovations 10.378
Date: Thu, 28 Nov 2024 AIC 5382.314
Time: 23:10:29 BIC 5450.835
Sample: 03-01-1960 HQIC 5408.781
- 06-01-2019
================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
const 1.3615 0.458 2.970 0.003 0.463 2.260
HOUSTNSA.L1 -0.2900 0.036 -8.161 0.000 -0.360 -0.220
HOUSTNSA.L2 -0.0828 0.031 -2.652 0.008 -0.144 -0.022
HOUSTNSA.L3 -0.0654 0.031 -2.106 0.035 -0.126 -0.005
HOUSTNSA.L4 -0.1596 0.031 -5.166 0.000 -0.220 -0.099
HOUSTNSA.L5 -0.0434 0.031 -1.387 0.165 -0.105 0.018
HOUSTNSA.L6 -0.0884 0.031 -2.867 0.004 -0.149 -0.028
HOUSTNSA.L7 -0.0556 0.031 -1.797 0.072 -0.116 0.005
HOUSTNSA.L8 -0.1482 0.031 -4.803 0.000 -0.209 -0.088
HOUSTNSA.L9 -0.0926 0.031 -2.960 0.003 -0.154 -0.031
HOUSTNSA.L10 -0.1133 0.031 -3.665 0.000 -0.174 -0.053
HOUSTNSA.L11 0.1151 0.031 3.699 0.000 0.054 0.176
HOUSTNSA.L12 0.5352 0.031 17.133 0.000 0.474 0.596
HOUSTNSA.L13 0.3178 0.036 8.937 0.000 0.248 0.388
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 1.0913 -0.0000j 1.0913 -0.0000
AR.2 0.8743 -0.5018j 1.0080 -0.0829
AR.3 0.8743 +0.5018j 1.0080 0.0829
AR.4 0.5041 -0.8765j 1.0111 -0.1669
AR.5 0.5041 +0.8765j 1.0111 0.1669
AR.6 0.0056 -1.0530j 1.0530 -0.2491
AR.7 0.0056 +1.0530j 1.0530 0.2491
AR.8 -0.5263 -0.9335j 1.0716 -0.3317
AR.9 -0.5263 +0.9335j 1.0716 0.3317
AR.10 -0.9525 -0.5880j 1.1194 -0.4120
AR.11 -0.9525 +0.5880j 1.1194 0.4120
AR.12 -1.2928 -0.2608j 1.3189 -0.4683
AR.13 -1.2928 +0.2608j 1.3189 0.4683
------------------------------------------------------------------------------
plot_predict は予測を可視化します。ここでは、モデルによって捉えられた強い季節性を示す多数の予測を生成します。
[7]:
fig = res.plot_predict(720, 840)
plot_diagnositcsは、モデルがデータの主要な特徴を捉えていることを示しています。
[8]:
fig = plt.figure(figsize=(16, 9))
fig = res.plot_diagnostics(fig=fig, lags=30)
季節ダミー¶
AutoReg は季節性ダミーをサポートしており、これは季節性をモデル化するための別の方法です。ダミーを含めることで、動的性質が AR(2) のみに短縮されます。
[9]:
sel = ar_select_order(housing, 13, seasonal=True, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: HOUSTNSA No. Observations: 725
Model: Seas. AutoReg(2) Log Likelihood -2652.556
Method: Conditional MLE S.D. of innovations 9.487
Date: Thu, 28 Nov 2024 AIC 5335.112
Time: 23:10:31 BIC 5403.863
Sample: 04-01-1959 HQIC 5361.648
- 06-01-2019
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.2726 1.373 0.927 0.354 -1.418 3.963
s(2,12) 32.6477 1.824 17.901 0.000 29.073 36.222
s(3,12) 23.0685 2.435 9.472 0.000 18.295 27.842
s(4,12) 10.7267 2.693 3.983 0.000 5.449 16.005
s(5,12) 1.6792 2.100 0.799 0.424 -2.437 5.796
s(6,12) -4.4229 1.896 -2.333 0.020 -8.138 -0.707
s(7,12) -4.2113 1.824 -2.309 0.021 -7.786 -0.636
s(8,12) -6.4124 1.791 -3.581 0.000 -9.922 -2.902
s(9,12) 0.1095 1.800 0.061 0.952 -3.419 3.638
s(10,12) -16.7511 1.814 -9.234 0.000 -20.307 -13.196
s(11,12) -20.7023 1.862 -11.117 0.000 -24.352 -17.053
s(12,12) -11.9554 1.778 -6.724 0.000 -15.440 -8.470
HOUSTNSA.L1 -0.2953 0.037 -7.994 0.000 -0.368 -0.223
HOUSTNSA.L2 -0.1148 0.037 -3.107 0.002 -0.187 -0.042
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 -1.2862 -2.6564j 2.9514 -0.3218
AR.2 -1.2862 +2.6564j 2.9514 0.3218
-----------------------------------------------------------------------------
季節ダミーは、今後 10 年間のすべての期間において重要な季節要素を持つ予測において明らかです。
[10]:
fig = res.plot_predict(720, 840)
[11]:
fig = plt.figure(figsize=(16, 9))
fig = res.plot_diagnostics(lags=30, fig=fig)
季節ダイナミクス¶
AutoReg はパラメータ推定に OLS を使用するため、季節成分を直接サポートしていませんが、季節ARにおける制約を課さない過剰にパラメータ化された季節 AR を使用することで、季節的な動態を捉えることは可能です。
[12]:
yoy_housing = data.HOUSTNSA.pct_change(12).resample("MS").last().dropna()
_, ax = plt.subplots()
ax = yoy_housing.plot(ax=ax)
まず、最大のラグのみを選択する簡単な方法を使用してモデルを選択します。それより小さいラグは自動的に含まれます。確認する最大のラグは13に設定されます。これは、モデルが短期的な AR(1) 成分と季節的な AR(1) 成分の両方を持つ季節的ARを次に選択するためです。したがって、
\[(1-\phi_s L^{12})(1-\phi_1 L)y_t = \epsilon_t\]
これは展開すると、
\[y_t = \phi_1 y_{t-1} +\phi_s Y_{t-12} - \phi_1\phi_s Y_{t-13} + \epsilon_t\]
となります。
AutoReg はこの構造を強制することはありませんが、ネストされたモデルを推定することができます。
\[y_t = \phi_1 y_{t-1} +\phi_{12} Y_{t-12} - \phi_{13} Y_{t-13} + \epsilon_t\]
ここで、すべての13のラグが選ばれたことがわかります。
[13]:
sel = ar_select_order(yoy_housing, 13, old_names=False)
sel.ar_lags
[13]:
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
すべての13のラグが必要であるとは考えにくいです。glob=Trueを設定することで、最大13のラグを含む \(2^{13}\) のモデルをすべて探索することができます。
ここでは、最初の三つと 7 番目、最後に 12 番目と 13 番目のラグが選択されているのがわかります。これは、上記で説明した構造に表面的に似ています。
モデルを適合させた後、診断プロットを確認しますが、この仕様がデータの動態を捉えるのに適切であることを示しています。
[14]:
sel = ar_select_order(yoy_housing, 13, glob=True, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: HOUSTNSA No. Observations: 714
Model: Restr. AutoReg(13) Log Likelihood 589.177
Method: Conditional MLE S.D. of innovations 0.104
Date: Thu, 28 Nov 2024 AIC -1162.353
Time: 23:10:39 BIC -1125.933
Sample: 02-01-1961 HQIC -1148.276
- 06-01-2019
================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
const 0.0035 0.004 0.875 0.382 -0.004 0.011
HOUSTNSA.L1 0.5640 0.035 16.167 0.000 0.496 0.632
HOUSTNSA.L2 0.2347 0.038 6.238 0.000 0.161 0.308
HOUSTNSA.L3 0.2051 0.037 5.560 0.000 0.133 0.277
HOUSTNSA.L7 -0.0903 0.030 -2.976 0.003 -0.150 -0.031
HOUSTNSA.L12 -0.3791 0.034 -11.075 0.000 -0.446 -0.312
HOUSTNSA.L13 0.3354 0.033 10.254 0.000 0.271 0.400
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 -1.0309 -0.2682j 1.0652 -0.4595
AR.2 -1.0309 +0.2682j 1.0652 0.4595
AR.3 -0.7454 -0.7417j 1.0515 -0.3754
AR.4 -0.7454 +0.7417j 1.0515 0.3754
AR.5 -0.3172 -1.0221j 1.0702 -0.2979
AR.6 -0.3172 +1.0221j 1.0702 0.2979
AR.7 0.2419 -1.0573j 1.0846 -0.2142
AR.8 0.2419 +1.0573j 1.0846 0.2142
AR.9 0.7840 -0.8303j 1.1420 -0.1296
AR.10 0.7840 +0.8303j 1.1420 0.1296
AR.11 1.0730 -0.2386j 1.0992 -0.0348
AR.12 1.0730 +0.2386j 1.0992 0.0348
AR.13 1.1193 -0.0000j 1.1193 -0.0000
------------------------------------------------------------------------------
[15]:
fig = plt.figure(figsize=(16, 9))
fig = res.plot_diagnostics(fig=fig, lags=30)
季節のダミーを含めることもできます。モデルは前年比の変化を使用しているため、これらはすべて重要ではありません。
[16]:
sel = ar_select_order(yoy_housing, 13, glob=True, seasonal=True, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())
AutoReg Model Results
====================================================================================
Dep. Variable: HOUSTNSA No. Observations: 714
Model: Restr. Seas. AutoReg(13) Log Likelihood 590.875
Method: Conditional MLE S.D. of innovations 0.104
Date: Thu, 28 Nov 2024 AIC -1143.751
Time: 23:10:54 BIC -1057.253
Sample: 02-01-1961 HQIC -1110.317
- 06-01-2019
================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
const 0.0167 0.014 1.215 0.224 -0.010 0.044
s(2,12) -0.0179 0.019 -0.931 0.352 -0.056 0.020
s(3,12) -0.0121 0.019 -0.630 0.528 -0.050 0.026
s(4,12) -0.0210 0.019 -1.089 0.276 -0.059 0.017
s(5,12) -0.0223 0.019 -1.157 0.247 -0.060 0.015
s(6,12) -0.0224 0.019 -1.160 0.246 -0.060 0.015
s(7,12) -0.0212 0.019 -1.096 0.273 -0.059 0.017
s(8,12) -0.0101 0.019 -0.520 0.603 -0.048 0.028
s(9,12) -0.0095 0.019 -0.491 0.623 -0.047 0.028
s(10,12) -0.0049 0.019 -0.252 0.801 -0.043 0.033
s(11,12) -0.0084 0.019 -0.435 0.664 -0.046 0.030
s(12,12) -0.0077 0.019 -0.400 0.689 -0.046 0.030
HOUSTNSA.L1 0.5630 0.035 16.160 0.000 0.495 0.631
HOUSTNSA.L2 0.2347 0.038 6.248 0.000 0.161 0.308
HOUSTNSA.L3 0.2075 0.037 5.634 0.000 0.135 0.280
HOUSTNSA.L7 -0.0916 0.030 -3.013 0.003 -0.151 -0.032
HOUSTNSA.L12 -0.3810 0.034 -11.149 0.000 -0.448 -0.314
HOUSTNSA.L13 0.3373 0.033 10.327 0.000 0.273 0.401
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 -1.0305 -0.2681j 1.0648 -0.4595
AR.2 -1.0305 +0.2681j 1.0648 0.4595
AR.3 -0.7447 -0.7414j 1.0509 -0.3754
AR.4 -0.7447 +0.7414j 1.0509 0.3754
AR.5 -0.3172 -1.0215j 1.0696 -0.2979
AR.6 -0.3172 +1.0215j 1.0696 0.2979
AR.7 0.2416 -1.0568j 1.0841 -0.2142
AR.8 0.2416 +1.0568j 1.0841 0.2142
AR.9 0.7837 -0.8304j 1.1418 -0.1296
AR.10 0.7837 +0.8304j 1.1418 0.1296
AR.11 1.0724 -0.2383j 1.0986 -0.0348
AR.12 1.0724 +0.2383j 1.0986 0.0348
AR.13 1.1192 -0.0000j 1.1192 -0.0000
------------------------------------------------------------------------------
工業生産¶
鉱工業生産指数データを用いて予測を検討します。
[17]:
data = pdr.get_data_fred("INDPRO", "1959-01-01", "2019-06-01")
ind_prod = data.INDPRO.pct_change(12).dropna().asfreq("MS")
_, ax = plt.subplots(figsize=(16, 9))
ind_prod.plot(ax=ax)
[17]:
<Axes: xlabel='DATE'>
まず、最大 12 のラグを使用するモデルを選択します。 AR(13) は、多くの係数が重要でないにもかかわらず、BIC 基準を最小化します。
[18]:
sel = ar_select_order(ind_prod, 13, "bic", old_names=False)
res = sel.model.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: INDPRO No. Observations: 714
Model: AutoReg(13) Log Likelihood 2321.382
Method: Conditional MLE S.D. of innovations 0.009
Date: Thu, 28 Nov 2024 AIC -4612.763
Time: 23:10:55 BIC -4544.476
Sample: 02-01-1961 HQIC -4586.368
- 06-01-2019
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0012 0.000 2.783 0.005 0.000 0.002
INDPRO.L1 1.1577 0.035 33.179 0.000 1.089 1.226
INDPRO.L2 -0.0822 0.053 -1.543 0.123 -0.187 0.022
INDPRO.L3 7.041e-05 0.053 0.001 0.999 -0.104 0.104
INDPRO.L4 0.0076 0.053 0.143 0.886 -0.096 0.111
INDPRO.L5 -0.1329 0.053 -2.529 0.011 -0.236 -0.030
INDPRO.L6 -0.0080 0.052 -0.153 0.879 -0.111 0.095
INDPRO.L7 0.0556 0.052 1.064 0.287 -0.047 0.158
INDPRO.L8 -0.0296 0.052 -0.568 0.570 -0.132 0.073
INDPRO.L9 0.0920 0.052 1.770 0.077 -0.010 0.194
INDPRO.L10 -0.0808 0.052 -1.552 0.121 -0.183 0.021
INDPRO.L11 -0.0003 0.052 -0.005 0.996 -0.102 0.102
INDPRO.L12 -0.3820 0.052 -7.364 0.000 -0.484 -0.280
INDPRO.L13 0.3613 0.033 10.996 0.000 0.297 0.426
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 -1.0408 -0.2913j 1.0808 -0.4566
AR.2 -1.0408 +0.2913j 1.0808 0.4566
AR.3 -0.7803 -0.8039j 1.1204 -0.3726
AR.4 -0.7803 +0.8039j 1.1204 0.3726
AR.5 -0.2728 -1.0533j 1.0880 -0.2903
AR.6 -0.2728 +1.0533j 1.0880 0.2903
AR.7 0.2716 -1.0507j 1.0852 -0.2097
AR.8 0.2716 +1.0507j 1.0852 0.2097
AR.9 0.8010 -0.7285j 1.0827 -0.1175
AR.10 0.8010 +0.7285j 1.0827 0.1175
AR.11 1.0219 -0.2219j 1.0457 -0.0340
AR.12 1.0219 +0.2219j 1.0457 0.0340
AR.13 1.0560 -0.0000j 1.0560 -0.0000
------------------------------------------------------------------------------
必要に応じて、短いラグを必要とせずに長いラグを入力できるグローバル検索を使用することもできます。ここでは、多くのラグが削除されていることがわかります。モデルは、データに季節性がある可能性があることを示しています。
[19]:
sel = ar_select_order(ind_prod, 13, "bic", glob=True, old_names=False)
sel.ar_lags
res_glob = sel.model.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: INDPRO No. Observations: 714
Model: AutoReg(13) Log Likelihood 2321.382
Method: Conditional MLE S.D. of innovations 0.009
Date: Thu, 28 Nov 2024 AIC -4612.763
Time: 23:11:01 BIC -4544.476
Sample: 02-01-1961 HQIC -4586.368
- 06-01-2019
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0012 0.000 2.783 0.005 0.000 0.002
INDPRO.L1 1.1577 0.035 33.179 0.000 1.089 1.226
INDPRO.L2 -0.0822 0.053 -1.543 0.123 -0.187 0.022
INDPRO.L3 7.041e-05 0.053 0.001 0.999 -0.104 0.104
INDPRO.L4 0.0076 0.053 0.143 0.886 -0.096 0.111
INDPRO.L5 -0.1329 0.053 -2.529 0.011 -0.236 -0.030
INDPRO.L6 -0.0080 0.052 -0.153 0.879 -0.111 0.095
INDPRO.L7 0.0556 0.052 1.064 0.287 -0.047 0.158
INDPRO.L8 -0.0296 0.052 -0.568 0.570 -0.132 0.073
INDPRO.L9 0.0920 0.052 1.770 0.077 -0.010 0.194
INDPRO.L10 -0.0808 0.052 -1.552 0.121 -0.183 0.021
INDPRO.L11 -0.0003 0.052 -0.005 0.996 -0.102 0.102
INDPRO.L12 -0.3820 0.052 -7.364 0.000 -0.484 -0.280
INDPRO.L13 0.3613 0.033 10.996 0.000 0.297 0.426
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 -1.0408 -0.2913j 1.0808 -0.4566
AR.2 -1.0408 +0.2913j 1.0808 0.4566
AR.3 -0.7803 -0.8039j 1.1204 -0.3726
AR.4 -0.7803 +0.8039j 1.1204 0.3726
AR.5 -0.2728 -1.0533j 1.0880 -0.2903
AR.6 -0.2728 +1.0533j 1.0880 0.2903
AR.7 0.2716 -1.0507j 1.0852 -0.2097
AR.8 0.2716 +1.0507j 1.0852 0.2097
AR.9 0.8010 -0.7285j 1.0827 -0.1175
AR.10 0.8010 +0.7285j 1.0827 0.1175
AR.11 1.0219 -0.2219j 1.0457 -0.0340
AR.12 1.0219 +0.2219j 1.0457 0.0340
AR.13 1.0560 -0.0000j 1.0560 -0.0000
------------------------------------------------------------------------------
plot_predict は、予測値とともに信頼区間を表示するために使用できます。ここでは、最後の観測値から始めて 18 か月間の予測を作成します。
[20]:
ind_prod.shape
[20]:
(714,)
[21]:
fig = res_glob.plot_predict(start=714, end=732)
フルモデルと制限モデルからの予測は非常に似ています。また、非常に異なるダイナミクスを持つ AR(5) も含まれています
[22]:
res_ar5 = AutoReg(ind_prod, 5, old_names=False).fit()
predictions = pd.DataFrame(
{
"AR(5)": res_ar5.predict(start=714, end=726),
"AR(13)": res.predict(start=714, end=726),
"Restr. AR(13)": res_glob.predict(start=714, end=726),
}
)
_, ax = plt.subplots()
ax = predictions.plot(ax=ax)
診断結果は、モデルがデータの大部分の動態を捉えていることを示しています。自己相関関数（ACF）は季節的な周波数でパターンを示しており、より完全な季節モデル（SARIMAX）が必要である可能性があります。
[23]:
fig = plt.figure(figsize=(16, 9))
fig = res_glob.plot_diagnostics(fig=fig, lags=30)
予測¶
予測は、結果インスタンスの predict メソッドを使用して生成されます。デフォルトでは、1 ステップの予測を行う静的予測が生成されます。複数ステップの予測を行うには、dynamic=True を使用する必要があります。
次のセルでは、標本の最後の 24 期間に対して12ステップの予測を行います。このためにはループが必要です。
注：これらは技術的には標本内予測です。なぜなら、予測に使用されるデータはパラメータを推定するために使用されたからです。標本外（アウト・オブ・標本、OOS）予測を行うには、二つのモデルが必要です。最初のモデルは OOS 期間を除外し、2 番目のモデルは、OOS 期間を除外した短い標本モデルのパラメータを使用して、フル標本モデルから predict メソッドを呼び出します。
[24]:
import numpy as np
start = ind_prod.index[-24]
forecast_index = pd.date_range(start, freq=ind_prod.index.freq, periods=36)
cols = ["-".join(str(val) for val in (idx.year, idx.month)) for idx in forecast_index]
forecasts = pd.DataFrame(index=forecast_index, columns=cols)
for i in range(1, 24):
fcast = res_glob.predict(
start=forecast_index[i], end=forecast_index[i + 12], dynamic=True
)
forecasts.loc[fcast.index, cols[i]] = fcast
_, ax = plt.subplots(figsize=(16, 10))
ind_prod.iloc[-24:].plot(ax=ax, color="black", linestyle="--")
ax = forecasts.plot(ax=ax)
SARIMAXとの比較¶
SARIMAXは、外生変数を含む季節的自己回帰統合移動平均モデル（Seasonal Autoregressive Integrated Moving Average with eXogenous regressors）の実装です。以下の機能をサポートしています：
季節的および非季節的な AR および MA コンポーネントの指定
外生変数の含有
カルマンフィルターを使用した最尤推定の完全実行
このモデルは AutoReg よりも機能が豊富です。 SARIMAX とは異なり、 AutoReg は OLS（最小二乗法）を用いてパラメータを推定します。これにより計算が高速で、問題はグローバルに凸であるため、局所的な最小値に関する問題はありません。閉形式推定量とその性能は、 AR(P) モデルを比較する際に AutoReg が SARIMAX に対して持つ主な利点です。 AutoReg は季節的ダミー変数もサポートしており、ユーザーがそれを外生変数として含めることで SARIMAX と一緒に使用することもできます。
[25]:
from statsmodels.tsa.api import SARIMAX
sarimax_mod = SARIMAX(ind_prod, order=((1, 5, 12, 13), 0, 0), trend="c")
sarimax_res = sarimax_mod.fit()
print(sarimax_res.summary())
SARIMAX Results
=========================================================================================
Dep. Variable: INDPRO No. Observations: 714
Model: SARIMAX([1, 5, 12, 13], 0, 0) Log Likelihood 2303.920
Date: Thu, 28 Nov 2024 AIC -4595.840
Time: 23:11:22 BIC -4568.415
Sample: 01-01-1960 HQIC -4585.248
- 06-01-2019
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
intercept 0.0011 0.000 2.534 0.011 0.000 0.002
ar.L1 1.0801 0.010 107.331 0.000 1.060 1.100
ar.L5 -0.0847 0.011 -7.592 0.000 -0.107 -0.063
ar.L12 -0.4428 0.026 -17.302 0.000 -0.493 -0.393
ar.L13 0.4073 0.025 16.213 0.000 0.358 0.457
sigma2 9.136e-05 3.08e-06 29.630 0.000 8.53e-05 9.74e-05
===================================================================================
Ljung-Box (L1) (Q): 21.77 Jarque-Bera (JB): 959.65
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.37 Skew: -0.63
Prob(H) (two-sided): 0.00 Kurtosis: 8.54
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[26]:
sarimax_params = sarimax_res.params.iloc[:-1].copy()
sarimax_params.index = res_glob.params.index
params = pd.concat([res_glob.params, sarimax_params], axis=1, sort=False)
params.columns = ["AutoReg", "SARIMAX"]
params
[26]:
| AutoReg | SARIMAX | |
|---|---|---|
| const | 0.001234 | 0.001085 |
| INDPRO.L1 | 1.088761 | 1.080116 |
| INDPRO.L5 | -0.105665 | -0.084738 |
| INDPRO.L12 | -0.388374 | -0.442793 |
| INDPRO.L13 | 0.362319 | 0.407338 |
カスタム決定論的プロセス¶
deterministic パラメータを使用すると、カスタムの DeterministicProcess を使用することができます。これにより、たとえば二つの期間を持つ季節成分を含むような、より複雑な決定論的項を構築することができます。また、次の例のように、季節的ダミーではなくフーリエ級数を使用することも可能です。
[27]:
from statsmodels.tsa.deterministic import DeterministicProcess
dp = DeterministicProcess(housing.index, constant=True, period=12, fourier=2)
mod = AutoReg(housing, 2, trend="n", seasonal=False, deterministic=dp)
res = mod.fit()
print(res.summary())
AutoReg Model Results
==============================================================================
Dep. Variable: HOUSTNSA No. Observations: 725
Model: AutoReg(2) Log Likelihood -2716.505
Method: Conditional MLE S.D. of innovations 10.364
Date: Thu, 28 Nov 2024 AIC 5449.010
Time: 23:11:22 BIC 5485.677
Sample: 04-01-1959 HQIC 5463.163
- 06-01-2019
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.7550 0.391 4.485 0.000 0.988 2.522
sin(1,12) 16.7443 0.860 19.478 0.000 15.059 18.429
cos(1,12) 4.9409 0.588 8.409 0.000 3.789 6.093
sin(2,12) 12.9364 0.619 20.889 0.000 11.723 14.150
cos(2,12) -0.4738 0.754 -0.628 0.530 -1.952 1.004
HOUSTNSA.L1 -0.3905 0.037 -10.664 0.000 -0.462 -0.319
HOUSTNSA.L2 -0.1746 0.037 -4.769 0.000 -0.246 -0.103
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 -1.1182 -2.1159j 2.3932 -0.3274
AR.2 -1.1182 +2.1159j 2.3932 0.3274
-----------------------------------------------------------------------------
[28]:
fig = res.plot_predict(720, 840)
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
最終更新日:
2025年01月28日