【気象データで時系列解析⑦】SARIMAモデル

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

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

SARIMA⁡(p,d,q)⁢(P,D,Q)sは後方シフト表記B（後述）を用いると次のように書けます[^1]。
[^1]:https://ir.nctu.edu.tw/bitstream/11536/28935/1/000174392400007.pdf

```math
\phi(B)\Phi(B^s)(1-B)^d(1-B^s)^Dx_t = c + \theta(B)\Theta(B^s)\epsilon_t
```
ただし、
```math
\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になります：

```math
Bx_t \equiv x_{t-1} 
```

SARIMAで新たに出てきたパラメータP,D,QはARIMAのパラメータp,q,rと同様で、周期成分に対するAR次数、I次数、MA次数を表します。sは周期です。

# 気温をSARIMAで予測してみる
前回の記事ではARIMAを試してみましたが、振幅が減衰していってしまうという問題がありました。

https://qiita.com/matsuda_tkm/items/2ed9b56165cd36c3c6c6

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/7b8d951d-6576-f996-4899-220eb174365a.png)

では、SARIMAではどうなるでしょうか。モデル構築方法は前回のARIMAとほぼ同じですので説明は割愛します。

```python
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).
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/a96c8f04-e33e-bf3c-429c-e39d38928787.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/4148c5b3-04da-fc02-8a8a-46ebea8ddf0b.png)

ARIMAの場合はどんどん振幅が減衰していってしまいましたが、SARIMAでは周期性を入れたことで振幅がほぼ一定に保たれています。

また、前回は紹介し損ねたのですが、以下のコードで信頼区間も取得することができます。

```python
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%信頼区間
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/f587911d-c362-2b2e-a998-f9f80c9aebb2.png)

グレーの部分が95％信頼区間です。

# ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています。

本記事のコードは`01-JMA_data_analytics`の中に入っています。

https://github.com/matsuda-tkm/qiita-calender-2023