概要
シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、季節性ARIMA（SARIMA）モデルについて扱います。
SARIMAモデル
SARIMA（Seasonal ARIMA） は、ARIMAモデルを周期方向にも適用したモデルです。
SARIMA(p,d,q)(P,D,Q)sは後方シフト表記B（後述）を用いると次のように書けます1。
\phi(B)\Phi(B^s)(1-B)^d(1-B^s)^Dx_t = c + \theta(B)\Theta(B^s)\epsilon_t
ただし、
\begin{align}
\phi(B) &= 1-\phi_1B-\cdots-\phi_pB^p\
\Phi(B^s) &= 1-\Phi_1B^s-\cdots-\Phi_PB^{Ps}\
\theta(B) &= 1-\theta_1B-\cdots-\theta_qB^q\
\Theta(B^s) &= 1-\theta_1B^s-\cdots-\theta_QB^{Qs}
\end{align}
です。
ここで使った後方シフトBは作用素で、例えばxtにBを作用させるとxt−1になります：
Bx_t \equiv x_{t-1}
SARIMAで新たに出てきたパラメータP,D,QはARIMAのパラメータp,q,rと同様で、周期成分に対するAR次数、I次数、MA次数を表します。sは周期です。
気温をSARIMAで予測してみる
前回の記事ではARIMAを試してみましたが、振幅が減衰していってしまうという問題がありました。
では、SARIMAではどうなるでしょうか。モデル構築方法は前回のARIMAとほぼ同じですので説明は割愛します。
from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(df_train['temp'], order=(1, 1, 1), seasonal_order=(1, 1, 1, 24))
result = model.fit()
print(result.summary())
SARIMAX Results
==========================================================================================
Dep. Variable: temp No. Observations: 26304
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 24) Log Likelihood -25011.586
Date: Sun, 24 Dec 2023 AIC 50033.172
Time: 22:46:21 BIC 50074.054
Sample: 01-01-2020 HQIC 50046.372
- 01-01-2023
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.5657 0.016 36.415 0.000 0.535 0.596
ma.L1 -0.3542 0.016 -21.557 0.000 -0.386 -0.322
ar.S.L24 0.0381 0.004 8.494 0.000 0.029 0.047
ma.S.L24 -0.9655 0.001 -679.107 0.000 -0.968 -0.963
sigma2 0.3919 0.001 267.251 0.000 0.389 0.395
===================================================================================
Ljung-Box (L1) (Q): 5.32 Jarque-Bera (JB): 97640.56
Prob(Q): 0.02 Prob(JB): 0.00
Heteroskedasticity (H): 1.03 Skew: -0.85
Prob(H) (two-sided): 0.14 Kurtosis: 12.29
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
ARIMAの場合はどんどん振幅が減衰していってしまいましたが、SARIMAでは周期性を入れたことで振幅がほぼ一定に保たれています。
また、前回は紹介し損ねたのですが、以下のコードで信頼区間も取得することができます。
pred = result.get_prediction(start=pd.to_datetime('2023-01-01'), end=pd.to_datetime('2023-01-31'))
pred_ci = pred.conf_int(alpha=0.05) # 95%信頼区間
グレーの部分が95％信頼区間です。
ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています。
本記事のコードは01-JMA_data_analyticsの中に入っています。