時系列解析　ことはじめ ②

概要

前回のポストでは, 時系列データの前処理について説明した. 適切な前処理を行なった時系列データは, 定常性, 正規性, 線形性と,とても解析しやすい性質をもつ.

本ポストでは, 定常時系列での主な量について説明する.定常時系列データの特徴をより正しく理解するうえでは必須のテクニックである. 主に1変量（１つの時系列データ）で話をすすめるが, 多変量の場合も少し扱うつもり.

https://zenn.dev/akitek/articles/f74753c2e95ce4

定常性時系列

時系列の定常性についてみていく.
まず定常の有無に関係なく, 次の3つの統計量を定義していく. ここで, 時系列データを$y_n$（$n$
は時刻, $N$を総データ数とし, ${y_1,y_2,\dots,y_N}$)とすると

1. 平均: $\mu_n = E[y_n]$
2. 分散: $Var(y_n) = E[(y_n - \mu_n)^2]$
3. 共分散: $Con(y_n, y_m) = E[(y_n - \mu_n)(y_m - \mu_m)]$

さて, 定常性とは「確率構造（性質）が時間変化しないこと」であり, その条件の厳しさに応じて, 強定常と弱定常にわけられる.

強定常と弱定常

時系列データの背景にある確率分布が時間変化しないもので, 任意の時点${y_{n_1},y_{n_2},\dots,y_{n_m}}$と, 任意の時間差（シフト量）$k$に対し, 以下の同時分布がいつでも等しいものをいう.

$$p(y_{n_1},y_{n_2},\dots,y_{n_m}) = p(y_{n_1 + k} ,y_{n_2 + k},\dots,y_{n_m + k})$$

実際はこの強定常を検証することは難しく, またそこまで強い条件が必要ないことも多いため, 次の弱定常性を考えることが多い.

弱定常とは, 平均・分散・共分散が時間によって変化しない（2次モーメントまで不変）ものをいう.つまり,任意の時間差(シフト量）$k \in {0, \pm 1, \pm 2, \dots}$に対し

$$\begin{align} E(y_n) &= E(y_{n+k}) \\ Var(y_n) &= Var(y_{n+k}) \\ Cov(y_n, y_m) &= Con(y_{n+k}, y_{m+k}) \ \end{align}$$

あきらかに強定常のほうが弱定常よりも厳しい条件である.そのため, 平均, 分散, 共分散が存在する強定常であれば, 同時に弱定常も満たす. 逆は必ずしもは成り立たない.

強定常なのに弱定常でない？

例えば, コーシー分布に従う時系列があげられる.

$$p(y_n) = \frac{1}{\pi(y_n^2 + \gamma^{2})}$$

コーシー分布は, 平均や分散が定義されない. そのため, 弱定常性の定義である「平均や分散が存在して...」が成り立たない. そのため, 弱定常性はないのだ.

Rでの実装例

x <- seq(10, -10, length = 101)

# dcauchy関数でコーシー分布の確率密度関数を生成
plot(x, dcauchy(x, location = 0, scale = 0.5, log = FALSE), type = 'l', ylim = c(0,0.7),col = 2)
par(new = T)
plot(x, dcauchy(x, location = 0, scale = 1, log = FALSE), type = 'l',  ylim = c(0,0.7),col = 3,ann=F)
par(new = T)
plot(x, dcauchy(x, location = 0, scale = 2, log = FALSE), type = 'l',  ylim = c(0,0.7),col = 4,ann=F)
par(new = T)
plot(x, dcauchy(x, location = 0, scale = 4, log = FALSE), type = 'l',  ylim = c(0,0.7),col = 5,ann=F)
par(new = T)
# 比較対象として, dnorm関数で（標準）正規分布の確率密度関数を生成
plot(x, dnorm(x, mean = 0, sd = 1),type = 'l',  ylim = c(0,0.7), col = 1 ,ann=F)
legend("topleft", legend=c("Gamma = 0.5", "Gamma = 1", "Gamma = 2", "Gamma = 4","Norm(0,1)"), lty=1, col=c(2,3,4,5,1))
title("Cauchy and Norm")

:::

弱定常時系列の自己共分散関数

改めて弱定常性の性質を整理すると,

1. 平均$E[y_n]$は時刻$n$に依存しない, つまり, $E[y_n] = \mu$とかける
2. 共分散$Cov(y_n, y_m)$は, 時間差$k$だけに依存する. つまり, $Cov(y_n, y_{n-k}) = E[(y_n - E[y_n])(y_{n-k} - E[y_{n-k}])] =E[(y_n - \mu)(y_{n-k} - \mu)] = C_k$

特に自己共分散$C_k$の性質として,1)偶関数, 2)最大は$C_0$,つまり$|C_k| \leq C_0$, 3)半正定値, がなりたつ. 大事な性質なので1つ１つ見ていこう.

---

自己共分散関数は偶関数

偶関数とは, $f(x) = f(-x)$となる性質をいい, 例えば$con(x)$はその一例である.
ここでは, $C_{-k} = C_{k}$となることを述べている.

証明：

$$\begin{align} C_k &= Cov(y_n, y_{n-k}) \\ C_{-k} &= Cov(y_n, y_{n+k}) \\ &= Cov(y_{n+k}, y_{n}) \\ &= Cov(y_{n}, y_{n-k}) = C_{-k} \\ \end{align}$$

式(6)では, 「共分散は要素を入れ替えても不変」, 式(7)では,「弱定常性より同じだけシフトしても問題ない」ことを利用した

---

$|C_k| \leq C_0$

$E(a^2)E(b^2) \geq E(ab)^2$が一般的になりたつことをつかう.

趣味の人向けの導出

任意の$\lambda$に対し, $E(\lambda a + b)^2 \geq 0$が成り立つ.
これをさらに分解すると

$$E(\lambda a + b)^2 \\ = \lambda^2E(a^2) + 2\lambda E(ab) + E(b^2) \\ = E(a^2)\Big( \lambda - \frac{E(ab)}{E(b^2)} \Big)^2 + E(b^2) - \Big( \frac{E(ab)}{E(a^2)} \Big)^2$$

さて, $a,b$の値を固定し, 上式を$\lambda$に関する関数とみなすと, $E(a^2) \geq 0$より, 上式は$\lambda$に関して下に凸な二次関数となる. そして, 最小値は$\lambda = \frac{E(ab)}{E(b^2)}$のときとなるが, 任意の$\lambda$に対し, $E(\lambda a + b)^2 \geq 0$が成り立つためには, 最小値に関し, $E(b^2) - \Big( \frac{E(ab)}{E(a^2)} \Big)^2 \geq 0$である必要があることから求まる.

ここで, $a = y_n - \mu, b = y_{n-k} - \mu$とおくと,

$$\begin{align} E(a^2) &= E(y_n - \mu)^2 = C_0 \\ E(b^2) &= E(y_{n-k} - \mu)^2 = C_0 \\ E(ab) &= E[(y_n - \mu)(y_{n-k} - \mu)] = C_h \\ \end{align}$$

となり, $C_0C_0 \geq C_h \rightarrow C_0 \geq |C_k|$で求まる.

---

半正定値

非負定値性とは以下が成り立つことをいう

$$\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i \alpha_j C_{i-j} \geq 0 \\$$

この性質より, Yule-Walker方程式と呼ばれる便利な解法がつかえる.

導出は,

$$\begin{align} E\Big( \sum_{i=1}^{m}\alpha_i y_{n-i}\Big)^2 &= \sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i \alpha_jE[y_{n-i}y_{n-j}] \\ &= \sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i \alpha_j C_{i-j} \geq 0 \\ \end{align}$$

式(11)から(12)へは,弱定常性よりラグ$i-j = (n - j) - (n - i)$に依存する自己共分散$C_{i-j}$への置き換えによる.

---

弱定常時系列の自己相関関数

自己共分散に関連してもうひとつの量として自己相関関数がある.
これは,

$$R_k \equiv Cor(y_n, y_{n-k}) = \frac{Cov(y_n, y_{n-k})}{\sqrt{Var(y_n}Var{y_{n-k}}}$$

でもとまる. しかし, 定常時系列の場合はとてもシンプルに書き直せる.

$$\begin{align} R_k &= \frac{Cov(y_n, y_{n-k})} { \sqrt{Var(y_n)Var(y_{n-k})} } \\ &= \frac{Cov(y_n, y_{n-k})} { \sqrt{Cov(y_n,y_n)Cov(y_{n-k},y_{n-k}) }} \\ &= \frac{C_k}{\sqrt{C_0C_0}} \\ &= \frac{C_k}{C_0} \end{align}$$

つまりは, 自己共分散関数$C_k$がわかれば自己相関関数$R_k$も計算できてしまう.

各種推定量の求め方

平均, 自己共分散関数, 自己相関関数を時系列データ$y_n$から求める.

• 標本平均

$$\hat{\mu} = \frac{1}{N}\sum_{n=1}^{N} y_n$$

• 標本自己共分散関数

$$\hat{C}_k = \frac{1}{N}\sum_{n=k+1}^{N} (y_n - \hat{\mu})(y_{n-k} - \hat{\mu})$$

• 標本自己相関関数

$$\hat{R}_k = \frac{\hat{C}_k}{\hat{C}_0}$$

注意をいれると, 標本自己共分散関数が求めているのは, 時刻$n = k + 1, ..., N$に渡って, ラグ$k$前のデータとの自己共分散の和をデータ総数$N$で割っている. 平均をとっているように見えるが, 平均だとすれば分母は$N-k$だといえる.

実際, この$\hat{C}_k$は不偏推定量ではない...が, 時系列解析では$N$のまま使われているようだ(『時系列解析』の講座より)

不偏推定量ではない証明

$$E[\hat{C}_k] = E\Big[ \frac{1}{N}\sum_{n=k+1}^{N} y_n y_{n-k} \Big] = \frac{1}{N} E[\sum_{n=k+1}^{N} y_n y_{n-k}] = \frac{N-k}{N}C_k$$

分母を$N$とすると, 推定量$\hat{C}_k$の期待値は, ラグ$k$やデータ数$N$の影響をうけるため, 不変推定量ではない. 逆に, 分母を$N-k$としていれば, 分母・分子が消えるので不変推定量となる.

とりあえず, 推定量の性質を求めておこう

$$\begin{align} E[\hat{C}_k] &= \frac{N-k}{N}C_k \\ Var(\hat{C}_k) &\sim \frac{1}{N}\sum_{j = -\infty}^{\infty} (C_j^2 + C_{j-k}C_{j+k}) \\ Var(\hat{R}_k) &\sim \frac{1}{N}\sum_{j = -\infty}^{\infty} ( R_j^2 + R_{j-k}R_{j+k} - 4R_j R_k R_{j-k} + 2R_j^2 R_k^2) \end{align}$$

白色雑音

平均, 自己共分散, 自己相関, それぞれの推定量の求め方がわかったところで, 時系列解析においてとても重要な時系列データである白色雑音を説明する.

白色雑音は, 期待値（平均）で$E[y_n] = 0$が成り立ち, 時間差$k > 0$で自己共分散$C_k = 0$という性質をもつ. つまりは, ある時点のデータは, その前後の任意の時点のデータとは無相間に変動する時系列である.以下は, 白色雑音の例である.

plot(as.ts(rnorm(1000))) # 白色雑音(標準正規分布から1000個乱数生成）

もう少し厳密に定義すると, 以下の式で示せる.

$$\begin{equation} C_k = \begin{cases} \sigma^2 & k = 0 \\ 0 & k \neq 0 \end{cases} \end{equation}$$

自己共分散が上のような性質のため, 自己共分散の比で求まる自己相関関数は,

$$\begin{equation} R_k = \begin{cases} 1 & k = 0 \\ 0 & k \neq 0 \end{cases} \end{equation}$$

となっている.

さて, 白色雑音はもちろん弱定常時系列であるため, 上述した平均, 自己共分散, 自己相関, それぞれの推定量を計算することができる. いま, 白色雑音の場合でこれらを計算すると,

$$\begin{align} \hat{C}_0 &\sim N(C_0, \frac{2\sigma^4}{N}) \\ \hat{C}_k &\sim N(0, \frac{\sigma^4}{N}) \\ \hat{R}_k &\sim N(0, \frac{1}{N}) \\ \end{align}$$

趣味の人向けの導出

まず推定値$\hat{C}_k$の期待値は, $\hat{C}_k = \frac{1}{N}\sum_{i=k+1}^{N} y_i y_{i-k}$より,

$$\begin{align} E[\hat{C}_k] &= \frac{N-k}{N}C_k = 0 \\ E[\hat{C}_0] &= \frac{N}{N}C_0 = C_0 (= \sigma^2)\\ \end{align}$$

そして, 推定値$\hat{C}_k$の分散は, 一般的に$Var(\hat{C}_k) = \frac{1}{N^2}\sum_{i = k+1}^{N} \sum_{j = k+1}^{N} (C_{i-j}^2 + C_{i-j-k}C_{i-j+k})$とかけることから,

$$\begin{align} Var(\hat{C}_0) & = \frac{1}{N^2}\sum_{i = 1}^{N}\sum_{j = 1}^{N}(C_{i-j}^2 + C_{i-j}^2) = \frac{2}{N^2} N C_0^2 = \frac{2\sigma^2}{N} \\ Var(\hat{C}_k) &= \frac{1}{N^2}\sum_{i = k+1}^{N}\sum_{j = k+1}^{N}(C_{i-j}^2 + 0) = \frac{N-k}{N^2} C_0^2 \cong \frac{\sigma^4}{N} \\ \end{align}$$

式(28)について, $k \neq 0$の場合, $(i - j - k)$か$(i - j - k)$のどちらかは必ず0以外の値をもつことから, $C_{i-j-k}C_{i-j+k} = 0$となることより

白色雑音での自己相関関数の推定量の分散は正規分布$N(0, \frac{1}{N})$となる.^[1]

逆にいえば, 扱う時系列データが白色雑音であれば自己相関関数は$\pm 1/\sqrt{N}$ に収まるはずである. この性質を使い, 簡単な白色性評価に用いることができる.

例えば, さきほどの白色雑音のデータに対し, 自己相関係数をプロットし, データサイズ$N$で定まる$\pm 1/\sqrt{N}$を図中に表示したものが以下である.

wn <- as.ts(rnorm(1000))
wn <- unicor(wn, plot = FALSE) # 一変量での平均, 自己共分散, 自己相関を同時に求める

# ホワイトノイズの自己相関関数: 
plot(wn$acov, ylim = c(-1,1),col=1, type = "l")
par(new = T)
# 1/√N のプロット（赤色）
plot(wn$acov.err, ylim = c(-1,1),col=2, type = "l", ann = F)
par(new = T)
plot(-wn$acov.err, ylim = c(-1,1),col=2, type = "l", ann = F)

この図は, 横軸を時間差$k$,縦軸が自己相関関数の推定量$\hat{C}_k$を示す.
赤色の基準線($\pm 1/\sqrt{N}$)の領域に対し, ある程度収まっていることがわかる.^[2]

次に, Sunspotデータで同じことを行うと,

data(Sunspot)
unicor(Sunspot)

まず自己相関関数だけをみると, とてもわかりやすい周期性があり, そして, 減衰が遅いことがわかる. つまりは, とても相関が強くなるラグ（周期）があり, そして, ある程度離れた時点でのデータとも影響が大きいという性質がみえてくる.

そして, 赤色の基準線($\pm 1/\sqrt{N}$)の領域に対し, ほとんど収まっていないことがわかる. つまりは, このデータは純粋な白色雑音ではない,,,と簡単な方法で確認できる^[3].

周期性の可視化

つぎに, 自己共分散関数を使ったスペクトル表現に表現する. スペクトル表現では,時系列データを時間軸から周波数軸に変換することで, 時間にとらわれず, 時系列データに含まれている周期的な動きをあぶりだすことができる.

...といっても, よくわからないという意見もあるだろう.(スペクトル表現に慣れない私もそのひとり）

なので, 正確性を欠く可能性を承知に, 次のように理解できる言葉に置き換えたい.

つまり, 時間に沿って変動する時系列データは, 実は複数の周期性のある関数（波; cos波やsin波みたいな）が合成されて実現しているとみなす. スペクトル表現では, この周期の大小から, 時系列データを構成する複数の波を分解する作業だといえる.

スペクトル表現にすることで, どのような周期をもつ波が, どれくらいの大きさ含まれているのか, が解析できるようになる.

スペクトル表現

さて, この（パワー）スペクトルは, 実は前述した自己共分散関数$C_k$のフーリエ変換で求めることができる. つまり,$f$を周波数として使うと,

$$p(f) = \sum_{k = -\infty}^{\infty} C_k e^{-2\pi ikf}, -\frac{1}{2} \leq f \leq \frac{1}{2}$$

と表せる. 一変量定常時系列では, 自己共分散$C_k$は偶関数であることから,

$$\begin{align} p(f) &= \sum_{k = -\infty}^{\infty} C_k e^{-2\pi ikf} \\ &= C_0 + \sum_{k = 1}^{\infty} C_k (e^{-2\pi ikf} + e^{2\pi ikf}) \\ &= C_0 + \sum_{k = 1}^{\infty} C_k (\cos{2\pi kf}-i\sin{2\pi kf} + \cos{2\pi kf}+i\sin{2\pi kf}) \\ &= C_0 + 2\sum_{k = 1}^{\infty} C_k \cos{2\pi kf} \end{align}$$

と書き直せる.式からわかるが, このパワースペクトルは偶関数である^[4]

また, 逆フーリエ変換により, パワースペクトルから自己共分散$C_k$に戻すこともできる.

$$\begin{align} C_k &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)e^{2\pi ikf}df \\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)\cos(2\pi kf)df \end{align}$$

趣味の人向け逆フーリエ変換の導出

$$\begin{align} C_k &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)e^{2\pi ikf}df \\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)(\cos(2\pi kf) + i\sin(2\pi kf) )df \\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)(\cos(2\pi kf) + i\sin(2\pi kf) )df \\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)\cos(2\pi kf)df + \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f) i\sin(2\pi kf)df \\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)\cos(2\pi kf)df + \int_{0}^{\frac{1}{2}} p(f) i\sin(2\pi kf)df + \int_{-\frac{1}{2}}^{0} p(f) i\sin(2\pi kf)df\\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)\cos(2\pi kf)df + \int_{0}^{\frac{1}{2}} p(f) i\sin(2\pi kf)df - \int_{\frac{1}{2}}^{0} p(-f) i\sin(-2\pi kf)df\\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)\cos(2\pi kf)df + \int_{0}^{\frac{1}{2}} p(f) i\sin(2\pi kf)df + \int_{0}^{\frac{1}{2}} - p(f) i\sin(2\pi kf)df\\ &= \int_{-\frac{1}{2}}^{\frac{1}{2}} p(f)\cos(2\pi kf)df \end{align}$$

ピリオドグラム

パワースペクトル$p(f)$が, 自己共分散$C_k$のフーリエ変換で求まることがわかった.
実際の運用では次のように,

1. まず時系列データに対し, 何らかの時系列モデルで学習し, 自己共分散の推定量である$\hat{C}_k$をえる.
2. $\hat{C}_k$をフーリエ変換して, パワースペクトルの推定量$\hat{p}(f)$をえる.

時系列データは有限数なので, 作成される自己共分散の推定量も有限個となる.そのため, パワースペクトルの推定量$\hat{p}(f)$は離散的に求まる.

このパワースペクトルの推定量$\hat{p}(f)$を, 周波数$f$を横軸にプロットしたものをピリオドグラムとよぶ. ピリオドグラムをみることで, どのような周波数成分（周期性の波）が含まれているのが把握できるようになる.

data(WHARD) # TSSSパッケージ内のWHARDデータ
period(WHARD, window = 0) # ピリオドグラムの作成：window = 0は平滑化なし(後述）

上記は, WHARDデータ（あるハードウェアの毎月の卸売り高を記録したデータ）のピリオドグラムを作成したものである. 一定間隔でスペクトルが高まっていることがわかる. 実は, 12個間隔で高まる点があることから, 12ヶ月の周期性があるようにみえる.

パワースペクトル, ピリオドグラムをみれば, どのような周期性の波が含まれているのか, どれくらいの強度なのかを分析できるようになるので, 時系列データの分析時にはぜひ作ってみるとよい.

白色雑音のピリオドグラム

例として, 白色雑音のピリオドグラムを求めてみる. 白色雑音の自己共分散関数$C_k$は, $k=0$のとき,$C_0 = \sigma^2$で, それ以外では, $C_k = 0$である.

これを, 自己共分散$C_k$のフーリエ変換に代入すると,

$$\begin{align} p(f) &= C_0 + 2\sum_{k = 1}^{\infty} C_k \cos{2\pi kf} \\ &= C_0 = \sigma^2 \end{align}$$

となり, 周波数に関係なく一定値をとる.^[5] スペクトル表現は時系列データを周期成分に分解するというものだが, 白色雑音の場合は, さまざまな周期性の波が均等に含まれているということになる. 白色光は, いろいろな周波数の可視光線の合成波だ, ということからもイメージできるかも

ピリオドグラムの性質

ここでピリオドグラムの性質として, 不偏推定量だが, 一致推定量ではないことを述べておく

$$\begin{align} \lim_{n \rightarrow \infty}E[\hat{p}(f)] &= p(f) \text(漸近的に不偏)\\ \lim_{n \rightarrow \infty}\hat{p}(f) &\neq p(f) \text(一致推定量ではない)\\ \end{align}$$

これはどういうことかというと, イメージとして, データ数が増えるほど,期待値的には真のスペクトル密度関数$p(f)$に近づくが, しかし, 真のスペクトル密度関数$p(f)$と同じ形状にはならないという.実際, パワースペクトルのスペクトル密度の値はカイ２乗分布に従い, データサイズによって分布は変化しない.

例として, 白色雑音から作成した時系列データを, 200個, 800個, 3200個に増やしていって, それぞれピリオドグラムを作成した.

r <- rnorm(3200) # 乱数発生
spec.pgram(as.ts(r[1:200]), log = "yes")
spec.pgram(as.ts(r[1:800]), log = "yes")
spec.pgram(as.ts(r[1:3200]), log = "yes")

図をみると,　データサイズが増えると, スペクトルの平均は1に近づくが, 分散は特に改善されていない.^[6]

より良いピリオドグラムを作るためには, 全体のデータを部分区間ごとにわけ, 各区間のデータだけでピリオドグラムをもとめ, 最後に, それらを平均化するといったテクニックがある.これを使えば, データ全体でピリオドグラムを作った場合のスペクトルの分散を$l$(部分区間の数)だけ小さくできる.

$$\hat{p}_j = \frac{1}{l} \sum_{k=1}^{l}p_j^{(k)}$$

これを使って, 3200個の白色雑音データを, 1)3200個フルで使った場合,2)1600個ごとに2分割した場合, 3)800個ごとに4分割した場合, 4)200個ごとに16分割した場合, を比較する.
(Rコードは, こちらのサイトのものをお借りした)

図をみるとわかるが, 分割数をふやすほど, 分散が小さくなっていることがわかる.

他にも, ピリオドグラムの平滑化も行うことができる.スペクトルは本来なめらかであるという前提から, 平滑化しようというものである.

詳細は省くが, Rであれば, TSSSパッケージのperiod関数のwindowオプションで平滑化したピリオドグラムを作れる.

多変量時系列

最後に, 多変量時系列の場合について簡単に説明する. 時系列解析における多変量($m$変量） とは, 同じ時刻で複数の観測点（センサー）で計測された$m$個の時系列データのことである.

$$y_n = \begin{pmatrix} y_n(1) \\ y_n(2) \\ \vdots \\ y_n(m) \\ \end{pmatrix} y_n \in R^{m}, n = 1, \dots, N$$

一例として, TSSSパッケージ内にあるHAKUSANデータ(船舶の方向角速度(yaw rate), 横揺れ（rolling), 縦揺れ(pitching)および舵角（rudder angle) を1秒ごとに記録したもの; 4変量）を扱う.

data(HAKUSAN)
plot(HAKUSAN[,1])

図には4変量の時系列データを,時間軸を合わせてプロットした.

1変量で求めた自己共分散関数, 自己相関関数から拡張し, 二つの時系列データ間での相互共分散関数, 相互相関関数を定義しよう.

定常時系列の相互共分散関数

$$\begin{align} C_k(i,j) &= Cov(y_n(i), y_{n-k}(j)) \\ &= E[(y_n(i) - \mu_i)(y_{n-k}(j)-\mu_j)] \end{align}$$

上の式は, $i$番目の時系列データと$j$番目の時系列データで時間差$k$での相互共分散関数を定義している. つまり, 別の時系列データにおける時間差$k$でのデータとを比べている.

そして, $C_k(i,j)$を$m$個全ての時系列データ間で求めた相互共分散関数行列を以下で定義.
ここで, $i = j$の場合は, $i$番目の時系列データの自己共分散関数である.

$$C_k = \begin{pmatrix} C_k(1,1) & \cdots & C_k(1,m) \\ \vdots & \ddots & \vdots \\ C_k(m,1) & \cdots & C_k(m,m) \end{pmatrix}$$

注意として, 多変量の場合には, 相互共分散関数は偶関数ではないことを覚えておく.
一般に,

$$C_{-k} = C_k^{\top}$$

とかける.

証明

$$\begin{align} C_{-k}(i,j) &= Cov(y_n(i), y_{n+k}(j)) \\ &= Cov(y_{n+k}(j), y_n(i)) \\ &= Cov(y_n (j), y_{n-k}(i)) \\ &= C_k(j,i) \end{align}$$

より
なお, 式(31)(32)では, 共分散の要素入れ替えによる結果不変性を利用. 式(33)(34)にかけては, 定常性時系列における時間差$k$を維持したシフトでは結果が不変であることを利用.

定常時系列の相互相関関数

$$R_k(i,j) = \frac{C_k(i,j)}{\sqrt{C_0(i,i)C_0(j,j)} }$$

上の式は, $i$番目の時系列データと$j$番目の時系列データで時間差$k$での相互相関関数を定義している. つまり, 別の時系列データにおける時間差$k$でのデータとを比べている.

相互共分散関数と同じく, $R_k(i,j)$を$m$個全ての時系列データ間で求めた相互相関関数行列を以下で定義.ここで, $i = j$の場合は, $i$番目の時系列データの自己相関関数である.

$$R_k = \begin{pmatrix} R_k(1,1) & \cdots & R_k(1,m) \\ \vdots & \ddots & \vdots \\ R_k(m,1) & \cdots & R_k(m,m) \end{pmatrix}$$

HAKUSANデータで相互相関関数を求める以下のようになる(TSSSパッケージのcrscor関数が便利）

data(HAKUSAN)
y <- as.matrix(HAKUSAN[, 1:4])
crscor(y)

可視化してみると, 時系列データ間である時間差$k$（の倍数）での相互相関係数が強まったり, 弱まったりしていることがわかる.例えば, 横揺れと舵角の相互相関は安定して大きい. また時間遅れがどちらか（どちらが先行して他のデータに影響を与えているのか,）もわかってくる

成分分解

どれも不規則な変動しているが, 実は時系列データ間で何らかの相互作用が働いている可能性がある.例えば, 方向角速度の変動が横揺れの変動に影響を与えている,,,かもしれない. そして面白いことに, その逆も考えられる.

多変量では,個々の時系列データはデータ間での複雑な相互作用の結果生じているため, うまく成分を分解する必要がある. 詳細な説明は『時系列モデルの推定』のポスト（予定）で説明しようと思うが, ここでは簡単に成分分解を行なった結果だけ示す.

# 多変量でのクロススペクトル(TSSSパッケージ）
data(HAKUSAN)
yy <- as.matrix(HAKUSAN[, 1:4])
nc <- dim(yy)[1]
n <- seq(1, nc, by =2)
y <- yy[n, ]
z <- marfit(y, lag = 20)
marspc(z$arcoef, v = z$v, col = c(1,2,3,4,5))

上図はHAKUSANデータに対し, 各成分(船舶の方向角速度(yaw rate), 横揺れ（rolling), 縦揺れ(pitching)および舵角（rudder angle))のパワースペクトルの計算から求めたパワー分解・パワー寄与率である.

１変量でのパワースペクトルをみれば, どの周波数成分の波が強く現れているかを判断できた. しかし, 多変量の場合, ある時系列データでのある周波数成分がそのデータ自身のみで構成されているわけではない. 実際, 他の成分（時系列データ）からの影響も含まれている.

そのため, パワー分解・パワー寄与率では, ある時系列データのスペクトルが, 自身あるいは他の時系列データからの影響をどれくらいの割合でうけているかを把握することができる.

図では, 例えば, 舵角（D)のパワースペクトルからは, 舵角（D)のデータは, 周波数 0.13-0.14の波が主として含まれていることがわかる. しかし, 右のパワー寄与率をみると, 方向角速度（A)の影響も大きくなっていることがわかる. つまりは, 舵角（D)のデータは, 方向角速度（A)の寄与が大きいということがパワー分解にからわかってきそうだ.

まとめ

本ポストでは（弱）定常時系列データの性質について整理した. 時系列データから自己共分散関数$C_k$が計算できれば, 自己相関関数$R_k$やピリオドグラムを求めることができる.

次は, いよいよ時系列モデルを導入し, パラメータの推定やモデル選択などに触れていこう

雑記

今回はやや理論的な話になってしまった. 元々理論や数式が得意ではないので, 間違っている箇所もあるかもしれないが, 優しく指摘してもらえると幸いである.
しかし, 理論から知ることで, 時系列解析の勘所みたいなのがいろいろ学べてよいなと思う.
（数式をMarkdownで書くのが大変だが💦）

参考

https://chaosmemo.com/entry/2022/07/11/232004

ピリオドグラムの平均化処理のコードをお借りした

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

北川先生による時系列解析の講座. 今回も説明内容をかなり参照させていただいた.

https://www.asakura.co.jp/detail.php?book_code=12247&srsltid=AfmBOoq4qLoIG3p4c67Xiy-ptZ7Z_bxJ4X6lwItROeO2O_nQjrpTOvJ5

赤池先生・北川先生編著の時系列解析の実応用事例集. 時系列解析の工夫や結果の解釈など, 実務応用でのエッセンスがたくさん得られる. いつかこの本を紹介できればと思う.

https://jasp.ism.ac.jp/ism/TSSS/

本ポストでのRコードやデータセットの提供パッケージTSSS

https://mathlandscape.com/cauchy-distrib/

コーシー分布の性質を知りたい人向けに

脚注

1. 本来は, 白色雑音は$R_k = 0$であるが,なんらかのノイズが含まれているため, 分散の値は0となっていない.

2. ここでは, データ数が1000なので,赤色の基準線はおおよそ $\pm 0.03$に引かれている.

3. 多くの時系列解析ソフト, パッケージにはこの白色性評価のため判定ラインが一緒にプロットされている

4. $C_k$および$cos(x)$は偶関数で, $p(f)$はそれらの線形和であるため

5. グラフで描画すると, $p(f) = C_0$となる, 水平線を引くようなイメージ

6. 白色雑音のスペクトルは, 周波数に関わらず$\sigma^2 = 1$(今回は標準正規分布なので)で一定より