時系列データへのクロスバリデーション法を用いて、予測精度の評価を行う方法を説明します。
R言語のforecastパッケージのtsCV関数を用いると、効率的な短いコードで実装が可能です。

この記事では、時系列データの前処理~モデル化~予測~予測の評価、といった一連の流れをすべて通して解説します。
今回は標準的な時系列モデルであるSARIMAモデルを用いますが、このモデル以外でもおおよその手順は変わらないと思います。

予測の評価における基本的な事項は『予測の評価方法:誤差の指標とナイーブな予測』も合わせて参照してください。

コードはGitHubから参照することができます。



スポンサードリンク

目次

  1. 分析の準備と前処理
  2. SARIMAモデルの構築
  3. SARIMAによる予測
  4. テストデータを使った予測の評価
  5. クロスバリデーション法による予測の評価
  6. スライド型のクロスバリデーション法の実行
  7. 後記

 

1.分析の準備と前処理

この記事では、時系列データの前処理~モデル化~予測~予測の評価、といった一連の流れをすべて通して解説します。
逆に言えば、時系列データを予測するというだけでも、まじめにやればこれだけの手順を踏まなければいけないということですね。

この記事の方針は以下の通りです。
1.前処理としては、データの対数変換を行います。
2.モデル化としては、SARIMAモデルを構築します。
3.予測は、forecastパッケージの関数を使えば素直なコードで実行できます。
4.予測の評価としては以下の3つを行います
4-1.訓練データとテストデータに分けて、テストデータでの予測精度を評価する
4-2.訓練データを1時点ずつ増やしていき、予測を何度も行い、その時の予測精度を評価する
4-3.訓練データをスライドさせていき、予測を何度も行い、その時の予測精度を評価する

 
まずは、分析に必要なパッケージを読み込みます。
必要であれば1~3行目のコメントを外して、パッケージをインストールしてください。
分析に使うのは主にforecastパッケージです。グラフを描くために残りのパッケージを読み込んでおきました。

 
続いて、データの読み込みと対数変換です。
今回は飛行機乗客数のデータを使いました。月単位のデータです。R言語における標準的な時系列データのクラスであるts型で保存されています。

 
データを読み込んだら、最初にやることは図示です。
以下の1行のコードで美麗なグラフが描けます。ACFは自己相関でPACFは偏自己相関と呼ばれる指標です。例えば『時系列解析_理論編』等を参照してください。

 
実際に数値を眺めて、どんなデータなのか確認をしておきます。

 
最後に、訓練データとテストデータに分割しておきます。
ts型のデータの場合は、window関数の引数にstartやendと指定することで簡単にデータの切り出しが可能です。

 

2.SARIMAモデルの構築

SARIMAモデルを構築します。
forecastパッケージのauto.arima関数を使うことで、ほぼ自動的に次数の選択やモデルの構築が完了します。
auto.arima関数の使い方に関しては『ARIMAモデルによる株価の予測』等も参照してください。

推定結果を出力します。
ARIMA(0,1,1)(0,1,1)が選択されたようです。

 
残差のチェックを行います。
checkresiduals関数を使うと、残差の自己相関の検定と残差の図示を同時に行ってくれます。

結果を見ると『p-value = 0.5618』となっているので有意な自己相関は見られず、残差のグラフを見ても大きな問題は見当たらなかったのでそのまま進めます。



スポンサードリンク

 

3.SARIMAによる予測

推定されたモデルを用いて、予測を行います。
予測はforecastパッケージのforecast関数を使います。
予測結果を『f_sarima』という変数に格納したうえで、図示を行います。
コードと結果をまとめて載せます。

 

4.テストデータを使った予測の評価

予測の評価はaccuracy関数を使うことで簡単に実行できます。
accuracy関数に関しては『予測の評価方法:誤差の指標とナイーブな予測』も合わせて参照してください

例えばRMSEの値を見ると、訓練データでは0.03539809、テストデータでは0.09593236となりました。
しかし、0.09といった数値を見て「なるほど!」と思う方はあまりいないのではないかと思います。
RMSEは誤差の大きさですので、小さければ小さい方が良いのですが、どれくらい小さければ十分なのでしょうか。目安がほしいところですね。

そこで「複雑な技術を使わなくても行うことができる予測」すなわち「ナイーブ予測」を用います。
ナイーブ予測のRMSEと、頑張って構築したSARIMAモデルのRMSEを比較して、SARIMAモデルの方が小さければ、十分に良い予測ができていると判断できます。

ナイーブ予測の種類に関してはやはり『予測の評価方法:誤差の指標とナイーブな予測』も参照していただければと思いますが、今回は以下の3つのタイプのナイーブ予測を用いることにします。

1.訓練データの最終時点を予測値とする
2.1にたいして、全体のトレンド(今回のデータでは右肩上がり)を組み込む
3.1周期前(今回は12時点前)を予測値とする

1.は平均に回帰しない、非定常な時系列データに対して用いられるナイーブ予測です。
2.は、さらに右肩上がりのトレンドを組みこみました。
3.は、周期性を持つデータに対してしばしば用いられる方法で、前年同期の値を予測値としてそのまま使っているわけです。

こういった単純な予測に負けてしまうようでは情けないですね。こいつらに勝てるかどうかを評価しようということです。
ナイーブ予測はforecastパッケージを使うと簡単に計算できます。

ナイーブ予測の結果も図示してみます。
autoplotの結果を各々p1~p3に格納して、3つまとめて表示します。

このグラフを見ると、3つ目の季節付きのナイーブ予測はうまく予測ができているように見えます。
単なる前年同期を使っただけなのですが、このような単純な予測でも十分であることはしばしばあります。
こいつに勝てるかどうかが1つのポイントですね。

全部のRMSEを比較します。

訓練データでもテストデータで見ても、共にSARIMAモデルのRMSEが最小となっていることがわかりました。
まずは、うまく予測できてよかったね、というところでしょうか。

 
最後におまけとして、予測誤差のグラフを描いてみます。
数値で直接見てもなかなかわかりにくいですが、棒グラフなどにすると視覚的に理解できます。

ggplot2というパッケージを使うことで美麗なグラフを描くことができます。
そのためには、まずはデータフレームの形式になるように結果を整形する必要があります。

訓練データとテストデータのRMSEを並べて比較したいので、まずは凡例を用意しておきましょう。

ここからの方針ですが、RMSEの計算結果を以下のように整形すると、ggplot2を使ってきれいなグラフを描くことができます。

このコードをあと3回繰り返して頑張って整形しても良いのですが、面倒なので、簡単に整形できる便利関数を作ってみました。
この記事の中でしか使いませんので、少々適当に作ってありますが、最低限の動作はしてくれます。

この関数は、名称と予測誤差のタイプ(今回は訓練データorテストデータ)とRMSEを渡すことで、整形されたデータフレームを出力してくれます。

動作を確認してみます。

うまくできていそうなので、これを使って4種類の予測のRMSEを一気に整形して、rbind関数を使って4種類の予測の結果を結合させます。

こんな結果になります。
きれいになりましたね。こういう形式を整然データと呼びます。整然データに関しては『整然データとは何か|Colorless Green Ideas』等も参照してください。

後はこの結果を使ってグラフを描きます。ggplot2関数を使います。

これなら結果が一目ですね。
SARIMAモデルのRMSEが最も小さいようです。

今回は標準的なやり方でデータを整形しましたが、dplyer等を使うともう少し短いコードで済むかもしれません。
もっといいやり方があれば、ご指摘いただけますと幸いです。

 

5.クロスバリデーション法による予測の評価

先ほどはデータを2分割して予測誤差を評価しました。
しかし、これにはいくつか不満な点があります。

1.予測を1回しか行っていない→複数回の予測を行い、予測誤差を正確に評価したい
2.何時点先の予測が当たりやすいor当たりにくいかがわからない→1時点先の予測誤差、2時点先の予測誤差と別々に求めたい

2番は地味に重要でして、例えば「2時点先(今回のデータでは2か月先)までは精度よく予測できるが、3時点先になると予測誤差が大きくなる」みたいなことがわかれば、予測の運用において大きな情報となります。予測を用いた意思決定をしたいならば、なるべく2時点先までの結果だけを使ってね、と運用側に注意を呼び掛けるわけです。

こういうのを簡単にやっつけてしまおうということで使うのが、forecastパッケージのtsCV関数です。

tsCV関数での予測の評価には大きく2つのやり方があります。
1つ目は訓練データの「開始時点」を1時点目に固定して、そこから1つずつ訓練データを増やしていくというものです。
■を訓練データ
□をテストデータ
とします。
また、2時点先までを予測して、その評価を行うとしましょう。

以下のように訓練データ、テストデータが用いられます。
1回ごとに、訓練データが増えていくことに注意してください。
データが少ない場合などにしばしば使われます。
■□□
■■□□
■■■□□
■■■■□□
■■■■■□□

実装してみましょう。
tsCV関数の引数には、予測結果を出力する関数を入れます。
そこで、訓練データを入れると、SARIMAモデルを推定して、予測結果を出力してくれる関数を自作するところから始めましょう。

まずは、最適な次数がいくつだったかを確認しておきます。

ARIMA(0,1,1)(0,1,1)でしたね。
この次数を使って、予測結果を出力してくれる関数を自作します。

動作確認をします。

 
うまくできているようなので、あとはこれをtsCV関数の引数に入れるだけです。
SARIMAモデルだけではなく、ナイーブな予測の予測の評価もあわせて行いました。
データが小さいので、数秒で終わります。

結果を確認してみます。SARIMAモデルの評価結果を出力します。

注目すべきはNAが大量にある点です(中略で消したところにNAはありません)。
SARIMAモデルは季節性を組み込んだモデルですので、最低でも1年以上の訓練データがなければモデルを構築することができません。そのため、最初の1年と少しはすべてNAです。
また、3時点先までを予測しようとしているので、3時点先の回答データがない末尾においてもやはりNAが出力されます。

出力の説明をしておきます。
『Mar 1950』の行は「Mar 1950までのすべてのデータを使ってモデルを構築した結果」が出力されています。
『h=1』列目は1時点先の予測残差(実測値-予測値)が出力されます。
『Mar 1950』行の『h=1』列目は、1950年4月の予測残差となり、
『Mar 1950』行の『h=2』列目は、1950年5月の予測残差となります。

 
次に移る前に、訓練データの大きさと予測誤差の対応を図示して確認してみます。
1955年以降は、安定してそうに見えます。

 
次に行きましょう。
やることは、ナイーブな予測とのRMSEの比較です。
NAが大量にあるので、まずはこれを削除します。na.omit関数を使います。
また、テスト期間は1959年以降ですので、この結果だけを使うことにします。

 
データの抽出が終わったので、あとはRMSEを計算するだけなんですが、残差からRMSEを計算する関数がforecastパッケージには用意されていないようだったので、自作します。
コメントが長いんですが、実質1行で終わる簡単な関数です。

この関数を使って、SARIMAモデルにおけるRMSEを計算します。
apply関数の2つ目の引数に「2」を、3つ目の引数にcalc_rmseを入れると、列を対象にしてcalc_rmseを一気に実行した結果を出力してくれます。

前回と同様に、整形した後に、ggplot2を用いて図示します。

やはりSARIMAモデルのRMSEが最も小さいようです。

続いて、何時点先までの予測であれば、精度よく予測できるのか、確かめてみましょう。
予測期間をh=12として長くしたうえで、tsCV関数を実行します。

1年先までを予測対象としているので、最終年は評価結果に表れないことに注意してください。
予測期間を変えると、評価対象となるデータまで変わってしまうのが難しいところです。

前回同様、データを整形したうえで図示します。
凡例が12パターンもあって作るのが大変なので、ちょっと工夫して作ります。
formatC関数で『width=2, flag=”0″』と指定すると、頭を0埋めにしてくれます。グラフの並び順をきれいにする(h=1の隣にh=10が来ないようにする)ために指定しました。

後はほぼ同じ方針で図示できます。
SARIMAモデルしか対象になっていないので、『fill = type』と指定するのを無くし、色分けをやめました。

グラフで見ると、2時点先までは精度が良くて、3時点先からはちょっとRMSEが大きくなりそうだなということがわかります。

 

6.スライド型のクロスバリデーション法の実行

最後に、スライド型のクロスバリデーション法を実行します。
訓練データの大きさをwindowで指定し、これをスライドさせることで、一定の大きさの訓練データを常に使います。

簡単な概略図を示します。
■:訓練データ
□:テストデータ
×:使わないデータ
window=3, h=2の時
■■■□□
×■■■□□
××■■■□□

tsCV関数に『window』という引数を入れることを除けば、ほとんど同じコードで実装できます。
予測残差を求めたうえで、NAを排除するところまでまとめてやってしまいます。

RMSEを求めます。

せっかくですので、2つのクロスバリデーション法の比較をしてみます。

グラフを見ると、大して変化はなさそうです。

 

7.後記

後記というよりかは、補足であったりいいわけであったりをここに記します。

tsCV関数に関しては、forecastパッケージの作成者の以下の記事も参照してください。
Cross-validation for time series

時系列データへのクロスバリデーション法を説明したわけですが、時系列データだからといって常にこの方法でなければならないというわけではありません。単純な自己回帰モデルで表現可能な場合は、通常のk-fold-CVを使っても支障ないようです。
A note on the validity of cross-validation for evaluating autoregressive time series predictionという論文も参照してください。

今回はtsCV関数を使ったわけですが、こういった関数を自作することも難しくはありません。tsCV関数は指定できる箇所がやや少ないです(例えば固定型のCVにおいて、データを2個ずつ増やしたい、と思っても、この関数では実行できない)。この記事の最後に参考文献としてあげた「前処理大全」などにtsCVを使わない方法が載っています。

後は、書いていて思ったことを雑に記します。
今回はやや長いコードとなりました。長いコードをRで書く場合はセクションを区切ると便利です。RStudioを使っている場合は「Ctrl+Shift+R」でセクションを挿入することができます。

関数の作成をする箇所も3カ所ありましたね。本来は関数を書くのは単一のセクションにまとめてしまうと見通しが良くなります。私の場合は実際の処理を行う前に「utility関数」みたいなセクションを作ってその中に突っ込むことがしばしばあります。今回は説明のためにあっちこっちで関数を定義しましたが、できるだけまとめましょう。別のファイルに関数を定義してsource関数を使って読み込むこともあります。

関数には、なるべくコメントを残しておいた方が、後で読みやすいです。calc_rmseみたいに、書くのが馬鹿らしくなるような短い関数もあります。しかし、引数が残差だったのか予測値と実測値両方入れるんだったか、忘れちゃった、みたいな問題は、コメントを書くことで防ぐことができます。この記事のコメントの書き方が最善であるとは考えにくいので、色々なサイトの書き方を参照してみてください。

命名規則はやや適当なので、もう少し考えてコードを書いた方が良かったですね。

ggplot2はbiopapyrusさんの記事を参考にさせていただきました『geom_bar|ggplot で棒グラフを描く方法』。

0埋めに関しては、『数値の頭に0を詰めて桁を揃える|My Life as a Mock Quant』を参照させてもらいました。

最後に、一瞬だけ出てきた整然データですが、この用語は小学校の「情報」の教科書に出てきてもいいくらいのものだと思います。宣伝のために、もう一度リンクを載せておきます『整然データとは何か|Colorless Green Ideas』。
論文も出ているようです。引用しやすくなったので助かります。
整然データとは何か

 

参考文献


時系列分析と状態空間モデルの基礎:RとStanで学ぶ理論と実装

 
管理人が書いた時系列分析の入門書です。SARIMAモデルの仕組みを知りたいという方や、他の時系列モデルについて学びたいという方は、こちらも読んでいただけると幸いです。
サポートページはこちらです
重版に合わせて、出版社サイトにて第1部の1~2章が立ち読みできるようになりました。
 

前処理大全

 
通称Awesome本と呼ばれる書籍です。Awesomeとは「素晴らしい」とか「驚くばかり」という意味らしいですね。
タイトルの通り、データの効率的な前処理の方法が記されていて、便利です。
 



スポンサードリンク

 
更新履歴
2018年4月24日:新規作成
2018年4月25日:GitHubへのリンクを追加