表記の本を読んでまとめた、個人勉強用メモ。
大事そうな表現・ことば、あるいは知っておくべき事項と思ったことを抜粋してまとめました。
※【2017-06-22更新】12章まで更新しました。
１章 時系列データの解析とその準備
定常時系列
一定の確率モデルで表現できる、一見不規則な現象。平均・分散などが時間的・周期的に変化しないもの。変化してしまうのは非定常時系列となる。
ガウス型時系列
分布がガウス分布に従うもの。変動の仕方に上下非対称性があるなど、ガウス分布に従わない場合は非ガウス型時系列という。
~時系列データの前処理~
変数変換
値の増加に伴って変動の大きさも増加する場合、分散が変動しているので非定常となる。その場合はzn=log(yn)として分散が一様になるように調整する。
また、値があることが起こる確率や割合のように(0,1)の範囲のみを取る場合はzn=log(yn/(1-yn))というロジット変換を行い(-∞,∞)の値に変換する。詳しくは４章で。
差分
トレンドや季節変化がある場合は変化の微分、二階微分、...をして一次、二次のトレンド成分を特定して除く。
前期比、前年同期比
１つ前のデータと１年前の同じ時期とのデータの日のこと。単調に増加していくようなトレンドをもつデータに対しては、これによって前者は成長率が、後者は周期成分を除去したデータが得られる。
移動平均
データの変動が激しすぎて見にくい場合は移動平均をとってグラフを滑らかにする。これをすると分散が変化するので注意。
２章 共分散関数
※定常な時系列に行う。
自己共分散
時系列の過去の変動との比較。
値が大きいほど、ラグkとの関係が強いことが分かる。
$$ 共分散：Cov(y_n,y_{n-k})=E\left[(y_n-µ)(y_{n-k}-µ)\right]=\frac{1}{N}\sum_{n=k+1}^{N}(y_n-µ)(y_{n-k}-µ) $$
$$ 分散：Var(y_n)=E\left[(y_n-µ)(y_{n}-µ)\right] $$
のとき、（標本）自己相関関数は
$$ R_k=\frac{Cov(y_n,y_{n-k})}{\sqrt{Var(y_n)Var(y_{n-k})}}\frac{1}{N}\sum_{n=k+1}^{N}(y_n-µ)(y_{n}-µ) $$
で求められる。
＊ホワイトノイズ
自己共分散関数$C_k$が$k=0$で$0$、それ以外で$σ^2$の場合、ホワイトノイズと呼ばれる。
相互共分散
他の時系列の変動との関連の強さ。
分散共分散行列と同じ。
時系列$i,j$はそれぞれ別の変数の時系列データとして、
$$ 相互共分散：C_k(i,j)=\frac{1}{N}\sum_{n=k+1}^{N}\left[y_n(i)-µ(i)\right]\left[y_{n-k}(j)-µ(j)\right] $$
より、（標本）相互相関関数は
$$R_k(i,j)=\frac{C_k(i,j)}{\sqrt{C_0(i,i)C_0(j,j)}}$$
で求められる。
３章 スペクトルとピリオドグラム
スペクトル
自己共分散関数を周期成分に分ける解析法。
$$ \sum^∞_{k=-∞}\left|C_k\right|<∞ $$
が成立するときフーリエ変換が実行でき、パワースペクトル密度関数（スペクトル）が計算できる。
自己共分散関数をフーリエ変換することで、どのぐらいの周期に相関があるのか（＝元のデータがどのぐらいの周期なのか）がわかる。
ピリオドグラム
スペクトルが自己共分散関数を使っている一方で、標本自己共分散関数を使うとピリオドグラムになる。
ピリオドグラムは非連続な標本値を使って作るので、定義域が離散的になる。この定義域を連続に拡張したものが標本スペクトル。
基本的に観測データは連続ではなく離散的な値となるはずなので、普通「スペクトルを出す」といえば、この「ピリオドグラム」を指しているはずである。
ピリオドグラム -------> 標本スペクトル -------> スペクトル
定義域が連続 データが標本でなく全部
【注意】ピリオドグラムの「期待値」はデータ数の増加でスペクトルの「期待値」に近づくが、ピリオドグラムがスペクトルに近づくわけではない。つまり、データが多いからといって真のスペクトルを推定できるとは限らない！
【注目】一方で標本自己相関関数はデータ数の増加で真の自己相関関数に近づく
ピリオドグラムの平均化と平滑化
前項の問題点であった「データが多いからといって真のスペクトルを推定できるとは限らない」を解決して、真のスペクトルを標本で表現する方法を考える。
ここでは$0〜N/2$までフーリエ変換するのではなく、適当な値$L$を定め、$f_j=\frac{j}{2L}$を用いて
$$p_j=C_0+2\sum^{L-1}_{k=1}C_kcos(2\pi kf_j)$$
とフーリエ級数展開する。
この$p_j$を生スペクトルとよぶ。
これはデータ数が増えれば増えるほど、スペクトルに近づく事が多い。
それは生スペクトルは各L個のデータのピリオドグラムの平均みたいなだからである。 （スペクトルの平均化）
ピリオドグラムを平滑化するためにはスペクトルウィンドウというものを使って、隣の値と平均のようなことをする。
そうすると急激な変動が消え、よい推定が可能となる一方、大事なスペクトルが消えることもあるので注意。
４章 モデリング
4.6 データの変換
平均値の増加とともに分散が増大 したり 分布が著しく偏ったり するデータは簡単なモデル（共分散とか）で分析するのはNG。
そういう場合はデータの変換が必要となる。
ここではデータの変換手法の一つとして、対数変換をふくんだ一般形であるBox-Cox変換を紹介する。
まずはベースとなる最尤法について説明する。
最大対数尤度
対数尤度の最大値。
尤度というのはもっともらしさの意。
パラメタリックモデルのパラメータ最適化に用いられる。
統計におけるモデル決定は、「今得られた有限個の観測データ（$y_1,y_2,...,y_N$）から、真の分布を推定する」ことである。
分布の形が既存のモデル（正規分布など）とわかっている場合に、最も尤もらしいパラメータをデータから逆算するのが、このパラメータの最適化である。
最適化は、以下で示される対数尤度関数を解析し、最大値を算出することで行われる。
対数尤度関数は「θの場合、yが起こる確率の和」のようなイメージ。
$$ l(θ)=\left \{
\begin{array}{ll}
\sum^N_{n=1}log f(y_n|θ) & （観測が独立の場合） \\
log f(y_1,y_2,...,y_N|θ) & （一般の場合）
\end{array}
\right. $$
ただし、対数尤度関数は一般には極めて複雑な式になり、解析的に最大値を算出するのは困難。
そこで一般的には疑似ニュートン法などの数値計算手法で近似的に計算する。
AIC（赤池情報量規準）
$$AIC=-2(最大対数尤度)+2(パラメータ数)$$
「次に示す対数尤度を大きくしつつもパラメータの数が少ない方が良い」という、最適なモデル決定のための指標。
小さいほどよいモデルであるといえる。
Box-Cox変換
Box-Cox変換を用いると データの分散が絶対値の増加と共に増加する場合 や 変動に上下非対称性があった場合 にそれらを軽減することができる。
$y_n$を元データとした場合、
$$ z_n = \left \{
\begin{array}{ll}
λ^{-1}(y_n^λ-1) & (λ \neq 0) \\
logy_n & (x = 0)
\end{array}
\right. $$
で示される$z_n$がBox-Cox変換である。
この式中の$λ$はAIC（赤池情報量規準）で最適化（AICが最小となる場合のλを採用する）して算出する。
AICで、対数尤度を算出するときには「θのときに$y_n$が生じる確率」を考えたが、
ここでは「λのときに変換前のデータの確率密度関数」を考える。
したがって、λを使って変換前のデータの確率密度関数を表現する必要がある。
元のデータを$y_n$、Box-Cox変換を$z_n=h(y_n)$、変換後の密度関数を$f(z)$とすると、
元のデータ$y_n$の密度関数は、
$$g(y)=\left|\frac{dh}{dy}f\left(h(y)\right)\right|$$
となる。
これを使ってAICないしはその中の最大対数尤度を計算する。
5章 最小二乗法
$$ y_n=\sum_{i=1}^ma_ix_{ni}+\epsilon_n $$
で示される 回帰モデル はパラメータの最尤推定を行うことで解けるが、それは以下の式で示される二乗誤差を最小化することと等しい。
$$ \sum_{n=1}^N\left(y_n-\sum_{i=1}^ma_ix_{ni}\right)^2 $$
(6~10章 省略)
11章 トレンドの推定
トレンドは多項式回帰モデルを最小二乗法で解くことで簡単に特定できる。
何次の多項式にするかは、AICを計算することで決定する。（詳しくは本書の５章を参照）
ただし、これを用いると全データのなめらかな傾向を示すので、
急な傾向の変化などを捉えることが出来ない。
これを改善するものにトレンド推定法がある。
トレンド成分モデル
以下のような差分方程式を トレンド成分モデル と呼ぶ。
$$Δ^kt_n=v_n$$
ただし、$Δt_n=t_n-t_{n-1}$で定義される差分、
$v_n$は平均$0$、分散$τ^2$の正規分布に従うノイズであるとする。
これは$k=1$の場合、
$$ t_n=t_{n-1}+v_n $$
で示され、$v_n$がたまたま0のときは$t_n=t_{n-1}$で一定となる値を取る。
また、$k=2$の場合は、
$$ t_n=2t_{n-1}-t_{n-2}+v_n $$
となり、$v_n=0$の時に直線となる。
一般に、
$$ t_n=\sum_{i=1}^kc_it_{n-1}+v_n $$
と表せるので
$$ x_n= \left[
\begin{array}{c}
t_n \\
t_{n-1} \\
... \\
t_{n-k+1}
\end{array}
\right],\space F= \left[
\begin{array}{c}
c_1 & c_2 & ... & c_k \\
1 & 0 & ... & 0 \\
& ... & \\
0 & 0 & 1 & 0
\end{array}
\right], G= \left[
\begin{array}{c}
1 \\
0 \\
... \\
0
\end{array}
\right], H=\left[
\begin{array}{c}
1 & 0 & ... & 0
\end{array}
\right] $$
とすると、
$$ x_n=Fx_{n-1}+Gv_n $$
$$ t_n=Hx_n $$
と表わせ、これはAR（自己回帰）モデルとみなせる。
トレンドモデル
実際の時系列はトレンド成分モデル＋観測ノイズモデルで
$$ y_n=t_n+w_n $$
となる。（$t_n$が単なる多項式で、それを最小二乗法で解くのが本章冒頭）
$t_n$を前節で述べた$Δ^kt_n=v_n$とすると先程の$x_n,F,G,H$を用いて、
$$ x_n=Fx_{n-1}+Gv_n $$
$$ y_n=Hx_n+w_n $$
と表せる。
これは$v_n,w_n$の分散を決めれば、カルマンフィルタとカルマンスムーザーのアルゴリズムを用いて、$x_n$の平滑値:$x_{1|N},...,x_{N|N}=x_{n|N}$を算出できる。
トレンドの平滑値$t_n$は$Hx_{n|N}$で算出できる。
$k,v_n,w_n$はAICを計算することで最適な値を特定できる。
12章 季節調整モデル
季節成分モデル
前章と同様に季節のトレンドを推定する。
以下の式を 季節成分モデル と呼ぶ。
$$ \left(\sum_{i=0}^{p-1}B^i\right)^ls_n=v_{n2}$$
$$ ただし\space v_{n2}〜N(0,τ2^2),\space B^ps_n=s{n-p} $$
これは$l=1$のとき、
$$\sum_{i=0}^{p-1}B^is_n=\sum_{i=0}^{p-1}s_{n-i}=v_{n2}$$
となり、季節変化一周期分の和が$v_{n2}$であることを示している。
これを
$$ \left(\sum_{i=0}^{p-1}B^i\right)=1-\sum_{i=1}^{l(p-1)}d_iB^i $$
とすると、
$$ x_n= \left[
\begin{array}{c}
s_n \\
s_{n-1} \\
... \\
s_{n-k+1}
\end{array}
\right],\space F= \left[
\begin{array}{c}
d_1 & d_2 & ... & d_k \\
1 & 0 & ... & 0 \\
& ... & \\
0 & 0 & 1 & 0
\end{array}
\right], G= \left[
\begin{array}{c}
1 \\
0 \\
... \\
0
\end{array}
\right], H=\left[
\begin{array}{c}
1 & 0 & ... & 0
\end{array}
\right] $$
として、前章のトレンド成分モデルと同様にARモデルとみなすことができる。
ただし、$l=１$では$d_i=-1$、$l=2$では$d_i=-i-1 (i\leq p-1), d_i=i+1-2p(p\leq i\leq 2(p-1))$となる。
標準的季節調整モデル
前節の季節成分モデル$s_n$を用いて、
$$y_n=t_n+s_n+w_n$$
の３つの成分に分解する。
このモデルは、
トレンド成分モデルの状態空間マトリックスを$F_1,G_1,H_1$、
季節成分モデルの状態空間マトリックスを$F_2,G_2,H_2$とすると、
$$F=\left[
\begin{array}{c}
F_1 & O \\
O & F_2
\end{array}
\right],\space G=\left[
\begin{array}{c}
G_1 & O \\
O & G_2
\end{array}
\right],\space H=\left[
\begin{array}{c}
H_1 & H_2
\end{array}
\right]$$
$$ x_n=Fx_{n-1}+Gv_n $$
$$ y_n=Hx_n+w_n $$
とすることで、前章と同様に解くことができる。
定常AR成分を含む分解
前節の方法で成分分解をすると、トレンドが上下に細かく変動してしまう場合がある。
この状態のまま将来予測を行うと、予測の分散が非常に大きくなる（＝信頼性が低い予測）。
トレンドの分散を小さくするためには トレンド成分モデル のノイズ$v_{n1}$の分散を最尤推定で求めた値より小さくすればよいが、
そうするとこんどは日々の推定を精度良く行えなくなる。
そこで、こうした場合は以下のような式に分解を行う。
$$y_n=t_n+s_n+p_n+w_n$$
ただし、$p_n$は
$$p_n=\sum_{i=1}^{m_3}a_ip_{n-i}+v_{n3}$$
に従う定常AR(Auto Regressive)成分モデルと呼ばれる。
これはこれまでの$t_n$や$s_n$のような長期的な傾向ではなく、短期的な循環変動成分を示す。
この時、$a_i<1$であることからこのARは定常であるといえる。
$$ x_n= \left[
\begin{array}{c}
p_n \\
p_{n-1} \\
... \\
p_{n-m_3+1}
\end{array}
\right],\space F_3= \left[
\begin{array}{c}
a_1 & a_2 & ... & a_{m3} \\
1 & 0 & ... & 0 \\
& ... & \\
0 & 0 & 1 & 0
\end{array}
\right], G= \left[
\begin{array}{c}
1 \\
0 \\
... \\
0
\end{array}
\right], H=\left[
\begin{array}{c}
1 & 0 & ... & 0
\end{array}
\right], Q_3=τ_3^2 $$
とする。
この時、$Q$はシステムノイズを示している。
カルマンフィルタの手法により、
$$F=\left[
\begin{array}{c}
F_1 & O & O \\
O & F_2 & O \\
O & O & F_3
\end{array}
\right],\space G=\left[
\begin{array}{c}
G_1 & O & O \\
O & G_2 & O \\
O & O & G_3
\end{array}
\right],\space H=\left[
\begin{array}{c}
H_1 & H_2 & H_3
\end{array}
\right],\space
Q=\left[
\begin{array}{c}
\tau^2_1 & O & O \\
O & \tau^2_2 & O \\
O & O & \tau^2_3
\end{array}
\right] $$
を解くことで$y_n$を分解することができる。
曜日調整項を含む分解
「平日は多いが土日は少ない」といった曜日に関わるトレンドを取り出すために、
以下の式で定義される曜日調整項を導入する。
$$td_n=\sum_{i=1}^7\beta_{ni}d^*_{ni}$$
ただし、$\beta_{ni}$ は $i$ 曜日の増減が $y_n$ に及ぼす影響の係数、
$d^*_{nk}$は$n$月に含まれる$k$曜日の日数を指す。
$\beta_{ni}$は月によらず一定の場合が多いため、
以下のように状態空間モデルを定義することでこれまでと同様に分解することができる。
$$
F_{n4}=G_{n4}=I_6,\space
H_{n4}=[d_{n1},...,d_{n6}],\space
x_{n4}=\left[
\begin{array}{c}
\beta_{n1}\\
...\\
\beta_{n6}\\
\end{array}
\right]
$$
ただし、$d_{ni}$は各曜日の日数から土曜日の日数を引いたもの。
これまでのをすべて含めて、分解して
$$y_n=t_n+s_n+p_n+td_n+w_n$$
とする。