医療職からデータサイエンティストへ

統計学、機械学習に関する記事をまとめています。

トップ > 時系列解析 > じっくり学ぶ時系列解析~ARIMAの予測と季節調整SARIMA編~

じっくり学ぶ時系列解析~ARIMAの予測と季節調整SARIMA編~

時系列解析 R

前回は時系列解析の流れ、そしてARIMAモデルをデータから同定する方法について書きました。

www.medi-08-data-06.work

今回は、ARIMAモデルの予測について、そして単純なARIMAモデルでは説明できない周期性を、季節調整を使ってモデリングする方法について書いていきます。

ARIMAモデルの予測について

ARIMAモデルの予測とはどのようになるのでしょうか。ARIMA過程の差分系列はARMA過程になるので、ここでは簡単なARMA過程の予測を考えてみます。

例えば、こんなARMA(2,1)過程があったとしましょう。

y_{t} = 0.2y_{t-1} + 0.5y_{t-2} +0.5 \epsilon_{t-1}+\epsilon_{t}

ここで、 y_{t} = -3 \,, y_{t-1} = 1\,,\epsilon_{t} = 1 \,と与えられている時、一期先の予測値y_{t+1}は逐次的に求めることができて、

y_{t+1} = 0.2 \times -3 + 0.5\times 1 +0.5 \times 1 = 0.4

となります。 \epsilon_{t+1}は誤差の期待値が0であることから、最適予測解は0になるため無視できます。

ここから2期、3期先の予測値は、

y_{t+2} = 0.2\times 0.4 + 0.5\times-3 = 0.8-1.5 = -0.7

y_{t+3} = 0.2\times -0.7 + 0.5\times0.4 = -0.14-0.2 = -0.34

これを繰り返していくと予測値はARMA過程の平均値に収束していきます。つまり、ARMA過程の予測値は、超短期的にはある程度予測として成り立ちますが、長期的にはほとんど予想できないのです。

実際にARMA(2,1)過程をシミュレーションで生成し、50期先まで予測をしてみると

時系列

青線が原系列、赤線が予測値になります。

初めは若干下降傾向にありますが、一定の値に収束しているのがわかるかと思います。ちなみに灰色は95%信頼区間を表し、データの±2SD（赤破線）に収束します。

グラフをみるとそりゃそうだという感じです。

ARMAモデルは過去の自分の情報しか持ち合わせていないため、予測はあまりパッとしないのです。これが一変量時系列モデルの限界ですね。

Rでの実践

library(forecast)
library(ggplot2)
#シミュレーションデータの作成
set.seed(123)
data <-  arima.sim(list(order = c(2,0,1), 
                        ar = c(0.2,0.5),
                        ma=c(0.5)), n = 300)

#ARMAモデルの作成
mod <- Arima(data,order = c(2,0,1))

作成したモデルから予測を行うにはforecast{forecast}かpredict{stat}を使います。どちらでもいいですが、相性はforecastの方がいいのでこちらを使います。

ちなみにグラフはplot{base}を使えば一発で終わりますが、綺麗に書きたいという欲望から{ggplot2}を使っています。

#50期先までを予測
#levelはc(50,95)などとすれば複数の信頼区間を出してくれる

pred <- forecast(mod,h = 50,level = 95)

#データフレームに加工し、可視化
data_frame(time=c(seq(1,300),seq(300,349)),
           value=c(data,pred$mean),
           lower = c(data,pred$lower),
           upper=c(data,pred$upper),
           type = rep(c("row","pred"),c(300,50))) %>% 
  ggplot(aes(x=time))+
  geom_line(aes(y=value,col=type))+
  geom_ribbon(aes(ymax=upper,ymin=lower),alpha=0.3)+
  geom_hline(yintercept = sd(data)*2,col="red",linetype=2)+
  geom_hline(yintercept = -sd(data)*2,col="red",linetype=2)

#plotを使えば一発でグラフにできる
plot(pred)

季節調整モデル

さてさて、次はある周期に従うデータのモデル推定と予測を行なっていきます。

今回扱うのは有名な'Airpassenger'のデータです。これは飛行機の乗客数を1949年から1960年まで記録したデータです。Rにもデフォルトで入っていますが、今回はcsvから読み込んで使います。データはこちらからダウンロードできます。※外部サイトへ行きます。Pythonをお使いの方はこちらも参考になるかと思います。

Complete guide to create a Time Series Forecast (with Codes in Python)

早速データをみてみましょう。グラフにしてみるとこんな感じになります。

時系列

右肩上がりのギザギザしたグラフですね。分散もだんだんと大きくなっているようです。

まずは対数差分で前処理をして、トレンドの除去と分散の均一化をしてみましょう。こうなります。

時系列

いい感じに定常過程っぽくなりました。次に自己相関とのコレログラムを書いてみます。

時系列

どうやら12期ごとに強い自己相関がありそうです。今回のように自己相関が、ある一定の周期で強くなる場合には季節調整モデル、通称SARIMAモデルを使います。

SARIMAモデルを使う前に、まずは通常のARIMAモデルで予測をしてみましょう。

前回同様にAICでモデル選択をしてみると、ドリフト（定数項）のあるARIMA(3,1,3)モデルが一番良さそうなのでこれでいきましょう。5年先までの予測値をグラフに書いてみると

時系列

こんな感じになります。周期がなくなり、先ほどと同様にARIMA過程の期待値へと収束していきます。ちなみにドリフトがあるARIMAモデルは、右肩上がりになります。

モデルの残差のコレログラムを書いてみると

時系列

正しくモデルが作成できていれば残差はホワイトノイズになり、自己相関がなくなるはずですが、周期性のある自己相関が残っています。

さて、ここで季節調整を加えていきましょう。季節といっても年周期だけでなく、曜日による週の周期、４半期周期など仮定する周期はなんでもOKです。

季節調整方法については、様々な議論があるので、詳しくは触れませんが、トラディショナルな移動平均を使った調整法や、x-11、x-12と呼ばれる方法があります。詳しくはこちらをご参考ください。

入門 季節調整―基礎知識の理解から「X‐12‐ARIMA」の活用法まで

入門 季節調整―基礎知識の理解から「X‐12‐ARIMA」の活用法まで

今回はloessアルゴリズムという局所重みづけ回帰の手法を使って、季節変動を捉えています。まずは季節調整済みの値をグラフに書いてみましょう。こんな感じになります。

時系列

季節の周期性がなくなり、ギザギザの小さい右肩上がりのグラフなっていますね。 この季節調整済みの値の対数差分自己相関コレログラムをみてみると

時系列

こちらも周期性がなくなっていることが分かります。

それでは、12期を一周期と過程した、SARIMAモデルを使って先ほどの同じように予測をしてみると

時系列

おお！うまく季節の周期性を捉えて長期的にもそれっぽい予測ができているようです！

モデルの残差を確認してみても、自己相関のないホワイトノイズになっていることが分かりますね。

時系列

Rでの実践

通常ARIMAでの予測

library(forecast)
library(ggplot2)
library(dplyr)

#データの読み込みと、変数名の加工
data <- read_csv("AirPassengers.csv",
                 col_types = cols(Month=col_date(format = "%Y-%m")))%>% 
              rename(passengers = `#Passengers`)

#可視化
data %>% 
  ggplot()+
  geom_line(aes(x=Month,y=passengers))

#対数差分系列の可視化
data %>% 
  mutate(log_diff = c(NA,diff(log(passengers)))) %>% 
  ggplot()+
  geom_line(aes(x=Month,y=log_diff))

#48期までの自己相関コレログラム
acf(diff(log(data$passengers)),lag.max = 48)

ここからauto.arima{forecast}でAICに基づいてモデルを選択し、予想までを行います。forecast{forecast}で予測した値は対数変換されているため、グラフに書くときにはexp{base}で戻してグラフを書きます。

#auto.arimaでARIMA(3,1,3)with driftが選ばれる
mod <- auto.arima(log(data$passengers))

#予測をする
pred <- forecast(mod,12*5,level = 95)

#exp()で予測値を戻してからグラフ化する
data_frame(Time = seq(1,12*17),
           value = c(data$passengers,exp(pred$mean)),
           lower = c(data$passengers,exp(pred$lower)),
           upper=c(data$passengers,exp(pred$upper)),
           type=rep(c("row","pred"),c(12*12,12*5))) %>% 
  ggplot(aes(x=Time))+
  geom_line(aes(y=value,col=type))+
  geom_ribbon(aes(ymax=upper,ymin=lower),alpha=0.3)

時系列

SARIMAモデル

季節調整を行うためには、stl{stats}を使いますが、tsオブジェクトしか受け付けてくれないため、まずはts{stats}を使って、tsオブジェクトに変換します。

#対数変換したデータを12期を一周期としたtsオブジェクトに変換する
s_data <- ts(log(data$passengers),frequency = 12)

#ここにも周期を指定して、季節調整を行う。
sy <- stl(s_data,s.window =12)

stlでは、s.windowを指定する必要があり、ここにも一周期を指定します。stlからは、季節調整値、トレンド、残差が返ってきます。

こちらのグラフもplotを使えば一発で書けますが、ggplotを使っています。好みの問題です。

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

#季節調整値、トレンド、残差をグラフ化する
ggplot(data=data_frame(time=data$Month,season=sy$time.series[,1]),
       aes(x=time,y=season))+
  geom_line()+
ggplot(data=data_frame(time=data$Month,trend=sy$time.series[,2]),
         aes(x=time,y=trend))+
  geom_line()+
ggplot(data=data_frame(time=data$Month,resid=sy$time.series[,3]),
         aes(x=time,y=resid))+
  geom_bar(stat="identity")+
  plot_layout(nrow=3)

#plotを使えば一発です。
plot(sy)

時系列

トレンド、季節調整値ともにうまく推定できているようですね。 季節調整値12年分の各月平均値をみてみます。

#時刻を扱うパッケージ
library(lubridate)

data %>% 
  mutate(season = month(Month),
         adj_s = exp(sy$time.series[,1])) %>% #対数変換しているため元に戻す
        group_by(season) %>% 
        summarise(mean=mean(adj_s)) %>% #月ごとに集計して可視化
　　ggplot(aes(x=season,y=mean,fill=season))+
　　geom_bar(stat="identity")

時系列

夏がプラスに調整され、冬はマイナス方向に調整されているようです。夏は旅行客が多いみたいですねー

ここから季節調整済の値をグラフに書くには以下のようにします。

#原系列から季節調整値を引いたものをexp()で元に戻す。
data %>% 
  mutate(non_s = exp(log(passengers)-sy$time.series[,1])) %>% 
  ggplot()+
  geom_line(aes(x=Month,y=passengers))+
  geom_line(aes(x=Month,y=non_s),col="red")

#季節調整済の自己相関コレログラム
acf(diff(log(data$passengers)-sy$time.series[,1]),lag.max = 48)

時系列

それでは、SARIMAモデルを作っていきましょう。SARIMAモデルを作るには、Arima{forecast}に季節調整項を加えて係数を推定してもらいます。また、auto.arima{forecast}に周期を指定したtsオブジェクトを渡すと、AICに基づいてモデルを推定してくれます。

#ArimaにARIMAの次数と季節調整項を指定してモデルを作る
mod <- Arima(log(data$passengers),
      order = c(0,1,1),
      seasonal = list(order=c(0,1,1),period=12))

#auto.arimaに周期を指定したtsオブジェクトを渡してモデルを推定してもらう。
s_data <- ts(log(data$passengers),frequency = 12)
mod <- auto.arima(s_data,seasonal = T)

モデルを作ってしまえば、あとは先ほどと同じように予測を行います。

pred <- forecast(mod,12*5,level = 95)

data_frame(Time = seq(1,12*17),
           value = c(data$passengers,exp(pred$mean)),
           lower = c(data$passengers,exp(pred$lower)),
           upper=c(data$passengers,exp(pred$upper)),
           type=rep(c("row","pred"),c(12*12,12*5))) %>% 
  ggplot(aes(x=Time))+
  geom_line(aes(y=value,col=type))+
  geom_ribbon(aes(ymax=upper,ymin=lower),alpha=0.3)

#残差の自己相関コレログラム
acf(mod$residuals,lag.max = 48)

時系列

まとめ

一変量時系列は過去の自分の情報から予測を行うため、どうしても長期の予測は難しくなりますが、季節変動やトレンドを捉えることで、ある程度精度の良い予測をすることができます。

今回で一変量時系列を一区切りにして、次回からはいよいよ多変量時系列に入っていきます。

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

参考書籍

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

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

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

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

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

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

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

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

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

その他

Pythonをお使いの方はこちらが参考になるかと思います。 Complete guide to create a Time Series Forecast (with Codes in Python)

Rの方はこちら A Complete Tutorial on Time Series Modeling in R