Pythonで時系列予測に使える機械学習モデルの実行例まとめ

# はじめに

時系列データの予測は、トレンドを把握し今後の見通しを立てるために必要な要素の一つです。この記事では、過去のデータから未来を予測する際に利用されるさまざまな機械学習モデルについてまとめて紹介します。各モデルの理論的な説明はそこそこに、モデルの名前と「動く」Pythonサンプルコードを通じて各モデルの特徴を直感的に概観することに重点を置いています。

本記事では一次元時系列データを想定したモデルのサンプルコードを紹介します。多次元データに適用可能なモデルの中には、一次元の手法の拡張として理解できるものも多くあります。この記事が複雑なモデルを理解するためのスモールステップになれば幸いです。


## 動作環境

サンプルコードは Jupyter Notebook のような環境で使用することを想定しています。というか、そもそもこの記事自体が手元の Notebook を移植してできています。

```python:import
import platform                  # Python 3.9.16
import sys

import numpy             as np   # Numpy 1.23.2
import pandas            as pd   # Pandas 2.0.3
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams.update({
    'font.size'     : 16,
    'axes.grid'     : True,
    'grid.linestyle': '--'
})

import statsmodels.api  as sm    # Statsmodels 0.14.0
import tensorflow       as tf    # Tensor Flow 2.12.0
import tensorflow.keras as keras # Keras 2.12.0

```


## 使用するデータセットについて

statsmodels に用意されているCO<sub>2</sub>濃度のデータセットを使用します。   
データ詳細: <https://www.statsmodels.org/dev/datasets/generated/co2.html>

生データのままでは扱いにくいので、ここでは欠損が少ない1965年以降のデータから月次平均を算出して使用します。


### データの読み込み
```python: load_data
import statsmodels.datasets.co2 as co2

co2_raw = co2.load().data
df = co2_raw.iloc[353:] # 1965以降のデータを抽出
df = df.resample('M').mean() # 月次データに変換 (月の最終日を取得)"
df.index.name = "YEAR"
```


### データの特徴
予測に移る前に、読み込んだデータの特徴を確認しておきましょう。
まず、STL分解でトレンドや季節性変動を見てみます。

- STL分解 <https://doc.arcgis.com/ja/insights/latest/analyze/stl.htm>

```python:STL_decomposition
# STL分解
stl = sm.tsa.STL(df["co2"]).fit()

# それぞれの成分を描画
fig, ax = plt.subplots(4, 1, figsize=(15, 12), sharex=True)

df["co2"].plot(ax=ax[0], c='black')
ax[0].set_title("Original Data")

stl.trend.plot(ax=ax[1], c='black')
ax[1].set_title("Trend")

stl.seasonal.plot(ax=ax[2], c='black')
ax[2].set_title("Seasonal")

stl.resid.plot(ax=ax[3], c='black')
ax[3].set_title("Residual")

plt.tight_layout()

plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/bcb01814-7c74-63a6-3546-234a60813c1f.png)

データ期間全体にわたる直線的な増加傾向と、一年周期の季節性変動が見られます。機械学習モデルがこれらの変動をとらえた予測を出せるかどうかが注目のポイントとなります。

続いて、自己相関を確認しておきます。

```python:correlogram
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# コレログラムを描画
fig, ax = plt.subplots(2, 1, figsize=(15,7))
plot_acf (df["co2"], lags=40, ax=ax[0]) # 自己相関係数
plot_pacf(df["co2"], lags=40, ax=ax[1]) # 偏自己相関係数
plt.xlabel('Lag [month]')
plt.tight_layout()
plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/55e38da2-a385-818f-2f34-740ffedca03d.png)

**自己相関係数**とは「過去の自分と今の自分との間の相関係数」です。自己相関がなければ、過去のデータから未来を予測することはできないということになります。今回のデータでは自己相関係数が十分大きいので、過去から学んで未来を予測することができます。

**偏自己相関係数**とは、自己相関係数から純粋に2時点間の関係を取り出したものになります。今回のデータでは、12ヶ月周期の変動があることがわかります。また、前月のCO<sub>2</sub>濃度が大きければ、当月も大きくなることもわかります。



# 準備

機械学習モデルを動かす前に、データの変換や関数の定義を行います。モデルを動かした結果だけ知りたい方は飛ばしていただいてOKです。

<details><summary>準備の内容を表示</summary>

### データの変換
機械学習モデルが扱いやすいように、数値データを [0,1] の範囲に正規化しておきます。

```python:scale_0to1
min_val, max_val = np.min(df["co2"].values), np.max(df["co2"].values)
scale = max_val - min_val
df["y"] = (df["co2"].values - min_val) / scale #正規化

display(df)
```


> 実行結果：
> <div>
> <style scoped>
>     .dataframe tbody tr th:only-of-type {
>         vertical-align: middle;
>     }
> 
>     .dataframe tbody tr th {
>         vertical-align: top;
>     }
> 
>     .dataframe thead th {
>         text-align: right;
>     }
> </style>
> <table border="1" class="dataframe">
>   <thead>
>     <tr style="text-align: right;">
>       <th></th>
>       <th>co2</th>
>       <th>y</th>
>     </tr>
>     <tr>
>       <th>YEAR</th>
>       <th></th>
>       <th></th>
>     </tr>
>   </thead>
>   <tbody>
>     <tr>
>       <th>1965-01-31</th>
>       <td>319.400</td>
>       <td>0.036486</td>
>     </tr>
>     <tr>
>       <th>1965-02-28</th>
>       <td>320.450</td>
>       <td>0.055083</td>
>     </tr>
>     <tr>
>       <th>1965-03-31</th>
>       <td>320.925</td>
>       <td>0.063496</td>
>     </tr>
>     <tr>
>       <th>1965-04-30</th>
>       <td>322.000</td>
>       <td>0.082536</td>
>     </tr>
>     <tr>
>       <th>1965-05-31</th>
>       <td>322.060</td>
>       <td>0.083599</td>
>     </tr>
>     <tr>
>       <th>...</th>
>       <td>...</td>
>       <td>...</td>
>     </tr>
>     <tr>
>       <th>2001-08-31</th>
>       <td>369.425</td>
>       <td>0.922512</td>
>     </tr>
>     <tr>
>       <th>2001-09-30</th>
>       <td>367.880</td>
>       <td>0.895147</td>
>     </tr>
>     <tr>
>       <th>2001-10-31</th>
>       <td>368.050</td>
>       <td>0.898158</td>
>     </tr>
>     <tr>
>       <th>2001-11-30</th>
>       <td>369.375</td>
>       <td>0.921626</td>
>     </tr>
>     <tr>
>       <th>2001-12-31</th>
>       <td>371.020</td>
>       <td>0.950762</td>
>     </tr>
>   </tbody>
> </table>
> <p>444 rows × 2 columns</p>
> </div>

作ったデータを、学習用データと試験用データに分割します。

```python:make_train_test
ratio = 0.8 #学習用データの割合（80%）
train_size = int(ratio * len(df)) #前半(学習用) と 後半(試験用)に分ける
test_size  = len(df) - train_size #試験用データのサイズ

train = df.iloc[:train_size] #全データから前半を抜粋
test  = df.iloc[train_size:] #後半を抜粋

# 正規化されたCO2濃度の推移を図示して確認する
fig = plt.subplots(figsize=(15, 7))
df["y"]   .plot(label='All data',   color='black')
train["y"].plot(label='Train data', color='C0')
test["y"] .plot(label='Test data',  color='C1')
plt.ylabel("Scaled")
plt.title(f"Scaled = ( Data - {min_val:.1f} ) / {scale:.1f}")
plt.legend()

plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/f160e2f2-91db-37cb-9b43-0b04e76c0435.png)


### 関数の定義
オレオレ関数で恐縮ですが、実行結果を表示するための関数を定義しておきます。

```python:my_get_rmse
def my_get_rmse(pred: pd.Series, true: pd.Series) -> pd.Series:
    """真値と予測値のRMSEを計算する関数
    Args:
        pred: 予測値
        true: 真値
    Returns:
        rmses: RMSEのリスト
    """
    y_pred = pred.values.reshape(-1, 1)
    y_true = true.values.reshape(-1, 1)
    y_error = y_true - y_pred
    rmses = np.sqrt(np.mean(y_error**2, axis=1))

    return pd.Series(rmses, index=pred.index)
```

```python:my_plot_results
def my_plot_results(pred_test: pd.Series) -> plt.figure:
    """真値と予測値の比較をプロットする関数
    Args:
        pred_test: 試験用データの予測
    Returns:
        fig: Figure object.
        ax: Axes object (左列, 元データ).
        bx: Axes object (右列, RMSE).
    """

    # STL分解
    stl_true = sm.tsa.STL(test["co2"]).fit()
    stl_test = sm.tsa.STL(pred_test)  .fit()

    # それぞれの成分を描画
    fig, axes = plt.subplots(4, 2, figsize=(15, 12), sharex=True)
    ax,  bx   = axes.transpose()[0], axes.transpose()[1]

    for i in range(4):
        ax[i].set_ylabel(["Original Data", "Trend", "Seasonal\n", "Residual"][i])

    # 元データ
    ax[0].set_title("ax[0].set_title(name)")
    test["co2"]      .plot(ax=ax[0], c='black', label='Original Data')
    pred_test        .plot(ax=ax[0], c='C1'   , label='Prediction')

    stl_true.trend   .plot(ax=ax[1], c='black')
    stl_test.trend   .plot(ax=ax[1], c='C1')

    stl_true.seasonal.plot(ax=ax[2], c='black')
    stl_test.seasonal.plot(ax=ax[2], c='C1')

    stl_true.resid   .plot(ax=ax[3], c='black')
    stl_test.resid   .plot(ax=ax[3], c='C1')

    # RMSE
    rmse_original = my_get_rmse(pred_test,         test["co2"])
    rmse_trend    = my_get_rmse(stl_test.trend,    stl_true.trend)
    rmse_seasonal = my_get_rmse(stl_test.seasonal, stl_true.seasonal)
    rmse_resid    = my_get_rmse(stl_test.resid,    stl_true.resid)

    bx[0].set_title("bx[0].set_title(name)")
    rmse_original.plot(ax=bx[0], c='C1')
    rmse_trend   .plot(ax=bx[1], c='C1')
    rmse_seasonal.plot(ax=bx[2], c='C1')
    rmse_resid   .plot(ax=bx[3], c='C1')

    ax[0].legend()
    plt.tight_layout()

    return fig, ax, bx
```
</details>



# 機械学習モデルまとめ

## 古典的統計モデル

### 自己回帰モデル (Auto Regressive Model, AR)
$$
\begin{aligned}
y_t &= c + φ_1 \ y_{t-1} + φ_2 \ y_{t-2} + ... + φ_p \ y_{t-p} + ε_t \cr
    &= c + \sum^p_{i=1} φ_i \ y_{t-i} + ε_t
\end{aligned}
$$

未来のデータは過去のデータの線型結合である、と予想するモデルです。
どれだけ過去に遡るかを表すパラメータ $p$ は、あらかじめ決定しておく必要があります (ハイパーパラメータ)。今回は $p=10$ としています。

```python:AR(10)
p = 10

ar = sm.tsa.SARIMAX(train["y"].values,
                    order=(p, 0, 0)
                    ).fit(maxiter=1000)

pred_val = ar.forecast(test_size)

pred_test = pd.Series(pred_val * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title(f"Auto Regressive Model (p={p})")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ```
> UserWarning: Non-stationary starting autoregressive parameters found.
> Using zeros as starting parameters.
> ```
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/c2c5c119-2b96-2c06-2436-89eed4340987.png)

季節性周期変動は捉えられていますが、トレンドは捉えられていません。また、実行すると警告 `UserWarning` が表示されます。この理由については後で触れます。

:::note info
チューニングパラメータ $φ_i$ は、予測したい時刻の $i$ ステップ前のデータが予測したい時刻のデータにどのような影響を与えているかを表しています。各係数の値や学習済みのモデルに関する情報は `summary()` メソッドで確認できます。
```python:model_summary
display(ar.summary())
```
```shell:result
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  355
Model:              SARIMAX(10, 0, 0)   Log Likelihood                1138.927
Date:                Mon, 18 Sep 2023   AIC                          -2255.853
Time:                        18:25:58   BIC                          -2213.260
Sample:                             0   HQIC                         -2238.909
                                - 355                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          1.2297      0.059     20.794      0.000       1.114       1.346
ar.L2         -0.1626      0.083     -1.955      0.051      -0.326       0.000
ar.L3         -0.6233      0.086     -7.273      0.000      -0.791      -0.455
...
ar.L9         -0.0172      0.087     -0.197      0.844      -0.188       0.154
ar.L10         0.2637      0.058      4.581      0.000       0.151       0.377
sigma2      9.041e-05   7.19e-06     12.566      0.000    7.63e-05       0.000
===================================================================================
Ljung-Box (L1) (Q):                  27.12   Jarque-Bera (JB):                 0.57
Prob(Q):                              0.00   Prob(JB):                         0.75
Heteroskedasticity (H):               1.06   Skew:                             0.09
Prob(H) (two-sided):                  0.74   Kurtosis:                         3.08
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
:::

#### さらに詳しい内容を知りたい方へ

- ARモデル | 東邦大学メディアネットセンター: 
  <https://www.mnc.toho-u.ac.jp/v-lab/yobology/AR_model/AR_model.htm>
- 9.3 自己回帰モデル | 予測: 原理と実践 (第3版):
  <https://otexts.com/fppjp/AR.html>


### 移動平均モデル (Moving Average Model, MA)
$$
\begin{aligned}
y_t &= c + ε_t + θ_1ε_{t - 1} + θ_2ε_{t - 2} + ... + θ_qε_{t - q}\cr
    &= c + ε_t + \sum^q_{i=1}θ_iε_{t - i}
\end{aligned}
$$

未来のデータは過去の誤差項（予測値と実際の値との差）の線型結合である、と予想するモデルです ($\varepsilon _i$ は分散 $\sigma ^2$ のホワイトノイズ)。
どれだけ過去に遡るかを表すパラメータ $q$ は、あらかじめ決定しておく必要があります (ハイパーパラメータ)。今回は $q=8$ としています。

```python:MA(8)
q = 8

ma = sm.tsa.SARIMAX(train["y"].values,
                    order=(0, 0, q)
                    ).fit(maxiter=1000)

pred_val = ma.forecast(test_size)

pred_test = pd.Series(pred_val * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title(f"Moving Average Model (q={q})")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ```
> UserWarning: Non-stationary starting MA parameters found.
> Using zeros as starting parameters.
> ```
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/6b4e95a0-e3d6-e066-a857-c997b1da4f43.png)

うまく予測できませんでした。また、実行すると警告 `UserWarning` が表示されます。この理由については後で触れます。

#### さらに詳しい内容を知りたい方へ
- 9.4 移動平均モデル | 予測: 原理と実践 (第3版):
  <https://otexts.com/fppjp/MA.html>



### 自己回帰移動平均モデル (ARMA)
$$
\begin{aligned}
y_t &= c + φ_1 y_{t-1} + φ_2y_{t-2} + ... +φ_py_{t-p} + ε_t + θ_1ε_{t - 1} + θ_2ε_{t - 2} + ... + θ_qε_{t - q}\cr
    &= c + ε_t + \sum^p_{i=1}φ_iy_{t-i} + \sum^q_{i=1}θ_iε_{t - i}
\end{aligned}
$$

自己回帰モデル と 移動平均モデル の単純な足し合わせになっています (AR + MA)。

```python:ARMA(10, 8)
p, q = 10, 8

arma = sm.tsa.SARIMAX(train["y"].values,
                      order=(p, 0, q)
                      ).fit(maxiter=1000)

pred_val = arma.forecast(test_size)

pred_test = pd.Series(pred_val * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title(f"ARMA{(p, q)}")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ```
> UserWarning: Non-stationary starting autoregressive parameters found.
> Using zeros as starting parameters.
> ...
> STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT 
> ConvergenceWarning: Maximum Likelihood optimization failed to 
> converge. Check mle_retvals
> ```
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/2ac40836-7374-9190-7ed0-922fc10e5aad.png)

ARモデルの場合と似た結果になりました。また、実行すると警告 `UserWarning` が表示されます。

#### 警告の意味 - 弱定常性
実は、AR、MA、ARMAモデルでは長期的なトレンドを捉えることができません。ARMAモデルを使用できるのは時間の推移に対して安定的に変動する時系列データに限られます。このような性質を **弱定常性** といいます。
今回のCO<sub>2</sub>濃度データにはデータ期間全体にわたる長期的な増加トレンドがあり、非定常なデータとなっています。このため、ARMAモデルではうまく将来を予測することができず、実行時に警告が表示されたのです。

:::note info
系列データの弱定常性を調べる方法として「系列に単位根があり、かつ系列が定常でない」という帰無仮説を検定する「**拡張 Dickey–Fuller (ADF) 検定**」があります。帰無仮説が棄却されることで、差分モデルで記述すると定常になることが確認されます。
:::

ところで、データが非定常であっても、その差分 (階差数列) は定常な場合があります。実際、今回のデータの二階差分は以下のように、増加トレンドが消えた定常なデータとなります。

```python:diff
from statsmodels.tsa.seasonal import STL

diff1 = df["co2"].diff().dropna()
diff2 = diff1    .diff().dropna()

# STL分解
stl = STL(diff2).fit()

# それぞれの成分を描画
fig, ax = plt.subplots(4, 1, figsize=(15, 12), sharex=True)

diff2.plot(ax=ax[0], c='black')
ax[0].set_title("(Data(t) - Data(t-1)) - (Data(t-1) - Data(t-2))")

stl.trend.plot(ax=ax[1], c='black')
ax[1].set_title("Trend")

stl.seasonal.plot(ax=ax[2], c='black')
ax[2].grid(True, linestyle='--', alpha=0.5)
ax[2].set_title("Seasonal")

stl.resid.plot(ax=ax[3], c='black')
ax[3].set_title("Residual")

plt.tight_layout()

plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/1bf8dba1-6158-79a3-b913-d1422a471383.png)

この差分データに対してARMAモデルを適用するのが、次に紹介するARIMAモデルとなります。

#### さらに詳しい内容を知りたい方へ
- Autoregressive Moving Average (ARMA): Sunspots data: 
<https://www.statsmodels.org/stable/examples/notebooks/generated/statespace_arma_0.html>


### 自己回帰和分移動平均モデル (ARIMA)
$$
\begin{aligned}
z_t &= \Delta ^d (y_t) = \sum ^d_{k=0} (-1)^k \ \_d \mathrm{C} \_k \ y_{t - k} \cr
z_t &= c + ε_t + \sum ^p_{i=1} φ_i z_{t-i} + \sum^q_{i=1} θ_i ε_{t - i}
\end{aligned}
$$

階差 $d$ の時系列を $p$ 次のAR、$q$ 次のMAモデルに適用したものです。

```python:ARIMA(10, 2, 6)
p, d, q = 10, 2, 6

arima = sm.tsa.SARIMAX(train["y"].values,
                       order=(p, d, q)
                       ).fit(maxiter=1000)

pred_val = arima.forecast(test_size)

pred_test = pd.Series(pred_val * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title(f"ARIMA{(p, d, q)}")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/f24d1d2d-39b9-a606-14f4-2edcbbd801cf.png)

エラーなく実行が完了しました。季節性周期変動だけでなく、長期的な増加トレンドも現れており、これまでのモデルと比較して精度が大きく上がっています。

:::note info
今回は $(p, d, q)$ を適当に決めましたが、より精度と汎化性能を向上させるためには **赤池情報量基準 (AIC)** などをもとに適切なハイパーパラメータを決定する必要があります。AICは、`summary`メソッドや`.aic`で確認できます。
:::

#### さらに詳しい内容を知りたい方へ
一階差分が定常となるようなARIMAモデル (単位根過程) を使用する際の理論的な話や注意点については、TJO氏のブログで詳細に議論されています。
- 見せかけの回帰について（そして単位根過程・共和分など）:
  <https://tjo.hatenablog.com/entry/2013/04/23/190417>
- Rで計量時系列分析：単位根過程、見せかけの回帰、共和分、ベクトル誤差修正モデル: 
  <https://tjo.hatenablog.com/entry/2013/08/16/095536>


### Seasonal ARIMAモデル (SARIMA)
$$
\begin{aligned}
y_t = c + ε_t + &\sum^p_{i=1}φ_i y_{t-i} + \sum^q_{i=1}θ_iε_{t - i} \cr + &\sum_{i=1}^{P_S} \phi_i y_{t-s i} + \sum_{i=1}^{Q_S} \eta_i \varepsilon_{t-s i}
\end{aligned}
$$

ARIMAモデルに季節成分を追加したモデルです。ARIMAモデルのパラメータ $(p, d, q)$ に4つのパラメータ $(P_S, D_S, Q_S, s)$ を加えた7つのハイパーパラメータを持ちます。このうち、$s$ が注目したい季節性変動の周期 (stride) であり、数式としては「$s$ 個ずつ飛ばしたデータを持ってきて線型結合」という形になっています。

```python:SARIMA(10, 2, 8) x (2, 0, 2)
p,  d,  q  = 10, 2, 8
ps, ds, qs =  2, 0, 2
s = 12 # 季節性変動の周期

sarima = sm.tsa.SARIMAX(train["y"].values,
                        order=(p, d, q),
                        seasonal_order=(ps, ds, qs, s)
                        ).fit(maxiter=1000)

pred_val = sarima.forecast(test_size)

pred_test = pd.Series(pred_val * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title("SARIMA Model")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/057fc5d2-fc4e-cdb2-95c1-8ac3bf4de1f9.png)


ARIMAモデルと比較して精度が上がっています。

#### さらに詳しい内容を知りたい方へ

- SARIMAX: Introduction - statsmodels 0.14.0: 
<https://www.statsmodels.org/stable/examples/notebooks/generated/statespace_sarimax_stata.html>



## ニューラルネットワークモデル

今回は「過去のデータから未来を予測する」ことを目的としているため、ニューラルネットワークモデルの概形は次の図のようになります。

![NN_architecture.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/3275c44d-2abd-8d8d-5dac-094205470c83.png)

これに合わせてデータの形式を変換しておきます。

```python:input_data
# 入力データの形状を変換
"""
[イメージ]
0, 1, 2,・・・ 11          +  正解データt=12
   1, 2,・・・ 11, 12      +  正解データt=13
      2,・・・ 11, 12, 13  +  正解データt=14
"""

p = 12

train_input   = np.array([df["y"][i : i + p] for i in range(0, train_size - p)])
train_input   = train_input.reshape(train_input.shape[0], p, -1)

train_correct = np.array([df["y"][i + p]     for i in range(0, train_size - p)])
train_correct = train_correct.reshape(train_correct.shape[0], -1)
```

### 再帰的ニューラルネットワーク (Recurrent Neural Network)
<img src=http://colah.github.io/posts/2015-08-Understanding-LSTMs/img/LSTM3-SimpleRNN.png>

(画像元: [Christopher Olah](http://colah.github.io/posts/2015-08-Understanding-LSTMs/))

隠れ層に再起的な構造を取り入れたニューラルネットワークです。系列データに対応したニューラルネットワークの中で最も基本的かつ素朴なモデルとなっています。RNN をもとに LSTM(後述) や Reservoir Computing、 Transformer といったモデルが開発されました。

```python:RNN
tf.keras.utils.set_random_seed(42)
tf.config.experimental.enable_op_determinism() # 演算の再現性を担保

n_hidden = 64
act_func = "relu"

rnn = keras.Sequential()
rnn.add(keras.layers.SimpleRNN(n_hidden,
                               activation=act_func,
                               input_shape=(p, 1)
                               )
        )

rnn.add(keras.layers.Dense(1))

rnn.compile(optimizer="adam", loss="mean_squared_error")

rnn_history = rnn.fit(train_input,
                      train_correct,
                      batch_size=4,
                      epochs=10
                      )

pred_val = np.zeros(test_size + 1)
snapshot = train_input[-1].copy()

for t in range(test_size + 1):
    step = rnn.predict(snapshot.reshape(1, -1, 1))
    
    pred_val[t] = step[0]
    
    snapshot[:-1] = snapshot[1:]
    snapshot[-1]  = step[0]

pred_test = pd.Series(pred_val[1:] * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title("RNN")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/21e95f0f-70bf-a446-1660-ff468837af48.png)

季節性変動っぽいものは出ましたが、トレンドは右下がりで、変動もしぼんでしまっています。隠れ層のニューロン数 `n_hidden` や活性化関数 `act_func` を調整するとうまくいくかもしれません。また、私の手元では入力のラグ `p` を倍の 24 にするとかなり良い結果が得られました。

:::note info

`rnn_history.history['loss']` から、学習の様子 (epoch が進むにつれて loss が低下していく様子) を確認することができます。

:::

#### さらに詳しい内容を知りたい方へ

- ニューラルネットワークで時系列データの予測を行う
<https://qiita.com/icoxfog417/items/2791ee878deee0d0fd9c>
- The Unreasonable Effectiveness of Recurrent Neural Networks: 
<http://karpathy.github.io/2015/05/21/rnn-effectiveness/>


### Long Short-Term Memory (LSTM)
<img src=http://colah.github.io/posts/2015-08-Understanding-LSTMs/img/LSTM3-chain.png>

(画像元: [Christopher Olah](http://colah.github.io/posts/2015-08-Understanding-LSTMs/))

RNN全体にわたって情報を横断させるためのメモリ(Cell State)を追加したモデルです。これによってRNNで指摘されていた勾配消失(爆発)問題が解決し、長い系列データを扱えるようになりました。強力なモデルですが、RNNと比較して収束が遅く、学習に時間がかかるのが欠点です。

```python:LSTM
tf.keras.utils.set_random_seed(42)
tf.config.experimental.enable_op_determinism() # 演算の再現性を担保

n_hidden = 64
act_func = "relu"

lstm = keras.Sequential()
lstm.add(keras.layers.LSTM(n_hidden,
                           activation=None,
                           input_shape=(p, 1)
                           )
        )

lstm.add(keras.layers.Dense(1))

lstm.compile(optimizer="adam", loss="mean_squared_error")

lstm_history = lstm.fit(train_input,
                        train_correct,
                        batch_size=4,
                        epochs=50 # LSTMは収束が遅いのでエポック数を増やす
                        )

pred_val = np.zeros(test_size + 1)
snapshot = train_input[-1].copy()

for t in range(test_size + 1):
    step = lstm.predict(snapshot.reshape(1, -1, 1))
    
    pred_val[t] = step[0]
    
    snapshot[:-1] = snapshot[1:]
    snapshot[-1] = step[0]

pred_test = pd.Series(pred_val[1:] * scale + min_val, 
                      index=test.index)

fig, ax, bx = my_plot_results(pred_test)

ax[0].set_title("LSTM")
bx[0].set_title("Root Mean Square Error")
plt.show()
```
> 実行結果：
> ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/684531/d67e0ca9-e0cf-202e-2176-e300c759340f.png)

完璧ではありませんが、増加トレンドと季節性変動が再現できました。RNNよりはマシな結果を出してくれましたね。

#### さらに詳しい内容を知りたい方へ
- 今度こそわかるぞRNN, LSTM編
<https://qiita.com/kazukiii/items/df809d6cd5d7d1f57be3>



# まとめ

古典的統計モデルには、主に以下の種類があります。いずれも、多くのハイパーパラメータを人力で調整する必要があります。

- AR(p), MA(q)
- ARMA(p, q)
- ARIMA(p, d, q)
- SARIMA(p, d, q) x (Ps, Ds, Qs, s)

ニューラルネットワークモデルは、RNNをもとに様々なモデルが提案されています。大規模かつ複雑なモデルゆえ「なぜその予測結果になるのか」を解析することは難しいですが、非線形なデータの予測には大きな力を発揮します。

以上です。ありがとうございました。

#### 参考URL
- Rで計量時系列分析：AR, MA, ARMA, ARIMAモデル, 予測:
  <https://tjo.hatenablog.com/entry/2013/07/12/184704>
- 時系列解析_理論編: 
  <https://logics-of-blue.com/%e6%99%82%e7%b3%bb%e5%88%97%e8%a7%a3%e6%9e%90_%e7%90%86%e8%ab%96%e7%b7%a8/>
- Time Series Forecasting with ARIMA , SARIMA and SARIMAX:
  <https://towardsdatascience.com/time-series-forecasting-with-arima-sarima-and-sarimax-ee61099e78f6>

- kerasで学習が再現できない人へ
<https://qiita.com/okotaku/items/8d682a11d8f2370684c9>

#### 参考文献

- Aileen Nielsen (著), 山崎 邦子 (翻訳), 山崎 康宏 (翻訳).
  「実践 時系列解析 ―統計と機械学習による予測」 オライリージャパン (2021).

- 沖本 竜義 (著).
  「経済・ファイナンスデータの軽量時系列分析」 朝倉書店 (2010).

- 白石 博 (著).
  「時系列データ解析」 森北出版 (2022).