じっくり学ぶ時系列解析~準備編~

時系列データとは、ある一時点ではなく、時間軸に沿ってデータが収集され、変数の並び方にも意味があるデータのことを指します。商品の売り上げを月ごとに集積したデータや、企業の株価データ、気温データ、電車の乗客数などが時系列データの例としてあげられますね。

しかし、世の中の多くのデータが時系列データとして集計されているにも関わらず、これを解析するには通常の線形回帰モデルでは火力不足で、正しい解釈を導き出すことが困難です。

そこで、今回から複数回に分けてそんな時系列データ解析の理論やその実践法についてまとめていきます。

記事中の内容は、Rで再現できるようにコードを記載していくすので、参考にしてみてください。（近々pythonでも再現できるようにいたします。）

時系列分析の目的

時系列解析の目的は、実は普通の線形回帰モデルと変わりません。例えば、アイスクリームの売り上げデータを解析するとします。アイスの売り上げや気温も時系列データですが、一つの例として、気温を説明変数に、売り上げを目的変数にした単回帰分析をしてみましょう。

すると

という統計モデルができるはずです。このモデルを使えば、気温が一度上がった場合にアイスがどれくらい売れるのか予測することもできますし、決定係数を使えば気温がどれぐらい売り上げに貢献しているのか推測することもできます。

これと同じように株価や乗客数などの時系列データに対して統計モデルを作成し、データ変動に与える要因の推測や将来の予測を行うこと、これが時系列解析の目的になります。

時系列データの必要性

さて、何の代わり映えもない単回帰分析でもアイスの売り上げデータの統計モデルを作成できてしまいましたが、これでは不十分です。

時系列は並べ順にも意味があるデータです。つまり、他の変数だけでなく過去の自分や過去の他の値からも影響を受けることがあります。

アイスの売り上げも過去にほとんど売れていない店舗が、気温が上がったからといって、急に売れるなんてことは考えにくいですよね。おそらく気温がどうであれ、今までの売り上げも現在の売り上げに影響を与えているはずです。

他にも秋のコメの収穫量は夏の降水量に影響を受けることもありますよね。

このように、時系列データはその並びも考慮して解析する必要があるので、通常の回帰モデルでは火力不足になってしまうのです。

また、単回帰モデルでは、未来の目的変数を予測するのに未来の説明変数の値を知らねばならないというのもこのモデルの限界です。

それではどんな統計モデルを作れば良いのでしょうか？

いにしえの時系列モデル

まずは解析のイメージを掴むため、いにしえから伝わる時系列解析手法を紹介していきましょう。

トレンドモデル

トレンドモデルは1番シンプルなモデルです。このモデルは目的変数が時間とともにある関数系に従って変化することを前提とします。

例えばイケイケ企業の株価を予測するために

とモデルを作ったとしましょう。

このモデルは、多少の誤差はあるものの企業がイケてるので、あるt時点の株価は時間が経つほど高くなることを意味します。 つまり、10年先の株価は

となり、20年先の株価は

となるということです。

また、関数は対数型にしたり指数系にしたり、前提とする仮定によって様々に設定できます。

しかし、トレンドモデルでは解析者が独自に関数を設定するため、恣意性が強く、モデルが時間のみに依存するため現実では考えにくい仮定が必要になるというデメリットがあります。

平滑化モデル

他のモデルを考えてみましょう。今度は時間ではなく過去の自分から影響を受けること仮定したモデルです。

これは

と表し、あるt時点の値は平均値にその時の誤差を足した値であることを意味します。

ここで、平均値が時間ともに変化するとして、自分の過去と現在の値から算出される移動平均値を使うモデルを平滑化モデルと呼びます。

例えば、1~5期までの株価がわかっているとして、その移動平均値を使えば6期目の株価の予測値は、

となります。

このように自分の過去の値からのみ影響を受けているモデルを一変量時系列モデルといい、他の過去の値（お米の収穫量に対する過去の降水量など）からも影響を受けるモデルを多変量時系列モデルといいます。

まずは一変量時系列モデルの理解から深めていきましょう。

時系列モデルのその前に

これら二つのモデルは現在ではほとんど用いられていません。しかし、根底となる考え方は脈々と現代の時系列解析においても受け継がれているのです。

ここからは現代でも用いられる解析手法を紹介していきましょう。

とその前に、今回は準備編ということで、今後の理解を深めるために時系列解析の考え方とちょっと特殊な用語を紹介しましょう。

確率過程（重要）

まずは確率過程です。実は時系列解析をするにあって、この確率過程の意味を理解するとはとても重要なのです。

例えば、単回帰モデルでは"未来の目的変数を予測するのに未来の説明変数の値を知らねばならない"という限界があることを書きました。

"そうはいってもデータの未来を予測したい...そうだ、タイムマシンを作ろう”と思うところですが、そうはいきません。

そこで、時系列のデータはある確率分布を元に、時間にそって生成されていると仮定します。これを確率過程と呼びます。（細かい定義は間違っている可能性がありますので、詳しくはこちらを...）

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

例えば、毎分1と0がランダムに生成される機械の値を10分間観察しました。すると

という値が出たとします。これはどのような確率分布を背景にした値でしょうか？1と0 がそれぞれ5回ずつでているので、おそらく確率が0.5のベルヌーイ分布から生成されているだろうと推測できます。

www.medi-08-data-06.work

これが確率過程です。このように一見ランダムに見える時系列データでも、ある確率分布から生成されていると仮定して、いかにもっともらしい確率過程をモデリングするかということが時系列解析のきもになり、この確率過程が統計モデルとなるのです。

具体的にいきましょう。ある店舗ではアイスクリームの売り上げ毎月記録しています。まずは過去2年間分(2023~2025)のデータを見てみてください。

ここからわかることは、一年を通して売り上げが10万ぐらい変化してはいるものの、前の月と次の月との違いは５万円前後であることが分かります。

ここで、来月の売り上げを予想してみます。タイムマシンに乗って未来を見に行く訳には行かないので、過去のデータからある仮説を立てます。

それは、次の月の売り上げは今月の売り上げにプラスマイナス５万円ぐらいの誤差を足したものに違いない、というものです。さらにこの誤差はプラス方向へもマイナス方向へも同じ確率で変動しているっぽいので、正規分布に従うだろう、と仮定します。

www.medi-08-data-06.work

すると、来月の売り上げは、

という平均値が今月の売り上げ、標準偏差が５万円の正規分布に従うことが予測されるので、

ぐらいだろうと推測できます。

これを一般化すると

となります。つまり、t時点のアイスの売り上げは、平均値がt-1時点の売り上げ、標準偏差が５万円の正規分布に従うという確率過程から生成されていると推測でます。

この

のようなモデルを作成することが時系列解析の目的となるのです。

平均・自己相関・自己共分散

さて、時系列解析の目的を何となく理解したところで、時系列データにおける平均、自己相関、自己共分散の紹介もしておきます。

通常の統計学で用いられる平均、相関、分散の時系列バージョンだと思ってください。平均値は簡単で、純粋にデータの平均値です。データが過去の自分とどれぐらい似ているかを表すのが自己相関、データがどれぐらいばらついているのかを表すのが自己共分散になります。

先ほどのアイスクリームの売り上げは、先月の売り上げに誤差を足したものが今月の売り上げになっていました。ということは、先月の売り上げと今月の売り上げには相関があるでしょう。

以下は縦軸に月の売り上げ、横軸にその一月前の売り上げをプロットしたものです。

この相関係数を計算してみると0.7ぐらいになります。 さらに、これを10ヶ月分算出して、横軸をずらした期間、縦軸を相関係数としてみるとこんなグラフが書けます。

これをコレログラムと言います。青の点線は相関係数の95%信頼区間を表していて、これより相関係数が小さければ相関しているとは言えません。

予想通り、１ヶ月前の売り上げとは相関がありますが、２ヶ月以上離れるとほぼ無相関であることが分かりますね。

定常過程（重要）

準備編の最後に時系列データにおいて重要な確率過程の一つである定常過程を紹介します。定常過程の定義は、時点によらず平均値が一定で、自己共分散・自己相関が時間間隔にのみ依存する確率過程のことです。（厳密には弱定常、強定常があり、これは弱定常の定義なのですがとりあえずこのまま進めます。）

よくわからないですね。何がともあれとりあえず定常過程とそうではない時系列データを見てみましょう。

上の二つのグラフは2020年から2040年まで毎月記録されたあるデータです。最初のグラフが定常過程ではないデータ、二つ目のグラフが定常過程に従うデータです。定常過程のデータはギザギザしていますが、値は10の周りに分布していることがわかるかと思います。

これが時点によらず平均値が一定であるということです。また、自己共分散が時間間隔のみに依存するということは、2030年から2040年のデータでの今月と先月の関係も2020年から2030年のデータでの今月と先月の関係も同じであるということです。

定常だと何が嬉しいのか？

定常過程というものがだいたい想像がついたかと思います。しかし、この定常過程は何が嬉しいのでしょうか？

先ほど、時系列解析の目的はもっともらしい確率過程を見つけることであると書きました。

もっともらしい確率過程とは過去から未来にかけて、同じような確率分布からデータが生成されていると仮定できるものです。つまり、時系列が定常過程であれば、平均値と分散を過去のデータから決めることができるので、未来のデータに対しても、もっともらしい確率過程を見つけ出すことが可能であるということです。

アイスクリームの売り上げも先月の売り上げから５万円ぐらいバラつくことがわかるよりも、平均的に10万円の売り上げから5万円ぐらいバラつくとわかった方が見通しが立てやすいですよね。

さらに平均と分散がわかるということは、統計的な検定ができたりもします。 これが未来の見えない時系列解析において定常過程が重要な役割をもつ理由です。

----- 3/13 追記 -----

ホワイトノイズ

最後に最もよく用いられる定常過程であるホワイトノイズを紹介します。ホワイトノイズとは平均値が0、分散が一定で自己相関を持たない定常過程のことを指します。

何に使われるかというと、時系列データにおける確率的変動部分（撹乱項）を表す時に用いられます。

よく分かりませんね。例えば、先ほど紹介したトレンドモデルで表される時系列データがあるとしましょう。この時系列データは時間が経つごとに値が大きくなっていきますが、年によって多少の誤差があるのが普通です。その誤差を表すのにホワイトノイズは用いられます。

言い換えると、もっともらしい確率過程の統計モデルを作成し、そのモデルの誤差を表現するのにホワイトノイズは利用されるのです。以下は一つ目がホワイトノイズ、二つ目がトレンドモデルにホワイトノイズを足したグラフです。

通常このホワイトノイズは、平均が0の正規分布を仮定することが多いです。

まとめ

今回は、準備編ということで、時系列解析における基本事項をまとめてみました。

• 時系列データは過去の値からも影響される
• 時系列解析の目的は、もっともらしい確率過程を見つけ出すこと
• 定常過程であれば、過去のデータからもっともらしい確率過程を見つけられること
• 確率過程の誤差にはホワイトノイズを使用すること

この基本事項を押さえておけば、時系列解析を学ぶ時に役立つと思います。

次回からは、実際に時系列データをモデリングする方法をまとめていきます。

www.medi-08-data-06.work

※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論、筆者の勘違い等で誤りがあった際はご指摘頂けると幸いです。

参考

最初に時系列の概略を学ぶにはオススメの書籍です。時系列解析はイメージが湧きにくく難解な書籍が多いですが、この書籍は非常に読みやすい内容になってます。

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

Rを使っている方は是非ともオススメしたい書籍です。少しだけ古いパッケージを使用している箇所もありますが、理論、実践ともに分かりやす、時系列解析のイメージが掴みやすいのです。

Rによる計量経済分析 (シリーズ〈統計科学のプラクティス〉)

言わずもがな時系列解析といえば一番有名な沖本本です。理論的な背景がしっかりとまとめてあり、時系列を学ぶための必読書といっても過言ではありません。少しレベルは高いので、上記２冊のような入門書の後に読むことをオススメします。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

Rでの実行方法

今回のデータセットの作成方法とコレログラムの書き方です。今回は主にコレログラムをplotするacf関数さえ押さえておけば問題ありません。

library(ggplot2)
library(dplyr)
library(lubridate)

#アイスクリームデータの作成
set.seed(123)
date <- seq(ymd("2023/1/1"),ymd("2025/1/1"),length.out  = 12*2 )
y <- rep(0,24) 
eps <- rnorm(24,0,3)
y[1] <- 50
for(i in 2:24){
  y[i] <- y[i-1]+eps[i]
}

#作図
data_frame(date=date,ice=y) %>% 
  ggplot(aes(x=date,y=ice))+
  geom_line()


#一月前との自己相関の出し方
#dplyrのlag関数は一つずらしたベクトルを返す
cor(y,lag(y),use = "pair")

#コレログラムを書く
acf(y)

#lag.maxで、ずらすlagの最大値を決められる
acf(y ,lag.max=10)

# 定常過程データの作り方
set.seed(123)
date <- seq(ymd("2020/1/1"),ymd("2040/1/1"),length.out  = 12*20 )
y <- rep(0,12*20) 
eps <- rnorm(12*20,0,1)
y[1] <- 10
for(i in 2:240){
  y[i] <- 5+0.5*y[i-1]+eps[i]
}

# 非定常過程データの作り方
date <- seq(ymd("2020/1/1"),ymd("2040/1/1"),length.out  = 12*20 )
set.seed(123)
eps <- rnorm(12*20,0,1)
y[1] <- 10
for(i in 2:240){
  y[i] <- y[i-1]+eps[i]
}

# 作図
data_frame(date=date,value=y) %>% 
  ggplot(aes(x=date,y=value))+
  geom_line()

#ホワイトノイズの作成
x <- seq(1,100)
wh <- rnorm(100,0,1)
y <- 0.1*x+wh

library(patchwork)#ggplotで複数のグラフを一枚にまとめるパッケージ

#作図
ggplot(data=data_frame(value=wh,time=x),aes(x=time,y=value))+
  geom_line()+
ggplot(data=data_frame(value=y,time=x),aes(x=time,y=value))+
  geom_line()+
  geom_abline(intercept = 0,slope = 0.1,col="red")+
plot_layout(nrow=2)