Excelで時系列予測

初めに

事務作業の必需品ともいえるExcelですが、意外とその機能の多彩さは知られていないように思います。
今回はExcelを使った、 時系列データの見える化と時系列予測 をご紹介します。また、オマケとしてRというフリーの統計ソフト上でprophetという時系列予測ツールを使った結果と、pythonというプログラミングソフトでLSTMを実行してみた結果も、お見せします。
なお、本記事で使用するExcelは Excel2016 です。時系列予測機能は、2016以降に搭載された機能です。

動画

この記事に関する動画がこちらのアドレスにあります。

時系列データ

時系列データとは

時系列データとは、時間などに沿って順番に測定された一続きのデータのことです。ここでは一定の時間間隔で測定された数字のデータを見てみます。

東山動植物園の入場者数

名古屋市の東山動植物園は長く地元に愛されています。
この動植物園の入場者数はどのように推移しているのでしょう。近年は増えてきているのでしょうか？

東山動植物園の入場者数のデータは名古屋市が市のWebサイトでオープンデータとして公開しています。具体的には名古屋市観光客・宿泊客動向調査というページから閲覧できるPDF文書の「資料編」という箇所に掲載されています。ここの「施設別・月別入込客数」という表（例えば平成30年分では58ページにあります）の数値のうち、今回は平成19年から平成30年までのものを利用します。なお、平成29年から年度単位が暦年単位に変更になっており平成18年は1～3月がありません。このため、平成18年12月までのデータを削除しました。データ数としては144ヶ月分です。

お手元で以下の操作を再現される際にはこれらのデータをご取得頂くか、総務省統計局などから同様のデータを取得するなどしてご用意ください。その場合のおすすめは「住民基本台帳人口移動報告」の転入や転出の月次データです。

データの見える化

折れ線グラフ

まずは 見える化（可視化） です。
「年月」と「入場者数」がセットになった表と、年月を「年」と「月」に分けて「入場者数」とセットにした表の二つを用意します。
なお、年と月は例えば「2007年1月」がセルA2に入っているとすると

=year(A2)

で年が求まり

=month(A2)

で月が求まります。

最初にデータをグラフで確認します。「年月」と「入場者数」がセットの表のいずれかのセルを選択した状態で、「挿入」リボンから「折れ線グラフ」を選んでください。

このグラフでは、毎年春と秋にピークがあるのがわかります。こうしたデータは、 周期性（季節性） があるといいます。

ピボットグラフ

より一層の見える化のための集計と分析には、ピボットグラフが適しています。

「年」「月」「入場者数」がセットになった表のいずれかのセルを選択して、「挿入」リボンから「ピボットテーブル」を選んで「OK」としてください。

まずは各年の月ごとの変化を確認しましょう。
右の「ピボットテーブルのフィールド」という欄で

• 「年」を「列」に
• 「月」を「行」に
• 「転入者数」を「値」に

それぞれつまんで、ドロップします。そして「挿入」リボンから「折れ線グラフ」を選んでください。
これで各月の変化が年ごとに表示されます。
毎年のピークは5月と10・11月ですね。

考えてみるとこの時期は遠足のシーズンです。やはり動植物園は遠足の定番スポットで、動植物園は子供たちに支えられているということでしょうか。

では、入場者数は年々増えているのでしょうか？減っているのでしょうか？
「ピボットテーブルのフィールド」欄で

• 「年」を「行」に
• 「月」を「フィルター」に

変更すると年ごとの集計結果のグラフになります。
2011年から2012年にかけて落ち込みがありますが全体としては増加しているようです。
この落ち込みは東日本大震災の影響でしょうか。

ここで、グラフの線の上で右クリックして「データ系列の書式設定」から線を消してマーカーを見やすく変更し、さらにグラフの線の上で右クリックして「近似曲線の追加」をしてください。
上昇傾向が直線で分かりやすく表示されます。
こうした上昇や下降の傾向を トレンド と言います。

時系列予測

このままいくと、来年の入場者数はどう推移するでしょうか？
こうした予測を行うのが 時系列予測 です。

もとのデータがあるシートに戻ります。
「年月」と「入場者数」がセットになった表のどこかのセルが選ばれた状態で、「データ」リボンから「予測シート」を選びます。
「オプション」という項目をクリックすると設定できる項目が増えます。

• 「予測開始」を最終測定時点よりも前（今回は「2017/12/1」）
• 「予測終了」を最終測定時点（今回は「2018/12/1」）

として、「作成」としてください。

新しいシートが作成されて、予測がグラフとして表示されます。

元の転入者数データが青線で、予測がオレンジ色の線です。また、細いオレンジ色の線は、信頼区間 といって、実際のデータがこの線より外にはみ出ることはめったにないだろう、という予測の線です。

もう少し詳しくいいますと。信頼区間に「95%」という数字が入っていますが、これは同じ測定を100回行ったら95回まではその範囲に入るだろうということを意味しています。大雑把にいえば、「この範囲の外に出ることは5%以下で滅多にない」とこの予測シートは主張しているわけです。

つまり、実際のデータがこの範囲の幅の中なら予測は当たり、範囲外なら外れです。逆に言うと、信頼区間の幅が広ければどんな値が来ても「当たり」になるので、この区間幅は予測結果にどの程度の信頼性があるのかも示してくれます。

予測開始時点以後では、実際のデータと予測が一緒に表示されています。

比較してみてどうでしょうか。割と予測できていますね。

予測シートのG134に

=ABS(B134-C134)/ABS(B134)

と入力してそれをG145までコピーし、G147に

=MEDIAN(G134:G145)*100

と入力してください。
これは、予測値などの評価に使われる MER (Median Error Rate)という値で、数値は何％外しているかを示します。
約6.6％となりました。9割がたは合っているわけですから、良い予測と言えます。

このまま行くと翌年はどうなるでしょうか。「予測終了」を伸ばせば見ることができます。

さて、この通りになったかは、今年度の名古屋市の発表を待ちたいと思います。

なお、名古屋市に限らず全国の多くの自治体がPDFなどの文書の表としてしかデータを公開していません。
そのため、今回のように分析を行う場合に入力し直すなどの無用の労力が発生するうえに、自治体が自らの立場に基づいた分析を行い文書に取りまとめてから公表するという手順が挟まるので、公表までに少なくとも数か月、長ければ一年以上ものタイムラグがあります。
こうした状況はタイムリーで円滑な情報利用を妨げており、「国民参加・官民協働の推進を通じた諸課題の解決」など目的に官民データ活用推進基本法が目指すオープンデータのあり方と合致しないのではないでしょうか。「データが整い次第まずはCSV形式で公表。その後、行政として結果を分析し発表する」といった運用が一般的になることを希望します。

まとめ

今回は、Excel2016を使った時系列データの見える化と時系列予測を行いました。
見える化にはピボットグラフを用いるのが早道でした。
また、時系列予測はExcel2016以降の予測シートを用いると簡単に出来ました。

ただ、時系列予測はあくまでも「このままの状態が続くこと」を前提としています。東山動植物園も新型コロナの影響で一時閉園していましたし、その後も入場制限が行われたりしていて、入場者数は大きく減少したはずです。時系列予測ではこうしたことまでは予測はできないのです。

オマケ

Rでprophet

R で prophet という時系列予測ツールを利用して予測を行ったコードと結果のグラフを載せます。（Rやprophetについての詳細は他のサイトに譲ります）

以下のコードを実行しました。

library(prophet) # prophetのライブラリ（インストールしておく）

# --- ファイル読み込み
df <- read.csv('***.csv')  # csvデータ読み込み
# 注：形式は
# 1行目: 年月,入場者数
# 2行目: 2007/1/1,104917

df$年月 <- as.Date(df$年月) # 日付型データに変換する
colnames(df) <- c("ds", "y")
# 注：prophetは列名が指定されているのでそれに合わせる

# --- prophet実行
m = prophet(df[1:132,])  
# 注：2017/12までのデータで学習（予測開始：2018年1月）

f = make_future_dataframe(m, periods = 12, freq = "month")
# 注：学習結果を用いて12ヶ月先まで予測（予測終了：2018年12月）

pred = predict(m, f)       # 予測計算を実行


# --- グラフ用のデータにする
# 実際のデータ
y <- ts(df$y, start = c(2007,1,1), frequency = 12)
# 予測
yhat <- ts(pred$yhat, start = c(2007,1,1), frequency = 12)
# 信頼区間
lcl <- ts(pred$yhat_lower, start = c(2007,1,1), frequency = 12)
ucl <- ts(pred$yhat_upper, start = c(2007,1,1), frequency = 12)

MER <- median( abs(y[133:144] - yhat[133:144])/abs(y[133:144]) )* 100. # MER計算

# --- グラフ化
# 実際のデータ
plot(y, col = "blue", yaxt = 'n', ylim = c(min(y)*0.8, max(y)*1.2), lwd = 2,
     xlim = c(2017.75,2019), xaxt = 'n')
axis(side = 1, at = 2018:2019)
# 予測
lines( window(yhat, start = c(2018,1)),
       type = "l", col = "orange1", lwd = 2, ylim = c(min(y)*0.8, max(y)*1.2))
# 信頼区間
lines(window(lcl, start = c(2018,1)),
      type = "l", col = "orange1", lwd = 1, ylim = c(min(y)*0.8, max(y)*1.2))
lines(window(ucl, start = c(2018,1)),
      type = "l", col = "orange1", lwd = 1, ylim = c(min(y)*0.8, max(y)*1.2))

# 2018年1月に区切り線
abline(v = 2018.0, col = "red")

# 【scientific（乗数表示）ではないｙ軸ラベルを書く】
ys <- par()$yaxp
y_seq <- seq(ys[1],ys[2],by = (ys[2]-ys[1])/ys[3])
y_seq <- as.integer(y_seq)
axis(2, y_seq, labels = format(y_seq, scientific = F))

# タイトルなど
title("Predicted by prophet")
mtext( paste("MER =", sprintf("%4.1f",MER), "%"), side = 3, line = 0, adj = 1 )

グラフ化周りのコードが多くなっていますが、prophetの実行と予測自体は3行です。

予測結果と実際のデータの比較箇所のグラフがこちら（オレンジ色の太い線が予測結果で、細い線が信頼区間です）。

信頼区間の幅が狭く、予測の信頼性は高いことが伺えます。
ただし、外している箇所もいくつか見られます。なによりMERが良くありません。

prophetがExcelの予測シートの結果を下回ることはしばしば起こります。

しかしこれはprophetの性能がExcelより劣っている、ということではありません。
prophetは、曜日や天候、イベントなどの外的要因を考慮した予測を行うことを目的として作られています。動植物園のデータでいうなら、幼稚園や小学校の行事日程や教員による児童の熱中症などへの配慮、人気の動物の体調などを考慮するべきでしょう。こうしたデータも学習して、その条件を与えられたときの入場者数を予想することを得意とするツールなのです。

つまり、「多少バラツキがある」「周期性が強い」「外部からの影響は考えない」という条件だけでは、手がかりが少ないので、複雑なことは考えずにシンプルに予測したほうが、正解率が高くなるというわけです。
こうした単純な予測では、prophetなどの高度な予測法の出番はあまりないかもしれません。

pythonでLSTM

ついでですので、 python の pytorch という深層学習用のシステムを使って、 LSTM による予測も行ってみました。
pythonやpytorch, LSTMについての解説も他のサイトに譲ります。

なお、 jupyter-Notebook で作業していますので、コードはセルごとに載せます。
バージョンはPytorchが1.5.0で、pythonは3.8.3です。

まずはライブラリ読み込み。
（datetimeはデフォルトで入っており、numpyはtorchインストール時に入りますが、pandas, matplotlib, scikit-learnはインストールしておく必要があります）

import numpy as np                # 数値演算用ライブラリ
import pandas as pd               # データフレーム用ライブラリ
import matplotlib.pyplot as plt   # グラフなどの画像表示ライブラリ
import datetime                   # 日付時刻用ライブラリ
from sklearn.preprocessing import MinMaxScaler # 正規化用関数
import random                     # 乱数ライブラリ

import torch                      # Pytorch本体
from torch import nn             # 深層学習ネットワーク構築用ライブラリ
import torch.optim as optim      # 深層学習用最適化ライブラリ

次にデータの読み込みと加工。

# データ準備
df = pd.read_csv('***.csv', encoding='cp932') # csvデータ読み込み（Rで使用したものと同じ）
nd = df.iloc[:,1].values                      # 数値列取り出し

scaler = MinMaxScaler(feature_range=(-1, 1))              # 正規化用スケーラ
nd_normalized = scaler.fit_transform(nd.reshape(-1, 1))   # 正規化実施
tc_normalized = torch.FloatTensor(nd_normalized).view(-1) # torchテンソル化

# 予測用時系列データ作成関数
def create_inout_sequences(input_data, windowlenght):
    out_set = []                                   # 結果の保存先
    L = len(input_data)                            # データ数
    for i in range(L-windowlenght):
        data_X = input_data[i:i+windowlenght]      # シーケンス長の分のデータをX
        data_y = input_data[i+windowlenght]        # その次をｙとする
        out_set.append((data_X ,data_y))           # Xとyをセットにして保存
    return out_set

# データ作成
seqlenght = 36                                    # 36ヶ月分で次の1月を予測
tc_setdata = create_inout_sequences(tc_normalized, seqlenght) 

# train/test分割
p = len(tc_setdata) - 12 # ラスト12ヶ月をtest用として分離する
train_set = tc_setdata[:p]
test_set = tc_setdata[p:]

そしてLSTMの本体。

# LSTMのモデル定義
class LSTM(nn.Module):
    def __init__(self, input_size=1, hidden_layer_size=64, output_size=1):
        super().__init__()
        self.hidden_layer_size = hidden_layer_size
        self.lstm = nn.LSTM(input_size, hidden_layer_size)              # LSTM層
        self.linear = nn.Linear(hidden_layer_size, output_size)         # 全結合層
        self.hidden_cell = (torch.zeros(1,1,self.hidden_layer_size),
                            torch.zeros(1,1,self.hidden_layer_size))    # メモリセルを設定

    def forward(self, input_seq):
        lstm_out, self.hidden_cell = self.lstm(input_seq.view(len(input_seq) ,1, -1), self.hidden_cell) # LSTM層
        predictions = self.linear(lstm_out.view(len(input_seq), -1))    # 全結合層
        return predictions[-1]

そして学習実行です。

# モデル実装
model = LSTM()
loss_function = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

# エポック数の指定
epochs = 100

# 保存先準備
val_result_list = []

# 学習実行
for i in range(epochs):
    train_loss = 0
    val_loss = 0

    # 【学習モード】
    model.train()
    trains = random.sample(train_set, len(train_set)) # シャッフル
    for X, y in train_set:
        optimizer.zero_grad()             # 勾配の初期化
        model.hidden_cell = (torch.zeros(1, 1, model.hidden_layer_size),
                        torch.zeros(1, 1, model.hidden_layer_size)) # メモリセルの初期化

        y_pred = model(X)[-1]             # LSTM計算

        loss = loss_function(y_pred, y)   # 評価関数計算
        train_loss += loss.item()         # 評価値を累計
        loss.backward()                   # 誤差逆伝播計算
        optimizer.step()                  # 最適化実行

    avg_train_loss = train_loss / len(train_set) # エポック毎の評価値の平均

    # 【評価モード】
    model.eval()
    L = len(test_set)
    val_ypred = np.zeros(L)    # 結果の保存先の準備
    k = 0
    with torch.no_grad():
        for X, y in test_set:
            y_pred = model(X)[-1]           # LSTMで予測値を計算
            loss = loss_function(y_pred, y) # 評価関数計算
            val_loss += loss.item()         # 評価値を累計
            val_ypred[k] = y_pred           # 予測値を保存
            k += 1
        avg_val_loss = val_loss / len(test_set) # エポック毎の評価値の平均

    if i%2 == 0:
        # ＜経過を表示＞
        dt = datetime.datetime.now().time()
        print(' {}: Epoch [{}/{}], Loss: {:.4f}, Val_loss: {:.4f}'.format(dt, i+1,epochs, avg_train_loss, avg_val_loss))

    # ＜数値を保存先に追加＞
    val_result_list.append(val_ypred)

print('{}: Epoch: {} Val_loss: {:.5f}'.format(dt,i+1,avg_val_loss)) # 最終の評価値の表示

あとはグラフ化です。

# エポックの指定と、そのエポックの結果の取り出し
epoch = 50 # 見たいエポックを指定
t = epoch -1
pred = scaler.inverse_transform(val_result_list[t].reshape(-1,1)) # 結果の取り出し（正規化したデータなので元のスケールに戻す）

# 表示範囲
x0 = len(tc_normalized)-len(test_set)
x1 = len(tc_normalized)
x = np.arange(x0, x1, 1)

# MER計算
MER = np.median( np.abs( (nd[x] - pred.reshape(1,-1)[0])/nd[x]) )* 100.

# グラフ化
fig = plt.figure()
ax =  fig.add_subplot(1,1,1)
ax.set_title('Real vs Predicted')          # タイトル
ax.set_ylabel('Value')                     # Y軸ラベル
ax.grid(True)                              # グラフにグリッド追加
ax.plot(nd, label='real')                  # 実際のデータをプロット
ax.plot(x,pred, label = 'predict')         # 予測結果をプロット
ax.set_xlim([x0,x1-1])                     # X軸表示範囲
ax.text(0.1, 0.9, 'epoch = {}'.format(epoch), transform=ax.transAxes)  # エポック番号書き込み
ax.text(0.1, 0.8, 'MER = {:.1f}%'.format(MER), transform=ax.transAxes)  # MER書き込み
plt.legend(loc = 'upper right')            # 凡例
plt.show()

そして、得られた中で一番成績が良かった予測結果はこうでした（predictが予測結果です）。

確かに近いのですが、苦労したわりにExcelとMERがあまり違いません。
これは、データ数が少ないことも原因ですが、やはりこうした手がかりの少ないデータではLSTMなどの高度な予測法は本領を発揮出来ないのだと思われます。

また、信頼区間がありません。
基本的にAIは信頼区間などの分かりやすい目安を計算してくれません。そこで、AIの出来はMERなどの値だけで判断することになります。
このため、テストの点数だけで人を判断できないのと同じで、大体は良さそうでも時々とんでもない間違いをする、というAIが出来上がることもあります。

＜以上＞