統計数理（2017）
第 65 巻 第 1 号 21–38
c 2017 統計数理研究所


特集「高頻度金融データに基づく統計的推測とモデリング」［研究ノート］Lévy 駆動型確率微分方程式の段階的推定について
上原 悠槙† ・増田 弘毅†（受付 2016 年 6 月 28 日；改訂 10 月 13 日；採択 10 月 17 日）

要

旨

非正規型 Lévy 過程で駆動される確率微分方程式モデルの推定を考える．指数的エルゴード
性とデータの高頻度性の下，正規型疑似スコア関数に基づいてスケール係数およびドリフト係
数をこの順に段階的に推定する方法を提案し，推定量の漸近正規性および裾確率評価を導出す
る．推定対象を分割することで最適化の計算負荷が削減され，より安定した推定精度が得られる．拡散過程の場合と異なり，特に両係数のパラメータが共通因子を持つ場合には，ドリフト係数の漸近共分散行列は同時推定の場合と異なる形をとる．
キーワード：エルゴード性，正規型疑似スコア関数，高頻度観測，段階的推定，Lévy
駆動型確率微分方程式．
1.

はじめに

連続時間確率過程で記述される統計モデルによって観測頻度に伴う推定精度の差異がモデリング可能となり，モデルの具体系のみならず推定方式によって漸近分布論も様々な形式を取り得る．本稿では，常微分方程式 dxt = a(xt , α)dt の時間発展に確率的なノイズを加えた，フィル
ター付き確率空間（Ω, F , (Ft )t∈R+ , P ）上に定義された次の d 次元確率微分方程式（SDE）モデル(1.1)

dXt = a(Xt , α, γ)dt + c(Xt− , γ)dJt .

を考える．ここで J は (Ft ) - 適合な r 次元非正規型 Lévy 過程であり，ドリフト係数 a: Rd × Θα ×Θγ → Rd とスケール係数 c: Rd × Θγ → Rd ⊗ Rr は有限次元未知パラメータθ := (α, γ) ∈ Θ := Θα × Θγ ⊂ Rpα × Rpγ

を除いて既知であるとし，初期変数 X0 は J と独立とする．簡単のためパラメータ空間 Θ は有（1.1）
の解過程から高頻度観測（Xt0 , Xt1 , . . . , Xtn ）が得られている
界凸領域とする．SDE モデルnと仮定する．ここで，tj = tj = jhn であり，正数列 (hn )n∈N は(1.2)

Tn := nhn → ∞,nh2n → 0

を満たしているとする．目的は θ の真値 θ0 ∈ Θ を推定することである．
Lévy 過程とは独立定常増分性を有する確率過程であり，離散時間ランダムウォークの自然なLévy
連続時間版と解釈される．これは互いに独立な Wiener 過程と非正規型（純粋ジャンプ型）
過程の畳み込みとして，また同時に複合 Poisson 過程のある種の極限として定義される．そのジャンプ構造・挙動の多様性により，Lévy 過程は非正規性が色濃い金融データや工学・生体信†

九州大学大学院 数理学研究院：〒 819–0395 福岡県福岡市西区元岡 744

統計数理

22

第 65 巻

第 1 号 2017号の表現として広く用いられている．Lévy 過程の理論を包括的に扱った標準的な教科書として
Sato（1999）
を挙げておく．
現在膨大な高頻度データの入手が可能となり，上述の現象記述精度の高さから Lévy 駆動型
SDE モデルの実現象への応用の機運が高まっている．しかしながらその高頻度統計理論は J に対応する無限分解可能分布族を（相当）限定した場合でさえ未だ十分に整備されていない．これは一つには，ジャンプ構造が過度に多様なため統一的な推定手法の構築が不可能であるという事実が災いしている．実際，最尤推定量の最適収束率も観測頻度とジャンプ構造の双方が本質的に関与して決まる．この問題への一つの対処法として，Masuda（2013）は（1.1）でドリフト係
数 a(x, α, γ) が γ によらない場合に，漸近効率を犠牲にした正規型疑似最尤推定量の漸近挙動を調べた．この手法により，漸近有効性を放棄することで陽な疑似尤度関数を用いつつ駆動ノイ
ズ J の具体的な分布構造を仮定せずに θ の推定を行える．
ところで単一の疑似尤度に基づいた θ の推定は α と γ の同時最適化を伴い，Masuda（2013）の推定方式も例外ではない．このため，係数が複雑な場合には計算負荷の高さおよび同時最適の場合にこの問題に対す化の不安定さが懸念されよう．近年，拡散過程（J が標準 Wiener 過程）
る種々のアプローチ（Uchida and Yoshida, 2012 や Kamatani and Uchida, 2015 とその参考文献を参照）が提案され，拡散係数とドリフト係数を個別に，しかも漸近有効性を保持しつつ推定する形で，数値統計の観点からも発展がなされた．これらの方法論においては α と γ の最適収束率の差異が有効活用されている．翻って Lévy 駆動型 SDE の場合についてはその限りではない．
√
実際，正規型疑似尤度は α, γ 共に Tn の収束率しか持たず，さらに一般に漸近共分散行列は拡
散過程の場合のようにブロック対角型にはならないため，段階的な推定手法の定式化は自明な話ではない．この点を踏まえ，本稿の主眼は θ の推定を α, γ 個々の推定へ分割する段階的推定法を定式化し，当該推定量の漸近分布と多項式型裾確率評価を導出することにある．提案手法
によって γ, α をこの順に別個に推定可能となり，計算負荷の大きな軽減の可能性が期待できる．
本稿の構成は以下の通りである．まず 2 節で段階的推定量 θ̂n 構成の概略を述べ，その漸近的性質を導出するための仮定を導入する．3 節において主結果である θ̂n の漸近正規性および裾確率評価を与える．特に，漸近分布の形を通じて拡散過程の場合との差異を浮き彫りにする．4 節
では段階的推定量のパフォーマンスを確認するために幾つかの数値実験結果を紹介し，最後に5 節において主結果の証明を与える．
2.

段階的係数推定法の段階的係数推定法について述べていく．
本稿の主眼である（1.1）
J が有限共分散をもつとき，標準化した J˜t := cov(J1 )−1/2 {Jt − E(J1 )t} を用いて
（1.1）を書き
換えると
(2.1)

dXt = (a(Xt , α, γ) + c(Xt , γ)E(J1 ))dt + c(Xt− , γ)cov(J1 )1/2 dJ˜tとなる．新たに E(J1 ), cov(J1 ) をパラメータとみなせば（2.1）は（1.1）での特別な場合となるため，（1.1）において一般性を失うことなく J は標準化されているとしてよい．ここではより強く，以．
下を仮定する（Ir は r 次元単位行列を表す）
仮定 2.1. E(J1 ) = 0, E(J1⊗2 ) = Ir , また任意の q > 0 に対し E[|J1 |q ] < ∞.
任意の確率過程 Y と Rd × Θ̄ 上の可測関数 f に対して
Yj = Yt j ,Δj Y = Yj − Yj−1 ,fj (θ) = f (Yj , θ)Lévy 駆動型確率微分方程式の段階的推定について

23と適宜略記する．
（1.1）の形式的な Euler - 丸山近似は
Xj ≈ Xtj−1 + hn aj−1 (α, γ) + cj−1 (γ)Δj J(2.2)

で与えられる．以下では
S := c⊗2と表す．ただし任意の行列に対し x⊗2 := xxT（T は転置を表す）である．形式的に局所正規近似
Δj J ≈ N (0, hn Ir ) を適用した場合の条件付き尤度の近似
L(Xj |Xtj−1 = x) ≈ N (x + hn a(x, α, γ), hn S(x, γ))(2.3)に基づき，Masuda（2013）
では，ドリフト係数とスケール係数内のパラメータの重複がない場
において，同時正規型疑似尤度合（ドリフト係数が a(x, α) の形）


n
1
1 −1
(2.4)
Mn (θ) := −Sj−1 (γ)[(Δj X − haj−1 (θ))⊗2 ]log |Sj−1 (γ)| +2 j=1hnから得られる疑似最尤推定量の漸近正規性，モーメント収束を導出した．局所正規近似（2.3）そのものは理論上正しくないが，条件付き平均・共分散構造に関するモーメント法の一種として推定に利用可能であることが分かる．このような正規型疑似尤度の有用性は時系列モデルにおいてもよく知られているが，Lévy 駆動型 SDE モデルの場合には収束率が真に下がる．
以下では ∂x で変数 x に関する偏微分を，また |A| で正方行列 A の行列式を表す．また，多重線形形式 M := {M (i1 ···iK ) : ik = 1, . . . , dk ; k = 1, . . . , K} ∈ Rd1 ⊗ · · · RdK と dk 次元ベクトル(j)uk = {uk } について
M [u1 , . . . , uK ] :=

d1

i1 =1

···

dK


(i )

(i

M (i1 ,...,iK ) u1 1 · · · uKK

)

iK =1

と表す．
正規型疑似スコア関数 θ → (G1,n (γ), G2,n (α, γ))，ただし

n 
1
1
−1
G1,n (γ) :=∂γ (Sj−1)(γ)[(Δj X)⊗2 ]∂γ log |Sj−1 (γ)| +
n j=1
hn


n
1
1
pγ
−1
−1
−1
=
∂γi Sj−1 )(γ)]i=1−(Sj−1(∂γ Sj−1 )Sj−1)(γ)[(Δj X)⊗2 ] ,[trace(Sj−1n j=1hn
1  −1
S (γ)[∂α aj−1 (θ), Δj X − hn aj−1 (θ)],
Tn j=1 j−1n

G2,n (α, γ) :=に対して以下で定義される段階的推定量 θ̂n := (α̂n , γ̂n ) を提案する．
（1）γ̂n を G1,n (γ) = 0 のランダムな解とする．
（2）上記 γ̂n を用いて，α̂n を G2,n (α, γ̂n ) = 0 のランダムな解とする．
このように定義された θ̂n は，元々の疑似尤度 Mn による (pα + pγ ) - 次元同時最適化を γ と α 個別の最適化へ分断可能とするため，推定値計算負荷の削減による時間短縮へつながる．
G1,n (γ) は，Mn の定義でドリフト項を無いものと勝手にみなした場合の γ に関する疑似スコ
ア関数をスケーリングしたものに相当する．このドリフト非依存正規型疑似尤度は拡散過程の拡散項の推定において今日標準的となったが，そこに至るまでは当該分野において数多の試行
錯誤があったと思われる．近年においては，計算負荷的観点から疑似尤度関数の同時最適化変

統計数理

24

第 65 巻

第 1 号 2017

数の削減や初期推定量の構築などの目的にも利用されている．関連する先行研究については，
Uchida and Yoshida（2012）
や Kamatani and Uchida（2015）とその参考文献を参照されたい．
ドリフト項，拡散項に共通パラメータを持つ拡散過程に上記の段階的推定法を適用すると，収√√束レートの違い（ドリフトパラメータは Tn ，拡散パラメータは n）から，その漸近分散は共通
パラメータを持たない場合と同様のものが得られる．しかしながら，Lévy 駆動型 SDE において正規型疑似尤度最尤推定量はドリフトパラメータ，スケールパラメータともに収束レートが√Tn であるため（Masuda, 2013），状況は本質的に異なる．実際，後述の定理 3.4 によって，疑似スコア関数 (G1,n (γ), G2,n (α, γ)) に基づく段階的推定量 θ̂n = (α̂n , γ̂n ) は，ドリフト・スケー
ル係数が共通パラメータを持つ場合には α̂n の漸近分散が変わることが明示的に示される．
3.

主結果

本節の定理 3.4 において，推定量の漸近正規性に加え，推定量依存型統計量の確率・モーメ
ント評価を保証する裾確率評価を与える．
仮定 3.1. ドリフト係数 a(·, θ0 ) とスケール係数 c(·, γ0 ) は Lipschitz 連続で，各 i ∈ {0, 1, 2},j ∈ {0, . . . , 5}, k ∈ {0, . . . , 5} に対して連続な導関数 ∂xi ∂αj ∂γk a および ∂xi ∂γk c が存在し，ある非負定数 C(i,j,k) に対して
(3.1)

sup
(x,α,γ)∈Rd ×Θ̄α ×Θ̄γ

1{|∂xi ∂αj ∂γk a(x, α, γ)| + |∂xi ∂γk c(x, γ)| + |S −1 (x, γ)|} < ∞1 + |x|C(i,j,k)が成り立つ．
仮定 3.2. ある確率測度 π0 と正定数 a, C に対して





(3.2)
sup  f (y){P (Xt ∈ dy|X0 = x) − π0 (dy)} ≤ Ce−at gq (x),f :|f |≤gq

x ∈ R, q > 0,が成立する（gq (x) := 1 + |x|q ）．さらに任意の q > 0 に対して
(3.3)

sup E[|Xt |q ] < ∞.
t∈R+

仮定 3.2 の下，エルゴード定理より任意の多項式増大関数 f についてn1Pf (Xtj−1 ) −→ f (x)π0 (dx)n j=1が成立する．さらに大数の法則の収束率を加味したモーメント評価を得ることができる（Masuda,．
2013, Lemma 4.3）∞
p
関数 G∞ (θ) := (G∞α (α), Gγ (γ)) ∈ R および定数行列 Iα ∈ Rpα ⊗ Rpα , Iγ ∈ Rpγ ⊗ Rpγ を以

pα

pγ
下で定義する: u, u ∈ R と v, v ∈ R に対して

G∞
(α)
:=
∂α a(x, θ0 )T S −1 (x, γ0 )[a(x, α0 , γ0 ) − a(x, α, γ0 )]π0 (dx),
α

(γ)
:=
∂γ S −1 (x, γ)[S(x, γ0 ) − S(x, γ)]π0 (dx),G∞γIα [u, u ] := − (S −1 (x, γ0 )[∂α a(x, θ0 )[u], ∂α a(x, θ0 )[u ]])π0 (dx),Iγ [v, v  ] := trace{(S −1 ∂γ S ⊗ S −1 ∂γ S)(x, γ0 )[v, v  ]}π0 (dx).

Lévy 駆動型確率微分方程式の段階的推定について

モデルの識別可能性・非負定値性のため，以下の仮定を導入する．
仮定 3.3. Iα , Iγ はそれぞれ逆行列を持ち，さらにある正定数 χα , χγ が存在して|G∞
α (α)| ≥ χα |α0 − α|,

(3.4)

|G∞
γ (γ)| ≥ χγ |γ − γ0 |,が任意の (α, γ) に対して成り立つ．
ν0 で J の Lévy 測度を表すとき，i1 , . . . , im ∈ {1, . . . , r} に対して

νi1 ,...,im (m) = z (i1 ) · · · z (im ) ν0 (dz)と書く（仮定 2.1 の下，これらは有限値として定まる）
．記号




Σα
Σ̂α,n Σ̂α,γ,n
Σα,γΣ=, Σ̂n =
,
Sym.
Σγ
Sym.
Σ̂γ,n




−Σ̂α,n B̂n
Iα B, În :=
I :=
0
Îγ,n
0 Iγ
を導入する．ここで u, u ∈ Rpα と v, v  ∈ Rpγ に対して，Σα [u, u ] := S −1 (x, γ0 )[∂α a(x, θ0 )[u], ∂α a(x, θ0 )[u ]]π0 (dx),Σγ [v, v  ] :=r

νs,t,s ,t (4){∂γ S −1 (x, γ0 )[v, c(·s) (x, γ0 ), c(·t) (x, γ0 )]}s,t,s ,t =1· {∂γ S −1 (x, γ0 )[v  , c(·s ) (x, γ0 ), c(·t ) (x, γ0 )]}π0 (dx),rΣα,γ [u, v] :=νk1 ,k2 ,k3 (3)S −1 (x, γ0 )[∂α a(x, θ0 )[u], c(·k1 ) (x, γ0 )]k1 ,k2 ,k3 =1· ∂γ S −1 (x, γ0 )[v, c(·k2 ) (x, γ0 ), c(·k3 ) (x, γ0 )]π0 (dx),B[u, v] := − S −1 (x, γ0 )[∂α a(x, θ0 )[u], ∂γ a(x, θ0 )[v]]π0 (dx),Σ̂α,n [u, u ] :=1  −1S (γ̂n )[∂α aj−1 (θ̂n )[u], ∂α aj−1 (θ̂n )[u ]],n j=1 j−1Σ̂γ,n [v, v  ] :=1 (∂γ S −1 ⊗ ∂γ S −1 )(γ̂n )[(v, (Δj X)⊗2 ), (v  , (Δj X)⊗2 )],Tn j=1Σ̂α,γ,n [u, v] :=1  −1−1(S⊗ ∂γ Sj−1)(γ̂n )
Tn j=1 j−1

n

n

n

· [(∂α aj−1 (θ̂n )[u], Δj X − hn aj−1 (θ̂n )), (v, (Δj X)⊗2 )],1  −1S (γ̂n )[∂α aj−1 (θ̂n )[u], ∂γ aj−1 (θ̂n )[v]],n j=1 j−1nB̂n [u, v] := −Îγ,n [v, v  ] :=

1
−1
−1
trace{(Sj−1∂γ Sj−1 ⊗ Sj−1
∂γ Sj−1 )(γ̂n )[v, v  ]}n j=1nである．最後に φ(·; 0, V ) で N (0, V ) の密度関数を表す．

25

統計数理

26

第 65 巻

第 1 号 2017

定理 3.4. 仮定 2.1, 3.1, 3.2, 3.3, さらに行列 Σ の正定値性の下で以下が成り立つ．
（1） 裾確率評価: 任意の K > 0 に対してある定数 CK > 0 が存在して
(3.5)

√
√
CK
sup P (| Tn (α̂n − α0 )| > r) + sup P (| Tn (γ̂n − γ0 )| > r) ≤ K ,rnnr > 0.

（2） 漸近標準正規性:
√

(3.6)

L

În (θ̂n − θ0 ) −→ N (0, Ip ).
Tn Σ̂−1/2n（3） モーメント収束: 高々多項式増大度を持つ任意の連続関数 f に対して

√
(3.7)
E[f ( Tn (θ̂n − θ0 ))] → f (y)φ(y; 0, I −1 Σ(I −1 ) )dy.
注意 3.5. 定理 3.4 より特に，正規型疑似尤度を適用する場合の Δj X の微小時間変動においてはスケール係数部分が支配的であることがわかる．すなわち，元々の α, γ の同時推定量の収√（独立である）とは限らない状束率は共に同じ Tn でしかもその同時推定量は漸近的に直交する況であるにも関わらず，ドリフト係数を完全に無視しても，γ の同時推定と段階的推定は漸近同等となる．
注意 3.6. 行列 B は共通パラメータの存在により生じる有限バイアスの補正的役割を果たしている．実際，係数間のパラメータの重複がない場合には B = 0 であり，α̂n と γ̂n の漸近分散のそれと等しく，推定分離による推定量の分散の変化は生じない．さらにこは Masuda（2013）（独立になる）．
の時，ν(3) = 0 であれば α̂n と γ̂n は漸近的に直交する
注意 3.7. d = r (dim(Xt ) = dim(Jt )) の場合には，駆動ノイズの増分の代替として Euler 残差と同様，駆δ̂j := cj−1 (γ̂n )−1 (Δj X − hn aj−1 (θ̂n )) を用いることで，Masuda and Uehara（2017）
動ノイズ J に関するモーメント法を定式化できる．証明の概略，推定関数の満たすべき条件等は同論文を参照されたい．
例 3.8. （2.1）の部分族である一次元 SDE モデル


 pγ

pα

1
(3.8)
dXt = a0 (Xt , γ) +αl a1l (Xt , γ) dt + exp
γk ck (Xt− ) dJt
2
l=1
k=1
は，係数の状態変数 Xt について多様な非線形性を表現可能であると同時に，γ および α の推定の各段階が凸最適化になるという点において特殊である．ここでa0 (x, γ),a1 (x, γ) = (a11 (x, γ), . . . , a1pα (x, γ)),c(x) = (c1 (x), . . . , cpγ (x)),

はその関数形が既知であるとしている．このとき，γ̂n は推定方程式
n


1−

j=1

の解であり，α̂n は
1
α̂n =
hn
·

 n


n

j=1

(Δj X)2exp(−γ · c(Xj−1 )) c(Xj−1 ) = 0hn−1

exp(−γ̂n · c(Xj−1 ))a⊗2
1 (Xj−1 , γ̂n )j=1

exp(−γ̂n · c(Xj−1 )){Δj X − hn a0 (Xj−1 , γ̂n )}a1 (Xj−1 , γ̂n )Lévy 駆動型確率微分方程式の段階的推定について

27と陽に与えられる．また，他の第一段階が凸最適化となる scale 係数の関数形の例として，
pγc(x, γ) = ( k=1
γk ck (x))−1/2 などが挙げられる．
4.

数値実験

先に述べた通り，駆動ノイズは基準化されたものを考えればよい．まず次の SDE モデルを取
り扱う．


γXt
α2
(4.1)
dXt =
−
α
X
+
X0 = 0.
dt + c(Xt− , γ)dJt ,1t
1 + Xt2(1 + Xt2 )γ
ここで，スケール係数および各パラメータの真値は以下の 2 種類を想定する．
γ(i) c(x, γ) = exp(− 1+x2 ) および θ0 = (α1,0 , α2,0 , γ0 ) = (1, 2, 3)
γ(ii) c(x, γ) = − 1+x
2 および θ0 = (α1,0 , α2,0 , γ0 ) = (1, 2, 1)

また，J1 の分布は，N IG(1, 0, t, 0), N IG(10, 0, 10t, 0), N IG(25/3, 20/3, 9/5t, −12/5t) の 3 通り
（Barndorﬀを対象とする．N IG 分布は IG 分布の正規平均尺度混合により定義される確率分布
であり，N IG(α, β, δ, μ) の特性関数は，Nielsen, 1998）
φX1 (u) = exp(−δ(α2 − (β + iu)2 ) −α2 − β 2 ) exp(iμu),で与えられる．この表現と Lévy-Kchintchin の公式から，
βE[Xt ] = t μ + δ, V ar[Xt ] = tδα2 − β 2(α2
α 2 − β 2 )3が得られる．これらを用いて，E[Jt ] = 0, E[Jt2 ] = t であることが確認できる．また，Masuda（2013, Proposition 5.4）と係数の形を考慮すれば，（4.1）は仮定 3.2 を満たしていることがわかる．
本 数 値 実 験 で は ，各 駆 動 ノ イ ズ に つ い て (Tn , hn , n) を (10, 0.05, 200), (10, 0.01, 1000),(50, 0.05, 1000), (50, 0.01, 5000), (100, 0.05, 2000), (100, 0.01, 10000) の 6 通りに分け，それぞれ
500 回ずつ独立に Euler - 丸山法に基づきパスを生成し段階的推定量 θ̂n := (α̂1,n , α̂2,n , γ̂n ) を計算した．(i) の γ についてはパラメータ空間は Θγ = [0, 10] とし，初期値は 2 として最適化を行った（他の推定量は陽に書ける）．また，パスの生成と γ̂n の最適化には R の YUIMA パッケージを用いた．表 1–3 にそれらの推定結果を示した．
（Brouste et al., 2014）• ドリフトパラメータの推定結果を見ると，表 1–3 全ての場合においてターミナル Tn の増加に応じて推定精度が向上している．また，ターミナル固定の下で観測幅が小さくなっても推定√精度の向上はほとんど見られない．これは α̂n の収束レートが Tn であること，また注意 3.5
で触れた通り X の微小時間変動においてはスケール項が支配的であることからも妥当な結果だと考えられる．
• スケールパラメータの推定精度は表 1–3 を見るとターミナル Tn の増加に応じて良くなっ√ており，これは拡散過程とは異なる傾向であるが，ここでは γ̂n の収束レートが Tn であることから自然である．観測幅 hn の観点からは，表 1, 表 3 ではターミナル固定の下で観測幅が小さくなる時，標準偏差の改善が見られない．しかし表 2 では，ターミナル固定の下では観測幅の小さい場合の推定精度が大きい場合の推定精度を優越している．この点は，N IG(δ, 0, δ, 0) は
δ → ∞ で N (0, 1) に全変動収束するため，ノイズの微小時間増分系列について正規近似が機能していると捉えれば整合的である．
• スケール係数の違いに関して見ると，表 1–3 いずれにおいても（i）より（ii）の推定精度が良

統計数理

28

第 65 巻

第 1 号 2017

表 1．L(Jt ) = N IG(1, 0, t, 0) の場合の段階的推定量 θ̂n := (α̂1,n , α̂2,n , γ̂n ) の標準平均・標本標準偏差（括弧内）．

表 2．L(Jt ) = N IG(10, 0, 10t, 0) の場合の段階的推定量 θ̂n := (α̂1,n , α̂2,n , γ̂n ) の標準平均・標本標準偏差（括弧内）．

いことが確かめられるが，これは当然パラメータの真値に依存して結果は異なる．γ̂n の漸近分散に着目すると，（ii）の設定（パラメータが線形に含まれている場合）だと漸近分散が不変測度
π0 の変化に加えて，真値の 2 乗に比例することが容易に確認できる（ケース（i）では γ̂n の漸近分（ii）の場合に真散は π0 (dy) を通じてのみ θ0 に依存する）．ここには記載していないが，実際に値を大きくすると標準偏差が悪化することを確認している．
また，図 1, 2 において L(J1 ) = N IG(25/3, 20/3, 9/5t, −12/5t), (Tn , n, hn ) = (100, 10000, 0.01)の場合（表 3）のヒストグラムおよび標準正規密度関数（実線）を示した．標準正規分布近似は概ねうまくいっている．
次に scale 項に複数のパラメータを持つ例として下記の SDE モデルを想定する:
(4.2)

dXt =

−α1 Xt +α2
1 + Xt2

dt + exp −γ0 + γ1 Xt
2(1 + Xt2 )dJt ,X0 = 0.

ここで，L(Jt ) = N IG(10, 0, 10t, 0), Θγ = [0, 50] × [0, 50], θ0 = (α1,0 , α2,0 , γ1,0 , γ2,0 ) = (1, 2, 3, 4)，

Lévy 駆動型確率微分方程式の段階的推定について

29

表 3．L(Jt ) = N IG(25/3, 20/3, 9/5t, −12/5t) の場合の段階的推定量 θ̂n := (α̂1,n , α̂2,n , γ̂n )の標準平均・標本標準偏差（括弧内）．

図 1．標準化（スチューデント化）された段階的推定量 θ̂n のヒストグラム（（i）の場合）．

図 2．標準化（スチューデント化）された段階的推定量 θ̂n のヒストグラム（（ii）の場合）．

初期値をパラメータ空間上の一様乱数とし，その他の設定は先程のものと同様とする．数値実験の実験結果も併結果を表 4 に示す．ここでは比較のため拡散過程（L(Jt ) = N (0, t) とした場合）記した．表 4 から，両モデルともにドリフトの推定精度は比較的良い一方で，scale パラメータの推定値の平均が真値とかけ離れており，今のモデルでは平均構造を無視したことで生じるバ
イアスの影響が顕著であると考えられる．この問題に対し，三段階目の推定として通常の正規

30

統計数理

第 65 巻

第 1 号 2017

表 4．L(Jt ) = N IG(10, 0, 10t, 0) および L(Jt ) = N (0, t) の場合の段階的推定量
θ̂n := (α̂1,n , α̂2,n , γ̂1,n , γ̂2,n ) の標準平均・標本標準偏差（括弧内）．

表 5．L(Jt ) = N IG(10, 0, 10t, 0) および L(Jt ) = N (0, t) の場合の適応型推定法に基づく推（括弧内）．
定量 γ̃n := (γ̃1,n , γ̃2,n ) の標準平均・標本標準偏差型疑似尤度（2.4）に第二段階で得られた α̂n をプラグインし，scale パラメータの推定値をアップ
デートすることを考える．すなわち，γ の推定値を以下の関数の最大点として新たに定義する:

n 
1 −1
1
⊗2 
S (γ) (Δj X − haj−1 (α̂n , γ))γ → −
log |Sj−1 (γ)| +.
2 j=1
hn j−1
ここで最適化の初期値には第一段階で得た推定量を用いることとする．このアイデアは拡散過（1.2）の下，程の適応型推定法（cf. Uchida and Yoshida, 2012）に基づくものである．観測頻度条件
本稿で扱っている Lévy 型 SDE の scale 推定に関しても，三段階目の推定量の漸近挙動は一段階目の推定量のそれと等しいことが示される．表 5 にアップデートした推定量 γ̃n を示した．拡
散過程と同様に，推定量の平均値に関する改善が見て取れる．標準偏差は悪化しているが，これは上述の通り漸近分散の改善が理論上見込めないことや，α̂n の推定値の揺らぎに起因するものであると考えられる．

Lévy 駆動型確率微分方程式の段階的推定について

5.

31

定理 3.4 の証明

まず本節で用いる記号を述べておく．E j−1 [·] で Ftj−1 に関する条件付き期待値を表す．任意（記号 fj−1 (θ) = f (Xtj−1 , θ) との重複の R × Θ 上の関数 f について fs (θ) = f (Xs , θ) と定義するはあるが，混乱は生じないであろう）．特に θ = θ0 の時は fs = f (Xs , θ0 ) と略記する．またM  (x, θ) := ∂α a(x, θ)T S −1 (x, γ) ∈ Rpα ⊗ Rd ,M  (x, γ) := −∂γ S −1 (x, γ) ∈ Rpγ ⊗ Rd ⊗ Rdとする．ある正定数 C が存在して十分大きな任意の n について xn ≤ Cyn となる時 xn  yn と．
書く（C は n に依存せず，現れるごとに異り得る）
5.1 裾確率評価（3.5）の証明(γ̂n , α̂n ) は，コントラスト関数 M1,n (γ) := −nhn |G1,n (γ)|2 , M2,n (α) := −nhn |G2,n (α, γ̂n )|2 に対応した M 推定量とみなすことができる．これらに対して以下の確率場を定義する:

1
J1,n (u1 ) := exp M1,n γ0 + √ u1 − M1,n (γ0 ) ,Tn1J2,n (u2 ) := exp M2,n α0 + √ u2 , γ̂n − M2,n (α0 , γ̂n ) ,Tn

ここで u1 , u2 はそれぞれ，集合 U1,n := {u1 ∈ Rpγ : γ0 + √1Tn u1 ∈ Θγ }, U2,n := {u2 ∈ Rpα :α0 + √1T u2 ∈ Θα } の元とする．J1,n (u1 ), J2,n (u2 ) の定義から，n√√Tn (γ̂n − γ0 ) ∈ argmax J1,n (u1 ),Tn (α̂n − α0 ) ∈ argmax J2,n (u2 )u1 ∈U1,nである．ゆえに，任意の r > 0 について√
P (| Tn (γ̂n − γ0 )| > r) ≤ P

u2 ∈U2,nsupu1 :|u1 |>r,u1 ∈U1,n

J1,n (u1 ) ≥ J1,n (0) ≡ 1√が成り立つ（ Tn (α̂n − α0 ) についても同様）．すなわち，本定理の主張は J1,n (·), J2,n (·) の確率評価に帰着する．その評価のため，統計的確率場の多項式型大偏差不等式（Yoshida, 2011）の理論を援用する．証明は複数のモーメント評価の検証を伴うが，今モーメントの存在は必要なだけ，したがって一般性を失うことなく d = pγ = pα = 1
仮定している状況であり（仮定 2.1, 3.1, 3.2）としてよいので，以下その場合で見ていく．ここでは疑似尤度 Mn が混合収束率を持つことにより対応する統計的確率場の直接評価はできないため，疑似スコア関数から構成される確率場を対象とした多項式型大偏差不等式を適用する．

仮定の下，G1,n (·) に関しては，および正定数 M, が存在して



√
K
k
K
max sup E[| Tn G1,n (γ0 )| ] + E sup |∂γ G1,n (γ)|<∞
k∈{1,2,3} n∈N

 

(5.1)
sup
n∈N

E

sup |

(K > 0),

γ∈Θγ

√

γ∈Θγ



M +
Tn (G1,n (γ) − G∞
γ (γ))|√+ E[| Tn (∂γ G1,n (γ0 ) − Iγ )| ]が成り立つことを示せば十分である．任意の k ∈ {1, . . . , 4} に対して1 
−1∂γ Sj−1{(Δj X)2 − hn Sj−1 }Tn j=1n

G1,n (γ0 ) =


M

<∞

統計数理

32

2 k
∂γ
n j=1
n

∂γk G1,n (γ) =第 65 巻

∂γ cj−1 (γ)
cj−1 (γ)

第 1 号 20171  k+1 −1∂γ Sj−1 (γ)(Δj X)2
Tn j=1n+であることに注意すれば，（5.1）
は Sobolev の不等式と Masuda and Uehara（2017, Lemma 5.3）の証明内の評価からしたがう．また，G1,n (γ) =nn2  ∂γ cj−1 (γ) 22  ∂γ cj−1 (γ)
2
(γ)
−
c
)
−
(c
{(Δj X)2 − hn c2j−1 }j−1
j−1
n j=1 c3j−1 (γ)
Tn j=1 c3j−1 (γ)

であり，Sobolev の不等式と Masuda（2013, Lemma 4.3），Masuda and Uehara（2017, Lemmaの証明内の評価より，任意の K > 0 について
5.3）


√
∞
K
sup E sup | Tn (G1,n (γ) − Gγ (γ))|
<∞
γ∈Θγ

n∈N

を得る．同様にして，nn 
4  (∂γ cj−1 )22  (∂γ2 cj−1 )cj−1 − 3(∂γ cj−1 )2
2
2
−
{(Δ
X)
−
h
c
}
j
n j−1
n j=1
c2j−1
Tn j=1
c4j−1∂γ G1,n (γ0 ) =なので√
sup E[| Tn (∂γ G1,n (γ0 ) − Iγ )|K ] < ∞n∈Nを得る．以上より，Masuda（2013, Theorem 3.5）を適用して√sup P (| Tn (γ̂n − γ0 )| > r) ≤n√Tn (γ̂n − γ0 ) の裾確率評価

CK
rK

を得る．
Taylor の公式から k ∈ {0, 1, 2, 3} について√√Tn ∂αk G2,n (α, γ̂n ) = Tn ∂αk G2,n (α, γ0 ) 1√+∂γ ∂αk G2,n (α, γ0 + u(γ̂n − γ0 ))du [ Tn (γ̂n − γ0 )]0√となる．すでに任意の K > 0 に対して Tn (γ̂n − γ0 ) の LK - 有界性を確保しているため，右辺
第 2 項は LK - 有界である．ゆえに G2,n (α, γ0 ) については，任意の K > 0 とある正定数 M, に
対して



√
max sup E[| Tn G2,n (α, γ0 )|K ] + E sup |∂αk G2,n (α, γ0 )|K< ∞,k∈{1,2,3} n∈N
α∈Θα

 

√
√
∞
M +
M
+ E[| Tn (∂α G2,n (α, γ0 ) − Iα )| ] < ∞sup E sup | Tn (G2,n (α, γ0 ) − Gα (α))|
n∈N

α∈Θα

を示せば十分である．今の場合
1  −1S ∂α aj−1 (α, γ0 )(Δj X − hn aj−1 )Tn j=1 j−1n

G2,n (α, γ0 ) =1  −1S ∂α aj−1 (α, γ0 )(aj−1 − aj−1 (α, γ0 )),n j=1 j−1n+Lévy 駆動型確率微分方程式の段階的推定について33

1  −1 k+1
S ∂α aj−1 (α, γ0 )(Δj X − hn aj−1 (α, γ0 ))Tn j=1 j−1n∂αk G2,n (α, γ0 ) =1  −1 k−1
S ∂α (∂α aj−1 (α, γ0 ))2
n j=1 j−1n−√であることに注意すれば， Tn (γ̂n − γ0 ) のときと同様にして裾確率評価√
CK
sup P (| Tn (α̂n − α0 )| > r) ≤ K
rnを得る．以上で（3.5）
が示された．
5.2 漸近標準正規性
（3.6）
の証明
√
P
から Tn (θ̂n − θ0 ) = Op (1)，特に一致性 θ̂n −→ θ0 がしたがうことに注意す
裾確率評価（3.5）
る．まず




n
1 
(5.2)
Gn (θ0 ) := G2,n (α0 , γ̂n ) +Mj−1 [∂γ aj−1 ] [γ̂n − γ0 ], G1,n (γ0 )
n j=1

について
√

(5.3)

L

Tn Gn (θ0 ) −→ N (0, Σ)√を示す． Tn Gn (θ0 ) の第一成分について，Taylor 展開から√√√
Tn G2,n (α0 , γ̂n ) = Tn G2,n (θ0 ) + ∂γ G2,n (θ0 )[ Tn (γ̂n − γ0 )]
 1
√
1
+√
∂γ⊗2 G2,n (α0 , γ0 + u(γ̂n − γ0 ))du [( Tn (γ̂n − γ0 ))⊗2 ]Tn0 1 ⊗2を得る．Sobolev の不等式より 0 ∂γ G2,n (α0 , γ0 + u(γ̂n − γ0 ))du = Op (1) が示され，右辺第三項は op (1) となる．また，モーメント条件の下，a の Lipschitz 連続性および Masuda (2013,
Lemma 4.5) から E j−1 [| j (as − aj−1 )ds|2 ]  h3n ，さらに j cs− dJs のマルチンゲール性により，
E j−1 [ j cs− dJs ] = 0, E j−1 [| j cs− dJs |2 ]  hn が従うため，Genon-Catalot and Jacod（1993,より
Lemma 9）∂γ G2,n (θ0 ) =

n
n
1 
1 

∂γ Mj−1[Δj X − hn aj−1 ] −Mj−1 [∂γ aj−1 ]Tn j=1n j=1

1 
Mj−1 [∂γ aj−1 ]n j=1n= op (1) −が成り立つ．これを踏まえて√√√
Tn Gn (θ0 ) = ( Tn G2,n (θ0 ), Tn G1,n (γ0 )) + op (1)

n
1  
= √
Mj−1 [cj−1 Δj J],
Tn j=1
n
1 
−1

√
{hn trace(Sj−1∂γ Sj−1 ) − Mj−1[(cj−1 Δj J)⊗2 ]}Tn j=1+ op (1)を得る．最右辺第一項は各 n について (Ftj ) - マルチンゲール差分になっていることに注意する．

統計数理

34

第 65 巻

第 1 号 2017

モーメント評価 E[|Jhn |q ]  hn により，最右辺第一項について Lyapunov 条件が容易に示される．以下 u ∈ Rpα , v ∈ Rpγ を任意に固定し，二次特性量を計算する．
• E[(Δj J)⊗2 ] = hn Ir よりn
1  j−1E[(Mj−1⊗ Mj−1)[cj−1 Δj J, cj−1 Δj J]] [u, u ]Tn j=11 Mj−1 [u, cj−1 ] × Mj−1[u , cj−1 ]
n j=1
n

=

−→ Σα [u, u ].
P• 第二要素について，まずn1  j−1⊗2 ⊗2E[{Mj−1 [hn Sj−1 − (cj−1 Δj J) ]} ] [v, v  ]Tn j=1


n
hn 


=−
(Mj−1 ⊗ Mj−1 )[Sj−1 , Sj−1 ] [v, v  ]
n j=1


n
1  j−1


⊗2
⊗2
+
E
[(Mj−1 ⊗ Mj−1 )[(cj−1 Δj J) , (cj−1 Δj J) ]] [v, v  ]Tn j=1n
1  j−1


⊗2
⊗2
=
E
[(Mj−1 ⊗ Mj−1 )[(cj−1 Δj J) , (cj−1 Δj J) ]] [v, v  ] + op (1)Tn j=1に注意する．Lévy-Khintchin 表現



T

T

log E[eiu Jt ] = t
Rr(eiu z − 1 − iuT z)ν0 (dz)を偏微分して
(i )

(i )

(i )

(i )

(i )

(i )

E[Jt 1 Jt 2 Jt 3 ] = νi1 ,i2 ,i3 (3)t(i )E[Jt 1 Jt 2 Jt 3 Jt 4 ] = νi1 ,i2 ,i3 ,i4 (4)t + O(t2 )を得る．これを適用することで，n1  j−1


⊗2
⊗2
E
[(Mj−1 ⊗ Mj−1 )[(cj−1 Δj J) , (cj−1 Δj J) ]] [v, v  ]Tn j=1
1  j−1 E[Mj−1 [v, (cj−1 Δj J)⊗2 ] × Mj−1
[v  , (cj−1 Δj J)⊗2 ]]Tn j=1
 d
 r
 r

n
 (i s) (i t)1  j−1 
(i1 i2 )
(s)
(t)
1
2
E
(Mj−1 [v])cj−1 (Δj J)
cj−1 (Δj J)=
Tn j=1s=1
t=1
i1 ,i2 =1
 r
 r
⎤
d



 
 
 


(i1 s )
(i2 t )
×
(Mj−1
[v  ])(i1 i2 )cj−1
(Δj J)(s )cj−1(Δj J)(t ) ⎦n=i1 ,i2 =1

=

⎧
n ⎨


1
Tn j=1 ⎩s =1r
s,t,s ,t =1

t =1




(E j−1 [(Δj J)(s) (Δj J)(t) (Δj J)(s ) (Δj J)(t ) ])Lévy 駆動型確率微分方程式の段階的推定について

⎛

d


×⎝

(i s) (i t)



1
2
(Mj−1
[v])(i1 i2 ) cj−1
cj−1
× (Mj−1[v  ])(i1 i2 )

35

⎞⎫
⎬

(i1 s ) (i2 t )cj−1
cj−1 ⎠

i1 ,i2 ,i1 ,i2 =1

1 ⎝
n j=1
n

=

⎛

r


(·s )⎭⎞(·t )νs,t,s ,t (4) × Mj−1
[v, cj−1 , cj−1 ] × Mj−1
[v  , cj−1 , cj−1 ]⎠
(·s)

(·t)

s,t,s ,t =1−→ Σγ [v, v  ].
P

• E[|Jhn |2 ] = Op (hn ) より
1

−1
(Mj−1
[cj−1 Δj J] ⊗ trace(Sj−1∂γ Sj−1 ))[u, v] = op (1)n j=1nであり，さらに−n#
1  j−1 " Mj−1 [cj−1 Δj J] ⊗ Mj−1
E[(cj−1 Δj J)⊗2 ] [u, v]
Tn j=1 dn
1  j−1E(Mj−1)(·s) (cj−1 Δj J)(s)=−Tn j=1
s=1


 d


(·tu) "
⊗2 #(tu)(cj−1 Δj J)(Mj−1 )[u, v]
⊗

⎡

=−

t,u=1

r
d
n

1  j−1 ⎣ 
E
Δj J (k1 ) Δj J (k2 ) Δj J (k3 )Tn j=1s,t,u=1
k1 ,k2 ,k3⎤


×cj−11 cj−12 cj−13 ((Mj−1)(·s) ⊗ (Mj−1)(·tu) )[u, v]⎦(sk ) (tk ) (uk )

=−

r
n
d

1 
(sk ) (tk ) (uk )νk1 ,k2 ,k3 (3)cj−11 cj−12 cj−13 ((Mj−1)(·s) ⊗ (Mj−1)(·tu) )[u, v]
n j=1
s,t,u=1
k1 ,k2 ,k3P−→ Σα,γ [u, v].

以上の下，マルチンゲール中心極限定理（e.g. Dvoretzky, 1972）から（5.3）がしたがう．
これまでの評価を踏まえて，行列"#Iα (θ) :=∂α M  (x, θ)[a(x, θ0 ) − a(x, θ)] − S −1 (x, γ)[(∂α a(x, θ))⊗2 ] π0 (dx),"Iγ (γ) := trace 2(S −1 ∂γ SS −1 ∂γ SS −1 )(x, γ)S(x, γ0 ) − (S −1 ∂γ SS −1 ∂γ S)(x, γ)#+(S −1 ∂γ⊗2 S)(x, γ) − (S −1 ∂γ⊗2 SS −1 )(x, γ)S(x, γ0 ) π0 (dx),'&
B(θ) := −∂γ M  (x, θ)[a(x, θ0 ) − a(x, θ)] + S −1 (x, γ0 )[∂α a(x, θ), ∂γ a(x, θ)] π0 (dx)に対して一様収束P

sup |∂γ G1,n (γ) − Iγ (γ)| −→ 0,

γ∈Θγ

P

sup |∂α G2,n (θ) − Iα (θ)| −→ 0,
θ∈Θ

統計数理

36

第 65 巻

第 1 号 2017P

sup |∂γ G2,n (θ) − B(θ)| −→ 0θ∈Θが成り立つことがわかる．これと Gn の定義（5.2）,（5.3）より√∂α G2,n (θ̂n ) ∂γ G2,n (θ̂n )
α̂n − α0
L
−1/2
(5.4)
Tn Σ
−→ N (0, Ip )γ̂n − γ00
∂γ G1,n (γ̂n )を得る．例えば Σ̂α,n について Taylor 展開を適用すると
Σ̂α,n [u, u ] =1  −1S [∂α aj−1 [u], ∂α aj−1 [u ]] + op
n j=1 j−1
n

1
√
Tn

が得られる．これにエルゴード定理と Genon-Catalot and Jacod（1993, Lemma 9）を適用すれば
PPP
Σ̂α,n −→ Σα が従う．同様にして Σ̂γ,n −→ Σγ および Σ̂α,γ,n −→ Σα,γ も示され，（3.5）,（5.4）とSlutsky の補題より
（3.6）
を得る．
5.3 モーメント収束
（3.7）
の証明
P
P
（3.6）により，すでに示した Σ̂n −→ Σ および În −→ I と√L
Tn (θ̂n − θ0 ) −→ N (0, I −1 Σ(I −1 ) )√が成り立つ．（3.5）
から任意の q > 0 に対する {| Tn (θ̂n − θ0 )|q }n∈N の一様可積分性が成り立つ
ので（3.7）が得られる．

謝辞

丁寧に査読していただき，多くの改善点の指摘および有益なコメントを下さった査読者に御礼申し上げます．本研究は JST CREST，JSPS 科研費 JP26400204（増田）の助成を受けたもの
です．

参

考

文

献

Barndorﬀ-Nielsen, O. E. (1998). Processes of normal inverse Gaussian type, Finance and Stochastics,
2(1), 41–68.
Brouste, A., Fukasawa, M., Hino, H., Iacus, S., Kamatani, K., Koike, Y., Masuda, H., Nomura, R.,
Ogihara, T., Shimizu, Y., Uchida, M. and Yoshida, N. (2014). The YUIMA project: A computational framework for simulation and inference of stochastic diﬀerential equations, Journal of
Statistical Software, 57(4), 1–51.
Dvoretzky, A. (1972). Asymptotic normality for sums of dependent random variables, Proceedings of
the Sixth Berkeley Symposium on Mathematical Statistics and Probability (University of California, Berkeley, California, 1970/1971 ), Vol. II : Probability Theory, 513–535, University ofCalifornia Press, Berkeley, California.
Genon-Catalot, V. and Jacod, J. (1993). On the estimation of the diﬀusion coeﬃcient for multidimensional diﬀusion processes, Annales de l’Institut Henri Poincaré (B ) Probabilités et Statistiques, 29, 119–151.
Kamatani, K. and Uchida, M. (2015). Hybrid multi-step estimators for stochastic diﬀerential equations
based on sampled data, Statistical Inference for Stochastic Processes, 18(2), 177–204.
Masuda, H. (2013). Convergence of Gaussian quasi-likelihood random ﬁelds for ergodic Lévy driven

Lévy 駆動型確率微分方程式の段階的推定について

37

SDE observed at high frequency, Annals of Statistics, 41, 1593–1641.
Masuda, H. and Uehara, Y. (2017). Two-step estimation of ergodic Lévy driven SDE, Statistical Inference for Stochastic Processes, 20(1), 105–137.
Sato, K. (1999). Lévy Processes and Infinitely Divisible Distributions, Cambridge University Press,
Cambridge, U.K.
Uchida, M. and Yoshida, N. (2012). Adaptive estimation of an ergodic diﬀusion process based on sampled data, Stochastic Processes and Their Applications, 122(8), 2885–2924.
Yoshida, N. (2011). Polynomial type large deviation inequalities and quasi-likelihood analysis forstochastic diﬀerential equations, Annals of the Institute of Statistical Mathematics, 63(3), 431–
479.

38

Proceedings of the Institute of Statistical Mathematics Vol. 65, No. 1, 21–38 (2017)

On Stepwise Estimation of Lévy Driven Stochastic Diﬀerential Equation
Yuma Uehara and Hiroki Masuda
Graduate School of Mathematics, Kyushu University

We consider estimation of a non-Gaussian Lévy driven stochastic diﬀerential equation. Under high-frequency data and exponential ergodicity, we propose the stepwise estimation procedure based on a Gaussian quasi-score function: ﬁrst we estimate the scaleparameter while ignoring the drift coeﬃcient, and then focus on the drift parameter byplugging in the estimated scale parameters, and derive the asymptotic normality andthe tail probability estimate of the proposed estimators. This stepwise strategy not only
reduces computational cost but may also stabilize estimate accuracy. Unlike the diﬀusion case, the asymptotic covariance matrix associated with the drift parameter takes adiﬀerent form when there is a common parameter in the coeﬃcients.

Key words: Ergodicity, Gaussian quasi-score function, high-frequency sampling, Lévy driven stochastic
diﬀerential equation, stepwise estimation.