概要

シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、季節性ARIMA（SARIMA）モデルについて扱います。

SARIMAモデル

SARIMA（Seasonal ARIMA） は、ARIMAモデルを周期方向にも適用したモデルです。

SARIMA⁡(p,d,q)⁢(P,D,Q)sは後方シフト表記B（後述）を用いると次のように書けます¹。

\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の中に入っています。

https://ir.nctu.edu.tw/bitstream/11536/28935/1/000174392400007.pdf ↩