SARIMAX: はじめに¶
このノートブックは、StataのARIMA時系列推定および事後分析のドキュメントからの例を再現します。
まず、以下の4つの推定例を再現します：http://www.stata.com/manuals13/tsarima.pdf
米国卸売物価指数（WPI）データセットに対するARIMA(1,1,1)モデル。
例1の変形で、ARIMA(1,1,1)仕様にMA(4)項を追加し、加法的な季節効果を考慮したもの。
月次航空会社データに対するARIMA(2,1,0) x (1,1,0,12)モデル。この例では、乗法的な季節効果を考慮しています。
外生変数を持つARMA(1,1)モデル。消費を自己回帰過程として記述し、マネーサプライも説明変数であると仮定しています。
次に、事後分析機能を使用してhttp://www.stata.com/manuals13/tsarimapostestimation.pdfを再現します。例4のモデルを使用して、以下を示します：
1ステップ先の標本内予測
nステップ先の標本外予測
nステップ先の標本内動的予測
[1]:
%matplotlib inline
[2]:
import numpy as np
import pandas as pd
from scipy.stats import norm
import statsmodels.api as sm
import matplotlib.pyplot as plt
from datetime import datetime
import requests
from io import BytesIO
# Register converters to avoid warnings
pd.plotting.register_matplotlib_converters()
plt.rc("figure", figsize=(16,8))
plt.rc("font", size=14)
ARIMAの例1: ARIMA¶
例2のグラフからわかるように、卸売物価指数（WPI）は時間とともに増加しており（すなわち定常ではない）、そのためARMAモデルは適切な仕様ではありません。この最初の例では、元の時系列が1次の統合過程であると仮定し、その差分が定常であると仮定して、1つの自己回帰ラグと1つの移動平均ラグ、および切片項を含むモデルを適合させます。
仮定されたデータ過程は次のようになります：
ここで、\(c\)はARMAモデルの切片、\(\Delta\)は1階差分演算子、\(\epsilon_{t} \sim N(0, \sigma^2)\)と仮定します。これを遅延多項式を強調する形で書き換えると（以下の例2で有用です）：
ここで、\(L\)は遅延演算子です。
Stataの出力と以下の出力との違いの1つに注目してください。Stataは次のモデルを推定します：
ここで、\(\beta_0\)は過程\(y_t\)の平均です。このモデルは、statsmodelsのSARIMAXクラスで推定されたモデルと同等ですが、解釈が異なります。等価性を確認するために、次のように注意してください：
したがって、\(c = (1 - \phi_1) \beta_0\)です。
[3]:
# Dataset
wpi1 = requests.get('https://www.stata-press.com/data/r12/wpi1.dta').content
data = pd.read_stata(BytesIO(wpi1))
data.index = data.t
# Set the frequency
data.index.freq="QS-OCT"
# Fit the model
mod = sm.tsa.statespace.SARIMAX(data['wpi'], trend='c', order=(1,1,1))
res = mod.fit(disp=False)
print(res.summary())
SARIMAX Results
==============================================================================
Dep. Variable: wpi No. Observations: 124
Model: SARIMAX(1, 1, 1) Log Likelihood -135.351
Date: Wed, 28 Aug 2024 AIC 278.703
Time: 15:00:51 BIC 289.951
Sample: 01-01-1960 HQIC 283.272
- 10-01-1990
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
intercept 0.0943 0.068 1.389 0.165 -0.039 0.227
ar.L1 0.8742 0.055 16.028 0.000 0.767 0.981
ma.L1 -0.4120 0.100 -4.119 0.000 -0.608 -0.216
sigma2 0.5257 0.053 9.849 0.000 0.421 0.630
===================================================================================
Ljung-Box (L1) (Q): 0.09 Jarque-Bera (JB): 9.78
Prob(Q): 0.77 Prob(JB): 0.01
Heteroskedasticity (H): 15.93 Skew: 0.28
Prob(H) (two-sided): 0.00 Kurtosis: 4.26
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
このように、最大尤度推定値は、上記のプロセスに対して次のような式を示唆しています：
ここで、\(\epsilon_{t} \sim N(0, 0.5257)\)です。最後に、\(c = (1 - \phi_1) \beta_0\)を思い出してください。ここでは\(c = 0.0943\)および\(\phi_1 = 0.8742\)です。Stataの出力と比較するために、平均を計算できます：
注：この値は、Stataのドキュメントに記載された値\(\beta_0 = 0.7498\)とほぼ同じです。わずかな違いは、丸め誤差や数値最適化アルゴリズムの停止基準の微妙な違いによるものと思われます。
ARIMA例2: 加法季節効果を持つARIMAモデル¶
このモデルは、例1のモデルの拡張です。ここでは、データが次のプロセスに従うと仮定されています：
このモデルの新しい部分は、年間の季節効果を許容することです（周期が4であっても、データセットが四半期ベースであるため、年間季節効果と見なされます）。2つ目の違いは、このモデルがデータのレベルではなく、ログを使用することです。
データセットを推定する前に、次のグラフを表示します：
時系列（ログで表示）
時系列の1階差分（ログで表示）
自己相関関数
部分自己相関関数
最初の2つのグラフから、元の時系列は定常でないように見えるのに対し、1階差分は定常であることがわかります。これにより、データの1階差分に基づいてARMAモデルを推定するか、または1階差分を持つARIMAモデルを推定するアプローチが支持されます（後者のアプローチを選択していることを思い出してください）。最後の2つのグラフは、ARMA(1,1,1)モデルの使用を支持しています。
[4]:
# Dataset
data = pd.read_stata(BytesIO(wpi1))
data.index = data.t
data.index.freq="QS-OCT"
data['ln_wpi'] = np.log(data['wpi'])
data['D.ln_wpi'] = data['ln_wpi'].diff()
[5]:
# Graph data
fig, axes = plt.subplots(1, 2, figsize=(15,4))
# Levels
axes[0].plot(data.index._mpl_repr(), data['wpi'], '-')
axes[0].set(title='US Wholesale Price Index')
# Log difference
axes[1].plot(data.index._mpl_repr(), data['D.ln_wpi'], '-')
axes[1].hlines(0, data.index[0], data.index[-1], 'r')
axes[1].set(title='US Wholesale Price Index - difference of logs');
[6]:
# Graph data
fig, axes = plt.subplots(1, 2, figsize=(15,4))
fig = sm.graphics.tsa.plot_acf(data.iloc[1:]['D.ln_wpi'], lags=40, ax=axes[0])
fig = sm.graphics.tsa.plot_pacf(data.iloc[1:]['D.ln_wpi'], lags=40, ax=axes[1])
このモデルをstatsmodelsでどのように指定するかを理解するために、まず例1で使用した以下のコードを思い出してください。これはARIMA(1,1,1)モデルを指定するためのコードです：
mod = sm.tsa.statespace.SARIMAX(data['wpi'], trend='c', order=(1,1,1))
order引数は、次の形式のタプルです：(AR指定、差分順序、MA指定)。差分順序は整数でなければなりません（例えば、ここでは1回の差分を仮定しているので1と指定しています。もしデータがすでに定常であれば、純粋なARMAモデルでは0になります）。
AR指定とMA指定の部分には2つの可能性があります。1つ目は、対応する遅延多項式の最大次数を指定する方法で、この場合、成分は整数になります。例えば、ARIMA(1,1,4)過程を指定したい場合、次のように書きます：
mod = sm.tsa.statespace.SARIMAX(data['wpi'], trend='c', order=(1,1,4))
そして、対応するデータ過程は次のようになります：
または
遅延多項式の最大次数として指定された場合、指定された次数までのすべての多項式項が含まれます。これは求めているモデルではないことに注意してください。なぜなら、\(\epsilon_{t-2}\)や\(\epsilon_{t-3}\)の項が含まれてしまうからです。
私たちが求めているのは、1次と4次の項を持ち、2次と3次の項を除外する多項式です。そのためには、指定パラメータとしてタプルを指定し、そのタプルで遅延多項式自体を記述する必要があります。具体的には、次のように指定します：
ar = 1 # これは最大次数指定
ma = (1,0,0,1) # これは遅延多項式指定
mod = sm.tsa.statespace.SARIMAX(data['wpi'], trend='c', order=(ar,1,ma)))
これにより、データ過程は次のような形になります：
これが私たちが求めているものです。
[7]:
# Fit the model
mod = sm.tsa.statespace.SARIMAX(data['ln_wpi'], trend='c', order=(1,1,(1,0,0,1)))
res = mod.fit(disp=False)
print(res.summary())
SARIMAX Results
=================================================================================
Dep. Variable: ln_wpi No. Observations: 124
Model: SARIMAX(1, 1, [1, 4]) Log Likelihood 386.034
Date: Wed, 28 Aug 2024 AIC -762.067
Time: 15:00:52 BIC -748.006
Sample: 01-01-1960 HQIC -756.355
- 10-01-1990
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
intercept 0.0024 0.002 1.482 0.138 -0.001 0.006
ar.L1 0.7811 0.094 8.274 0.000 0.596 0.966
ma.L1 -0.3996 0.126 -3.178 0.001 -0.646 -0.153
ma.L4 0.3089 0.120 2.574 0.010 0.074 0.544
sigma2 0.0001 9.81e-06 11.108 0.000 8.97e-05 0.000
===================================================================================
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 45.11
Prob(Q): 0.90 Prob(JB): 0.00
Heteroskedasticity (H): 2.58 Skew: 0.29
Prob(H) (two-sided): 0.00 Kurtosis: 5.91
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
ARIMAの例 3: 航空会社モデル¶
前の例では、季節効果を加法的な方法で含めました。これは、プロセスが4つ目のMA遅延に依存するような項を追加することを意味しています。しかし、季節効果を乗法的な方法でモデル化したい場合もあります。この場合、モデルは通常、ARIMA \((p,d,q) \times (P,D,Q)_s\) の形で記述されます。ここで、小文字の文字は非季節成分の仕様を、同じく大文字の文字は季節成分の仕様を示し、\(s\)は季節の周期（例えば、四半期データの場合は4、月次データの場合は12など）を示します。データプロセスは以下のように一般的に表現できます：
ここで：
\(\phi_p (L)\) は非季節的な自己回帰遅延多項式
\(\tilde \phi_P (L^s)\) は季節的な自己回帰遅延多項式
\(\Delta^d \Delta_s^D y_t\) は時系列で、\(d\)回差分をとり、\(D\)回季節差分をとったもの
\(A(t)\) はトレンド多項式（切片を含む）
\(\theta_q (L)\) は非季節的な移動平均遅延多項式
\(\tilde \theta_Q (L^s)\) は季節的な移動平均遅延多項式
時には、次のように書き換えることもあります：
ここで、\(y_t^* = \Delta^d \Delta_s^D y_t\) です。これは、単純な場合と同様に、差分（ここでは非季節的および季節的差分）をとってデータを定常化した後、結果として得られるモデルがARMAモデルであることを強調しています。
例として、航空会社モデルARIMA \((2,1,0) \times (1,1,0)_{12}\)（切片あり）を考えます。データプロセスは次のように表現できます：
ここでは：
\(\phi_p (L) = (1 - \phi_1 L - \phi_2 L^2)\)
\(\tilde \phi_P (L^s) = (1 - \phi_1 L^{12})\)
\(d = 1, D = 1, s=12\) は、\(y_t^*\) が \(y_t\) から1回の差分と12回の季節差分をとって導出されることを示しています。
\(A(t) = c\) は定数トレンド多項式（つまり、切片のみ）
\(\theta_q (L) = \tilde \theta_Q (L^s) = 1\)（つまり、移動平均効果はない）
時間系列変数の前に2つの遅延多項式があるのはまだ混乱するかもしれませんが、遅延多項式を掛け合わせて次のモデルを得ることができることに注目してください：
これは次のように書き換えることができます：
これは例2の加法的季節モデルに似ていますが、自己回帰の遅延の前の係数は、実際には基礎となる季節的および非季節的なパラメータの組み合わせです。
statsmodelsでモデルを指定するには、seasonal_order 引数を追加するだけです。この引数は、(季節的AR仕様、季節的差分順序、季節的MA、季節的周期性) の形式でタプルを受け入れます。季節的ARとMAの仕様は、前と同じように、最大多項式次数で表現することも、遅延多項式そのもので表現することもできます。季節的周期性は整数で表されます。
航空会社モデルARIMA \((2,1,0) \times (1,1,0)_{12}\)（切片あり）の場合、コマンドは次のようになります：
mod = sm.tsa.statespace.SARIMAX(data['lnair'], order=(2,1,0), seasonal_order=(1,1,0,12))
[8]:
# Dataset
air2 = requests.get('https://www.stata-press.com/data/r12/air2.dta').content
data = pd.read_stata(BytesIO(air2))
data.index = pd.date_range(start=datetime(int(data.time[0]), 1, 1), periods=len(data), freq='MS')
data['lnair'] = np.log(data['air'])
# Fit the model
mod = sm.tsa.statespace.SARIMAX(data['lnair'], order=(2,1,0), seasonal_order=(1,1,0,12), simple_differencing=True)
res = mod.fit(disp=False)
print(res.summary())
SARIMAX Results
==========================================================================================
Dep. Variable: D.DS12.lnair No. Observations: 131
Model: SARIMAX(2, 0, 0)x(1, 0, 0, 12) Log Likelihood 240.821
Date: Wed, 28 Aug 2024 AIC -473.643
Time: 15:00:53 BIC -462.142
Sample: 02-01-1950 HQIC -468.970
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.4057 0.080 -5.045 0.000 -0.563 -0.248
ar.L2 -0.0799 0.099 -0.809 0.419 -0.274 0.114
ar.S.L12 -0.4723 0.072 -6.592 0.000 -0.613 -0.332
sigma2 0.0014 0.000 8.403 0.000 0.001 0.002
===================================================================================
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 0.72
Prob(Q): 0.91 Prob(JB): 0.70
Heteroskedasticity (H): 0.54 Skew: 0.14
Prob(H) (two-sided): 0.04 Kurtosis: 3.23
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
ここで、追加の引数 simple_differencing=True を使用したことに注意してください。これはARIMAモデルにおける積分の順序がどのように処理されるかを制御します。simple_differencing=True の場合、endog として提供された時系列は文字通り差分を取られ、その結果得られた新しい時系列にARMAモデルが適合されます。これにより、差分処理によりいくつかの初期の期間が失われますが、他のパッケージ（例えば、Stataの arima
は常に単純差分を使用）の結果と比較する場合や、季節的周期が大きい場合には必要となることがあります。
デフォルトは simple_differencing=False で、この場合、積分成分は状態空間の定式化の一部として実装され、元のデータ全体が推定に使用されます。
ARIMA例4: ARMAX（フリードマン）¶
このモデルは説明変数（ARMAXのX部分）の使用を示しています。外生的回帰変数が含まれる場合、SARIMAXモジュールは「SARIMA誤差を伴う回帰」の概念を使用します（ARIMA誤差を伴う回帰と他の仕様の詳細については、http://robjhyndman.com/hyndsight/arimax/ を参照）。したがって、モデルは以下のように指定されます：
最初の方程式は単なる線形回帰であり、2番目の方程式は誤差成分がSARIMAに従う過程を説明しています（例3で説明したように）。この仕様が選ばれている理由の1つは、推定されたパラメータが自然な解釈を持つからです。
この仕様は多くの簡単な仕様を包含しています。例えば、AR(2)誤差を伴う回帰は次のように記述されます：
この例で考慮されるモデルは、ARMA(1,1)誤差を伴う回帰です。その過程は次のように書かれます：
ここで、\(\beta_0\) は上記の例1で説明したように、trend='c' で指定された切片とは異なることに注意してください。上記の例では、モデルの切片をトレンド多項式を通じて推定しましたが、ここでは外生的データセットに定数を追加することで\(\beta_0\)自体を推定する方法を示しています。出力では、\(beta_0\) は const と呼ばれ、上記では切片\(c\)が出力で intercept と呼ばれていたことに注意してください。
[9]:
# Dataset
friedman2 = requests.get('https://www.stata-press.com/data/r12/friedman2.dta').content
data = pd.read_stata(BytesIO(friedman2))
data.index = data.time
data.index.freq = "QS-OCT"
# Variables
endog = data.loc['1959':'1981', 'consump']
exog = sm.add_constant(data.loc['1959':'1981', 'm2'])
# Fit the model
mod = sm.tsa.statespace.SARIMAX(endog, exog, order=(1,0,1))
res = mod.fit(disp=False)
print(res.summary())
SARIMAX Results
==============================================================================
Dep. Variable: consump No. Observations: 92
Model: SARIMAX(1, 0, 1) Log Likelihood -340.508
Date: Wed, 28 Aug 2024 AIC 691.015
Time: 15:00:53 BIC 703.624
Sample: 01-01-1959 HQIC 696.105
- 10-01-1981
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -36.0612 56.642 -0.637 0.524 -147.078 74.956
m2 1.1220 0.036 30.825 0.000 1.051 1.193
ar.L1 0.9348 0.041 22.717 0.000 0.854 1.016
ma.L1 0.3091 0.089 3.488 0.000 0.135 0.483
sigma2 93.2556 10.889 8.565 0.000 71.914 114.597
===================================================================================
Ljung-Box (L1) (Q): 0.04 Jarque-Bera (JB): 23.49
Prob(Q): 0.84 Prob(JB): 0.00
Heteroskedasticity (H): 22.51 Skew: 0.17
Prob(H) (two-sided): 0.00 Kurtosis: 5.45
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
ARIMA後処理: 例1 - 動的予測¶
ここでは、statsmodelsのSARIMAXのいくつかの後処理機能について説明します。
まず、例で使用したモデルを基に、最後の数観測値を除外したデータを使用してパラメータを推定します（これは例として少し人工的ですが、標本外予測のパフォーマンスを考慮することができ、Stataのドキュメントとの比較を容易にします）。
[10]:
# Dataset
raw = pd.read_stata(BytesIO(friedman2))
raw.index = raw.time
raw.index.freq = "QS-OCT"
data = raw.loc[:'1981']
# Variables
endog = data.loc['1959':, 'consump']
exog = sm.add_constant(data.loc['1959':, 'm2'])
nobs = endog.shape[0]
# Fit the model
mod = sm.tsa.statespace.SARIMAX(endog.loc[:'1978-01-01'], exog=exog.loc[:'1978-01-01'], order=(1,0,1))
fit_res = mod.fit(disp=False, maxiter=250)
print(fit_res.summary())
SARIMAX Results
==============================================================================
Dep. Variable: consump No. Observations: 77
Model: SARIMAX(1, 0, 1) Log Likelihood -243.316
Date: Wed, 28 Aug 2024 AIC 496.633
Time: 15:00:54 BIC 508.352
Sample: 01-01-1959 HQIC 501.320
- 01-01-1978
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.6779 18.492 0.037 0.971 -35.565 36.921
m2 1.0379 0.021 50.329 0.000 0.997 1.078
ar.L1 0.8775 0.059 14.859 0.000 0.762 0.993
ma.L1 0.2771 0.108 2.572 0.010 0.066 0.488
sigma2 31.6978 4.683 6.769 0.000 22.520 40.876
===================================================================================
Ljung-Box (L1) (Q): 0.32 Jarque-Bera (JB): 6.05
Prob(Q): 0.57 Prob(JB): 0.05
Heteroskedasticity (H): 6.09 Skew: 0.57
Prob(H) (two-sided): 0.00 Kurtosis: 3.76
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
次に、推定されたパラメータ（データのサブセット上）を使用して、データセット全体の結果を取得します。
[11]:
mod = sm.tsa.statespace.SARIMAX(endog, exog=exog, order=(1,0,1))
res = mod.filter(fit_res.params)
predict コマンドは、まずここで標本内予測を得るために適用されます。full_results=True 引数を使用することで、信頼区間を計算できるようにします（predict のデフォルトの出力は予測値のみです）。
他の引数を指定しない場合、predict は標本全体に対する1ステップ先の標本内予測を返します。
[12]:
# In-sample one-step-ahead predictions
predict = res.get_prediction()
predict_ci = predict.conf_int()
私たちはまた、動的予測を行うことができます。一歩先の予測は、各ステップで内生変数の真の値を使用して次の標本内の値を予測します。動的予測は、データセットのある時点（dynamic引数で指定）まで一歩先の予測を使用し、その後は予測された前回の内生変数の値を使用して、新しい予測値を計算します。
dynamic引数は、start引数に対するオフセットとして指定されます。startが指定されていない場合は、0として扱われます。
ここでは、1978年第1四半期から動的予測を実行します。
[13]:
# Dynamic predictions
predict_dy = res.get_prediction(dynamic='1978-01-01')
predict_dy_ci = predict_dy.conf_int()
1 ステップ先の予測と動的予測 (および対応する信頼区間) をグラフ化して、それらの相対的なパフォーマンスを確認できます。動的予測が開始される時点 (1978:Q1) までは、2 つが同じであることに注意してください。
[14]:
# Graph
fig, ax = plt.subplots(figsize=(9,4))
npre = 4
ax.set(title='Personal consumption', xlabel='Date', ylabel='Billions of dollars')
# Plot data points
data.loc['1977-07-01':, 'consump'].plot(ax=ax, style='o', label='Observed')
# Plot predictions
predict.predicted_mean.loc['1977-07-01':].plot(ax=ax, style='r--', label='One-step-ahead forecast')
ci = predict_ci.loc['1977-07-01':]
ax.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color='r', alpha=0.1)
predict_dy.predicted_mean.loc['1977-07-01':].plot(ax=ax, style='g', label='Dynamic forecast (1978)')
ci = predict_dy_ci.loc['1977-07-01':]
ax.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color='g', alpha=0.1)
legend = ax.legend(loc='lower right')
最後に、予測 error .をグラフ化します。予想どおり、1 ステップ先の予測の方がかなり優れていることは明らかです。
[15]:
# Prediction error
# Graph
fig, ax = plt.subplots(figsize=(9,4))
npre = 4
ax.set(title='Forecast error', xlabel='Date', ylabel='Forecast - Actual')
# In-sample one-step-ahead predictions and 95% confidence intervals
predict_error = predict.predicted_mean - endog
predict_error.loc['1977-10-01':].plot(ax=ax, label='One-step-ahead forecast')
ci = predict_ci.loc['1977-10-01':].copy()
ci.iloc[:,0] -= endog.loc['1977-10-01':]
ci.iloc[:,1] -= endog.loc['1977-10-01':]
ax.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], alpha=0.1)
# Dynamic predictions and 95% confidence intervals
predict_dy_error = predict_dy.predicted_mean - endog
predict_dy_error.loc['1977-10-01':].plot(ax=ax, style='r', label='Dynamic forecast (1978)')
ci = predict_dy_ci.loc['1977-10-01':].copy()
ci.iloc[:,0] -= endog.loc['1977-10-01':]
ci.iloc[:,1] -= endog.loc['1977-10-01':]
ax.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color='r', alpha=0.1)
legend = ax.legend(loc='lower left');
legend.get_frame().set_facecolor('w')
