はじめに
前提として下記の知識があると良いです。
以下の知識は必要とされないですが、時系列に潜在変数を持たせるときには下記の知識の応用が用いられます。
今回の記事では時系列データ $X _ t$ を扱いますが、こいつらを直接観測できている前提のモデルを利用していきます。すなわち、問題設定として、時系列データ $X _ t \ \ \ (t = 1, ..., T)$ を観測できるとして、このデータから $X _ t$ の時間変化はどのようになっているのかを推論しておき、以後、推論結果を用いて未来予測を行いたいというケースを想定しています。
自己回帰モデル（ARモデル）
自己回帰モデル（ARモデル）とは、$X$ 自身で $X$ を回帰するようなモデルとなります。 ただし、その際には $X _ t$ は時刻が進む毎に値が変化するのだが、そこには何らかの規則があるはずだと想定しています。 従って、$X _ t$ を回帰モデルで表現したい時に過去のデータ $X _ {t - 1}, X _ {t - 2}, ..., X _ {t - N}$ を使って回帰ができるのではないかということが考えられるわけです。
このように過去の自分自身の値に依存して現在の値が決まるモデルをマルコフモデルと呼びます。特に現在の $X _ t$ が１つ前の $X _ {t-1}$ のみに依存する場合を１次マルコフモデルと呼びます（当然 N次マルコフモデルと言えば、過去 N 個に依存する）。今回は１次マルコフモデルを例に使っていきます。1次マルコフモデルにおいて、$X _ t$ と $X _ {t-1}$ の間にある変換（あるいは回帰関数）がどのようなものであるかに応じて、いろいろな自己回帰モデルが考えられることになります。
AR(1)
最も単純な例である１次マルコフモデルでの応用例として AR(1)モデルがあります。これは $\theta$ をパラメータに持つ関数で
$$ X _ t \simeq f _ {\theta} (X _ {t-1}) $$
と表現できます。特に一次関数
$$ X _ {t} \simeq \beta + \alpha X _ {t - 1} $$
が使われるケースが多いです（これは結局モデリングしている側が、データを時系列に並べることで未来予測のための特徴が取れると考え、過去のデータを特徴量とした単回帰を実施しているだけに過ぎない…！といえば気が楽になるだろう。ただし、気を緩めると、訓練データセットと検証データセットをランダムに選択してしまったりする。この点は、自分で未来予測をしたいと問題設定を決めたことをしっかり思い出し、過去のデータだけで学習をし、検証には未来のデータを使わなければ意味が無いことには注意が必要だ）。
さて、 $X _ t$ と $X _ {t-1}$ を回帰で結ぼうとしても、当然そこには何らかの誤差 $\epsilon$ が乗っているはずであり、この $\epsilon$ を考慮して初めて等式が成り立ち
$$ X _ {t} = \beta + \alpha X _ {t - 1} + \epsilon _ t $$
という表記になります。 $\epsilon _ t$ は毎時刻、同一の分布 $\mathcal N (0, \sigma)$ から独立に生起していると考える場合が多いです。何はともわれ、登場人物が確定したので、この"一時刻間の出来事だけに着目した"同時分布は
$$ p(X _ {t}, X _ {t-1}, \alpha, \beta, \sigma) = p(\alpha)p(\beta)p(\sigma)p(X _ {t-1})p(X _ t \mid X _ {t-1}, \alpha, \beta,\sigma) $$
と表せます。さて $X _ {t-1}, \alpha, \beta, \sigma$ に個々に事前分布があり、それらから生起した値に条件付けられて $X _ t$ が生成されるようなモデルになっています。 実際には系列が $X _ 1, ..., X _ T$ と得られるのであれば、全体は下記のようになっているはずです。ここでは紙面の都合上、 $\theta = (\alpha, \beta, \sigma)$ とまとめてしまいました。
これらを踏まえて、同時分布を書き直すと （$\theta = (\alpha, \beta, \sigma)$ として）
$$ \begin{align} p(X _ {1}, X _ {2},..., X _ {T}, \theta) &= p(\theta)p(X _ {1})p(X _ 2 \mid X _ {1}, \theta)...p(X _ T \mid X _ {T-1}, \theta)\\ & = p(\theta)p(X _ {1})\prod _ {t=2} ^ Tp(X _ t \mid X _ {t-1}, \theta) \end{align} $$
と書き表すことができます。容易すべきは事前分布 $p(X _ 1)$ と $p(\theta)$ となります。特に $p(X _ 1)$ を初期分布と言います（今回は時刻のインデックスは $t=1$ から取っていますが、文献、特にコード付きに文献では $t=0$ のインデックスが多いです）。こうしてAR(1)モデルの同時分布が書きくだせました。
さて、ここで $p(\theta) = p(\alpha)p(\beta)p(\sigma)$ とまとめて書いているだけのものでした。この事前分布と、初期分布 $p(X _ 1)$ とデータ生成モデル $p(X _ t \mid X _ {t-1}, \theta)$ に具体的な確率分布を仮定してみます。
$$ \begin{align} p(\alpha)& = \mathcal N(\alpha \mid \mu _ {\alpha}, \sigma _ {\alpha})\\ p(\beta) &= \mathcal N(\beta \mid \mu _ {\beta}, \sigma _ {\beta})\\ p(\sigma) &= {\rm LogNormal}(\sigma \mid \mu _ {\sigma}, \sigma _ {\sigma})\\ p(X _ 1) &= {\rm Delta}(X _ 1 \mid c) \\ p(X _ t \mid X _ {t-1}, \alpha, \beta, \sigma)& = \mathcal N (X _ t \mid \beta + \alpha X _ {t-1}, \sigma) \end{align} $$
としてみます。幾つかポイントを説明します。
事前分布は任意ですが、特に $\sigma$ は正の実数が生成されるような分布にする必要があるので、ここでは対数正規分布を選択しました（もちろん適当な正の実数を生成する分布から値を生成し、それを $\exp$ や $\rm softplus$ 関数などで正の値に変換する場合もある。これは再パラメータ化と呼ばれる方法で、変換は全単射である必要がある）。
また初期分布にはデルタ分布を置いており、${\rm Delta}(X \mid c)$ は$c$ にのみ$1$の値を持つ確率分布になります。なぜこんな分布を置くかというと、初期値として $X _ 1$ は得られているケースが多いからです。無論、ここに $X _ 0$ というダミーの確率変数を置いて初期分布を配置し、$X _ 1$ から観測が開始できた、という設定で推論しても構わないです。
データの生成モデル $\mathcal N (X _ t \mid \beta + \alpha X _ {t-1}, \sigma)$ は正規分布で、平均が一次関数 $\beta + \alpha X _ {t-1}$ で回帰されており、標準偏差がすべての時刻共通で $\sigma$ を用いるというものになっています。ここで、$X _ t \ \ \ (t = 1, 2, 3, ..., T)$ を観測できているとして、事後分布を推論する場合、グラフィカルモデルで観測済の変数をグレーで塗りつぶし
観測済の変数 $X _ 1, ..., X _ T$ だけの分布 $p(X _ 1, ..., X _ T)$ なるもので、これまで設計してきたすべての変数の同時分布を割ってやれば（乗法定理、ベイズの定理利用）
$$ \begin{align} p(\alpha, \beta, \sigma \mid X _ 1, ..., X _ T)&= \frac{p(X _ 1 , ..., X _ {T}, \alpha, \beta, \sigma)}{p(X _ 1, ..., X _ T)} \\ & = \frac{p(\theta)p(X _ {1})\prod _ {t=2} ^ Tp(X _ t \mid X _ {t-1}, \theta)}{p(X _ 1, ..., X _ T)} \end{align} $$
とかけて、あとは先程仮定した分布を具体的に代入していくだけになります。
寄り道：時系列データと独立同分布データ
完全に余談ですが、上記の時系列モデル、実はデータ間の関係性をより深く考察する上では常に念頭においておくべきモデルになります。 なぜならば、時系列データであることを考慮して上記のような$AR(1)$ モデルを考慮したとして、仮に時系列に一切の意味がなかった場合、 $\alpha = 0$ とすることで回帰モデルは
$$ X _ {t} = \beta + \alpha X _ {t - 1} + \epsilon _ t = \beta + \epsilon _ t $$
と表され、ここでは
$$ \begin{align} X _ t &= \beta + \epsilon _ t \\ \epsilon &\sim \mathcal N (\epsilon \mid 0, \sigma) \end{align} $$
となることによって、結局のところ
$$ p(X _ t\mid \beta, \sigma) = \mathcal N (X _ t \mid \beta, \sigma) $$
という正規分布での問題に落とし込まれるからです。$t$ のインデックスはもはや時刻としての意味を持たなくなり、単にデータをカウントするインデックスに成り代わります。こういったモデルに対しては、$X$ は互いに独立であるというような形式で推論が実施され、これをいわゆる独立同分布（i.i.d.） であると述べているのです。
とは言っても "AR(1) モデルを使えば、時系列に意味がない場合に $\alpha = 0$ といつでも推定されてくれるなんてことは無い"わけなので、推論アルゴリズムを楽に実行するために、時系列に意味がなさそうなのであれば、最初から i.i.d.でモデルを考えたほうが良いです。
ちなみに、グラフィカルモデルでは
と$X _ t$ と $X _ {t-1}$ 間の矢印をすべて撤廃してしまえば、同時分布が
$$ \begin{align} p(X _ {1}, X _ {2},..., X _ {T}, \theta) &= p(\theta)p(X _ {1})p(X _ 2 \mid \theta)...p(X _ T \mid, \theta)\\ & = p(\theta)p(X _ {1})\prod _ {t=2} ^ Tp(X _ t \mid \theta) \end{align} $$
となります。時系列を考慮して書いていた初期分布 $p(X _ {1})$ はここまで具体的にどんな分布であるかを明記していませんでしたが、この初期分布（もはや初期分布でも何でも無いが）を $p(X _ {1} \mid \theta)$ と選んでやれば
紛らわしい $T$ のインデックスの与え方を、データ数のイメージが強い $N$ に置き換えて
$$ p(X _ {1}, X _ {2},..., X _ {N}, \theta) = p(\theta) \prod _ {i=1} ^ Np(X _ i \mid \theta) $$
単に i.i.d の時の同時分布と同じものになりました。いずれにしても手に入っているデータをすべて確率変数とし、それに関わっているパラメータや潜在状態もすべて列挙する作業を通していれば、手に入れたデータが時間的に変化してそうか（データ間を矢印でつなぐ必要がありそうか）を考察する姿勢が自然と入るはずです。考えてみれば、時間のインデックスだけでなく、データの入手場所なども番地インデックスとして付与されているとしたら、番地間の関係性がありそうなのかも考慮に入れる必要が出たりと、単にデータ数 $N$ と思っていたそのインデックスにはより重要な情報が含まれているかもしれないのです。
例えば、性質が同等とされる実験対象を複数用意し、それぞれにセンサを幾つか配置し、そのセンサ値を一定時間計測した場合は、性質が同等とされる個体の数$M$ とセンサの個数 $K$ と計測時間 $T$ の $M\times K \times T$ のデータ数が得られるわけであるが、まさか、これが i.i.d. であると考えるわけにはいかないはずである。もしかしたら、「性質が同等とされる実験対象」を用意しているのであるから個体の数 $M$ の方向に関しては互いに（条件付き）独立同分布だと思っても良いかもしれない。いずれにしても、「独立なデータが○○個得られました」と言って良いのは、これまで見てきたとおり、かなり関連性を単純化できたケースに限るし、そう言えるようにデータを集めなければならないことも把握しておくべきでしょう。
ARCH(1)
さて次はARCH(1)というモデルを考えます。これもARモデルの一種でありマルコフモデルであるため、グラフィカルモデルは
となっています。ただし、具体的に当てはめていく分布が大きく変更されています。 具体的には、"時刻的に変動するのは各時刻の値の標準偏差である"ということを想定しています。
まずAR(1)モデルでは $X _ t$ と $X _ {t-1}$ の関係性を記している生成モデル $p(X _ t \mid X _ {t-1}, \theta)$ に関して
$$ p(X _ t \mid X _ {t-1}, \theta) = \mathcal N (X_ t \mid \beta + \alpha X _ {t-1}, \sigma) $$
という平均が過去の値で回帰されている例を書いていました。これは、割と思い浮かびやすい発想ですね。過去の値に応じて次の値の大体の値が決まってくるのだが、一定のバラつきを持って確率的に値が決まるという形です。しかしARCH(1)モデルでは、平均はずっと同じ値であるのだが、標準偏差が過去の値に依存して変動するような分布を書くことにして
$$ p(X _ t \mid X _ {t-1}, \theta) = \mathcal N \left(X_ t \mid \mu, \sqrt {\beta + \alpha (X _ {t-1} - \mu) ^ 2} \right) $$
と記されます（ここで $\theta = (\alpha, \beta, \mu)$）。あとは $p(\theta) = p(\alpha)p(\beta)p(\mu) $ などとしておいて、各パラメータに事前分布を設定してあげれば同時分布を設計してやることができます（グラフィカルモデル自体は同じなので同時分布の分解も同様に書ける。単に各パーツに代入する具体的分布が違うということである）。
ところで、こんな奇妙なモデルがどのように出てくるのかは自明ではありませんね。とは言いつつ、分散 $\sigma ^ 2$ という立場で考えれば、 $(X - \mu) ^ 2$ というのは分散の次元を持っているわけですから、１個前の時刻の平均からのずれの二乗（分散の次元）を一次関数で回帰していると思えば、それほど飛躍しているわけでもなさそうです。すなわち各時刻の分散 $V _ t$ に関して
$$ V _ t = \beta + \alpha V _ {t - 1} $$
みたいな一次関数を考えたいなと思った時に、１個前の分散の値を $ V _ {t - 1} = (X _ {t-1} - \mu) ^2$ と見積もる（あるいは見積もれるような $\beta, \alpha$を決める）と設定しただけの話です。
より詳しく進むと、
$$ X _ t = \mu + \eta _ t $$
という $\mu$ に対して $\eta _ t$ というバラつきが付与されるのだと考え、この $\eta _ t$ は
$$ \eta _ t = \sigma _ t \epsilon _ t $$
と分解できると考えます。 $\epsilon _ t$ は毎時刻 $\mathcal N (0, 1)$ の標準正規分布から独立に生成されるとしておきます。この分解を考えると、我々の興味は$X _ t$ にとってバラつきである $\eta _ t$ の標準偏差だけに興味があるのだと明確化できます（平均は $0$ であると決めている）。その標準偏差 $\sigma _ t$ のことをボラティリティなどと金融の分野では呼びます（なぜ固有名詞を与えているのだろう。それはこの時系列での標準偏差が非常に重要だからである。もう少し後で解説する）。
このボラティリティの二乗、つまり分散の次元の $\sigma _ t ^ 2$ を
$$ \sigma _ t ^ 2 = \beta + \alpha \sigma _ {t-1} ^ 2 $$
と一次関数で回帰してやろうと考えたわけです。ボラティリティが金融で重要なのは、仮に株価である $X _ t$ が仮に適正価格にいて平均的には $\mu$ にいるとしても、毎日のトレードで値は変動しているわけです。この時、変動の具合が大きい（＝ボラティリティが大きい）ときはかなり頻繁に（あるいは高額の？）取引が行われている状態であり、株価の平均価格自体は長い目で見れば何も変わっていないのに、運が良ければとても儲かる、そして運が悪ければかなり損をする状態にあるということになります。すなわちリスクやリターンを見積もるのに非常に重要な指標になるというわけです。例えば、アルゴリズムトレードでとりあえず購入し、次の売却を価格があたった瞬間に即座に売却できるシステムをあなたが持っているのだとすれば、無論、ボラティリティが高い状況を利用したいわけですね（そんなことは当然夢のような話なわけで、実際には既に述べたとおり、リスクやリターンという概念を数理的に扱うために利用する…はず。っていうか知ったかぶりの想像で書きました。金融のことなんてなんも知りません…）。
そうして、間にある不要な変数たちをひたすら代入消去していき、$X _ t$ の分布がどうなるのかを調べると先ほどのモデル
$$ p(X _ t \mid X _ {t-1}, \theta) = \mathcal N \left(X_ t \mid \mu, \sqrt {\beta + \alpha (X _ {t-1} - \mu) ^ 2} \right) $$
に至るというわけです。
最後に
とりあえず、確率モデリングにおける、時系列の基本的なことは抑えたように思います。（これを確率分布として考えずに、関数推定を点推定で実施していくことを考えると、信号処理や制御ではお馴染みの伝達関数的な概念が出てくる。AR(1) なんかは丸っきり１階差分方程式の同定なのだから自ずとZ変換周りの話と繋がるのである。無論差分方程式の安定性議論などもできるし、これから同定しようという差分方程式がもしも推定解付近で不安定だとしたらどうしようとかが心配になってくる。意外とこの辺を真面目に取り組むのは難しい。今回は確率モデリングのパートの一部なので踏み込みはしないが、機械学習や信号処理に慣れ親しんでいれば馴染みはあるので突っ込んでみても良いかもしれない。また、自己回帰成分だけでなく何らかの入力の移動平均成分があるなどの場合も考えることができる。）
今回ボラティリティの回帰モデルがあり、不要な変数を消去していくと $X _ t$ の生成モデルの標準偏差を回帰しているケースを考えることができました。しかし、実は観測している $X _ t$ の背後に別の確率変数が時系列的に変化しており、その結果として $X _ t$ が時刻的に変化していくことを直接表現できるモデルもあります。それがいわゆる状態空間モデルになります。
