【実践時系列解析】電力需要をSARIMAモデルで予測する

前回の記事「電力需要をARIMAモデルで予測する」に引き続き、電力需要だったら季節性があるだろうということで、季節性を考慮できるSARIMAモデルを利用して電力需要を予測する

前回記事：

https://qiita.com/icchi-ioo/items/bcbb8f439122887b5e54

当記事の内容は
・SARIMAモデルを適用し予測結果を可視化
・SARIMAモデルの季節性のパラメータをどのように調整したか
の2点をメインに記述していく

# import
```python
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
warnings.filterwarnings("ignore")
import statsmodels.api as sm
import pmdarima as pm
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
```

# パラメータの調整

季節性を発見するための方法ですが、1年単位や半年単位、1ヶ月単位、１週間単位など人間の生活リズムなどと関係があるという想定でその辺りの数字を探索します。自己相関や偏自己相関のグラフで可視化して、季節性の数字が妥当かどうか確認します。

```python
for i in [1,2]:
    for j in [3,7,14,30]:
        sns.lineplot(group_train.POWER.diff(i).diff(j))
        plt.show()

for i in [1,2]:
    for j in [3,7,14,30]:
        #自己相関係数のプロット
        sm.graphics.tsa.plot_acf(group_train.POWER.diff(i).diff(j).dropna(), lags=150)
        plt.xlabel(f'lags,i:{i},j:{j}')
        plt.ylabel('acf')
        

        #偏自己相関係数のプロット
        sm.graphics.tsa.plot_pacf(group_train.POWER.diff(i).diff(j).dropna(), lags=150)
        plt.xlabel('lags')
        plt.ylabel('pacf')
```

季節性を除いたからといって、定常性は発見できなかった。
原因として、電力需要は毎夏や毎冬などの年単位の季節性が存在すると考えられるが、データが1年分しかなく季節性を除外できなかったのかもしれない。

# 予測モデルの構築

### 手動でモデルを構築

```python
sarima_model_log = SARIMAX(np.log(df_train), order=(2,1,2), seasonal_order=(0,1,0,30))
sarima_fit_log = sarima_model_log.fit()

fig, ax = plt.subplots(figsize=(20, 6))

min_lag = get_key_with_min_value(results_list)[0] + 1

forecast = np.exp(sarima_fit_log.forecast(len_test))
forecast.index = df_test.index

ax.plot(forecast, ls="-", color="r", label="predicted")
ax.plot(group_train.POWER, ls="-")
```
前半の方はうまくいっているが、後半の予測はうまくいっていない。ただし、ARIMAモデルやARモデルよりも予測はうまくいっているように思える。

![スクリーンショット 2023-06-28 16.01.22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/e021f156-c87a-f659-459a-a2ecdf3cbc47.png)


### Auto_ARIMAを利用してモデル構築
Auto_arimaはおそらく、AICの最小化を基準にパラメータを選択しているように思われる。

```python
pm.auto_arima(np.log(df_train),
             seasonal=True,
             m=30,
             n_jobs=1,
             maxiter=100)
```
結論は下写真の通り

![スクリーンショット 2023-06-28 16.06.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/e919972c-3e52-0d08-e079-7c5cb888ecf0.png)

```python
sarima_model_log = SARIMAX(np.log(df_train), order=(1,1,2), seasonal_order=(0,0,0,30))
sarima_fit_log = sarima_model_log.fit()

fig, ax = plt.subplots(figsize=(20, 6))

min_lag = get_key_with_min_value(results_list)[0] + 1

forecast = np.exp(sarima_fit_log.forecast(len_test))
forecast.index = df_test.index

ax.plot(forecast, ls="-", color="r", label="predicted")
ax.plot(group_train.POWER, ls="-")
```

予測が全くうまくいっていない

![スクリーンショット 2023-06-28 16.08.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/38cd312d-86c4-e0c8-296a-98674314a682.png)

# モデルの妥当性
妥当性の確認方法は
・残差をプロット後定常性の確認
・残差のQQプロット
・残差の偏自己相関の確認
などの方法がある。

今回は残差の偏自己相関を見て、モデルの妥当性を評価する。

```python
# 残差の偏自己相関を確認
sm.graphics.tsa.plot_pacf(sarima_fit.resid, lags=50)
plt.xlabel('lags')
plt.ylabel('pacf')
```

ラグ３０あたりで偏自己相関係数の値が有意になっている。モデル自体はうまくいっていない可能性が高い。

![スクリーンショット 2023-06-28 16.13.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/f7793fb0-9906-16c2-d3ab-3b45d996d680.png)


# 参考文献

https://www.salesanalytics.co.jp/datascience/datascience087/

モデルの妥当性

https://bigdata-tools.com/arima-sarima-model/#SARIMA

https://qiita.com/atsuken/items/0c108d098ef54b69fadd#3sarima%E3%83%A2%E3%83%87%E3%83%AB%E3%81%A8%E3%81%AF

https://ai-trend.jp/basic-study/time-series-analysis/sarima_model/