時系列解析　ことはじめ ①

概要

統計手法でも特に難解な「時系列解析」の理論的・実装的な面での知識の整理を行いたい.

きっかけは 『時系列解析』(共立出版, 2017)で理論的な部分を勉強しようとした際, 時系列解析特有の考え方に慣れていないせいか, あまり理解が進まなかった. とはいえ, 解析さえできればよい, とライブラリやコードの書き方だけを知っていても意味がない^[1]

最近, 東京大学による公開講座で　名著をたくさん出されている北川先生による『時系列解析』の講座をみて, 理解がかなり進んだ.

本記事では, 公開講座や他の文献を参考に, 時系列データの分析で重要な「前処理・統計量の分析/可視化」を整理する. その際には, データセットを用いながら, プログラム例も提示する.

時系列データとは？

まず, 時系列データとは, 時間の経過とともに変動する現象を記録したデータ(例：降雨量, 株価, 感染者数, 地震動)である. そして, 時系列データに基づき, その背景にある複雑な現象を理解し, 予測, 制御や意思決定を行うための方法が時系列解析である.

時系列解析では, 時系列データの特徴に合わせた時系列モデルを持ち込むことが重要で, 解析からの結果内容や解釈に大きな影響を与える.特に,時系列データが, 線形で,定常性,正規性が仮定されるならば, とても扱いやすい強力なモデルを使用することができる.

しかし, 世の中の時系列データは, そこまで理想的なデータの性質をもっていることは少ない. そんな複雑なデータをそのまま扱う手もあるだろうが, やはり, モデリングを単純化するためにも, データそのものを単純なものに（極力）なおすことも検討すべきだろう.

中でも, 定常性という性質は時系列解析では特に重要である. 「定常」とは, 時間が経っても時系列データの性質（平均値や分散、共分散）が変化しないというものだ. こうした性質ゆえに, いくつかの時系列モデルはとても有効に働く(ARモデル, MAモデル, ARMAモデルなど）

なぜ定常性が大事？

時系列データの定常性が重要な理由として, 先に述べたモデリングのしやすさがある. それに関連して, 時系列データを統計手法を用いて解析できるために必要な前提に合わせる,という狙いもある

時系列データではない一般的なデータでは, 手元にある標本から, 母集団の特徴を推定するために, 大数の法則をうまく活用している. これは, 十分な標本サイズであれば, 標本に対する推定量が, 母集団の実際の量と一致する（近づいていく）というものだ.

時系列データも,時間に沿ってたくさんのデータが収集されているから大数の法則が効くのではないかと思うかもしれないが, じつはそうではない.

時系列解析の本で散見される次の文章が時系列データの扱いの難しさを説明している.

時系列解析が, ほかの統計解析ともっとも異なるのは, 観測値に順序があること, つまり時間の流れに沿って出現した値を扱う点にある. 時間は先へ進むだけで後戻りできないため, 実時間では繰り返し観測できず, 1本の時系列, つまり１本のサンプルパスだけから統計的推測を行う必要がある.
(『時系列解析』(共立出版, 2017)より)

もう少し具体的に説明すると, 本来, 時点$n$での観測値$y_n$を予測するモデルや解析量を得るには, この時点$n$での観測値$y_n$をいくつも収集しないといけない...が, 時系列データではまずできない作業である^[2].

そんな厄介さがある時系列データであるが, 定常性を仮定すれば, 時点$n$での観測値$y_n$の解析を, 別の時点でのデータ(例: $y_{n-3}$, $y_{n+1}$)を使ってもよくなる. なぜならどの時点のデータをとっても, 時間を通して時系列データの生成モデルの性質に変化がないためであるからだ.

そのため, 定常な時系列データを時間平均（エルゴード性）として扱うことが解析上重要である.

時系列データの可視化

できれば時系列データには”きれいな”定常性をみせてほしいが, 現実のデータで理想的な定常性をもつことは稀である. そのため, 何らかの前処理を経て, 解析しやすいデータに変換する必要がある.

前処理を施すにも, 元の時系列データがどのような性質をもっているかを知っておくことが大事である. なので, 時系列データが手に入ったら, 一度グラフなどに可視化し, 目視で確認することをお勧めする.

いろいろなデータセットをのぞく

説明のためのデータセットは, RパッケージTSSS内のものを利用する.

install.package("TSSS")

library("TSSS") # TSSSパッケージの読み込み

船舶の方向角速度データ: HAKUSAN

data(HAKUSAN)
hakusan1 <- HAKUSAN[:,1] # 方向角速度データのみを取得
plot(hakusan1)

データをみると,

1. $y \in [-2, 0]$のどこかに, 均一な平均線を引けそう（平均定常性）
2. 時間全体を通じて似たような変動をしている（分散・共分散定常性）
3. 変動のパターンに周期性がありそう

太陽黒点数データ: Sunspot

data(Sunspot) # 太陽黒点数データの読み取り
plot(Sunspot)

データをみると,

1. カウントデータなので, 波が上半分だけ(正値, 上下非対称)
2. 定期的にカウントが極大になっている, でも安定した周期かは不明（擬似周期性)
3. カウントの増加・減少はきれいな対称性にはなっていなさそう,,,どちらかといえば減少するときは少し遅い(前後非対称性)

日最高気温データ: Temperature

data(Temperature) # 日最高気温データ
maxtemp <- Temperature
plot(maxtemp, ylim = c(0, 35)) # y軸を 0 - 35 の目盛りに指定

データをみると(1日ごとの収集であることを考慮),

1. 年内で気温の変化にきれいなトレンドがありそう(長期の周期性)
2. トレンドを周りでは, その他の成分はほとんど定常な感じがする

MYE1F(東西方向地震波)データ: MYE1F

data(MYE1F)
plot(MYE1F)

データをみると,

1. トレンドはなさそう
2. 平均は定常（時間を通して同じ？）だが, 分散や共分散は非定常
3. しかし, 一定区間に限定すれば定常な感じ(区分的定常)

榛原地下水位・気圧データ: Haibara

data(Haibara)
plot(Haibara)

データをみると,

1. 同じ観測時点における2変量時系列データである(地下水位・気圧データ）
2. 二つの時系列データは正反対の変動をしている(逆相関)
3. ところどころデータが抜けたり, 飛び出したりしている（異常値, 欠測値）

時系列データの性質分類

時系列データといってもいろいろな性質をもっている.
一般に時系列の分類は以下のような分類ができる.

連続時間時系列と離散時間時系列

時間的に連続に記録されたデータは連続時間時系列(continuous time series)と呼び, ある時間間隔(例： 1日ごと）で収集されたデータを離散時間時系列(discrete time series)と呼ぶ.

実際には, データの収集は何らかのタイミングで行われると思われるので, 時系列データとしては離散的なものを扱うことがほとんどである.そして, さらに離散時間時系列の場合は, 等間隔時系列(例：10分おき）や, 逆に収集タイミングがバラバラな不等間隔時系列にわけられる.

１変量時系列と多変量時系列

各観測時点で１種類の情報のみを得るものが１変量(univariate)時系列で, ２つ以上の情報を同時に取得したものが多変量(multivariate)時系列である.

特に, 多変量時系列の場合, 時系列データ間での相関を考慮することで, 1変量時系列だけの分析よりも予測精度が上がったりする.

定常時系列と非定常時系列

時系列データをみると, どれも不規則な変動となっているが, 時系列解析では, この"不規則な変動"を確率モデルで表現する.
このとき, 一見したところ不規則な変動しているデータでも, その背景にある確率モデルの性質が時間に対して変化しないものを定常時系列(stationary time series) と呼ぶ.
一方で, 時間が変化すると, 確率モデルの性質が変化するようなものを非定常時系列(non-stationary time series) と呼ぶ.

非定常な時系列データを扱うにはいろいろと注意が必要である. しかし, 非定常な時系列データのなかには, 何らかの前処理(例：対数変換, 階差など)を施すことで, 定常にみなせる場合がある.

線形時系列と非線形時系列

時系列モデルが線形和で表される場合(例：過去のデータの重みづけ和で求まるARモデルや,元データが複数の成分の和で構成されている）は,線形時系列とよぶ. 他方, より複雑な非線形モデルが必要なものを非線形時系列という.

ガウス型時系列と非ガウス型時系列

時系列データとして観測された値の分布が正規分布に従うものがガウス型時系列で, そうでない場合を非ガウス型時系列とよぶ.

ここで時系列データがガウス型かどうかを判別する,...というイメージが難しい場合, 時系列データを観測値に関してカウントする, つまりヒストグラムを作ってみるとわかりやすい.

data(HAKUSAN)
hakusan1<- HAKUSAN[, 1]
plot(hakusan1, main = "")
hist(hakusan1, col = "#539952")

これは, 船舶の方向角速度データに関するヒストグラムを作成したものだ. ヒストグラムをよく眺めると観測値$y_n$は左右対称に広がっていることがわかる. これがガウス型時系列の様相である.

となると, 正のトレンドがある場合（全体的に量が増える）や, 上下非対称なデータでは, ヒストグラムを作成してもきれいな正規分布にはならない. これが非ガウス型時系列である.

非ガウス型時系列であっても, 何らかの前処理（例：対数変換）を行えば, 正規分布に近づけることができることもある.

欠測値と異常値

時系列データのなかには, 何らかの理由で観測値が記録されなかったり（欠測値）, 観測の異常やデータ送信時の問題で明らかな異常なデータ(異常値, 外れ値)が取得されることがある.

状態空間モデルを使えば, そうした異常性を意識せずに時系列解析できる.

時系列データの前処理

時系列データといってもさまざまな分類ができることを説明した. いろいろな性質のデータのなかでも, 正規性, 定常性, 線形といった時系列データは非常に扱いやすい. しかし,そのような理想的な時系列データが少ないことも事実である.

しかし, 一見, 非定常な, 非正規性な, 非線形なデータでも, 何らかの前処理を施すことで, 正規性, 定常性, 線形といった扱いやすいデータに変換することができる.

変数変換

まず基本的な要素として, 変数変換のやりかたを紹介する.
例として, 手元にあるデータ$X$が, 何らかの分布$f(x)$から生成されているとする(この$f(x)$は未知とする）.

そして, データ$x$は,変換関数$h$を用いて$y$に変換できるものとする. このとき, 逆変換も定義されているとする.($y = h(x), x = h(y)^{-1}$))　変換後の$y$は, 何からの既知の分布$g(y)$から生成されていることがわかっている.

すると, データ$X$を生成する分布$f(x)$は, 次のように求めることができる.

$$f(x) = g(h(x))\Big|\frac{dh}{dx}\Big|$$

ここで, $\frac{df}{dx}$はヤコビ行列(ヤコビアン） とよばれ以下の式で求まる.

$$\frac{\partial {h}}{\partial {x}} = \begin{pmatrix} \frac{\partial {h_1}}{\partial {x_1}} & \cdots & \frac{\partial {h_1}}{\partial {x_m}} \\ \vdots & \ddots & \vdots \\ \frac{\partial {h_m}}{\partial {x_1}} & \cdots & \frac{\partial {h_m}}{\partial {x_m}} \end{pmatrix}$$

変数変換によって確率分布も変形され, 同時に定義域が変わることもある

対数変換

シンプルな例として対数変換を説明する. 変換関数$h(x)$は, $h(x) = \log(x)$と定義されるため, ヤコビアンは, $\frac{dh}{dx} = \frac{1}{x}$と求まる. データ$y$が正規分布$g(y)$に従うとすると, $g(y) = \frac{1}{\sqrt{2\pi\sigma^{2}}}\exp^{-\frac{(y-\mu)^2}{2\sigma^{2}}}$となる.

すると, データ$x$を生成する$f(x)$の形状は次の対数正規分布となる.

$$f(x) = g(h(x))\Big|\frac{dh}{dx}\Big| = \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp^{-\frac{(\log{x} - \mu)^2}{2\sigma^{2}}} \frac{1}{x}$$

例として, TSSSパッケージ内のWHARDデータで対数変換を施す.

data(WHARD)
plot(WHARD) # 変換前のWHARDデータ
plot(log(WHARD)) #対数変換後のWHARDデータ

上の図はWHARDデータ(あるハードウェアの毎月の卸売り高を記録したデータ)を変換前（オリジナル）と対数変換後のデータをプロットしたものである.

オリジナルに見られる分散の非定常性が緩和されていること, 変換後もトレンドは維持されている.

Box-Cox変換

対数変換を含むさまざまなべき乗型変換をつくれる, 便利なBox-cox変換を紹介する.

まず定義であるが,

$$\begin{equation} y_{\lambda} = \begin{cases} \frac{x^{\lambda}-1}{\lambda} & \lambda \neq 0 \\ \log{x} & \lambda = 0 \\ \end{cases} \end{equation}$$

と, $\lambda$の値でさまざまな変換を表現できる. 例えば,先ほどの対数変換は$\lambda = 0$の場合のBox-Cox変換に対応する.

すると, $\lambda$はどのように決めるのか?という疑問がでてくる. これは, さまざまな値の$\lambda$の値でのBox-Cox変換を行い,変換後のデータが正規分布に近いかどうかを評価することで, ベストな$\lambda$の値を探すことができる.

Box-Cox変換のパラメータ選択

Box-Cox変換のパラメータ$\lambda$を, AIC(赤池情報量基準）を用いて決定する方法を説明する.
ざっくりいえば, 候補となる$\lambda$のそれぞれで元データをBox-Cox変換し, 変換後のデータの分布が正規分布に近いかどうかをAICで評価するというものだ.

TSSSパッケージには, 自動で$\lambda$の調整と, その値を使ったBox-Cox変換をしてくれるboxcox関数がある.

data(Sunspot)
boxcox(Sunspot) # Box-Cox変換（パラメータ自動調整）

Sunspotデータの変換前, Box-Cox変換($\lambda: 0.4)$後のプロット

Sunspotデータに対しては, 最適なパラメータ$\lambda$の値は, 0.4と決められた.
結果をみると, 上下非対称であった元データが, 変換後は, 上下対称なデータになっていることがわかる.つまりは, 時系列データがガウス型に変換されている.

実は間違っている変換後のAIC比較

boxcox関数では自動でAICに基づき最適な$\lambda$の値を出してくれるが, もし0から処理を実装したり, 別のソフトウェアでAICを計算するときには注意が必要である.

AICに基づき最適な$\lambda$を考える場合, 変換後のデータ$z$を用いて計算した$AIC_z$は, 実は比較してはいけないのだ. 比較するならば, 元データでの$AIC`_{z}$　を求める必要がある.

この値は, 変換関数のヤコビアンで$AIC_{z}$を補正することで求まる.

$$AIC`_{z} = AIC_{z} - 2\log{|\frac{df}{dx}\Big|}$$

以下に, 変換後のデータでの$AIC_{z}$, および, 元データでの$AIC`_{z}$を載せている. $AIC_{z}$では, $\lambda = 0.4$では最小になっていないことに注意が必要だ.

AICの結果

lambda	aic'	LL'	aic	LL	mean	variance
1.00	2360.37	-1178.19	2360.37	-1178.19	4.910866e+01	1.576345e+03
0.90	2335.48	-1165.74	2174.59	-1085.29	3.546734e+01	7.052906e+02
0.80	2313.88	-1154.94	1992.09	-994.04	2.591715e+01	3.200762e+02
0.70	2295.85	-1145.93	1813.16	-904.58	1.917797e+01	1.475255e+02
0.60	2281.75	-1138.87	1638.16	-817.08	1.438208e+01	6.915929e+01
0.50	2271.96	-1133.98	1467.48	-731.74	1.093834e+01	3.303407e+01
0.40	2267.00	-1131.50	1301.62	-648.81	8.441855e+00	1.611155e+01
0.30	2267.51	-1131.76	1141.24	-568.62	6.613792e+00	8.046373e+00
0.20	2274.40	-1135.20	987.23	-491.61	5.260971e+00	4.131024e+00
0.10	2289.01	-1142.51	840.94	-418.47	4.248735e+00	2.192978e+00
0.00	2313.40	-1154.70	704.44	-350.22	3.482608e+00	1.214490e+00
-	-	-	(以下略)	-	-	-

lambda = 0.40 AIC' minimum = 2267.00

ロジット変換

比率データ$x (0 < x < 1)$を, オッズ$\frac{x}{1-x}$の自然対数に変換するものをロジット変換とよぶ. 元のデータが確率（割合）ゆえに, 定義域が$[0,1]$であるが, 変換後は$[-\infty, \infty]$となるのが特徴的である.

$$y = \log{\frac{x}{1-x}}$$

階差（差分）

何らかの影響で平均や分散が一定ではない場合, ある時点からの変動を抽出すると, 実は定常な変化である場合がある.そのため, 階差（差分）を求める前処理も有効な場合がある.

時系列$y_n$の$k$階の階差を$\Delta^{k}y_{n}$と定義すると,

$$\begin{align} \Delta^{1}y_{n} &= y_n - y_{n-1} (\text{１階階差の場合}) \\ \Delta^{2}y_{n} &= \Delta^{1}y_{n} - \Delta^{1}y_{n-1} = y_n - 2y_{n-1} + y_{n-2} (\text{2階階差の場合}) \\ \end{align}$$

のように計算できる.
実際, TSSSパッケージ内のNikkei225のデータの（対数）差分をとると,

data(Nikkei225)
plot(Nikkei225) # 差分前のNikkei225データ

plot(diff(log(Nikkei225)))  # 対数差分後のNikkei225データ

上の図は, 元のNikkei225（日経225平均株価データ）であり,下の図が対数差分をとったものである. もともと平均非定常なデータが, ０を中心としたきれいな平均定常なデータに変換されている.

ここで, 一時点前のデータではなく, さらに前のデータからの差分をとることもできる. 特に,季節性のような周期のあるデータであれば, 季節の長さ$s$分の前のデータと比較することで季節的な変動が除去できる.

$$\Delta_{s} y_{n} = y_n - y_{n-s}$$

WHARDデータでs = 12とした季節階差

data(WHARD)
y <- log(WHARD)
n <- length(WHARD)
period <- 12 # 12ヶ月ごとの周期性
z <- rep(NA, n)
for (i in period + 1:n) {
    z[i] <- y[i] - y[i-period] # 季節階差
}
plot(as.ts(z))

他に, 前年比 $x_n = y_n/y_{n-1}$や, 前年同期比$x_n = y_n / y_{n-p}$という差分のとりかたもある.

移動平均フィルタと移動メディアン・フィルタ

不規則な変動をする時系列データには, トレンドや季節変動などの解釈上意味のある成分のほか, 誤差変動（ノイズ）が含まれている. この誤差変動（ノイズ）が大きすぎると, 元の時系列データから意味のあるパターンを読み取ることは困難である.

そのため, 誤差変動（ノイズ）を除去し, 変化のパターンをきれいにとりだす方法として, 平滑化(スムージング; smoothing)　が有効である.

有名な方法として, 移動平均フィルタと移動メディアン・フィルタがある.それぞれ,ある時点の前後のデータの平均, あるいは中央値で置き変える平滑化方法である.

移動平均フィルタは,

$$t_n = \frac{1}{2k+1}(y_{n-k} + \dots + y_{n} + \dots + y_{n-k})$$

移動メディアン・フィルタは,

$$t_n = median(y_{n-k},\dots,y_{n},\dots,y_{n-k})$$

と表せる.

data(Temperature)
maxtemp <- Temperature

# 移動平均フィルタの場合(SMA関数の使用）
plot(maxtemp, ylim = c(0,40))
x <- SMA(maxtemp, 34) # TTRパッケージのSMA関数を使用する; ウィンドウサイズは34
lines(x, col = 2, lwd = 2) # col: 色, lwd: 太さ

# 移動平均フィルタの場合(0から実装）
plot(maxtemp, ylim = c(0,40))
y <- maxtemp
ndata <- length(maxtemp)
y[1: ndata] <- NA
kfilter <- 17
n0 <- kfilter + 1
n1 <- ndata - kfilter
for (i in n0:n1) {
    i0 <- i - kfilter
    i1 <- i + kfilter
    y[i] <- mean(maxtemp[i0:i1])
}
lines(y, col=2, ylim = c(0, 40), lwd = 2)

# 移動メディアン・フィルタの場合(TTRのrunMedian関数の使用)
data(Temperature)
maxtemp <- Temperature
plot(maxtemp, ylim = c(0,40))
x <- runMedian(maxtemp, 34) # TTRパッケージのrunMedian関数を使用する; ウィンドウサイズは34 (17 + 17)
lines(x, col = 3, lwd = 2) # col: 色, lwd: 太さ

# 移動メディアン・フィルタの場合(0から実装）
plot(maxtemp, ylim = c(0,40))
y <- maxtemp
ndata <- length(maxtemp)
y[1: ndata] <- NA
kfilter <- 17
n0 <- kfilter + 1
n1 <- ndata - kfilter
for (i in n0:n1) {
    i0 <- i - kfilter
    i1 <- i + kfilter
    y[i] <- median(maxtemp[i0:i1])
}
lines(y, col=3, ylim = c(0, 40), lwd = 2)

移動平均フィルタの場合

移動メディアン・フィルタの場合

どちらの方法もウィンドウサイズ（考慮する時系列データの数）によって結果がかわるので, 適切なウィンドウサイズによる平滑化を考える必要がある.

また, 移動平均フィルタと移動メディアン・フィルタの各比較は次となる.

	移動平均フィルタ	移動メディアン・フィルタ
長所	滑らかな推定値	構造変化を正確に検出できる 異常値に頑健
短所	構造変化を正確に検出できない（無視して滑らかにしてしまう） 異常値に敏感	かえって変動を大きくする可能性も

まとめ

本ポストでは, 時系列データとはなにか？, 実際にデータを可視化しながら, 時系列データの性質を整理した. また, モデリングのしやすさのため, 前処理として, 階差や変数変換, 平滑化などのテクニックを確認した. 次回は, 時系列データの重要な量を紹介する.

雑記

時系列分析の方法についてはいろいろな記事が既にあがっているが, 実務的な視点で, しっかりと時系列分析の癖や手順を説明しているものは少ない. なぜ定常なデータにするのか, という疑問がなかなか解消されなかったことを思い出した.

北川先生による時系列解析の講座をみて, 時系列モデルの理論はもちろん, 実務的な解析手順や結果の見方までも丁寧に学べたことは幸せだったと思う.

ただ, 本ポストのボリュームで, まだ講座の1回目程度であるというのは恐ろしい...😱

参考

https://ocw.u-tokyo.ac.jp/course_11477/

北川先生による時系列解析の講座. とてもわかりやすい, かつ, あまり世間で見られる記事では知り得ない高度な話も含めている.

https://www.kyoritsu-pub.co.jp/book/b10003870.html

定常的な時系列モデルの説明を理論的な内容で紹介している. （わたしのような）初心者にはまだ難しく感じたが, いろいろ学んだあとに読み直すと, 理論がきれいかつ簡潔にまとまっていると気づけた.

https://www.kyoritsu-pub.co.jp/book/b10003204.html

Pythonでのコード事例が紹介されている.

脚注

1. データサイエンスとしては, 目の前のデータの性質を深く理解したうえで, 適切なモデルの選択や結果の解釈ができるようになりたい！

2. 例えば, 1980年のアイスクリームの値段についてはその一時点で一度しか得られない. 無理やり大数の法則を働かせるのであれば, "別の世界線"で1980年のアイスクリームの値段データを収集しないといけない