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


特集「高頻度金融データに基づく統計的推測とモデリング」［研究詳解］

高頻度データに基づく確率微分方程式モデルのハイブリッド推定
内田 雅之1,2,3（受付 2016 年 6 月 28 日；改訂 10 月 14 日；採択 10 月 17 日）

要旨

本論文では，高頻度データを用いた確率微分方程式モデルの未知パラメータにおけるハイブ
リッド推定について解説する．最適な収束率より遅いベイズ型推定量を初期推定量とした，マルチステップ推定法や適応的最尤型推定法を説明し，それらの推定量の漸近的性質について述
べる．エルゴード的拡散過程，非エルゴード的拡散過程，微小拡散過程の 3 種類の拡散モデル
を取り扱い，数値実験により提案した推定量の漸近的挙動を考察する．
キーワード：拡散過程，確率微分方程式，最尤型推定量，ハイブリッド推定量，ベイズ型推定量，マルチステップ推定量．
1.

はじめに

情報技術の発展により，金融などの時系列データが高頻度に観測可能となり，連続時間確率過程モデルの統計解析の需要が増してきている．特に，確率微分方程式で定義された拡散型確率過程モデルの統計解析は今までに盛んに研究されている（Kutoyants, 1984, 1994, 2004; Prakasa
Rao, 1988; Yoshida, 1992a, 1992c; Genon-Catalot and Jacod, 1993; Bibby and Sørensen, 1995;．非線形で非定常な連続時間確率
Kessler, 1997; Prakasa Rao, 1999; Iacus, 2008; Sørensen, 2012）過程である拡散過程モデルは表現力豊かな統計モデルである反面，尤度関数を明示的に求めることができない難点があるが，マリアバン解析や疑似尤度解析を援用することにより，統計的漸近理論を展開することが可能である．高頻度データにおける統計的モデリングやモデル選択のための情報量規準についても精力的に研究されていて，実際にデータ解析を行う上での問題点が指摘されている．例えば，エルゴード的拡散過程モデルでは，ドリフトパラメータとボラ
ティリティパラメータという 2 種類の未知パラメータを推測する必要があるが，疑似尤度関数が複雑で，パラメータ空間が高次元である場合，疑似尤度関数の最大化によって疑似最尤推定量（または最尤型推定量）を計算することが困難となることが多々ある．最尤型推定量の導出に失敗して，誤った値を算出した場合，情報量規準に基づくモデル選択に影響が出ることも知られている．このことから，理論的に保障された推定量であっても，実際にデータ解析において正しい推定量を算出できるかは別の問題となる．本論文では，確率微分方程式モデルのパラメー
タを効率よく推定するための統計的手法の一つであるハイブリッド推定法について解説する．
1

大阪大学大学院 基礎工学研究科：〒 560–8531 大阪府豊中市待兼山町 1–3
大阪大学 数理・データ科学教育研究センター（MMDS）：〒 560–8531 大阪府豊中市待兼山町 1–3
3
CREST
2

統計数理

40

第 65 巻

第 1 号 2017

最初に，I.I.D. モデルにおけるハイブリッドマルチステップ（HMS）推定量について考察する．
ln (θ) は θ ∈ R に関して滑らかな対数尤度関数とする．q ∈ (0, 1/2], J = [− log2 q] とすると，
2J−1 q ≤ 1/2 < 2J q となる．初期推定量 θ̂(0) は，M > 0 に対し，
sup Eθ∗ [|nq (θ̂(0) − θ∗ )|M ] < ∞

(1.1)

n

∗

を満たすとする．ここで，θ はパラメータ θ の真値で，Eθ∗ は真のモデル（分布）Pθ∗ の下での期待値を表す．k = 1, . . . , J に対し，k - ステップ推定量 θ̂(k) をθ̂(k) = θ̂(k−1) − [∂θ2 ln (θ̂(k−1) )]−1 [∂θ ln (θ̂(k−1) )]∂と定義する．ただし，∂θ = ∂θ, ∂θ2 = ∂θ∂ 2 で，∂θ2 ln (θ̂(k−1) ) = 0 とする．θ は 1 次元パラメータであることに注意する．テイラーの定理より，∂θ ln (θ∗ ) = ∂θ ln (θ̂(k−1) ) + ∂θ2 ln (θ̂(k−1) )[θ∗ − θ̂(k−1) ] + Rn [(θ∗ − θ̂(k−1) )2 ], 1(1 − t)∂θ3 ln (θ̂(k−1) + t(θ∗ − θ̂(k−1) ))dt
Rn =
0

だから，
θ̂(k) = θ̂(k−1) − [∂θ2 ln (θ̂(k−1) )]−1 [∂θ ln (θ∗ ) − ∂θ2 ln (θ̂(k−1) )[θ∗ − θ̂(k−1) ] − Rn [(θ∗ − θ̂(k−1) )2 ]]= θ̂(k−1) − [∂θ2 ln (θ̂(k−1) )]−1 [∂θ ln (θ∗ )] + (θ∗ − θ̂(k−1) ) + [∂θ2 ln (θ̂(k−1) )]−1 Rn [(θ∗ − θ̂(k−1) )2 ]となる．ゆえに，θ̂(k) − θ∗ = −[∂θ2 ln (θ̂(k−1) )]−1 [∂θ ln (θ∗ )] + [∂θ2 ln (θ̂(k−1) )]−1 Rn [(θ̂(k−1) − θ∗ )2 ].
√特に，2q ≤ 1/2 のとき（初期推定量 θ̂(0) が n - 一致性が保証されていない場合），
−1  2q
 
−1
1 2
1 2
n
1
n2q (θ̂(1) − θ∗ ) = −∂θ ln (θ̂(0) )∂θ ln (θ∗ ) +∂θ ln (θ̂(0) )Rn [(nq (θ̂(0) − θ∗ ))2 ].
nnnn

したがって， √1n ∂θ ln (θ∗ ) や n1 Rn の可積分性などの正則条件の下で，sup Eθ∗ [|nq (θ̂(0) − θ∗ )|M ] < ∞=⇒nsup Eθ∗ [|n2q (θ̂(1) − θ∗ )|M ] < ∞nが言える．これを繰り返すことによって，マルチステップ推定量は，ある正則条件の下でsup Eθ∗ [|n2q (θ̂(1) − θ∗ )|M ] < ∞

=⇒

n

2

sup Eθ∗ [|n2 q (θ̂(2) − θ∗ )|M ] < ∞,n..
.
sup Eθ∗ [|n2J−2

sup Eθ∗ [|n2

J−1

q

n

n

q

(θ̂(J−2) − θ∗ )|M ] < ∞=⇒(θ̂(J−1) − θ∗ )|M ] < ∞

=⇒

(J)

sup Eθ∗ [|n2
n

J−1

q

(θ̂(J−1) − θ∗ )|M ] < ∞,√sup Eθ∗ [| n(θ̂(J) − θ∗ )|M ] < ∞nとなり，J - ステップ推定量 θ̂ は漸近有効になる．注目すべきことは，（1.1）からわかるように，この初期推定量は最適な収束率をもつことを仮定していないが，マルチステップ推定量は√最適な収束率をもつことである．よく知られているワンステップ推定量は初期推定量が n - 一
致性をもつことを仮定している．このことからも，初期推定量の条件が緩いマルチステップ推定量がより有用であることがわかる．
（1.1）を満足する初期推定量として，最尤型推定量とベイズ型推定量を考えることができるが，

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

41

最尤型推定量の導出には最適化が伴い，初期値の設定が重要な鍵となる．一方で，ベイズ型推定量は，マルコフチェーンモンテカルロ（MCMC）法を用いた数値計算により導出される．最尤型推定量の数値計算に比べて，ベイズ型推定量の数値計算は初期値の影響が少ない．本論文では，ベイズ型推定量を初期推定量として用いて，最尤型推定（マルチステップ推定）
により，漸近有効推定量を導出する方法をハイブリッド推定法と呼ぶことにする．特に，初期ベイズ推定量を用いたマルチステップ推定法をハイブリッドマルチステップ推定法と表現することにする．
本論文の構成は次の通りである．第 2 節では，エルゴード的拡散過程について 3 種類のハイ
ブリッドマルチステップ推定量を提案し，モーメント収束性などの漸近的性質について論じる．
また，具体例とシミュレーション結果を示す．第 3 節では，非エルゴード的拡散過程について
ハイブリッドマルチステップ推定量を提案し，モーメント収束性を含めた漸近的性質について説明する．ハイブリッドマルチステップ推定量の具体例とシミュレーション結果を述べる．第4 節では，微小拡散過程について，適応的最尤型推定量と適応的ベイズ型推定量について考察
し，それらを混合したハイブリッド推定量を提案し，その漸近的性質について概説する．ハイブリッド推定量の具体例と数値実験結果について言及する．
言うまでもなく，推定量の極限分布は統計的漸近理論には必須であり，モデル選択問題における情報量規準を正当化するために，推定量のモーメントの収束性を示すことは不可欠である．幾，Kutoyants（1984,
多の統計的漸近理論の研究があるが，特に Ibragimov and Has’minskii（1981），Yoshida（2011）の精読を薦める．拡散過程のモデル選択に関しては，Uchida and Yoshida2004）（2001, 2004a, 2006, 2016），Uchida（2010）を参照．
2.

エルゴード的拡散過程の結果に従って，エルゴード的拡散過程のドリフトパ
本節は，Kamatani and Uchida（2015）
ラメータおよびボラティリティパラメータのハイブリッドマルチステップ推定量を導出し，その漸近的性質を考察する．
2.1 モデルと仮定
次の確率微分方程式で定義される d - 次元拡散過程を考える．すなわち，(2.1)

dXt = b(Xt , β)dt + a(Xt , α)dwt ,t ≥ 0,

X 0 = x0 .

ここで，w は r - 次元標準ウィーナー過程，x0 は初期値，θ = (α, β) ∈ Θ1 × Θ2 = Θ で，Θ1 ，Θ2はそれぞれ，Rm1 ，Rm2 のコンパクトで凸な部分集合とする．
さらに，b : Rd × Θ1 → Rd , a : Rd × Θ2 → Rd ⊗ Rr とする．θ∗ = (α∗ , β ∗ ) は θ の真値とする．θ∗ ∈ Int(Θ) であり，パラメータ空間は局所的にリプシッツ境界をもつと仮定する．Adams
を参照．データ Xn = (Xtni )0≤i≤n は離散観測され，tni = ihn とする．p はand Fournier（2003）
整数で p ≥ 2 とする．n → ∞ のとき，hn → 0，nhpn → 0 であるとし，十分大きな n に対し
n0 ≤ nhn を満たす 0 ∈ (0, (p − 1)/p) が存在するとする．hn = 1/n1/4 の例を考える．n → ∞のとき nhn → ∞，nh5n → 0 となるが，nh4n は 0 に収束しない．この場合 p ≥ 5 となる p を選択する必要がある．
Kutoyants（1984, 2004）に代表されるように連続観測データに基づくエルゴード的拡散過程モデルによる統計推測はこれまで精力的に研究されてきており，離散観測データからのエルゴー
ド的拡散過程モデルによる統計解析についても，今なお多くの研究者によって研究され続けて，Yoshida（1992c, 2011），Bibby and Sørensen（1995），Kessler
いる．Prakasa Rao（1983, 1988）（1995, 1997），Gobet（2002），Uchida（2010），Uchida and Yoshida（2001, 2011, 2012, 2014），

統計数理

42

第 65 巻

第 1 号 2017

Masuda（2013a, 2013b）
を参照．また，ジャンプ付拡散過程とレヴィ過程については Shimizu and，Shimizu（2006），Ogihara and Yoshida（2011），Masuda（2013a, 2013b）を参照．
Yoshida（2006）
特に Kessler（1995, 1997）
は，nhpn → 0 の仮定の下での 1 次元拡散過程の疑似対数尤度関数
Up,n (α, β) を提案している．k0 = [ p2 ]，A(x, α) = aa (x, α) とおく．ここで  は転置とする．Lθ


は拡散過程
（2.1）
の生成作用素，すなわち，Lθ = di=1 bi (x, β)∂i + 12 di,j=1 Aij (x, α)∂i ∂j とする．
pd

ΔXi = Xtni − Xtni−1 , Ai−1 (α) = A(Xtni−1 , α), bi−1 (β) = b(Xtni−1 , β) とする．→ と → はそれぞれ，確率収束と分布収束を表す．同じサイズの行列 A，B に対し，A⊗2 = AA ，B[A] = tr(BA )を定義する．| · | はユークリッド距離を表し，行列 A に対し，||A||2 = tr(AA ) とする．
疑似対数尤度関数 Up,n (θ) は次のように定義される．



k
k0
n
0


1
(j)
(k0 )
(j)
−1
j
⊗2
j
Up,n (θ) = −hn Di−1 (θ) [(Xtni − ri−1 (hn , θ)) ] +hn Ei−1 (θ) .
hn
2 i=1
j=0
j=0
k 0

hjn j
n
j=0 j! Lθ fm (Xti−1 ), fm (x) = xm である．さらに j = 0, 1, . . . に対
(j)
(j)
(j)
し，Di−1 (θ) := D (Xtni−1 , θ)，Ei−1 (θ) := E (j) (Xtni−1 , θ) を次のように定義する．q, r = 1, . . . , d(k0 )ここで (ri−1(hn , θ))m =に対し，(Ξ(l) (hn , x, θ∗ ))qr =

l


hvn

v=0

l−v

hw
n
w=0

w!

(v)

Lw
θ ḡqr,x,θ (x),(0)

ḡqr,x,θ (y) = (yq − xq )(yr − xr ),(j)

Ljθ fr (x)Lj fq (x)− (yr − xr ) θj!
j!
Lrθ fq (x) Lsθ fr (x), (1 ≤ j ≤ l).
r!
s!

ḡqr,x,θ (y) = −(yq − xq )+
r≥1,s≥1

r+s=j

(l)

(Ξ (hn , x, θ))qr =

l

j (j)
j=0 hn γqr (x, θ) とおけば，(0)
γqr(x, θ) = 0,(1)
γqr(x, θ) = Aqr (x, α),


d

1
(2)
γqr (x, θ) ={(∂xj bq (x, β))Ajr (x, α) + (∂xj br (x, β))Ajq (x, α)}Lθ Aqr (x, α) +
2
j=1
(l)
−1
となる．オーダー l までの (h−1のテイラー展開はn Ξ (hn , x, θ))

l

j
(j)
(x, θ) であり，
j=0 hn D

D(0) (x, θ) = A−1 (x, α),D(1) (x, θ) = −A−1 (x, α)γ (2) (x, θ)A−1 (x, α),D(2) (x, θ) = {(A−1 (x, α)γ (2) (x, θ))2 − A−1 (x, α)γ (3) (x, θ)}A−1 (x, α)
l
(l)
j
(j)
となる．同様に，オーダー l までの log det(h−1(x, θ)n Ξ (hn , x, θ)) のテイラー展開は
j=0 hn Eであり，E (0) (x, θ) = log det(A(x, α)),E (1) (x, θ) = tr(A−1 (x, α)γ (2) (x, θ)),1
E (2) (x, θ) = tr(2A−1 (x, α)γ (3) (x, θ) − {(A−1 (x, α)γ (2) (x, θ))2 })2

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

43となる．疑似対数尤度関数 Up,n (α, β) の詳細については Kessler（1997）や Uchida and Yoshida（2012）を参照．Kessler（1997）
は同時最尤型推定量 θ̂p,n = (α̂p,n , β̂p,n ) を Up,n (α̂p,n , β̂p,n ) =supα,β Up,n (α, β) と定義し，n → ∞，nhpn → 0 のとき，次のような漸近正規性を持つこと
を示した．
√
√
d
( n(α̂p,n − α∗ ), nhn (β̂p,n − β ∗ )) → (ζ1 , ζ2 ) ∼ Nm1 +m2 (0, Γ(θ∗ )−1 ).
ここで，Γ(θ∗ ) =
1
2


∗
Γij
1 (θ ) =∗Γij
2 (θ ) =∗(Γij
1 (θ ))1≤i,j≤m1

0

0

∗
(Γij
2 (θ ))1≤i,j≤m2


Rd

Rd

,

tr{A−1 (∂αi A)A−1 (∂αj A)(x, α∗ )}μθ∗ (dx),(∂βi b(x, β ∗ )) A(x, α∗ )−1 ∂βj b(x, β ∗ )μθ∗ (dx).

しかし，先述した通り，Θ の次元が大きいとき，同時最尤型推定量の導出は困難となる．Yoshida（1992c）は nh3n → 0 という仮定での適応的最尤型推定を考察した．適応的最尤型推定とはドリ
フトパラメータと拡散係数パラメータを別々に推定する有効な方法である．nhpn → 0 についてを参照．Uchida and Yoshida（2012, 2014）は Ibragimov-Has’minskii-Kutoyants
は Kessler（1995）の手法（Ibragimov and Has’minskii, 1981; Kutoyants, 1984, 2004）と Yoshida（2011）の結果を応用して，モーメント収束性を有する 3 種類の適応的最尤型推定量と適応的ベイズ型推定量を提] とする．ここで l0 ≤ k0 ≤ l0 + 1 であることに注案した．k0 ，l0 は整数で，k0 = [ p2 ]，l0 = [ p−1
2
(l0 )
(k )
, β̂p,n0 は k = 1, 2, . . . , k0 に対し，意する．タイプ I の適応的最尤型推定量 α̂p,n
(k−1)
(k)
(k−1)
Up,n (α̂p,n, β̂p,n) = sup Up,n (α̂p,n
, β),
β∈Θ2

(k)
(k)
, β̂p,n) =Up,n (α̂p,n(k)
sup Up,n (α, β̂p,n
)

α∈Θ1

(0)
(0)
(0)
(0)
(0)
で定義される．ここで，α̂p,n= α̂n は Un (α̂n ) = supα∈Θ1 Un (α) を満たし，

Un(0) (α) = −1  −1 −1{hn Ai−1 (α)[(ΔXi )⊗2 ] + log det(Ai−1 (α))}
2 i=1
n

である．
次に適応的ベイズ型推定について考える．事前分布 π1 (α)，π2 (β) は連続でinf π1 (α) > 0,α∈Θ1inf π2 (β) > 0

β∈Θ2

(l0 )
(k0 )
を満たすとする．タイプ I の適応的ベイズ型推定量 α̃p,n，β̃p,nは次のように定義される．
k = 1, 2, . . . , l0 に対し，
(k−1)

(k−1)
α̃p,n
=

Θ1

α exp{Hp,n(k−1)(α, β̃p,n)}π1 (α)dα

(k−1)
(k−1)
exp{Hp,n (α, β̃p,n )}π1 (α)dα
Θ1

(k)

,

(l )

(l0 )
=
α̃p,n

Θ1

0
α exp{Up,n (α, β̃p,n)}π1 (α)dα

Θ1

(l )

0
exp{Up,n (α, β̃p,n)}π1 (α)dα(l +1)

0k0 = l0 のとき，β̃p,n

Θ2

(k)
β̃p,n
=

(k−1)

β exp{H̃p,n (α̃p,n, β)}π2 (β)dβ

(k)
(k−1)
exp{H̃p,n (α̃p,n , β)}π2 (β)dβ
Θ2
(l )

,

(l0 +1)
β̃p,n
=

Θ2

0
β exp{Up,n (α̃p,n, β)}π2 (β)dβ

Θ2

(l )

0
exp{Up,n (α̃p,n, β)}π2 (β)dβは不要であることに注意する．さらに k = 1, 2, . . . , l0 に対し

.

,

統計数理

44

第 65 巻H(0)p,n (α, β) =1nUn(0) (α),

2
1− p

1

H(k)
p,n (α, β) =nH̃(k)p,n (α, β) =第 1 号 2017

1−

2(k+1)
p

Up,n (α, β),1(nhn )

2k1− p−1

Up,n (α, β).

は n → ∞，nhpn → 0 のとき，すべての M > 0 に対し，
p ≥ 3 の場合，Uchida and Yoshida（2014）
1(0)
sup Eθ∗ [|n p (α̃p,n− α∗ )|M ] < ∞,
n

1

(1)
sup Eθ∗ [|(nhn ) p−1 (β̃p,n− β ∗ )|M ] < ∞nが成り立つことを示した．α̃n(0) ，β̃n(1) の収束率は疑似対数尤度関数 Un(0) ，Up,n の正規化に依存することに注意する．Uchida and Yoshida（2012, 2014）は適応的最尤型推定量と適応的ベイズ型推定量が漸近正規性およびモーメント収束性を有することを示した．最尤型推定量を求める
ためには最適化が必要であり，ベイズ型推定量の計算には多大な時間を要することはよく知られている．
次にワンステップ推定量を考えてみる．すべての M > 0 に対し n → ∞，nhpn → 0 のとき
√
√
(0)
(0)
(2.2)
sup Eθ∗ [|( n(α̌p,n− α∗ ), nhn (β̌p,n− β ∗ ))|M ] < ∞
n∈N

(0)
(0)
(0)
(1)
(1)
(1)
を満たす初期推定量 θ̌p,n= (α̌p,n , β̌p,n ) を用いて，ワンステップ推定量 θ̌p,n = (α̌p,n , β̌p,n ) は(1)(0)(0) −1(0)= θ̌p,n− [∂θ2 Up,n (θ̌p,n)] ∂θ Up,n (θ̌p,n)θ̌p,nで定義される．ある正則条件の下で，n → ∞，nhpn → 0 のとき
√
√
(1)
(1)
Eθ∗ [f ( n(α̌p,n− α∗ ), nhn (β̌p,n− β ∗ ))] → E[f (ζ1 , ζ2 )]がすべての多項式増大な連続関数 f に対して成り立つ．ワンステップ推定量についての詳細はLehmann（1999）などを参照．しかし確率微分方程式モデルの場合，モーメント条件（2.2）を満(0)
たす初期推定量 θ̌p,n の導出は容易ではない．
以上の理由から，本節では，nhpn → 0 という一般的な仮定の下で，初期推定量の条件を緩めたハイブリッドマルチステップ推定量を提案し，その推定量が漸近正規性とモーメント収束性をもつことを示す．
これ以降に用いる記号を定義しておく．C↑k,l (Rd × Θ; Rd ) は次の条件を満たす関数 f の空間とする：(i) f (x, θ) は Rd × Θ 上の Rd - 値関数である．(ii) f (x, θ) は x に関して k 回連続微分可能
で，それらの導関数は θ について一様に，x に関する多項式増大である．(iii) |n| = 0, 1, . . . , k に対し，∂ n f (x, θ) は θ に関して l 回連続微分可能である．さらに，|ν| = 1, . . . , l，|n| = 0, 1, . . . , k に対し δ ν ∂ n f (x, θ) は θ について一様に，x に関する多項式増大である．ここで，n = (n1 , . . . , nd )，ν = (ν1 , . . . , νm ) はマルチインデックスであり，m = dim(Θ), |n| = n1 +. . .+nd , |ν| = ν1 +. . .+νm ,nm, δθi = ∂/∂θi である．F↑ (Rd ) は，f (x) が x に関
∂ n = ∂1n1 · · · ∂d d , ∂i = ∂/∂xi , δ ν = δθν11 · · · δθνmd
して多項式増大な R 上の R - 値関数となるような可測関数の空間とする．
次を仮定する．
[A1] (i) ある定数 K > 0 が存在して，すべての x, y ∈ Rd に対し，sup ||a(x, α) − a(y, α)|| + sup |b(x, β) − b(y, β)| ≤ K|x − y|.

α∈Θ1(ii) inf x,α det(A(x, α)) > 0.

β∈Θ2

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

45(iii) Xt の不変測度 μθ∗ がただ一つ存在し， Rd |f (x)|μθ∗ (dx) < ∞ を満たす任意の f ∈ F↑ (Rd )に対し，T → ∞ のとき，


1 T
p
f (Xt )dt →f (x)μθ∗ (dx)T 0Rdが成り立つ．
(iv) すべての M > 0 に対して，supt E[|Xt |M ] < ∞.
(v) Rd g(x)μθ∗ (dx) = 0 を満たす任意の g ∈ F↑ (Rd ) に対し，ある関数 G(x)，∂xi G(x) ∈ F↑ (Rd )(i = 1, . . . , d) が存在し，すべての x に対し，
Lθ∗ G(x) = −g(x).
∗(vi) Γ(θ ) は正則．
[A2] (k, l)a ∈ C↑k,4 (Rd × Θ1 ; Rd ). b ∈ C↑l,4 (Rd × Θ2 ; Rd ⊗ Rr ).

注 1. [A1] の十分条件については，Uchida and Yoshida（2012）を参照．
2.2 ハイブリッドマルチステップ推定量
初期推定量に対して次を仮定する．
[B] p ≥ 2，q = max{p, 4} とする．
(0)(i) α の初期推定量 α̂n は n → ∞, nhpn → 0 のとき，すべての M1 > 0 に対し
1(0)
n q (α̂n− α∗ ) →p 0,

1(0)
sup Eθ∗ [|n q (α̂n− α∗ )|M1 ] < ∞
n

を満たす．
(0)
(ii) β の初期推定量 β̂n は n → ∞，nhpn → 0 のとき，すべての M2 > 0 に対し
1(nhn ) q−1 (β̂n(0) − β ∗ ) →p 0,

1

sup Eθ∗ [|(nhn ) q−1 (β̂n(0) − β ∗ )|M2 ] < ∞nを満たす．
[B] を満たす初期推定量は次のように得ることができる．

det(A(x, α))1Y(α) = −tr[A(x, α)−1 A(x, α∗ ) − Id ] + log
2 Rd
det(A(x, α∗ ))μθ∗ (dx)とおき，以下を仮定する．
[A3] すべての α ∈ Θ1 に対し，Y(α) ≤ −χ|α − α∗ |2 となる正定数 χ が存在する．

初期最尤型推定量 α̌n(0) は Un(0) (α̌n(0) ) = supα∈Θ1 Un(0) (α) で定義される．ここでUn(0) (α) = −1  −1 −1{hn Ai−1 (α)[(ΔXi )⊗2 ] + log det(Ai−1 (α))}.
2 i=1n

初期ベイズ型推定量 α̃n(0) は
(0)
α̃n

=

Θ1

Θ1

で定義される．

1

α exp
n

1− 2
q

1

exp
n

1− 2
q

(0)

Un (α) π1 (α)dα(0)

Un (α) π1 (α)dα

統計数理

46

第 65 巻

第 1 号 2017

命題 1.（Uchida and Yoshida, 2012, 2014）p ≥ 2, q = max{p, 4} とする．[A1], [A2] (2, 2), [A3]を仮定する．n → ∞, nhpn → 0 のとき，すべての M > 0 に対して，次が成り立つ．
1(0)n q (α̌n− α∗ ) →p 0,

n

1
q

(0)
(α̃n
− α∗ ) →p 0,

さらに，Ỹ(β) = −

1
2


Rd

1

(0)
sup Eθ∗ [|n q (α̌n− α∗ )|M ] < ∞,
n

1

(0)
sup Eθ∗ [|n q (α̃n− α∗ )|M ] < ∞.
nA(x, α∗ )−1 [(b(x, β) − b(x, β ∗ ))⊗2 ]μθ∗ (dx)とし，もう一つの仮定をおく．
[A4] すべての β ∈ Θ2 に対し，Ỹ(β) ≤ −χ̃|β − β ∗ |2 となる正定数 χ̃ が存在する．

注 2. [A3] と [A4] の十分条件については，Fujii and Uchida（2014）の Remark 2 を参照．
初期最尤型推定量 β̌n(0) は Un(1) (α̌n(0) , β̌n(0) ) = supβ∈Θ2 Un(1) (α̌n(0) , β) で定義される．ここでUn(1) (α, β) = −1  −1 −1{hn Ai−1 (α)[(ΔXi − hn bi−1 (β))⊗2 ]}.
2 i=1n

初期ベイズ型推定量 β̃n(0) はβ̃n(0)

Θ2

1

β exp
(nhn )

=
Θ2

1−

1

exp
(nhn )

1−

2
q−1

2
q−1

(1)

(0)

(1)

(0)

Un (α̃n , β) π2 (β)dβUn (α̃n , β) π2 (β)dβで定義される．
命題 2.（Uchida and Yoshida, 2012, 2014）p ≥ 2, q = max{p, 4} とする．[A1], [A2] (2, 3), [A3],[A4] を仮定する．n → ∞, nhpn → 0 のとき，すべての M > 0 に対して，
1(nhn ) q−1 (β̌n(0) − β ∗ ) →p 0,
1(nhn ) q−1 (β̃n(0) − β ∗ ) →p 0,

1

sup Eθ∗ [|(nhn ) q−1 (β̌n(0) − β ∗ )|M ] < ∞,n1

sup Eθ∗ [|(nhn ) q−1 (β̃n(0) − β ∗ )|M ] < ∞nが成り立つ．
次に，マルチステップ推定量について考える．p ≥ 2, q = max{p, 4}, k0 = [ p2 ], l0 = [ q−1],2]とする．疑似対数尤度関数
U
(θ)
は，
p
=
2
のとき，
m0 = [ q−2p,n21  −1 −1{hn Ai−1 (α)[(ΔXi − hn bi−1 (β))⊗2 ] + log det(Ai−1 (α))}
2 i=1n

U2,n (θ) = −で定義される．p ≥ 3 に対しては，



k
k0
n
0


1
(j)
(k0 )
(j)
−1
j
⊗2
j
Up,n (θ) = −hn Di−1 (θ) [(Xtni − ri−1 (hn , θ)) ] +hn Ei−1 (θ)
hn
2 i=1
j=0
j=0
となる．
Jp,n (α, β) :=1 2∂α Up,n (α, β) は正則 ,n高頻度データに基づく確率微分方程式モデルのハイブリッド推定

Kp,n (α, β) :=47

1 2∂β Up,n (α, β) は正則 ,nhn

1 2
c (α,β) ,∂α Up,n (α, β)1Jp,n (α,β) + Em1 1Jp,nn
1 2c (α,β) ,Ξp,n (α, β) :=∂β Up,n (α, β)1Kp,n (α,β) + Em2 1Kp,nnhn
Γp,n (α, β) :=とおく．ここで，Em は m 次単位行列，1K (ω) は ω ∈ K のとき 1K (ω) = 1，ω ∈ K c のとき
(0)
(0)
(0)
(0)
1K (ω) = 0 であるとする．[B] において (α̂I,n , β̂I,n ) = (α̂n , β̂n ) とする．タイプ I のマルチス
(l0 )
(m )
, β̂I,n0 は，k = 1, . . . , l0 に対し，テップ推定量 α̂I,n
1
(k−1)
(k−1)
) ∂α Up,n (α̂I,n , β̂I,n ),
n
1
(k)
(k−1)
(k)
(k−1)
(k)
(k−1)
β̂I,n = β̂I,n − Ξ−1)∂β Up,n (α̂I,n , β̂I,n )p,n (α̂I,n , β̂I,n
nhn
(k)

(k−1)

α̂I,n = α̂I,n

(k−1)

− Γ−1
p,n (α̂I,n(k−1), β̂I,nで定義される．
(0)
(0)
(0)
(0)
, β̂II,n ) = (α̂n , β̂n )次に，タイプ II のマルチステップ推定量を考える．[B] において (α̂II,n(0)(1)とする．U0,n (θ) = Un (α), U1,n (θ) = Un (θ) とする．p = 2 (l0 = m0 = 1) のとき，
1
(1)
(0)
(0)
(0)
(0)
(0)
α̂II,n = α̂II,n − Γ−10,n (α̂II,n , β̂II,n ) ∂α U0,n (α̂II,n , β̂II,n ),
n
1
(1)
(0)
(1)
(0)
(1)
(0)
∂β U1,n (α̂II,n , β̂II,n )β̂II,n = β̂II,n − Ξ−1
1,n (α̂II,n , β̂II,n )nhnと定義する．p = 3 (l0 = m0 = 1) のとき，
(1)

(0)

(0)

(0)

β̂II,n = β̂II,n − Ξ−1
1,n (α̂II,n , β̂II,n )

1
(0)
(0)
∂β U1,n (α̂II,n , β̂II,n ),
nhn

1
(1)
(0)
(0)
(1)
(0)
(1)
α̂II,n = α̂II,n − Γ−13,n (α̂II,n , β̂II,n ) ∂α U3,n (α̂II,n , β̂II,n )nと定義する．
(l0 )(m )p ≥ 4 のとき，タイプ II のマルチステップ推定量 α̂II,n, β̂II,n0 は，k = 1, . . . , l0 に対し，
(k)
(k−1)
(k−1)
(k−1) 1(k−1)(k−1)α̂II,n = α̂II,n − Γ−1
2k+1,n (α̂II,n , β̂II,n ) ∂α U2k+1,n (α̂II,n , β̂II,n ),
n
1
(k)
(k−1)
(k)
(k−1)
(k)
(k−1)
β̂II,n = β̂II,n − Ξ−1∂β U2k+2,n (α̂II,n , β̂II,n )2k+2,n (α̂II,n , β̂II,n )nhnで定義される．ここで，U2k+2,n は p = 2k + 2 のときの Up,n と等しいことに注意する．例えば，
p = 4 (l0 = m0 = 1) のとき，
1
(1)
(0)
(0)
(0)
(0)
(0)
α̂II,n = α̂II,n − Γ−13,n (α̂II,n , β̂II,n ) ∂α U3,n (α̂II,n , β̂II,n ),
n
1
(1)
(0)
(1)
(0)
(1)
(0)
β̂II,n = β̂II,n − Ξ−1∂β U4,n (α̂II,n , β̂II,n ).
4,n (α̂II,n , β̂II,n )nhn

さらにタイプ III のマルチステップ推定量について考察する．疑似対数尤度関数は，k = 1, . . . , l0
に対し，




k+1
n
 j (j)1
−1 −1⊗2
V2k+1,n (α | θ̄) = −hn D̄i−1 (θ̄) + log det Ai−1 (α) ,hn Ai−1 (α) (Xtni − Xtni−1 ) −
2 i=1
j=2
⎡
⎤
⊗2
k+1
n
 j (j)1  −1 −1⎦
hn Ai−1 (ᾱ) ⎣ Xtni − Xtni−1 − hn bi−1 (β) −hn r̄i−1 (θ̄)
V2k+2,n (β | θ̄) = −
2 i=1
j=2

統計数理

48

第 65 巻

第 1 号 2017

で定義される．ここで，θ̄ = (ᾱ, β̄) である．さらに，l, m = 1, . . . , d に対し，fl (x) = xl ,
(j)
(j)
1 j
1 j
hlm (x) = (x − Xtni−1 )l (x − Xtni−1 )m , D̄i−1 (θ̄)lm = j!
Lθ̄ hlm (Xtni−1 ), r̄i−1 (θ̄)l = j!
Lθ̄ fl (Xtni−1 )である．V2k+1,n (α | θ̄) は V2k+1,n (α | α, β) でないことに注意する．例えば，∂α V2k+1,n (α | θ̄) =∂α V2k+1,n (α | α, β) である．
k = 1, 2, . . . , l0 に対し，θ = (α, β) とし，
1 2J2k+1,n,III (α, β) :=は正則 ,∂α V2k+1,n (α | θ̄)
n
θ̄=θ

1 2

∂β V2k+2,n (β | θ̄)は正則 ,K2k+2,n,III (α, β) :=
nhn
θ̄=θ

1

c
1J2k+1,n,III (α,β) + Em1 1J2k+1,n,III
Γ2k+1,n,III (α, β) := ∂α2 V2k+1,n (α | θ̄)(α,β) ,
n
θ̄=θ

1 2

c
Ξ2k+2,n,III (α, β) :=∂β V2k+2,n (β | θ̄)1K2k+2,n,III (α,β) + Em2 1K2k+2,n,III
(α,β)
nhn
θ̄=θ
(0)
(0)
(0)
(0)
(l0 )
(m0 )
とする．[B] において (α̂III,n, β̂III,n ) = (α̂n , β̂n ) とおく．p = 2, 3 のとき，(α̂III,n, β̂III,n
)=
(l0 )
(m0 )
(l0 )
(m0 )
(α̂II,n , β̂II,n ) とおく．p ≥ 4 のとき，タイプ III のマルチステップ推定量 α̂III,n , β̂III,n は，
k = 1, . . . , l0 に対して
(k)
(k−1)
(k−1)
(k−1) 1
(k−1)
(k−1)
(k−1)
α̂III,n = α̂III,n − Γ−1
2k+1,n,III (α̂III,n , β̂III,n ) ∂α V2k+1,n (α̂III,n | α̂III,n , β̂III,n ),
n
1
(k)
(k−1)
(k)
(k−1)
(k−1)
(k)
(k−1)
β̂III,n = β̂III,n − Ξ−1∂β V2k+2,n (β̂III,n | α̂III,n , β̂III,n )2k+2,n,III (α̂III,n , β̂III,n )nhnで定義される．
補題 1. （Kamatani and Uchida, 2015）k ∈ N , p ≥ 2k + 2, j = I, II, III とする．[A1],[A2] (2k0 , 2k0 + 1), [B] を仮定する．n → ∞, nhp → 0 のとき，すべての M > 0 に対して，
k

(k−1)

− α∗ )|M ] < ∞,

k

(k−1)

− β ∗ )|M ] < ∞.

sup Eθ∗ [|n p (α̂j,nn

sup Eθ∗ [|(nhn ) p−1 (β̂j,nn

定理 1.（Kamatani and Uchida, 2015）p ≥ 2, q = max{p, 4}, k0 = [ p2 ], l0 = [ q−1], m0 = [ q−2],22j = I, II, III とする．[A1], [A2] (2k0 , 2k0 + 1), [B] を仮定する．このとき，n → ∞, nhp → 0
の下，
√
√
d
(l0 )
(m )
( n(α̂j,n− α∗ ), nhn (β̂j,n0 − β ∗ )) → (ζ1 , ζ2 ) ∼ Nm1 +m2 (0, diag[Γ1 (θ∗ )−1 , Γ2 (θ∗ )−1 ])が成り立つ．さらに，すべての多項式増大な連続関数 f に対して
√
√
(l0 )
(m )
Eθ∗ [f ( n(α̂j,n− α∗ ), nhn (β̂j,n0 − β ∗ ))] → E[f (ζ1 , ζ2 )]が成り立つ．
3 種類のマルチステップ推定量は，Uchida and Yoshida（2012）における 3 種類の適応的最尤型推定量にそれぞれ対応している．以下，p ≥ 4 の場合について説明する．タイプ I のマルチステップ推定量は，最適な疑似尤度関数 Up,n に対して，必要回数だけニュートン・ラフソン
法を繰り返して導出される．タイプ II のマルチステップ推定量は，疑似尤度関数の精度を徐々に上げながら (U3,n , U4,n , . . ., U2k+1,n , U2k+2,n ), ニュートン・ラフソン法によって求める．タ
イプ III のマルチステップ推定量は，タイプ II のマルチステップ推定量で用いた疑似尤度関

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

49

数よりもシンプルな疑似尤度関数 (V3,n , V4,n , . . ., V2k+1,n , V2k+2,n ) を用いて導出される．理論的観点から言えば 3 種類のマルチステップ推定量は同じ漸近的性質を有するが，計算時間はタイプ I > タイプ II > タイプ III となっており，タイプ III のマルチステップ推定量の計算時間が一番短い．一方，タイプ I やタイプ II のマルチステップ推定量はタイプ III のマルチステッ
プ推定量よりも精密な疑似尤度関数を使用しているので，離散観測（サンプルサイズが有限）におけるタイプ I やタイプ II のマルチステップ推定量のパフォーマンスはタイプ III のマルチステップ推定量と同等もしくはそれよりも良いことが期待される．
2.3 例とシミュレーション結果
1 次元拡散過程(2.3)

β3
β1 − β 2 X t − 
1 + β2 Xt2

dXt =

dt +

α2 + Xt2
dwt ,1 + α1 Xt2

X0 = 2,を考える．ここで，θ = (α1 , α2 , β1 , β2 , β3 ) は未知パラメータ，パラメータの真値は θ∗ =(α1∗ , α2∗ , β1∗ , β2∗ , β3∗ ) = (0.3, 0.5, 3.0, 0.5, 2.0)，パラメータ空間は [0.1, 50]5 である．
p = 4 すなわち l0 = 1, m0 = 1 とする．バランス条件は nh4n → 0 となることに注意する．
(1)
(1)
(1)
(1)
(1)
(1)
(1)
θ̂B,III = (α̂III , β̂III ) := (α̂III,1 , α̂III,2 , β̂III,1 , β̂III,2 , β̂III,3 ) は，ベイズ型推定量を初期推定量としたタイプ III のマルチステップ推定量で，次のように表される．
(0)
α̂III

=

Θ1

n

=

Θ2

n

(0)

Un (α) π1 (α)dα

p

1
n

1−

2
p−1

1

exp
n

（初期ベイズ型推定量）,(0)

Un (α) π1 (α)dα
1− 2

β exp

Θ2

1− 2
p

1

exp
Θ1
(0)
β̂III

1

α exp

1−

2
p−1

(1)

(0)

(1)

(0)

Un (α̂III , β) π2 (β)dβ（初期ベイズ型推定量）,

Un (α̂III , β) π2 (β)dβ

(1)
α̂III

(0)
(0)
(0) 1
(0)
(0)
(0)
= α̂III − Γ−1
3,n,III (α̂III , β̂III ) ∂α V3,n (α̂III | α̂III , β̂III ),

(1)
β̂III

(0)
(1)
(0)
= β̂III − Ξ−1
4,n,III (α̂III , β̂III )

(1)

(1)

(1)

(1)

n
1
(0)
(1)
(0)
∂β V4,n (β̂III | α̂III , β̂III ).
nhn

(1)

(1)

(1)

θ̂M,III = (α̃III , β̃III ) := (α̃III,1 , α̃III,2 , β̃III,1 , β̃III,2 , β̃III,3 ) は，最尤型推定量を初期推定量としたタイプ III のマルチステップ推定量で，次のように表される．
(0)

Un(0) (α̃III ) = sup Un(0) (α) （初期最尤型推定量）
,
α∈Θ1

(0)
(0)
Un(1) (α̃III , β̃III ) =(0)

sup Un(1) (α̃III , β) （初期最尤型推定量）
,

β∈Θ2

(1)
(0)
(0)
(0) 1
(0)
(0)
(0)
α̃III = α̃III − Γ−1
3,n,III (α̃III , β̃III ) ∂α V3,n (α̃III | α̃III , β̃III ),
n
1
(1)
(0)
(1)
(0)
(0)
(1)
(0)
β̃III = β̃III − Ξ−1∂β V4,n (β̃III | α̃III , β̃III ).
4,n,III (α̃III , β̃III )
nhn
(1)

(1)

(1)

(1)

(1)

(1)

(1)

θ̂A,III = (α̌III , β̌III ) := (α̌III,1 , α̌III,2 , β̌III,1 , β̌III,2 , β̌III,3 ) は，タイプ III の適応的最尤型推定量で，次のように表される．
(0)

Un(0) (α̌III ) = sup Un(0) (α) （初期最尤型推定量）
,
α∈Θ1

統計数理

50

第 65 巻

第 1 号 2017

表 1．推定量の平均（上段）と標準偏差（下段）T = 250, h = 1/390, p = 4．

(0)

(0)

(0)

Un(1) (α̌III , β̌III ) = sup Un(1) (α̌III , β) （初期最尤型推定量）
,
β∈Θ2

(1)

(0)

(0)

(0)

(0)

(1)

(0)

V3,n (α̌III | α̌III , β̌III ) = sup V3,n (α | α̌III , β̌III ),
α∈Θ1

(1)

(1)

(0)

V4,n (β̌III | α̌III , β̌III ) = sup V4,n (β | α̌III , β̌III ).
β∈Θ2

(1)
初期推定量をベイズ型推定量としたタイプ III のマルチステップ推定量 θ̂B,III，初期推定量(1)を最尤型推定量としたタイプ III のマルチステップ推定量 θ̂M,III , そして，タイプ III の適応的
(1)の漸近挙動を検証する．シミュレーションでは，T = 250，hn = 1/390
最尤型推定量 θ̂A,III（1 年の取引日数が 250 日で 1 日の取引時間が 390 分）と設定し，真のモデルから 1000 本の独立なサンプルパスを発生させた．実行には R（R Development Core Team 2013）のパッケージ
Yuima（Brouste et al., 2014）と組込み関数 optim() を用い，初期ベイズ型推定量の計算のために
Kamatani（2014）
が提案したマルコフ連鎖モンテカルロ（MCMC）法のアルゴリズムを用いた．
3
MCMC 法におけるバーンイン回数 Bi は 10 とし，マルコフチェーンの生成数 M を 104 とした．推定値の平均と標準偏差を表 1 に示す．
初期最尤型推定量とタイプ III の適応的最尤型推定量は optim() を用いて計算した．その際に必要となる初期値は，真値に近いものを選び，θ̄ = (ᾱ1 , ᾱ2 , β̄1 , β̄2 , β̄3 ) = (0.5, 1.0, 4.0, 1.0, 2.0) とした．
3 つの推定量はすべて，理論上は定理 1 と同等の漸近的性質をもつが，数値計算上では違いが生じる．まず，θ̂B,III および θ̂M,III の計算は，関数 V3,n および V4,n の最適化を必要としない．これは数値計算する上では魅力的である．ただし，初期ベイズ推定量は MCMC 法により計算されるために θ̂B,III の導出には時間を要する．一方，θ̂M,III と θ̂A,III の計算に要した平均
CPU 時間は 9.772 秒と 14.802 秒にすぎない（PC Intel 2.8 GHz を使用）
．これはニュートン・ラ
フソン法および最適化 optim() が高速であることを意味する．初期ベイズ型推定量は正規化項
1− 2
1− 2
n p , (nhn ) p−1 により頑健である．正規化項はモンテカルロの用語では温度と言われており，計算速度を上げるために使われる．（Robert and Casella, 2004 を参照）
．これは，純粋な理論的結果とモンテカルロシミュレーション技法の興味深い関係を表していると言える．
表 1 の結果が示すように，3 種類の推定量には著しい違いは見られないが，これは、最適化のための初期値を真値に近いものにしたためである．θ̂M,III と θ̂A,III に使われる初期最尤型推定量の導出には，最適化が必要になるのでそれに伴い初期値の選択が要求される．一方，ベイズ型推定量の計算は初期値の影響を受けづらいので，最適化が困難な場合には θ̂B,III の方が好ましいと言える．初期値の効果を示すため，次の拡散過程を考える．

(2.4)

dXt = (β1 − β2 Xt − 2 sin(β3 Xt ))dt +α2 + Xt2
dwt ,1 + α1 Xt2

X0 = 2,高頻度データに基づく確率微分方程式モデルのハイブリッド推定

51

図 1．β1 , β2 の初期ベイズ型推定量（左）と初期最尤型推定量（右）．異なる初期値を用い，最適化には optim() を使用した．右図には 5 つのクラスターがあり，下から上に向かってそれぞれ 36, 1, 32, 11, 20 の点が集まっている．真値は (β1 , β2 ) = (3.0, 7.0) で，下の方にあるクラスターは真値に近いところに出来ている．描画の重なりを避けるために，
ggplot2 の関数 position jitter(w=0.02,h=0.02) を用いて実際の値の周囲に点を分散させている．

ここで，θ = (α1 , α2 , β1 , β2 , β3 ) の真値は，θ∗ = (α1∗ , α2∗ , β1∗ , β2∗ , β3∗ ) = (0.3, 0.5, 3.0, 7.0, 5.0) である．T, h, p およびパラメータ空間は前述のモデルと同じとする．図 1 は一回の実験で得られた
β（，β（の推定値の散布図である．左図は初期ベイズ型推定量，右図は初期最尤型
1 横軸）
2 縦軸）
推定量である．いずれもパラメータ空間から一様に選択された 100 個の異なる初期値から計算されている．MCMC 法においてはバーンイン 104 回，反復 105 回とした．右図の左上隅の 4 つのクラスターには 64% の点が存在しており局所最大になっている．左図では真値の近傍でただ一つのクラスターが存在するのみである．
（図 2）．
この初期ベイズ型推定量と初期最尤型推定量を用いて，θ̂B,III , θ̂M,III , θ̂A,III を計算するは約 50% の点が真値の近くに集まっている．θ̂A,III は θ̂M,III（右上）よりも安定しθ̂A,III（左下）は最も安定していることが分かる．
ており，θ̂B,III（左上）3.

非エルゴード的拡散過程

本節は，非エルゴード的拡散過程を含む一般の確率回帰モデルのボラティリティパラメータのハイブリッドマルチステップ推定量の構成およびその漸近的性質について解説する．詳細は，Kamatani et al.（2016）
を参照．
3.1 モデルと仮定
次の確率積分方程式で表される確率回帰モデルのボラティリティのパラメータ推定を考える．
 t
 t
(3.1)
Yt = Y0 +bs ds +σ(Xs , θ)dws , t ∈ [0, T ].
00

ここで w は確率空間（Ω, F , (Ft )t∈[0,T ] , P ）上の r 次元標準ウィーナー過程，b, X はそれぞれ Rm ,Rd 上の値をとる発展的可測過程，Y0 は Rm - 値の初期状態，σ は Rd × Θ 上で定義される Rm ⊗ Rr
値関数，θ∗ を θ の真値とする．Θ は局所リプシッツ境界をもつ Rp での有界領域とする．F - 安
定分布収束は →ds (F ) と表記する．F - 安定分布収束については，Jacod and Shiryaev（2003）の

52

統計数理

第 65 巻

第 1 号 2017

図 2．β1 , β2 の 3 つの推定量，θ̂B,III (左上), θ̂M,III (右上), θ̂A,III (下) の散布図．図 1 に示した初期ベイズ型推定量と初期最尤型推定量を用いている．右上図には下から上に向かってそれぞれ 36, 32, 32 の点が集まった 3 つのクラスターがあり，下図には 12, 36, 32, 20
個の点が集まった 4 つのクラスターがある．実際の値の周囲に点を描画させている．

VIII 章を参照．S(x, θ) = σ(x, θ)⊗2 , Δk Y = Ytk − Ytk−1 とおく．σ は Rd × Θ̄ 上の連続関数に拡張可能とし，それを σ と表記する．f ∈ Lp (P ) と p > 1 に対し，||f ||p = (E[|f |p ])1/p とする．
データ Zn = (Xtk , Ytk )0≤k≤n , tk = kh, h = hn = T /n は離散観測される．b は未知であることに注意する．極限は n → ∞ を考える．すなわち Zn は高頻度データである．
離散観測に基づく非エルゴード的拡散型確率過程モデルによる統計推測は多くの研究者に，Florens-Zmirou（1989），Genon-Catalot andよって発展してきている．例えば，Dohnal（1987），Gobet（2001）
を参照．Uchida and Yoshida（2013）は確率回帰モデルの最Jacod（1993, 1994）
尤型推定量もベイズ型推定量も漸近混合正規性とモーメント収束をもつことを示した．しかし，最尤型推定量の導出は数値最適化を必要とし，適切な初期値の選択が重要であり，ベイズ型推定量の計算には多大な時間を要する．先述した通り，ワンステップ推定量は非常に有効である√
が， n - 一致性をもつ初期推定量を見つけることは容易ではないため，非エルゴード的拡散過程のワンステップ推定量の導出は困難になる．
本節では，先述した Kamatani and Uchida（2015）の手法に基づいて，初期ベイズ型推定量を用いた確率回帰モデルのハイブリッドマルチステップ推定量を提案し，マルチステップ推定量
が漸近混合正規性とモーメント収束性をもつことを示す．
次を仮定する．

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

53[C1] (i) すべての p > 1 に対して，sup0≤t≤T bt p < ∞．
(ii) inf x,θ det S(x, θ) > 0, σ ∈ C↑2,4 (Rd × Θ; Rm ⊗ Rr ).
[C2] X は次のように表される:

 t

Xt = X0 +0 tb̃s ds +0

 t
as dws +0ãs dw̃s .

ここで，b̃, a, ã は，それぞれ Rd , Rd ⊗ Rr , Rd ⊗ Rr1 - 値の発展的可測過程で，すべての p > 1 に対して
X0 p + sup ( b̃t p + at p + ãt p ) < ∞t∈[0,T ]を満たす．w̃ は w と独立な r1 次元ウィーナー過程である．
疑似対数尤度関数 Hn (θ) は1
{log det S(Xtk−1 , θ) + h−1 S −1 (Xtk−1 , θ)[(Δk Y )⊗2 ]}2nHn (θ) = −k=1で与えられる．Yn (θ) = n1 {Hn (θ) − Hn (θ∗ )} とおくと，仮定 [C1] [C2] の下で，θ ∈ Θ について一様に，Y(θ) = −

1
2T



 T
log
0

det S(Xt , θ)
det S(Xt , θ∗ )へ確率収束する．
χ0 = inf∗θ=θ+ tr(S −1 (Xt , θ)S(Xt , θ∗ ) − Id ) dt

−Y(θ)|θ − θ∗ |2

とおく．次はインデックス χ0 の非退化性に関する条件である．
[C3] すべての L > 0 に対して，cL > 0 が存在して，すべての r > 0 に対し，P [χ0 ≤ r−1 ] ≤ rL
Lが成り立つ．
c[C3] は 1/χ0 がすべてのオーダーで有限なモーメントをもつことと同等であることに注意する．
[C3] の十分条件については Uchida and Yoshida（2013）
を参照．
3.2 ハイブリッドマルチステップ推定量
初期推定量についての仮定をおく．
(0)[D] q ∈ (0, 1/2] とする．θ の初期推定量 θ̂n は n → ∞ のとき，すべての M1 > 0 に対し，
sup Eθ∗ [|nq (θ̂n(0) − θ∗ )|M1 ] < ∞nを満たす．
[D] を満たす初期推定量は次のように求められる．q ∈ (0, 1/2] とする．事前分布 π : Θ → R+
(0)
に対する初期ベイズ型推定量 θ̃q,nは
 1
θ exp n1−2q
Hn (θ) π(θ)dθ
(0)
Θ
 1

θ̃q,n =
exp n1−2q
Hn (θ) π(θ)dθΘで定義される．π は連続で 0 < inf θ∈Θ π(θ) ≤ supθ∈Θ π(θ) < ∞ であるとする．

統計数理

54

第 65 巻

第 1 号 2017

q = 0.5 の場合が通常のベイズ型推定量である．q < 0.5 のときのベイズ型推定量を初期推定量としたマルチステップ推定量の方が通常のベイズ型推定量や通常のベイズ型推定量を初期推定量としたマルチステップ推定量よりもパフォーマンスが良くなる数値シミュレーションの結果を 3.3 節で述べる．
Uq,n = {u ∈ Rp ; θ∗ + n1q u ∈ Θ}, Vq,n (r) = {u ∈ Uq,n ; r ≤ |u|} とおく．Uq,n 上の統計的確
率場 Zq,n を


1
1
1
∗
(3.2)
Zq,n (u) = exp
H+u− 1−2q Hn (θ∗ )
θ
n
n1−2q
nq
n

と定義する．
命題 3.（Uchida and Yoshida, 2013）q ∈ (0, 1/2] とする．[C1], [C2], [C3] を仮定する．このとき，すべての L > 0 に対して，正定数 CL が存在して，すべての r > 0, n ∈ N について
CLP
sup Zq,n (u) ≥ e−r ≤ L
ru∈Vq,n (r)が成り立つ．
命題 4. （Kamatani et al., 2016）q ∈ (0, 1/2] とおく．[C1], [C2], [C3] を仮定する．このとき，
n → ∞ の下で，すべての M > 0 に対して，(0)
sup Eθ∗ [|nq (θ̃q,n− θ∗ )|M ] < ∞.
n

マルチステップ推定量について考察する．
1
Γn (θ) := ∂θ2 Hn (θ),
n
Kn (θ) := {Γn (θ) は正則},Γ̄n (θ) := Γn (θ)1Kn (θ) + Ep 1Knc (θ) ,とおく．ここで Ep は p 次単位行列である．
(J)q ∈ (0, 1/2] に対し，J = [− log2 q] とおく．マルチステップ推定量 θ̂n は，k = 1, . . . , J に対し(k−1) 1) ∂θ Hn (θ̂n(k−1) )θ̂n(k) = θ̂n(k−1) − Γ̄−1n (θ̂nnで定義される．
補題 2. （Kamatani et al., 2016）q ∈ (0, 1/2], J = [− log2 q] とおく．[C1], [C2], [C3], [D] を仮定する．このとき，k = 0, 1, . . . , J − 1 に対し，すべての M > 0 についてk

sup Eθ∗ [|n2 q (θ̂n(k) − θ∗ )|M ] < ∞nが成り立つ．
Γ(θ∗ ) = (Γij (θ∗ ))i,j=1,...,p , ただし
Γij (θ∗ ) =

1
2T

 T
0

tr((∂θi S)S −1 (∂θj S)S −1 (Xt , θ∗ ))dtとし，ζ は Γ(θ∗ ) とは独立な p 次元標準正規確率変数とする．
定理 2. （Kamatani et al., 2016）q ∈ (0, 1/2], J = [− log2 q] とおく．[C1], [C2], [C3], [D] を仮高頻度データに基づく確率微分方程式モデルのハイブリッド推定

55

定する．このとき，n → ∞ の下，√n(θ̂n(J) − θ∗ ) →ds (F ) Γ(θ∗ )−1/2 ζが成り立つ．さらに，高々多項式増大のすべての連続関数 f に対して，√
E[f ( n(θ̂n(J) − θ∗ ))] → E[f (Γ(θ∗ )−1/2 ζ)]が成り立つ．
3.3 例とシミュレーション結果
次式で定義される 1 次元拡散過程を考える．
dXt = −(Xt − 1)dt + [θ1 + θ2 {1 + sin(θ3 Xt )}]dwt ,X0 = 1.

t ∈ [0, 1],

ここで，真値は θ∗ = (1, 4, 8) とし，パラメータ空間は Θ = [0.01, 20] × [0, 20] × [0, 20] であるとする．データ (Xti )i=0,1,...,n は ti = ih, h = 1/104 , tn = nh = T = 1, サンプル数 n は 104 とする．ここでは，最尤型推定量 θ̂M,n（Genon-Catalot and Jacod, 1993），ベイズ型推定量 θ̂B,n（Uchida and Yoshida, 2013）
そして本節で提案されたハイブリッドマルチステップ（HMS）推定
量のシミュレーションを行う．最尤型推定量 θ̂M,n はHn (θ̂M,n ) = sup Hn (θ)θ∈Θ

で定義され，一様分布を事前分布とするベイズ型推定量 θ̂B,n はθ̂B,n :=

θ exp(Hn (θ))dθexp(Hn (θ))dθ
Θ

Θ

(0)
で定義される．q ∈ (0, 1/2] とすると，一様事前分布に対する初期ベイズ型推定量 θ̃q,nは
 1
θ exp n1−2q Hn (θ) dθ
(0)
 1

= Θ
θ̃q,n
exp n1−2q
Hn (θ) dθΘで定義される．Hn (θ) の最大化には，R 言語の組込み関数 optim() の L-BFGS-B 法を用いた．ベ
イズ型推定量は Kamatani（2014）
が提案した MCMC 法のアルゴリズムを使用して計算した．
真のモデルから 1000 本の独立なサンプルパスを発生させ，推定値の平均および標準偏差を計算し，結果を表 2–4 に示す．表 2 は 2 つの異なる初期値から得られた最尤型推定量 θ̂M,n である．最尤型推定量は組込み関数 optim() を用いて導出した．真値を初期値としたものは良いパ
フォーマンスをするが，Θ 上の一様分布から得られた初期値は真値から遠く離れた値をとり得るため，最適化に失敗している．
表 3 は一様事前分布を持つベイズ型推定量 θ̂B,n の結果である．マルコフチェーン生成数 M
を 5 × 104 , 5 × 105 , 107 とし，これに対するバーンイン回数 Bi をそれぞれ 5 × 103 , 5 × 104 ,
106 とした．(M, Bi) = (107 , 106 ) のとき，ベイズ型推定量は良いふるまいをしているが，(M, Bi) = (5 × 104 , 5 × 103 )，(5 × 105 , 5 × 104 ) のときは，MCMC 法によって生成されるマ
ルコフチェーンが定常分布に収束しきれていないため，ベイズ型推定量の計算に失敗している．
(0)と
表 4 は M = 5 × 104 , Bi = 5 × 103 の場合の，一様事前分布をもつベイズ型推定量 θ̃q,n(J)
q = 0.5, 0.45, 0.4, . . . , 0.1, 0.05, J = [− log2 q] の場合の HMS 推定量 θ̂q,n のシミュレーション結果(0)− θ∗ )|M ] < ∞ となることに注意されたい．
である．すべての M > 0 に対し supn Eθ∗ [|nq (θ̃q,n
この例では，HMS 推定量の中で q = 0.2 とした場合が一番良いことがわかる．理論的に最適な

56

統計数理

第 65 巻

第 1 号 2017

表 2．最尤型推定量 n = 1 × 104 ．

表 3．ベイズ型推定量 n = 1 × 104 ．

表 4．HMS 推定量 n = 1 × 104 , M = 5 × 104 , Bi = 5 × 103 ．

q を選び出すことは困難であるが，実際には q を様々な値に変えて，競合する HMS 推定量の中(J))
から最も良い推定量を得ることができる．最も良い推定量 θ̂n∗ とは Hn (θ̂n∗ ) = maxq∈K Hn (θ̂q,nを満たすものであり，K は q の値の組
（例えば K = {0.05, 0.1, 0.15, . . . , 0.5}）である．
次に，推定量を得るために要した計算時間に焦点をあてる．シミュレーションに使用したパ（3.4 GHz base clock/3.9 GHz Turbo, 12 MB cache）である．
ソコンのスペックは，Intel i7 4930K(J)
最尤型推定量 θ̂M,n , ベイズ型推定量 θ̂B,n , HMS 推定量 θ̂0.2,nの計算時間はそれぞれ 0.9 秒，1733
秒，1667 秒であった．θ̂M,n は Θ 上の一様分布からの乱数を初期値とした．θ̂B,n は一様分布を(J)
は，q = 0.2，一様事前分布，M = 5 × 104 とした初期ベ事前分布，M = 5 × 104 とした．θ̂0.2,n高頻度データに基づく確率微分方程式モデルのハイブリッド推定

57(0)
イズ型推定量 θ̃q,nを用いて計算した．事前分布を一様分布とし，M = 5 × 105 とした場合のベイズ型推定量 θ̂B,n の平均計算時間は 281 分である．事前分布を一様分布とし，M = 5 × 104 とした場合の初期ベイズ型推定量に基づいた HMS 推定量は，計算時間に関してはほとんど違いがみられない．q を変えた結果を表 4 に示す．数値計算上の見地から，短時間で推定値を得ることは極めて重要であるが，正確な推定値を得ることの方がより重要である．この意味で，optim()を使って求める最尤型推定量と比較して，初期ベイズ型推定量を用いて HMS 推定量を得るには多くの時間を要するが，このモデルでは最尤型推定量より HMS 推定量の方が優れていると
言える．

4.

微小拡散過程本節は，Nomura and Uchida（2016a, 2016b）の結果に従って，微小拡散過程のドリフトパラ
メータおよびボラティリティパラメータの適応的最尤型推定，適応的ベイズ型推定そしてハイブリッド推定について概説する．
4.1 モデルと仮定
次の確率微分方程式で定義される d 次元拡散過程を考える．
dXt = a(Xt , α)dt + b(Xt , β)dwt , t ∈ [0, T ],(4.1)
X 0 = x0 .

 ∈ (0, 1],

ここで，a : Rd × Θα → Rd , b : Rd × Θβ → Rd ⊗ Rr , w は r 次元標準ウィーナー過程， と T は既
知な定数，x0 は確定的初期条件，θ = (α, β) ∈ Θ = Θα × Θβ ，θ∗ = (α∗ , β ∗ ) ∈ Θ は θ の真値，Θα ,Θβ はそれぞれ Rp , Rq 上の有界で凸な開部分集合である．パラメータ空間は局所的リプシッツ
境界をもつとする．本節の微小拡散過程のドリフト係数 a とボラティリティ係数 b は 2 節のエ
ルゴード的拡散過程のドリフト係数 b とボラティリティ係数 a とは記号が異なっていることに注意する．データは離散観測され，Xn = (Xti )0≤i≤n , ti = ihn , hn = T /n であるとする． → 0,√n → ∞, √1 n = O(1) の場合を考える．また，ある γ ∈ (0, 1] が存在して，( n)γ = O(1) を満たすとする．（4.1）で定義される拡散過程は，小さな摂動をもつダイナミカルシステムあるいは微小拡．微散過程と呼ばれている（Azencott, 1982; Freidlin and Wentzell, 1998; Yoshida, 1992a を参照）
，Kunitomo and小拡散過程の数理ファイナンスや生命科学への応用については Yoshida（1992b），Takahashi and Yoshida（2004），Uchida and Yoshida（2004b），Fuchs（2013）
Takahashi（2001）とそれらの参考文献を参照．B(x, β) = bb (x, β), ΔXi = Xti − Xti−1 , ai−1 (α) = a(Xti−1 , α),Bi−1 (β) = B(Xti−1 , β) とする．Xt0 を  = 0 としたときの常微分方程式の解とする．すなわち
dXt0 = a(Xt0 , α∗ )dt, X00 = x0 である．
連続観測される微小拡散過程の統計推測は著しく発展している．例えば，Kutoyants（1984,，Yoshida（1992a, 1993, 2003），Prakasa Rao（1999），Iacus（2000），Iacus and Kutoyants
1994）（2001），Uchida and Yoshida（2004a）とそれらの参考文献を参照．離散観測された微小拡散過程
については，最尤型推定量や M 推定量に関する研究は数多く挙げられる（Genon-Catalot, 1990;
Laredo, 1990; Sørensen, 2000, 2012; Sørensen and Uchida, 2003; Uchida, 2003, 2004, 2006, 2008;．しかし，適応的推測やハイブリッド型推測
Gloter and Sørensen, 2009; Guy et al., 2014 など）の研究はほとんどない．
本節では，Nomura and Uchida（2016a, 2016b）で提案された適応的最尤型推定量，適応的ベイズ型推定量そしてハイブリッド推定量の漸近的性質（漸近正規性とモーメント収束性）について説明する．

統計数理

58

第 65 巻

第 1 号 2017

以下の仮定をおく．
[E1] (i) すべての x, y ∈ Rd に対して
sup |a(x, α) − a(y, α)| + sup b(x, β) − b(y, β) ≤ K|x − y|

α∈Θα

β∈Θβ

となる K > 0 が存在する．
(ii) inf det B(x, β) > 0.
x,β[E2] a(x, α) ∈ C↑6,4 (Rd × Θα ; Rd ), b(x, β) ∈ C↑6,4 (Rd × Θβ ; Rd ⊗ Rr ).

次に，Y(1) (α) = −

1
2


Rd

|a(Xt0 , α) − a(Xt0 , α∗ )|2 dt,

det B(Xt0 , β)1tr[B(Xt0 , β)−1 B(Xt0 , β ∗ ) − Id ] + log
2 Rd
det B(Xt0 , β ∗ )
1
B(Xt0 , β ∗ )−1 [(a(Xt0 , α) − a(Xt0 , α∗ ))⊗2 ]dt
Y(3) (α) = −2 Rd
Y(2) (β) = −

dt,

とおく．
[E3] すべての α ∈ Θα , β ∈ Θβ に対しY(1) (α) ≤ −χ(1) |α − α∗ |2 ,Y(2) (β) ≤ −χ(2) |β − β ∗ |2 ,Y(3) (α) ≤ −χ(3) |α − α∗ |2となる正定数 χ(1) , χ(2) , χ(3) が存在する．
フィッシャー情報行列 I(θ∗ ) をI(θ∗ ) =Iaij (θ∗ ) =(Iaij (θ∗ ))1≤i,j≤p

0

0(Ibij (β ∗ ))1≤i,j≤q,(∂αi a(Xt0 , α∗ )) B(Xt0 , β ∗ )∂αj a(Xt0 , α∗ )dt,
1tr{B −1 (∂βi B)B −1 (∂βj B)(Xt0 , β ∗ )}dt
Ibij (β ∗ ) =
2 Rd
Rd

とおく．
4.2 適応的推定量
(1)
(2)
(3)
(α), U,n (α, β), U,n (α, β) は次のように表される．
疑似対数尤度関数 U,n(1)U,n(α) = −

1
|ΔXi − hn ai−1 (α)|2 (2 hn )−1 ,
2 i=1
n

(2)
(α, β) = −
U,n

1
−1
{log det Bi−1 (β) + (2 hn )−1 Bi−1(β)[(ΔXi − hn ai−1 (α))⊗2 ]},
2 i=1

(3)
(α, β) = −
U,n

1 2
−1
( hn )−1 Bi−1(β)[(ΔXi − hn ai−1 (α))⊗2 ].
2 i=1

n

n

高頻度データに基づく確率微分方程式モデルのハイブリッド推定59(1)
初期推定量 α̂,nと適応的最尤型推定量 α̂,n と β̂,n は次で定義される．
(1)
(1)
(1)
U,n
(α̂,n) = sup U,n
(α),
α∈Θα

(2)
(1)
U,n
(α̂,n, β̂,n ) =(2)(1)
sup U,n(α̂,n
, β),

β∈Θβ

(3)
(3)
(α̂,n , β̂,n ) = sup U,n(α, β̂,n ).
U,n
α∈Θα

定理 3. （Nomura and Uchida, 2016a）[E1], [E2], [E3] を仮定する．このとき， √1 n = O(1)の下，√d(−1 (α̂,n − α∗ ), n(β̂,n − β ∗ )) → (ζ1 , ζ2 ) ∼ Np+q (0, I(θ∗ )−1 )が成り立つ．さらに，高々多項式増大のすべての連続関数 f に対して，√Eθ∗ [f (−1 (α̂,n − α∗ ), n(β̂,n − β ∗ ))] → E[f (ζ1 , ζ2 )]が成り立つ．
事前分布 π1 (α), π2 (β) は連続で，0 < inf α∈Θα π1 (α) ≤ supα∈Θα π1 (α) < ∞, 0 < inf β∈Θβ π2 (β) ≤(1)supβ∈Θβ π2 (β) < ∞ を満たすと仮定する．初期ベイズ型推定量 α̃,n および適応的ベイズ型推定量 α̃,n , β̃,n は次で定義される:
(1)

(1)
α̃,n
=

α exp{U,n (α)}π1 (α)dα

Θα

(1)

Θα

β̃,n =

Θβ

exp{U,n (α)}π1 (α)dα
(2)

(1)

(2)

(1)

,

β exp{U,n (α̃,n , β)}π2 (β)dβ
exp{U,n (α̃,n , β)}π2 (β)dβ

Θβ

,

(3)

α̃,n =

Θαα exp{U,n (α, β̃,n )}π1 (α)dα

Θα

(3)

exp{U,n (α, β̃,n )}π1 (α)dα.

定理 4. （Nomura and Uchida, 2016b）[E1], [E2], [E3] を仮定する．このとき， √1 n = O(1)の下，√d(−1 (α̃,n − α∗ ), n(β̃,n − β ∗ )) → (ζ1 , ζ2 ) ∼ Np+q (0, I(θ∗ )−1 )が成り立つ．さらに，高々多項式増大のすべての連続関数 f に対して，√Eθ∗ [f (−1 (α̃,n − α∗ ), n(β̃,n − β ∗ ))] → E[f (ζ1 , ζ2 )]が成り立つ．
次にハイブリッド推定量を定義する．r1 , r2 ∈ (0, 1], r2 ≤ 2r1 γ とする．n → ∞,  → 0 のと√き， √1 n = O(1), ( n)γ = O(1) であることに注意する．
2−2r1 (1)H(1)
U,n (α),
,n,r1 (α) = 
1
(2)
(2)
H,n,r2 (α, β) = √ 2−2r U,n(α, β)
2
( n)
(1)
(2)
とする．初期ベイズ型推定量 α̃,n,r
1 , β̃,n,r2 は

統計数理

60

第 65 巻

第 1 号 2017

(1)

(1)
α̃,n,r=1

α exp{H,n,r1 (α)}π1 (α)dα

Θα

(1)

Θα
(2)
β̃,n,r=2

Θβexp{H,n,r1 (α)}π1 (α)dα
(2)

(1)

(2)

(1)

,

β exp{H,n,r2 (α̃,n,r1 , β)}π2 (β)dβΘβexp{H,n,r2 (α̃,n,r1 , β)}π2 (β)dβにより定義され，ハイブリッド推定量 α̌,n , β̌,n は
(3)
(2)
(3)
(2)
U,n
(α̌,n , β̃,n,r) = sup U,n(α, β̃,n,r
),
2
2
α∈Θα

(2)
U,n
(α̌,n , β̌,n ) =(2)
sup U,n(α̌,n , β)β∈Θβにより定義される．
注 3. r2 ≤ 2r1 γ の条件は，α に関する a(x, α) の滑らかさによって緩和することが可能である．
命題 5. （Nomura and Uchida, 2016b）r1 , r2 ∈ (0, 1] とする．[E1], [E2], [E3] を仮定する．このとき，すべての M > 0 に対し，次が成り立つ．
(1)(i) sup Eθ∗ [|−r1 (α̃,n,r1 − α∗ )|M ] < ∞.
,n√(ii) sup Eθ∗ [|( n)r2 (β̃,n,r2 − β ∗ )|M ] < ∞.
,n(iii) sup Eθ∗ [|−1 (α̌,n − α∗ )|M ] < ∞.
,n√(iv) sup Eθ∗ [| n(β̌,n − β ∗ )|M ] < ∞.
,n定理 5. （Nomura and Uchida, 2016b）[E1], [E2], [E3] を仮定する．このとき， √1 n = O(1)の下，√d(−1 (α̌,n − α∗ ), n(β̌,n − β ∗ )) → (ζ1 , ζ2 )が成り立つ．さらに，高々多項式増大のすべての連続関数 f に対して，√Eθ∗ [f (−1 (α̌,n − α∗ ), n(β̌,n − β ∗ ))] → E[f (ζ1 , ζ2 )]が成り立つ．
注 4. r1 < 1, r2 = 1 のとき，高々多項式増大なすべての連続関数 f に対し，√
(2)
Eθ∗ [f (−1 (α̌,n − α∗ ), n(β̃,n,r− β ∗ ))] → E[f (ζ1 , ζ2 )]2となる．
4.3 例とシミュレーション結果
次式で与えられる 1 次元拡散過程を考える．
⎧
 2⎨ dX  = (α − α X  − 2 sin (α X  ))dt +  β2 + (Xt ) dw
1
2 t
3 t
t
t
1 + β1 (Xt )2⎩
X0 = 2.

t ∈ [0, 1],

 ∈ (0, 1],

ここで，θ = (α1 , α2 , α3 , β1 , β2 ) は未知パラメータ，真値は θ∗ = (α1∗ , α2∗ , α3∗ , β1∗ , β2∗ ) =(3, 7, 5, 0.5, 5)，パラメータ空間は Θ = [0.01, 20]5 であるとする．データ (Xti ) は  = 0.01, 0.05, 0.1,
n = 100, 1000 に対して発生させる．定理 5 のハイブリッド推定量 θ̌ = (α̌,n , β̌,n ) と定理 3 の適
応的最尤型推定量 θ̂ = (α̂,n , β̂,n ) をシミュレーションにより検証する．適応的最尤型推定量は，

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

61

表 5．適応的最尤型推定量 θ̂（初期値は真値）．

表 6．初期ベイズ型推定量 α̃(1) , β̃ (2)（初期値は真値）．

(1)
(1)
(1)
U,n
(α̂,n) = sup U,n
(α),
α∈Θα

(2)
(1)
(α̂,n, β̂,n ) =
U,n

(2)
(1)
sup U,n(α̂,n
, β),

β∈Θβ

(3)
(3)
(α̂,n , β̂,n ) = sup U,n(α, β̂,n )
U,n
α∈Θα

を満たすものである．
真のモデルから 1000 本の独立なサンプルパスを発生させ，各々の推定量の平均と標準偏差を計算する．適応的最尤型推定量 θ̂ とハイブリッド推定量 θ̌ は R 言語の組込み関数 optim() を使って計算し，初期ベイズ型推定量 α̃, β̃ は Kamatani（2014）が提案した MCMC 法のアルゴリ
ズムを用いて計算される．マルコフチェーンの生成数は 5 × 104 回，バーンイン回数は 5 × 103とする．初期値 (α10 , α20 , α30 , β10 , β20 ) には，真値あるいは Θ 上の一様乱数を用いる．
表 5–6 では真値，表 8–9 では一様乱数を初期値としている．表 7 と表 10 は，初期ベイズ型推

62

統計数理

第 65 巻

第 1 号 2017

表 7．ハイブリッド推定量 θ̌（初期値は初期ベイズ型推定量 α̃(1) , β̃ (2) ）．

表 8．適応的最尤型推定量 θ̂（初期値は乱数）．

定量 α̃, β̃ に対し r1 = r2 = 1 としている．表 11–12 では， = 0.05, n = 1000，初期ベイズ型推定量に対し r1 = 0.1, 0.2, . . . , 1, r2 = 1 としている．
とハイブリッド推定量（表 7）はほぼ同じ結果を示しており，いずれ
適応的最尤型推定量（表 5）も  = 0.1 のときに不安定となる．しかし表 8 に見られるように，初期値に Θ 上の一様乱数を用いた場合，適応的最尤型推定量は真値を用いた場合のようにうまく推定できていない．一方，とハイブリッド推定量（表 10）は初期値に Θ 上の一様乱数を用いたと
初期ベイズ型推定量（表 9）
しても良い結果を出している．
 = 0.05, n = 1000 としたときの推定量 α̃3 , α̌3 の標準偏差に着目すると，表 6 と表 9 の α̃3 の標準偏差に明らかな違いが生じており，表 7 と表 10 の α̌3 の標準偏差についても同様である．
これは表 9 の初期ベイズ型推定量がパラメータ α3 の推定に数回失敗していることに起因している．表 10 からハイブリッド推定量は初期ベイズ型推定量の影響を受けていることがわかる．
一方，表 11 の初期ベイズ型推定量 α̃3 の標準偏差は r1 = 0.2 または 0.3 のとき，表 6 の値に近

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

63

表 9．初期ベイズ型推定量 α̃(1) , β̃ (2)（初期値は乱数）．

表 10．ハイブリッド推定量 θ̌（初期値は初期ベイズ型推定量 α̃(1) , β̃ (2) ）．

づいており，表 12 のハイブリッド推定量 α̌3 の標準偏差は，r1 = 0.2 または 0.3 のとき，表 7 の値と非常に近い値を取っている．これは表 11 において r1 を調整することによって得られる初期ベイズ型推定量が良いふるまいをすること，表 12 では r1 = 0.2 または 0.3, r2 = 1 の場合にハイブリッド推定量が，真値を初期値として求められた適応的最尤型推定量（表 5）と同程度によいことが分かる．
5.

おわりに

本論文では，確率微分方程式モデルの未知パラメータに対して，高頻度データを用いたハイブ
リッド推定法について解説した．確率微分方程式（SDE）モデルの統計解析および数値シミュレー
ションの統計パッケージ Yuima を開発するための Yuima プロジェクト（Brouste et al., 2014）が立ち上がって 10 年以上の歳月が経っている．著者はほとんど貢献できていないが，実際に計算機で SDE の統計解析を実行しようとすると，SDE モデル特有の問題点も多々あるが，高次元パ

64

統計数理

第 65 巻

第 1 号 2017

表 11．初期ベイズ型推定量 α̃(1) , β̃ (2)（r2 = 1）．

表 12．ハイブリッド推定量 θ̌（r2 = 1, 初期値は初期ベイズ型推定量 α̃(1) , β̃ (2) ）．

高頻度データに基づく確率微分方程式モデルのハイブリッド推定

65

ラメータ空間における疑似尤度関数の最適化やベイズ型推定量の MCMC 法による統計計算など普遍的な問題が浮かび上がる．エルゴード的拡散過程モデルのパラメータ推定における顕著な性質は，ドリフトパラメータとボラティリティパラメータの推定量の収束率が異なる点である．この性質を用いて効率よく推定量の計算を実行したいというのが，ハイブリッド推定の最初の動機付けである．これは Yuima への実装という具体的な問題から生まれたものである．確，Prakasa Rao
率微分方程式モデルの統計解析については，Kutoyants（1984, 1994, 2004, 2015）（1999），Iacus（2008），Sørensen（2012），吉田（2003, 2010, 2011），林・吉田（2008），内田（2008a,，増田（2015）などを参照．SDE モデルの最尤型推定量やベイズ推定量の理論
2008b, 2009, 2013）
的性質に基づいて，これら推定量を効率よく計算する方法の開発が今後益々重要視されると考
えられる．
謝

辞

査読者から有益なコメントを頂いた．ここに記して感謝の意を表したい．本研究は JSPS 科
研費 24300107 および CREST, JST の助成を受けたものである．

参

考

文

献

Adams, R. A. and Fournier, J. J. F. (2003). Sobolev Spaces, 2nd ed., Pure and Applied Mathematics,
140, Elsevier/Academic Press, Amsterdam.
Azencott, R. (1982). Formule de Taylor stochastique et développement asymptotique d’intégrales deFeynmann, Séminaire de Probabilités XVI ; Supplément: Géométrie Diﬀérentielle Stochastique,
Lecture Notes in Mathematics, 921, 237–285, Springer Verlag, Berlin.
Bibby, B. M. and Sørensen, M. (1995). Martingale estimating functions for discretely observed diﬀusion
processes, Bernoulli, 1, 17–39.
Brouste, A., Fukasawa, M., Hino, H., Iacus, S., Kamatani, K., Koike, Y., Masuda, H., Nomura, R.,
Shimuzu, 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.
Dohnal, G. (1987). On estimating the diﬀusion coeﬃcient, Journal of Applied Probability, 24, 105–114.
Florens-Zmirou, D. (1989). Approximate discrete time schemes for statistics of diﬀusion processes,
Statistics, 20, 547–557.
Freidlin, M. I. and Wentzell, A. D. (1998). Random Perturbations of Dynamical Systems, 2nd ed.,Springer-Verlag, New York.
Fuchs, C. (2013). Inference for Diﬀusion Processes: With Applications in Life Sciences, Springer,
Berlin.
Fujii, T. and Uchida, M. (2014). AIC type statistics for discretely observed ergodic diﬀusion processes,
Statistical Inference for Stochastic Processes, 17, 267–282.
Genon-Catalot, V. (1990). Maximum contrast estimation for diﬀusion processes from discrete observations, Statistics, 21, 99–116.
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.
Genon-Catalot, V. and Jacod, J. (1994). Estimation of the diﬀusion coeﬃcient for diﬀusion processes:
Random sampling, Scandinavian Journal of Statistics, 21, 193–221.
Gloter, A. and Sørensen, M. (2009). Estimation for stochastic diﬀerential equations with a small diﬀu-

66

統計数理

第 65 巻

第 1 号 2017sion coeﬃcient, Stochastic Processes and Their Applications, 119, 679–699.
Gobet, E. (2001). Local asymptotic mixed normality property for elliptic diﬀusion: A Malliavin calculus
approach, Bernoulli, 7, 899–912.
Gobet, E. (2002). LAN property for ergodic diﬀusions with discrete observations, Annales de l’Institut
Henri Poincaré Probabilités et Statistiques, 38, 711–737.
Guy, R., Laredo, C. and Vergu, E. (2014). Parametric inference for discretely observed multidimensional diﬀusions with small diﬀusion coeﬃcient, Stochastic Processes and Their Applications,
124, 51–80.
林高樹, 吉田朋広 (2008). 高頻度金融データと統計科学,『21 世紀の統計科学 I：社会・経済の統計科学』,

第 10 章, 東京大学出版会, 東京.
Iacus, S. (2000). Semiparametric estimation of the state of a dynamical system with small noise, Statistical Inference for Stochastic Processes, 3, 277–288.
Iacus, S. M. (2008). Simulation and Inference for Stochastic Diﬀerential Equations: With R Examples,
Springer, New York.
Iacus, S. and Kutoyants Yu. A. (2001). Semiparametric hypotheses testing for dynamical systems withsmall noise, Mathematical Methods of Statistics, 10, 105–120.
Ibragimov, I. A. and Has’minskii, R. Z. (1981). Statistical Estimation, Springer Verlag, New York.
Jacod, J. and Shiryaev, A. N. (2003). Limit Theorems for Stochastic Processes, 2nd ed., Springer-Verlag,
NewYork.
Kamatani, K. (2014). Eﬃcient strategy for the Markov chain Monte Carlo in high-dimension withheavy-tailed target probability distribution, arXiv:1412.6231.
Kamatani, K. and Uchida, M. (2015). Hybrid multi-step estimators for ergodic diﬀusion processes fromdiscrete observations, Statistical Inference for Stochastic Processes, 18, 177–204.
Kamatani, K., Nogita, A. and Uchida, M. (2016). Hybrid multi-step estimation of the volatility forstochastic regression models, Bulletin of Informatics and Cybernetics, 48, 19–35.
Kessler, M. (1995). Estimation des paramètres d’une diﬀusion par des contrastes corrigés, Comptes
Rendus de l’Académie des Sciences, Series I, Mathematics, 320, 359–362.
Kessler, M. (1997). Estimation of an ergodic diﬀusion from discrete observations, Scandinavian Journalof Statistics, 24, 211–229.
Kunitomo, N. and Takahashi, A. (2001). The asymptotic expansion approach to the valuation of interest
rate contingent claims, Mathematical Finance, 11, 117–151.
Kutoyants, Yu. A. (1984). Parameter estimation for stochastic processes (ed. B.L.S. Prakasa Rao),
Heldermann, Berlin.
Kutoyants, Yu. A. (1994). Identiﬁcation of Dynamical Systems with Small Noise, Kluwer, Dordrecht.
Kutoyants, Yu. A. (2004). Statistical Inference for Ergodic Diﬀusion Processes, Springer-Verlag, London.
Kutoyants Yu. A. (2015). On multi-step MLE-process for ergodic diﬀusion, arXiv 1504.01869, Preprint.
Laredo, C. F. (1990). A suﬃcient condition for asymptotic suﬃciency of incomplete observations of adiﬀusion process, The Annals of Statistics, 18, 1158–1171.
Lehmann, E. L. (1999). Elements of Large-Sample Theory, Springer-Verlag, New York.
Masuda, H. (2013a). Asymptotics for functionals of self-normalized residuals of discretely observed
stochastic processes, Stochastic Processes and Their Applications, 123, 2752–2778.
Masuda, H. (2013b). Convergence of Gaussian quasi-likelihood random ﬁelds for ergodic Levy driven
SDE observed at high frequency, The Annals of Statistics, 41, 1593–1641.
増田弘毅 (2015). 非正規ノイズ型エルゴード過程の推定, 日本統計学会誌, 44(2), 471–495.
Nomura, R. and Uchida, M. (2016a). Adaptive maximum likelihood type estimation for small diﬀusion
processes from discrete observations, Preprint, http://www.sigmath.es.osaka-u.ac.jp/˜uchida/高頻度データに基づく確率微分方程式モデルのハイブリッド推定

67

Paper/AdaML small.pdf.
Nomura, R. and Uchida, M. (2016b). Adaptive Bayes estimators and hybrid estimators for small diﬀusion processes based on sampled data, Journal of the Japan Statistical Society, 46, 129–154.
Ogihara, T. and Yoshida, N. (2011). Quasi-likelihood analysis for the stochastic diﬀerential equation
with jumps, Statistical Inference for Stochastic Processes, 14, 189–229.
Prakasa Rao, B. L. S. (1983). Asymptotic theory for nonlinear least squares estimator for diﬀusion processes, Mathematische Operationsforschung und Statistik, Series Optimization, 14, 195–209.
Prakasa Rao, B. L. S. (1988). Statistical Inference from Sampled Data for Stochastic Processes, Contemporary Mathematics, 80, 249–284, American Mathematical Society, Providence, Rhode Island.
Prakasa Rao, B. L. S. (1999). Statistical Inference for Diﬀusion Type Processes, Edward Arnold, London.
R Development Core Team (2013). R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.org/.
Robert, C. P. and Casella, G. (2004). Monte Carlo Statistical Methods, 3rd ed., Springer Verlag, NewYork.
Shimizu, Y. (2006). M-estimation for discretely observed ergodic diﬀusion processes with inﬁnitely many
jumps, Statistical Inference for Stochastic Processes, 9, 179–225.
Shimizu, Y. and Yoshida, N. (2006). Estimation of parameters for diﬀusion processes with jumps fromdiscrete observations, Statistical Inference for Stochastic Processes, 9, 227–277.
Sørensen, M. (2000). Small dispersion asymptotics for diﬀusion martingale estimating functions,
Preprint No. 2000-2, Department of Statistics and Operations Research, University of Copenhagen, http://www.math.ku.dk/˜michael/smalld.pdf.
Sørensen, M. (2012). Estimating functions for diﬀusion-type processes, Statistical Methods for Stochastic Diﬀerential Equations (eds. M. Kessler, A. Lindner and M. Sørensen), 1–107, CRC Press,
Chapmann and Hall, Boca Raton.
Sørensen, M. and Uchida, M. (2003). Small diﬀusion asymptotics for discretely sampled stochastic
diﬀerential equations, Bernoulli, 9, 1051–1069.
Takahashi, A. and Yoshida, N. (2004). An asymptotic expansion scheme for optimal investment problems, Statistical Inference for Stochastic Processes, 7, 153–188.
Uchida, M. (2003). Estimation for dynamical systems with small noise from discrete observations,
Journal of the Japan Statistical Society, 33, 157–167.
Uchida, M. (2004). Estimation for discretely observed small diﬀusions based on approximate martingale
estimating functions, Scandinavian Journal of Statistics, 31, 553–566.
Uchida, M. (2006). Martingale estimating functions based on eigenfunctions for discretely observed
small diﬀusions, Bulletin of Informatics and Cybernetics, 38, 1–13.
Uchida, M. (2008). Approximate martingale estimating functions for stochastic diﬀerential equations
with small noises, Stochastic Processes and Their Applications, 118, 1706–1721.

内田雅之 (2008a). 確率微分方程式の母数推定,『21 世紀の統計科学 III：数理・計算の統計科学』, 第 7 章,東京大学出版会, 東京.
内田雅之 (2008b). 離散観測に基づく拡散過程の統計的推測, 数学, 60(4), 368–379.
内田雅之 (2009). 小さな拡散過程のドリフトパラメータの推定, 統計数理, 57(1), 67–81.
Uchida, M. (2010). Contrast-based information criterion for ergodic diﬀusion processes from discreteobservations, Annals of the Institute of Statistical Mathematics, 62, 161–187.
内田雅之 (2013). 確率微分方程式の統計的モデリング, 日本統計学会誌, 43(2), 335–358.
Uchida, M. and Yoshida, N. (2001). Information criteria in model selection for mixing processes, Statistical Inference for Stochastic Processes, 4, 73–98.
Uchida, M. and Yoshida, N. (2004a). Information criteria for small diﬀusions via the theory of Malliavin-

68

統計数理

第 65 巻

第 1 号 2017

Watanabe, Statistical Inference for Stochastic Processes, 7, 35–67.
Uchida, M. and Yoshida, N. (2004b). Asymptotic expansion for small diﬀusions applied to option pricing, Statistical Inference for Stochastic Processes, 7, 189–223.
Uchida, M. and Yoshida, N. (2006). Asymptotic expansion and information criteria. Dedicated to Professor Minoru Siotani on his 80th birthday, SUT Journal of Mathematics, 42, 31–58.
Uchida, M. and Yoshida, N. (2011). Estimation for misspeciﬁed ergodic diﬀusion processes from discrete observations, European Series in Applied and Industrial Mathematics: Probability andStatistics, 15, 270–290.
Uchida, M. and Yoshida, N. (2012). Adaptive estimation of an ergodic diﬀusion process based on sampled data, Stochastic Processes and Their Applications, 122, 2885–2924.
Uchida, M. and Yoshida, N. (2013). Quasi likelihood analysis of volatility and nondegeneracy of statistical random ﬁeld, Stochastic Processes and Their Applications, 123, 2851–2876.
Uchida, M. and Yoshida, N. (2014). Adaptive Bayes type estimators of ergodic diﬀusion processes fromdiscrete observations, Statistical Inference for Stochastic Processes, 17, 181–219.
Uchida, M. and Yoshida, N. (2016). Model selection for volatility prediction, The Fascination of Probability, Statistics and Their Applications, In Honour of Ole, E. Barndorﬀ-Nielsen (eds. Mark
Podolskij, Robert Stelzer, Steen Thorbjørnsen and Almut E. D. Veraart), 343–360, Springer,
Switzerland.
Yoshida, N. (1992a). Asymptotic expansion of maximum likelihood estimators for small diﬀusions via
the theory of Malliavin-Watanabe, Probability Theory and Related Fields, 92, 275–311.
Yoshida, N. (1992b). Asymptotic expansion for statistics related to small diﬀusions, Journal of theJapan Statistical Society, 22, 139–159.
Yoshida, N. (1992c). Estimation for diﬀusion processes from discrete observation, Journal of Multivariate Analysis, 41, 220–242.
Yoshida, N. (1993). Asymptotic expansion of Bayes estimators for small diﬀusions, Probability Theory
and Related Fields, 95, 429–450.
Yoshida, N. (2003). Conditional expansions and their applications, Stochastic Processes and Their
Applications, 107, 53–81.

吉田朋広 (2003). Malliavin 解析と数理統計, 数学, 55, 225–244.
吉田朋広 (2010). 確率過程の統計学：概観と展望, 日本統計学会誌, 40(1), 47–60.
Yoshida, N. (2011). Polynomial type large deviation inequalities and quasi-likelihood analysis forstochastic diﬀerential equations, Annals of the Institute of Statistical Mathematics, 63, 431–479.
吉田朋広 (2011). 拡散型確率過程の推定における極限定理, 統計数理, 59(1), 125–140.

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

Hybrid Estimation for Stochastic Diﬀerential Equations
Based on High-frequency Data
Masayuki Uchida1,2,3
1 Graduate School of Engineering Science, Osaka University
2 MMDS, Osaka University
3 CRESTIn this paper, we survey previous researches on hybrid estimation for unknown parameters of stochastic diﬀerential equations based on high-frequency data. Using a Bayes
type estimator with a non-optimal rate of convergence as the initial estimator, we obtain amulti-step estimator and an adaptive maximum likelihood type estimator, and show their
asymptotic properties. For three kinds of diﬀusion models, ergodic diﬀusions, non-ergodic
diﬀusions and small diﬀusions, we give some examples and simulation results.

Key words: Bayes type estimator, diﬀusion process, hybrid estimator, maximum likelihood type estimator, multi-step estimator, stochastic diﬀerential equation.