じっくり学ぶ時系列解析~基礎編~ - 医療職からデータサイエンティストへ

前回は時系列解析の考え方や基本用語についてまとめました。今回はもう一歩進んで、実践に近づいた内容にしていきます。

www.medi-08-data-06.work

今回扱う範囲は、時系列データの前処理方法、AR、MA、ARMA、ARIMAです。それではいきましょう。

時系列解析の概略
データの可視化
定常過程の確認
前処理
統計モデルの同定
Rでの実践
- 前処理編
- モデル同定編
まとめ
- 前処理方法
- 定常過程の特徴
参考書籍
その他

時系列解析の概略

時系列解析は以下のステップで行うのが基本的な流れになります。

f:id:h-wadsworth02:20190316173746j:plain

データの可視化

まずはデータを可視化します。時系列における可視化は、データがどんな動きをしているのか、トレンドがあるのか、定常過程に従いそうか、周期性がありそうか、などデータのイメージを掴むためにも重要です。

<可視化イメージ>

時系列

定常過程の確認

可視化でイメージを掴んだらデータが定常過程に従っているかどうかを確認します。定常過程に従っていれば統計モデルを作りやすいんでしたね。定常過程かどうかの検定方法については次回以降に扱っていきますが、見た目でもある程度は分かります。

www.medi-08-data-06.work

前処理

もし、定常過程に従っていなさそうであれば、以下の前処理を行っていきましょう。よく使われる手法を紹介しています。

対数変換

対数変換は時系列データの分散が一定ではない時に行うことで、分散を一定と見なせる場合があります。例えばこんなデータ

時系列

値が大きくなっていくごとに分散が広がっているのが分かります。これを対数変換すると

時系列

データのバラツキが少なくなりましたね。対数変換をさらに一般化した方法としては、Box-Cox変換というものもありますが、使用頻度はあまり多くありません。

差分（階差）

差分はデータのトレンドを除くために行われます。先ほどのデータは右肩上がりのトレンドを持つデータになっていました。一つ前の値から次の値を引く一階差分にするとトレンドを除去することができます。一階差分を行ったデータが以下になります。

時系列

トレンドが除かれていますね。ちなみに数学的には、時系列データを $y_{t}$ とし、トレンドが直線であれば、

$y_{t} = a + bt + \epsilon$

と表すことができ、ここから一階差分をとると

$\Delta y_{t} = (a + bn + \epsilon_{t}) - (a + b(t-1) + \epsilon_{t-1})\\ = b + \Delta \epsilon$

となります。定数に誤差を足した値になっていることから、トレンドが除去できていることが分かります。また、トレンドが除去できない場合は、一階差分をさらに一階差分にする二階差分を行う場合もあります。

対数差分

対数差分は時系列の前処理において、データの変動を明らかにしたい時に用いられます。例えば、企業収益やお店の売り上げなど一つ前の時点からどれぐらい変化したかを知りたい場合などが当てはまります。

対数差分は、対数変換した系列の一階差分ですが、一つ前の値と次の値の比をとったものを対数変換していることと同じなので、変動を表すことができるのです。

$\Delta(log(y_{t})) = log(y_{t}) - log(y_{t-1})\\ = log \dfrac{y_{t}}{y_{t-1}}$

先ほどのデータの対数差分をとってみましょう。

時系列

y軸をみて頂くとわかる通り、変動のバラツキがだいぶ小さくなりました。

以上のような前処理行うことでデータを定常過程と見なせる場合があります。

統計モデルの同定

さて、データが定常過程に従っていることが分かれば、いよいよもっともらしい確率過程を見つけていきます。しかし、確率過程と言っても無数存在します。その中からどれがもっともらしいか一つずつ試すわけにも行かないので、定常過程に従う統計モデルの型が必要になります。

それが今回扱うMA過程、AR過程、ARMA過程です。この中からもっともらしそうな型を決める作業が統計モデルの型を決めることになります。この作業をモデルの同定と呼びます。

これらの型は時系列データ特有の過去の自分と相関があるという特徴を持っている確率過程です。

MA過程（Moving Average process ）

MA過程は過去の自分の誤差と相関を持つ確率過程になります。一期前の自分と相関のあるMA過程をMA(1)とし、以下のように表されます。

$y_{t} = \mu + \theta_{1}\epsilon_{t-1}+\epsilon_{t}$

$\epsilon$ は正規ホワイトノイズと仮定することが多いです。ここで、一期前の値は

$y_{t-1} = \mu + \theta_{1}\epsilon_{t-2}+\epsilon_{t-1}$

と表すことができ、 $\epsilon_{t-1}$ を共通項に持つので、一期前と相関する確率過程を表すことができます。

$\mu = 10 \,, \theta = 0.5$ としたMA(1)をグラフにしてみるとこのような感じになります。

時系列

平均値10の周りを一定の分散でばらついているので、定常過程に従っていますね。

このMA(1)過程のコレログラムを書いてみると

時系列

見ての通り2期以上離れると自己相関が0になっています。

ここで偏自己相関というものも紹介しておきます。偏自己相関とは、ある時点tをt+k時点との自己相関を計算する時に、t+1、t+2、t+3.....t+k-1時点との自己相関を除いた自己相関と言えます。イメージしにくいかもしれませんが、今はそんなものがあるんだぐらいで構いません。

その偏自己相関のコレログラムを書いてみるとこんな感じになります。

時系列

少し分かりにくいですが、偏自己相関が少しずつ減衰しているように見えませんか？実は、MA過程の特徴として偏自己相関は段々と減衰することが分かっています。

ここまででMA過程の特徴を一般化してまとめておきます。

MA(q)過程

$y_{t} = \mu + \theta_{1}\epsilon_{t-1} + \theta_{2}\epsilon_{t-2} \dots \theta_{q}\epsilon_{t-q} + \epsilon_{t}$
q時点までと自己相関を持ち、q時点以降は０になる
偏自己相関は少しずつ減衰する。

AR過程(AutoRegressive process)

AR過程は過去の自分の値そのものと自己相関を持つ確率過程です。一期前の自分と相関を持つAR(1)過程は以下のように表されます。

$y_{t} = \phi_{1}y_{t-1} + \epsilon_{t}$

実は前回のアイスクリームの売り上げモデルである、今月の売り上げは先月の売り上げに誤差を足した値であるという

$y_{t} = y_{t-1} + \epsilon_{t}$

この式もAR過程だったのです。

ここで $y_{0} = 5 \,, \phi = 0.5$ としたAR(1)過程は以下のような感じになります。

時系列

先ほどと同じように10の周りを一定の分散で散らばっているので、定常過程に従っていそうです。

このAR(1)過程の自己相関のコレログラムを書いてみるとこんな感じになります。

時系列

AR(1)過程では、自己相関が少しずつ減衰することが分かります。

また、偏自己相関のコレログラムを書いてみると

時系列

1期までは偏自己相関がありますが、それ以上先になると0になってなっていますね。

AR過程の特徴を一般化すると

AR(p)過程

$y_{t} = \phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\dots\phi_{p}y_{t-p} + \epsilon_{t}$
自己相関は少しずつ減衰する。
p時点までと偏自己相関を持ち、p時点以降は0になる。
$| \phi | \lt 1$ の時に定常となる

一番下の特徴はAR過程において重要で、AR過程の全ての係数 $\phi$ の絶対値を足した値が1より小さい時に限りAR過程は定常になります。従って、先ほどのアイスクリームの売り上げモデルは係数が1であるため定常とはならないんですね。

ARMA・ARIMA過程

さて最後はAR過程、MA過程の特徴を合わせ持ったARMA過程です。ARMA過程はARとMAの二つのパラメーターを持つため、一次のARMA過程をARMA(1 , 1)とすると

$y_{t} =\mu + \phi_{1}y_{t-1} + \theta_{1}\epsilon_{t-1}+ \epsilon_{t}$ と表すことができます。

$\mu = 5 \,, \phi = 0.5 \,, \theta=0.5$ としたARMA(1,1)過程はこんな感じになります。

時系列

こちらも定常過程に従っていそうですね。 ARMAの自己相関、偏自己相関のコレログラムを書いてみると

時系列

どちらも少しずつ減衰しているように見えます。実はARMA過程の場合は自己相関、偏自己相関ともに指数的に減衰するという特徴があります。

特徴を一般化してまとめると

ARMA(p,q)過程

$y_{t} = \phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\dots\phi_{p}y_{t-p} + \theta_{1}\epsilon_{t-1} + \theta_{2}\epsilon_{t-2} \dots \theta_{q}\epsilon_{t-q} +\epsilon_{t}$
自己相関、偏自己相関は少しずつ減衰する。
p時点までと偏自己相関を持ち、p時点以降は0になる。
$| \phi | \lt 1$ の時に定常となる

MA過程は常に定常になります。ARMA過程が定常になる条件としては、AR過程が定常、つまり全ての係数の絶対値 $| \phi|$ を足した値が1より小さい時にARMA過程は定常となります。

さらに一階差分をとった系列がARMA過程に従う時、この系列をARIMA過程(Auto Regressive Integrated Moving Average process)と呼びます。

ARIMA過程はARIMA(p,d,q)と表し、dは差分値、p,qはAR,MA過程の次数を表します。

モデルの同定

さて、それではここから統計モデルを作成する方法について紹介してきましょう。

モデルの同定には以下の二つのステップが一般的です。

step1 : 自己相関、偏自己相関コレログラムからモデル候補を見つける。
step2 : 情報量基準を用いてモデルを同定する。

step1では、データを前処理した後の自己相関、偏自己相関のコレログラムの特徴からそれらしいモデル候補を見つけます。step2では、それらのモデル候補から情報量基準(一般的にはAIC)を使って、一番もっともらしいモデルを同定します。

www.medi-08-data-06.work

実際にやってみましょう。今回のデータはこんなデータです。

時系列

まずは前処理として差分をとってデータのトレンドを無くします。

時系列

定常過程っぽくなりましたね。このデータの自己相関（上）、偏自己相関（下）のコレログラムを書いてみましょう。

時系列

自己相関は徐々に減衰しているのでMA過程ではなさそうです。次に偏自己相関をみてみると少しずつ減衰しているようにも見えますが、4期目以降の自己相関が3期目までと比べて小さいのでAR過程の可能性もありそうです。

ここから、AR(3)過程とARMA過程を候補のモデルとします。ARMA過程はコレログラムから次数を判断することが難しいためARMA(2,2)までのARMA(1,1),ARMA(1,2),ARMA(2,1),ARMA(2,2)を網羅的にモデル候補とします。

これら5個のモデルのAICを求めて、AICが一番小さくなったモデルがもっともらしいモデルと言えそうです。

実際にやってみると

model	AIC
AR(3)	1431.832
ARMA(1,1)	1519.899
ARMA(1,2)	1407.173
ARMA(2,1)	1484.817
ARMA(2,2)	1396.538

となります。ここからARMA(2,2)過程、差分をとっているのでARIMA(2,1,2)過程が一番当てはまりのよいモデルであることが分かります。

Rでの実践

ARIMA過程を扱うには{forecast}パッケージをインストールします。コードは以下になります。

前処理編

library(forecast)
library(ggplot2)

#徐々に分散が大きくなるデータの作成
mu <- seq(100,5000,length.out = 100)
sd <- mu*0.2
data <- mu+rnorm(100,0,sd=sd)

#可視化
ggplot(data=data_frame(Time=seq(1,100),Value = data))+
  geom_line(aes(x=Time,y=Value))

# 対数変換したデータの可視化
ggplot(data=data_frame(Time=seq(1,100),Value = log(data)))+
  geom_line(aes(x=Time,y=Value))

#差分をとったデータの可視化
ggplot(data=data_frame(Time=seq(1,99),Value = diff(data)))+
  geom_line(aes(x=Time,y=Value))


#対数差分データの可視化
ggplot(data=data_frame(Time=seq(1,239),Value = diff(log(data))))+
  geom_line(aes(x=Time,y=Value))

モデル同定編

MA過程、AR過程、ARMA過程はパッケージを使わずに作成しましたが、最後の練習データであるARIMA過程は{forecast}パッケージを使用しました。

#MA過程
set.seed(123)
eps <- rnorm(500,0,10)
ma <- rep(0,500) 
ma[1] <- 10
for(i in 2:500){
  ma[i] <- 10+0.5*eps[i-1]+eps[i]
}
ggplot(data=data_frame(Time=seq(1,500),Value = ma))+
  geom_line(aes(x=Time,y=Value))

#自己相関のコレログラム
acf(ma)

#偏自己相関のコレログラム
acf(ma,type="partial")

#AR過程
set.seed(123)
eps <- rnorm(500,0,10)
ar <- rep(0,500) 
ar[1] <- 10
for(i in 2:500){
  ar[i] <- 5+0.5*ar[i-1]+eps[i]
}

ggplot(data=data_frame(Time=seq(1,500),Value = ar))+
  geom_line(aes(x=Time,y=Value))

acf(ar)
acf(ar,type="partial")

#ARMA過程
set.seed(123)
eps <- rnorm(500,0,10)
arma <- rep(0,500) 
arma[1] <- 10
for(i in 2:500){
  arma[i] <- 5+0.5*arma[i-1]+0.5*eps[i-1]+eps[i]
}

ggplot(data=data_frame(Time=seq(1,500),Value = arma))+
  geom_line(aes(x=Time,y=Value))

acf(arma)
acf(arma,type="partial")

{forecast}のarima.sim関数はARIMA過程に従うシミュレーションデータを作成できます。今回は一階差分系列が

$y_{t} = 0.2y_{t-1}+0.3y_{t-2}+0.6\epsilon_{t-1} + 0.6\epsilon_{t-2} +\epsilon_{t}$

のARMA(2,2)過程に従うデータを作成しました。

#ARIMA過程に従うデータの作成

set.seed(123)
arima <- arima.sim(list(order = c(2,1,2), 
                        ar = c(0.2,0.3),
                        ma=c(0.6,0.6)), n = 500)

#可視化
ggplot(data=data_frame(Time=seq(1,501),Value = as.numeric(arima)))+
  geom_line(aes(x=Time,y=Value))

#一階差分
arma <- diff(arima)

#コレログラム
acf(arma)
acf(arma,type="partial")

ARIMA過程に従うモデルを推定するには{forecast}のArimaを使います。次数をorderの中に(p,d,q)の順に指定するとその係数の推定やAICの計算などを自動でやってくれます。

>  Arima(arma,order=c(1,0,2))
Series: arma 
ARIMA(1,0,2) with non-zero mean 

Coefficients:
         ar1     ma1     ma2    mean
      0.3415  0.4307  0.6240  0.0762
s.e.  0.0717  0.0689  0.0439  0.1360

sigma^2 estimated as 0.9621:  log likelihood=-698.59
AIC=1407.17   AICc=1407.29   BIC=1428.25

これを使って5つのモデルでAICを計算します。

#5つのモデルでAICを計算

model_lis <- list(c(3,0,0),c(1,0,1),c(1,0,2),
                  c(2,0,1),c(2,0,2))
AIC <- rep(0,5)
i <- 1
for(mod in model_lis){
  res <-Arima(arma,order = c(mod))
  AIC[i] <- res$aic
  i <- i+1 
}

また、各モデルに対してそれぞれAICを計算しましたが、実は{forecast}のauto.arimaを使うとAICに基づいて勝手にモデルを推定してくれます。

> auto.arima(arma,max.d=0,max.p = 3,max.q = 3)
Series: arma 
ARIMA(2,0,2) with zero mean 

Coefficients:
         ar1     ar2     ma1     ma2
      0.1700  0.2346  0.5757  0.5922
s.e.  0.0668  0.0619  0.0567  0.0438

sigma^2 estimated as 0.9382:  log likelihood=-692.38
AIC=1394.76   AICc=1394.88   BIC=1415.83

同じようにARMA(2,2)が選ばれました。係数をみても、シミュレーションで設定した値に近い係数が推定されていますね。

まとめ

今回はデータの前処理からモデルの同定方法までをまとめました。次回からはいよいよ予測の方法や多変量時系列モデルについて書いていきます。

以下今回のまとめです。

前処理方法

対数変換：分散が一定ではない時に用いる
差分：トレンドを除去するために用いる
対数差分：変動を確認するために用いる

定常過程の特徴

MA(q)過程

$y_{t} = \mu + \theta_{1}\epsilon_{t-1} + \theta_{2}\epsilon_{t-2} \dots \theta_{q}\epsilon_{t-q} + \epsilon_{t}$
q時点までとの自己相関を持ち、q時点以降は０になる
偏自己相関は少しずつ減衰する。

AR(p)過程

$y_{t} = \phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\dots\phi_{p}y_{t-p} + \epsilon_{t}$
自己相関は少しずつ減衰する。
p時点までとの偏自己相関を持ち、p時点以降は0になる。
$| \theta | \lt 1$ の時に定常となる。

ARMA(p,q)過程

$y_{t} = \phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\dots\phi_{p}y_{t-p} + \theta_{1}\epsilon_{t-1} + \theta_{2}\epsilon_{t-2} \dots \theta_{q}\epsilon_{t-q} +\epsilon_{t}$
自己相関、偏自己相関は少しずつ減衰する。
$| \theta | \lt 1$ の時に定常となる。

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

www.medi-08-data-06.work