はじめに

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

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

動作環境

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

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₂濃度のデータセットを使用します。
データ詳細: https://www.statsmodels.org/dev/datasets/generated/co2.html

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

データの読み込み

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

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()

実行結果：

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

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

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()

実行結果：

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

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

準備

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

準備の内容を表示

データの変換

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

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)

実行結果：

co2 y

YEAR

1965-01-31 319.400 0.036486

1965-02-28 320.450 0.055083

1965-03-31 320.925 0.063496

1965-04-30 322.000 0.082536

1965-05-31 322.060 0.083599

... ... ...

2001-08-31 369.425 0.922512

2001-09-30 367.880 0.895147

2001-10-31 368.050 0.898158

2001-11-30 369.375 0.921626

2001-12-31 371.020 0.950762

444 rows × 2 columns

	co2	y
1965-01-31	319.400	0.036486
1965-02-28	320.450	0.055083
1965-03-31	320.925	0.063496
1965-04-30	322.000	0.082536
1965-05-31	322.060	0.083599
...	...	...
2001-08-31	369.425	0.922512
2001-09-30	367.880	0.895147
2001-10-31	368.050	0.898158
2001-11-30	369.375	0.921626
2001-12-31	371.020	0.950762

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

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()

実行結果：

関数の定義

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

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)

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

機械学習モデルまとめ

古典的統計モデル

自己回帰モデル (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$ としています。

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.

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

チューニングパラメータ $φ_i$ は、予測したい時刻の $i$ ステップ前のデータが予測したい時刻のデータにどのような影響を与えているかを表しています。各係数の値や学習済みのモデルに関する情報は summary() メソッドで確認できます。

model_summary

display(ar.summary())

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$ としています。

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.

うまく予測できませんでした。また、実行すると警告 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)。

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

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

警告の意味 - 弱定常性

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

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

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

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()

実行結果：

この差分データに対して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モデルに適用したものです。

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()

実行結果：

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

今回は $(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$ が注目したい季節性変動の周期であり、数式としては「$s$ 個ずつ飛ばしたデータを持ってきて線型結合」という形になっています。

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()

実行結果：

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

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

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

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

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

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

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)

(画像元: Christopher Olah)

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

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()

実行結果：

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

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)

(画像元: Christopher Olah)

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

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()

実行結果：

完璧ではありませんが、増加トレンドと季節性変動が再現できました。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).

	co2	y
YEAR
1965-01-31	319.400	0.036486
1965-02-28	320.450	0.055083
1965-03-31	320.925	0.063496
1965-04-30	322.000	0.082536
1965-05-31	322.060	0.083599
...	...	...
2001-08-31	369.425	0.922512
2001-09-30	367.880	0.895147
2001-10-31	368.050	0.898158
2001-11-30	369.375	0.921626
2001-12-31	371.020	0.950762