論文読みメモ： Empirical Likelihood Estimation of Levy Processes（その1）

以下のワーキングペーパーを読みます。

Naoto Kunitomo, Takashi Owada. Empirical Likelihood Estimation of Levy Processes. CIRJE-F-272, Graduate School of Economics, University of Tokyo, 2004.
https://www.carf.e.u-tokyo.ac.jp/research/1313/

※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。おかしい点がありましたらご指摘ください。

次回：まだ

f:id:cookie-box:20190101155733p:plain:w60

レヴィ過程と無限分解可能分布（安定分布を特別なケースとして含む）に対する、経験尤度法を用いた新しいパラメータ推定方法を提案したということですが…レヴィ過程と無限分解可能分布とは何でしょう？最大経験尤度推定量は、経験特性関数に対する制約の数が多いとき一致性、漸近正規性、漸近有効性を満たすということですが…制約の数が問題なんですか…？

f:id:cookie-box:20190101160814p:plain:w60

レヴィ過程と無限分解可能分布の解説まではこのペーパーにはかいてないかな？イントロにあるように Feller (1971), Zolotarev (1986), Sato (1999) を参照ってことみたいだね。たださすがに書籍だからなあ…3つ目は最初だけちょっとみたいな PDF は見つかるね。2つ目も一部抜粋っぽい PDF があるな。1つ目はちょいちょい PDF が落ちてるけど素性がわかんないなあ…どのみちかなり古いし…。

http://assets.cambridge.org/97811076/56499/frontmatter/9781107656499_frontmatter.pdf

とりあえずウィキペディアにかいてあることをざっとまとめておくね。レヴィ過程っていうのは、

Lévy process - Wikipedia

添え字 0 でとる値はほとんど確実に 0。
真に増大する添え字列をとると、隣り合う添え字の値の差は互いに独立（独立増分性）。
異なる2つの添え字をとって、大きい添字の値から小さい添え字の値を引いた値は、添え字の差のみに依存する（定常増分性）。
任意の添え字 $t$ 、任意の正数 $\epsilon$ について $\lim_{h \to 0} P(| X_{t+h} - X_t| > \epsilon) = 0$ （確率連続性）。

これらに加えて「添え字から値への写像が、ほとんど確実に右連続かつ左極限をもつ」というのを定義に追加する場合もあるっぽい？ Sato (1999) の Preface をみると、レヴィ過程の例としてブラウン運動、ポアソン過程、安定過程などが挙げられているね。さらに、レヴィ過程の定常増分性を落とした確率過程のことを加法過程というけど、加法過程は任意の添え字に止めた値が無限分解可能分布にしたがう。無限分解可能分布というのは、それににしたがう確率変数が「任意の個数の i.i.d. な確率変数の和で表せる」がその定義みたいだね。

Infinite divisibility (probability) - Wikipedia

レヴィ過程については、なぜ任意の添え字で止めた値が無限分解可能分布にしたがうかは上のウィキペディアにかいてあるね。簡単な話で、ある添え字 $t$ までを $n$ 等分すれば独立な $n$ 個の i.i.d. な確率変数の和に分解できるからね。この無限分解可能分布の例としては、ポアソン分布、負の二項分布、ガンマ分布、退化分布もそうだし、正規分布、コーシー分布といった任意の安定分布族もそうらしい。

Stable distribution - Wikipedia

安定分布っていうのは、その分布から生成した2つの独立な確率変数 $X_1, X_2$ を任意の重み $a, b > 0$ で足した和 $aX_1 + bX_2$ がしたがう分布が、ある $c>0, d$ を選べば $cX + d$ がしたがう分布と同じになる（ $X$ はその分布から生成した1つの確率変数）。その確率密度関数は一般的に解析的にかけないけど、特性関数はパラメトリックにかけるみたい。ファットテールな現象のモデリングに用いられるみたいだね。

えっと…今回分布のパラメータ推定をするということなのでパラメトリックな表式で「レヴィ過程はこういう分布」「無限分解可能分布はこういう分布」というのを知りたかったのですが、その説明はだいぶ期待とちがいますね…話を整理するとこうでしょうか。

レヴィ過程はほとんど確実にゼロからはじまり、時間を区切ったときの各区間の増分が独立で時間差のみに依存する、確率連続性をもつ確率過程である。
レヴィ過程（もっというと加法過程）は添え字を止めると無限分解可能分布にしたがう。
無限分解可能分布の例に安定分布がある。
安定分布は一般に確率密度関数はかけない。特性関数のみがパラメトリックにかける。

ということですか…最後これ、特性関数しか特定できないって、そんな状態でそもそもどうやって確率密度関数を推定するんです？「新しい手法を提案した」以前に、既存の手法がわかりません…。

そこはイントロの2段落目かな…経験分布のパーセンタイルや経験特性関数に基づいた手法が提案されているみたい。これらの詳細はこのイントロだけ読んでもわかんないな。

結局何もわかりませんね…続きに進むと、このペーパーでは経験尤度法に基づいたアプローチを提案するということですが、経験尤度法とは？

Owen (1988, 1990) によって提案され Qin and Lawless (1994) が拡張した手法らしい。

最後のアブストを読むと、未知の分布 $F$ の何らかのパラメータの経験尤度比統計量 $\theta(F)$ のしたがう分布はカイ2乗分布に近づくので、パラメトリックな尤度を用いる場合と同様に信頼区間を推定できるとか。このペーパーでは経験尤度法が安定分布のパラメータ推定に応用できることを示すらしい。そしてさらにレヴィ過程と無限分解可能分布のパラメータ推定などに拡張するみたいだけど、他の先行研究とは別のベクトルに Qin and Lawless (1994) の研究を拡張したものになっているっていっているね。

とにかく2節に進みますか。安定分布の経験尤度推定ということですが… $X_i \, (i=1, \cdots, n)$ は独立に同一の安定分布にしたがう確率変数ということですね。それで、 $X_i$ の特性関数を $\phi_\theta(t) \equiv \mathbb{E}[e^{itX_i}]$ とすると以下が成り立つということです。成り立つというか、どちらかというとこれが安定分布の定義ですかね。

$\phi_\theta(t) = \phi_\theta^R(t) + i \phi_\theta^I(t)$
$\phi_\theta^R(t) = e^{-|\gamma t|^\alpha} \cos \Bigl[ \delta t + \beta \gamma t \bigl( |\gamma t|^{\alpha - 1} - 1 \bigr) \tan \frac{\pi \alpha}{2}\Bigr]$
$\phi_\theta^I(t) = e^{-|\gamma t|^\alpha} \sin \Bigl[ \delta t + \beta \gamma t \bigl( |\gamma t|^{\alpha - 1} - 1 \bigr) \tan \frac{\pi \alpha}{2}\Bigr]$ パラメータ空間は $\Theta = \{0<\alpha \leqq 2, \; -1<\beta \leqq 1, \; 0 < \gamma, \; \delta \in \mathbb{R}\}$ ということですが…Wikipedia だと $\alpha=1$ は場合分けされているようにみえますが、整合的なんでしょうか…まあとりあえずパラメータ $\theta = (\alpha, \beta, \gamma, \delta)^\top$ に対応する安定分布を $F_\theta(\cdot)$ とかくということです。それで $\theta$ を推定したいわけですが、確率密度関数の形が明示的にわからないので通常の最尤推定ができません。また、安定分布は 1, 2 次のモーメントをもつと限らないので漸近理論でスタンダードな手法も適用できないと…？

モーメントは一致性をもつからモーメントを通じてパラメータを探索できるってことなのかな…？コーシー分布は平均、分散、歪度、尖度が未定義だね。

それで経験尤度法を導入するわけですが、まず、 $L_n(F_\theta) = \prod_{k=1}^n p_k$ なる経験尤度関数を定義します。 $p_k$ の制約は $0 \leqq p_k$ と $\sum_{k=1}^n p_k = 1$ のみということです。これで単に $L_n(F_\theta)$ を最大化するならば任意の $k$ に対して $p_k =1/n$ ということでそれはそうですが…解釈しづらいですね。例えばコインを100回投げたとして、それぞれに 0.01 という確率を割り振って何が楽しいのか…それぞれの結果は同等に起こりやすかったのだということになりますが…。

読み進めていこうか。経験尤度比関数を $R_n(F_\theta) = \prod_{k=1}^n np_k$ と定義するらしい。これは最適な $1/n$ からの比にしているようにみえるね。それで最大経験尤度推定量 $\hat{\theta}_n$ は、 $R_n(F_\theta)$ を以下の制約下で最大化したものだって。この $m$ が経験特性関数の制約の数らしい。 $\sum_{k=1}^n p_k \cos(t_l X_k)$ と $\sum_{k=1}^n p_k \sin(t_l X_k)$ は $m$ 点の $t_1 < t_2 < \cdots < t_m$ で評価した経験特性関数の実部と虚部だって。
$\mathcal{P}_n = \Bigl\{ p_k \geqq 0 \, (k = 1, \cdots, n), \; \sum_{k=1}^n p_k = 1, \Bigr.$
$\qquad \quad \, \sum_{k=1}^n p_k \bigl( \cos(t_l X_k) - \phi_\theta^R(t_l) \bigr) = 0 \, (l = 1, \cdots, m),$
$\qquad \quad \, \Bigl. \sum_{k=1}^n p_k \bigl( \sin(t_l X_k) - \phi_\theta^I(t_l) \bigr) = 0 \, (l = 1, \cdots, m)\Bigr\}$

経験特性関数？特性関数は $\phi_\theta(t) \equiv \mathbb{E}[e^{itX_i}]$ ですよね？とするとその経験特性関数はこの期待値 $\mathbb{E}$ を分布 $\{p_k\}$ 上での期待値で代えたものであって、ある点 $t_l$ で評価したものですね？そしてどの点 $t_l$ でも真の特性関数とのずれの $\{p_k\}$ 上での期待値がゼロになるようにする…？ $X_k$ が都合の悪い点だったら $p_k = 0$ とすることもできてしまいそうですね…？

それで以下のように $g(X_k, \theta)$ を定義すると $\mathbb{E}_{\theta_0} \bigl[ g(X, \theta_0) \bigr] = 0$ が満たされなければならないって。 $\theta_0$ は真のパラメータね。部長がこれまでかいていた $\mathbb{E}$ もこの表記なら $\mathbb{E}_{\theta_0}$ だね（まあ実際には分布 $\{p_k\}$ 上での期待値で代えるはずだけど）。

$\displaystyle g(X_k, \theta) = \begin{pmatrix} \cos(t_1 X_k) - \phi_\theta^R(t_1) \\ \vdots \\ \cos(t_m X_k) - \phi_\theta^R(t_m) \\ \sin(t_1 X_k) - \phi_\theta^I(t_1) \\ \vdots \\ \sin(t_m X_k) - \phi_\theta^I(t_m) \end{pmatrix}$ $\displaystyle g^R(X_k, \theta) = \begin{pmatrix} \cos(t_1 X_k) - \phi_\theta^R(t_1) \\ \vdots \\ \cos(t_m X_k) - \phi_\theta^R(t_m) \end{pmatrix}, \quad g^I(X_k, \theta) = \begin{pmatrix} \sin(t_1 X_k) - \phi_\theta^I(t_1) \\ \vdots \\ \sin(t_m X_k) - \phi_\theta^I(t_m) \end{pmatrix}$ そして、ラグランジュ関数が以下のように定義できる。 $\mu$ と $\lambda$ はラグランジュ乗数だね。 $\displaystyle \mathcal{L}_n(\theta) = \sum_{k=1}^n \log (n p_k) - \mu \left( \sum_{k=1}^n p_k - 1 \right) - n \lambda^\top \left( \sum_{k=1}^n p_k g(X_k, \theta) \right)$ これを $p_k$ で微分すると $p_k^{-1} = - \mu + n \lambda^\top g(X_k, \theta)$ を得るから、この両辺に $p_k$ をかけて全ての $k$ について足し上げれば $\hat{\mu} = n$ を得る。したがって $\hat{p}_k = (1/n)[1 + \lambda^\top g(X_k, \theta)]^{-1}$ とかける。これを用いた $0 = \sum_{k=1}^n \hat{p}_k g(X_k, \theta)$ の解が $\lambda(\theta)$ だ。

え、えっと？いま $p_k, \, \theta, \, \mu, \, \lambda$ が最適化の対象で、 $\hat{\mu}$ はすぐわかって、 $\hat{p}_k$ は $\lambda$ と $\theta$ の関数でかけて、さらに $\lambda$ も $\theta$ によって決まるんですね？

うん、ある $\theta$ に対してペーパーの (2.6) 式によって $\lambda$ が決まるね。さらにこの (2.6) 式自体を一番大きくする $\theta$ が最適な $\hat{\theta}$ だね。この最適化は通常2段階でやるらしいね。以下の書籍をみてって。

Amazon | Empirical Likelihood (Chapman & Hall/CRC Monographs on Statistics & Applied Probability) | Art B. Owen | Applied

最後ぶん投げましたね…。最尤法と経験尤度法をざっくり対比すると、

最尤法では尤度を最大化するようなパラメータを求める。あるパラメータでの各データの尤度はパラメトリックなモデルにそのパラメータを当てはめればわかる。
経験尤度法では「経験尤度」を最大化するようなパラメータを求める。あるパラメータでのデータ $X_k$ の尤度 $p_k$ は、経験特性関数と真の特性関数との点 $t_1, \, t_2, \, \cdots, \, t_m$ でのずれの $\{p_k\}$ 上での期待値がゼロになるように割り当てる。

…点 $t_1, \, t_2, \, \cdots, \, t_m$ をどこに取るべきかが示されていませんね？

後で議論されるのかな？

（その2があれば）つづく