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


特集「高頻度金融データに基づく統計的推測とモデリング」［原著論文］

高頻度データに対する Whittle 推定
深澤 正彰1,2（受付 2016 年 6 月 24 日；改訂 9 月 28 日；採択 10 月 17 日）

要旨連続伊藤過程の拡散項に対する，高頻度データに基づく統計的推定問題を考察する．対象の伊藤過程そのものの値が高頻度データとして与えられている場合には，既にかなりの理論が整備されている．とくに Euler - 丸山近似に基づく尤度の近似が漸近有効な推定を与えることはよく知られている．ここでは対象とする伊藤過程そのものの値は直接観測されず，その積分値のみが高頻度に観測される状況を扱う．データの数値微分による近似を安直に用いると，推定の一致性さえ崩れてしまう．本稿ではとくに伊藤過程の二次変分をデータで近似する際の中心極限定理を与える．定常 Gauss 過程に対する Whittle 尤度のアイデアを用いて，漸近分散が小さ
い推定量を構成する．
キーワード：高頻度データ，Whittle 推定，中心極限定理，安定収束，Langevin モデル．
1.

研究の背景

溶媒中の微粒子運動に対する Langevin モデルは以下のように記述される：
mŸt = −∇q(Yt ) − γ Ẏt + σ Ẇtまたは同値な表現として(1.1)

dYt = Xt dt,
mdXt = −∇q(Yt )dt − γXt dt + σdWt .

ここで X, Y はそれぞれ粒子の速度と座標を表し，m は質量，q はポテンシャル，γ > 0 は抵抗係数，W は標準 Brown 運動である．以下，物理学者が言う Brown 運動（Brown が発見した微粒子の不規則な運動）と数学的に定義された Brown 運動を区別する必要があるため，後者を指すときには常に，標準 Brown 運動と呼ぶことにする．上の式で標準 Brown 運動は W だが，物理学者にとっては Y が Brown 運動である．拡散係数 σ が 0 のとき，上式は Newton の運動方
程式である．そうでないときは q に対する適当な条件の下で (X, Y ) はエルゴード的拡散過程となり，その不変分布は
222G(dxdy) = Ce−mγx /σ e−2γq(y)/σ dxdy,

で与えられる．ここで C は規格化定数である（例えば Mattingly et al., 2002 参照）．この X, Y の不変分布がそれぞれ Maxwell 分布，Boltzmann-Gibbs 分布に等しいことを要請すると，Einstein
12

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

統計数理

72

第 65 巻

第 1 号 2017

関係 σ 2 = 2γκB T が従う．ただし κB は Boltzmann 定数，T は温度である．
Einstein は 1905 年前後の論文において，微粒子の Brown 運動を，目に見えない溶媒分子との衝突の結果だと仮定し，そこから理論的に予言される Brown 運動の性質を実験的に検証することで，目に見えない溶媒分子の仮定そのものを検証しようと提唱した．当時，溶媒分子の存在を直接確認する方法はなく，全ての物質が原子からなるという原子論は単なる仮説であったことに注意する．Langevin モデルは，この Einstein のアイデアを表現するものとして Langevinが 1908 年に導入したもの（の現代的な記述）である．同時期に Perrin が，Brown 運動する微粒子の軌跡を顕微鏡を使って 30 秒ごとに記録し，その実験データが Einstein の予言通りだった
ことで，初めて原子論が決定的となった．
この歴史的研究のもう少し詳細を見てみよう．Perrin は 2 次元データをとったが，ここでは（q = 0）とすると，（1.1）は線形なシステム
簡単のため 1 次元で考える．今ポテンシャルはないだから Y は Gauss 過程である．定常性の仮定の下，


2κB T
2κB T
m
E[(Yt+h − Yt )2 ] =h + (e−γh/m − 1) ≈hγγγと計算できる．ここで実験対象となる微粒子の抵抗と質量の比は例えば γ/m ≈ 105 程度である．したがって微粒子の変位の二乗 (Yt+h − Yt )2 の平均は h に比例するはずだと言ったのがEinstein であり，これを実験で検証したのが Perrin であった．

上の近似は，Brown 運動 Y を標準 Brown 運動の 2κB T /γ 倍とみなすことに相当するが，この近似を考えた Einstein 本人も述べたように（Nerburgh et al., 2006 参照），時間幅 h が小さい
ときには明らかな問題が生じる．実際計算してみれば(1.2)

E[(Yt+h − Yt )2 ]
= 0,
h→0
h
lim

E[(Yt+h − Yt )2 ]
κB T
=
h→0
h2
m
lim

である．しかし Perrin の h = 30 秒では全く問題なかった．以来，本来なら微分可能なはずのBrown 運動 Y を，標準 Brown 運動の定数倍と同一視することが物理その他の分野の標準となっている．ちなみに上の近似は γ が小さいときにもよくない．実際，気体中では抵抗が小さいので，粒子の軌跡は標準 Brown 運動とはかけ離れている．実際に観察した物理学者の言葉を（図 1 参照）．
借りれば，気体中の Brown 運動は日本舞踊，液体中の Brown 運動は盆踊りであると呼ばれる技術で，1 分子をトラップしてその動き
近年では，光ピンセット（optical tweezer）
を直接観察できるようになった（例えば Gittes and Schmidt, 1998 参照）．これはレーザーで人為的にポテンシャル q を与え，分子を一定期間顕微鏡の視野に止めておく技術である．時間解像度も進歩し，現在の観測間隔はサブミリ秒である．このような 1 分子計測の技術発展により，
特に分子生物学分野で新たな知見が期待されるが，この新しい高頻度データの解析法に関する研究は未発展である．特に上述の議論から，高頻度データに対しては Einstein 近似ではなく，Langevin モデル（1.1）
そのものによる解析が必要と思われる．実際 Li et al.（2010）では，気体中のガラスビーズの Brown 運動をマイクロ秒単位で観測し，（1.2）の成立を確かめている．本稿では，非平衡状態を含むよう一般化したモデル(1.3)

dYt = Xt dt
dXt = Ut dt +√Vt dWtにおいて，Y の高頻度データ {Yjh ; j = 0, 1, 2, . . . , [1/h]} に基づく，X の二次変分
 1
(1.4)
V̄ :=
Vt dt
0

高頻度データに対する Whittle 推定73

図 1．γ/m = 10（上），γ/m = 105（下）の場合の（1）のサンプルパス（ただし q = 0）．

の推定を考察する．ここで U , V は適合過程であり，Y は 1 次元とする．尚，このモデルは計
量ファイナンス分野におけるボラティリティ推定の問題にも現れる．
2.

関連する先行研究

上述の Langevin モデルで，調和振動，つまりポテンシャル q が二次関数の場合には，観測系
列 {Yjh } は Gauss 過程となる． 隠れ Markov 過程 Y は連続時間自己回帰過程と呼ばれる，統において既にこの連続過程の離散観測が議論され計でも歴史あるモデルとなり，Bartlett（1946）ている．離散観測系列は ARMA の構造を持つことが示され，観測時間間隔 h を固定，観測期間が無限に長くなる設定で，標準的な定常 Gauss 過程に対する時系列解析が適用できる．パラ

統計数理

74

第 65 巻

第 1 号 2017

メータの同定可能性が失われないかに関する議論が Pandit and Wu（1975）にある．ノイズをフラクショナルに拡張したモデルが Tsai and Chan（2005）で扱われている．ポテンシャルが 0 の場合に，Gloter（2001）
は差分をとることで，やはり定常 Gauss 過程に帰着させた．定常 Gauss
過程である限り，モデルは共分散構造で決定されるから，これら時系列解析（スペクトル解析）の方法は最尤法と同等な有効推定量を与える．
ポテンシャルが一般の非線形関数の場合にはもはや正規性は失われるが，もし Y が連続的に観測できるなら，Girsanov - 丸山の定理によって尤度が陽に記述できる．観測期間が長くなると同時に，間隔 h も 0 に収束する設定で，この連続観測の尤度の近似を利用する方法が Brockwell
et al.（2007）
で扱われている．また Pokern et al.（2009）で Gibbs サンプリングによる数値的な推定法が提案されている．またモデルの特異摂動極限の拡散過程でフィッティングする方法
が Papavasiliou et al.（2009），Pavliotis and Stuart（2007）で議論されている．観測されない過程 X が一次元拡散過程のとき，間隔 h を固定，観測期間を無限にする設定で，Ditlevsen and，間隔 h も 0 に近づける設定で Gloter（2006）の研究がある．
Sørensen（2004）
以上の研究はすべて観測期間を無限にする漸近論であり，このときエルゴード性の仮定の下，
X のドリフトに現れるパラメータの推定が可能である．一方，我々の問題のように観測期間を固定，間隔 h のみに関する漸近論ではドリフトの一致推定は不可能である．それは Girsanov 丸山の定理を思い出せば，たとえ連続観測しても有限のパスから，そのドリフトを同定できないことから明らかである．それでも拡散係数に現れるパラメータの一致推定は可能である．拡散過程
dXt = b(Xt )dt + σ(Xt , θ)dWt
に対し，その積分値 Y が時刻 0, 1/n, 2/n, . . . , (n − 1)/n, 1 で観測されるとき，Gloter（2000）は漸近混合正規性

 
−1
2
1
√
∂θ σ(Xt , θ)9n(θ̂n − θ) → MN 0,
dt
16
σ(Xt , θ)
0を持つ推定量 θ̂n を与えている．ここで MN は混合正規分布を表す．その推定量は実は漸近有の LAMN の結果からわかる．そこで彼らは，推定
効でないことが，Gloter and Gobet（2008）√
量の誤差 n(θ̂n − θ) の漸近分散の下界は
 
−1
2
1
∂θ σ(Xt , θ)
1
(2.1)
dt
2
σ(Xt , θ)0であることを示した．これは積分値 Y ではなく X の値そのものが観測された場合と同じ漸近下．LAMN は示されたが，下界を達成する推定量は知られていない．
界である（Gobet, 2001 参照）を扱う．仮に Yjh の代わりに粒子の速度 Xjh ,我々はノンパラメトリックなモデル（1.3）
j = 0, 1, 2, . . . , [1/h] が観測できたとすると，これは高頻度データ解析の最も単純な例題となり，まず離散二次変分の一致性（確率収束）
V̂h0 :=

(2.2)

[1/h]

(Xjh − X(j−1)h) )2 → V̄ (h → 0)

j=1

さらに漸近混合正規性（分布収束）
(2.3)

h

−1/2

(V̂h0 − V̄ ) → MN



 1
0, 2
0

Vt2 dt(h → 0)が成立する（Rootzén, 1980; Jacod and Protter, 1998）．ここで V̄ は（1.4）で定義したものである．

高頻度データに対する Whittle 推定75

この離散二次変分は V̄ の漸近有効推定量であると予想されているが，その証明は高頻度データ解析の重要な未解決問題である．例えば V が W と独立な場合などを含む，限定的な状況においては Clement et al.（2013）
による肯定的な結果がある．
我々の問題では粒子の速度 X を直接観測することはできない．自然な発想で，Xjh の代わりに，計算可能な数値微分 Xjh = (Yjh − Y(j−1)h )/h を用いたくなる．初期値 Y0 を固定すれば{Yjh } → {Xjh } の変換で情報は失われていないことに注意．しかし Gloter（2000）
が発見したことは[1/h]h(Xjh − Xj−1
)2 →

j=2

2
V̄ (h → 0),3

つまり，数値微分による安直な近似では推定の一致性すら崩れてしまうという事実である．一
致推定量として彼は
[1/h]

(2.4)

3
h
(Xjh − Xj−1
)2
2 j=2

を提案し，安定収束
⎫
⎧



⎬
⎨ [1/h]
9 1 2
3
h
h−1/2
(Xjh − Xj−1
)2 − V̄ → MN 0,
Vt dt , (h → 0)
⎭
⎩2
4 0
j=2
を示した．先行研究においてこの推定量（2.4）
より漸近分散が小さいものは得られていない．
において，
（2.4）
より漸近分散が小さい V̄ の推定量，とくにその値が本稿ではこのモデル
（1.3）
 1
(2.5)
2
Vt2 dt0であるものを構成する．高頻度データ解析に初めて Whittle 推定のアイデアを導入することが2
における考察（3 節参照），また上述の本研究の特色である．特殊ケース U = 0, V = σ（定数）Gloter and Gobet（2008）のパラメトリックモデルに対する結果から，この値（2.5）が V̄ の推定√量の漸近分散の下限であることを予想している．たとえば（2.1）において σ = θ（定数）とする（2.5）に対応する．
と，下界は 2θ2 と計算でき，この値が3.

Gauss 過程に対する最尤法について

平均 0 の Gauss 過程 (x1 , . . . , xn ) の尤度は，Σ を自己共分散行列とすれば
1x Σ−1 xexp −, x = (x1 , x2 , . . . , xn )
2(2π)n |Σ|となる．未知パラメータ θ と既知の行列 R に対して Σ = θR なる構造があるとき，このパラ
メータ θ に対する最尤推定量は
n

(3.1)

θ̂n =

1
aij xi xj .
n i,j=1

ここで A = [aij ] は R = [rij ] の逆行列である．この推定量の誤差分散は

統計数理

76

第 65 巻1
E[(θ̂n − θ) ] = 2 E
n
1
n2

=

θ
n2

aij (xi xj − θrij )i,j=1
naij akl (E[xi xj xk xl ] − θ2 rij rkl )i,j,k,l=1
n

2

=

2

n

2

=

第 1 号 2017aij akl (rik rjl + ril rjk )i,j,k,l=12θ2.
n

ここで Wick の公式
E[xi xj xk xl ] = E[xi xj ]E[xk xl ] + E[xi xk ]E[xj xl ] + E[xi xl ]E[xj xk ]を用いた．
で X0 = 0, U = 0, V = （θ 定数），h = 1/n のとき，{h−1/2 (Xjh − X(j−1)h )},我々のモデル
（1.3）
−1/2
h
h
(Xj − X(j−1) )}, いずれをデータとしても，上の設定に当てはまり，θ の有効推定量の分{h
散は 2θ2 /n となる．前者の場合，最尤推定量は離散二次変分 V̂h0 に一致する．後者の場合，まず自己共分散を計算する必要がある：
 t


√
√
1 ih
1 ih
h
h
Xih − Xi−1=
(Xt − Xt−h )dt =Vs dWs dt = h(Ui−1+ Uih ),h (i−1)hh (i−1)h t−h
ここで
h
Ui−1
= h−3/2

Uih = h−3/2 (i−1)h(i−2)h ih(i−1)h√(s − (i − 2)h) Vs dWs ,√(ih − s) Vs dWs .

今，V = θ だから，
E[|Uih |2 ] = E[|Uih |2 ] =となり，i, j ≥ 2 に対して

⎧
2
⎪
⎪
⎨3

if i = jrij = 16⎪⎪⎩

θθ
, E[Uih Uih ] =3
6

if |i − j| = 10otherwise.

一般に未知パラメータに依存する共分散行列 C(θ) に対して，その逆行列を陽に求めるのはは定常過程の共分散関数の逆行列そのものではなく，その近似難しい．そこで Whittle（1953）
を用いた擬似尤度を提案した．共分散関数 c(h, θ) = Eθ [x0 xh ] に対してスペクトル密度はf (λ, θ) =

1
2π

∞

√

c(h, θ)e− −1hλh=−∞で定義される．共分散行列 C(θ) = [cij (θ)], cij (θ) = c(i − j, θ) の逆行列に対する Whittle 近似は行列 W (θ) = [wij (θ)], wij (θ) = w(i − j, θ), π √−1hλe1w(h, θ) =dλ(2π)2 −π f (λ, θ)高頻度データに対する Whittle 推定77で定義される．適当な条件の下で Whittle 推定量は漸近有効となる（例えば Dzhaparidze, 1986
参照）．
我々の行列 R の逆行列に対する Whittle 近似は∞

f (λ) =√
1
1 2 + cos λr0j e− −1jλ =2π j=−∞
2π
3w(k) =

1
2π

より
(3.2)

 π

√

√ √
3e −1kλdλ = 3( 3 − 2)|k|−π 2 + cos λと計算できる．よって Whittle 推定量は
（3.1）
の aij を w(i − j) で置き換えたV̂hw :=

(3.3)

[1/h]

√ √
h
h
3( 3 − 2)|i−j| (Xih − Xi−1)(Xjh − Xj−1
)

i,j=2

となる．
データ {h−1/2 (Xjh −X(j−1)h )} に対する最尤推定量 V̂h0 は，V が定数とは限らない一般の場合に
h
，
（2.3）
を持つのであった．ならばデータ {h−1/2 (Xjh −Xj−1)}
も，V̄ の推定量として良い性質（2.2）
wに対する Whittle 推定量 V̂h も，一般の場合に同様の良い性質を持つのではないだろうか．これが本研究のアイデアである．
4.

高頻度データ解析のツール

フィルター付き確率空間 (Ω, F , P, {Ft }t≥0 ) は通常の仮定を満たすとする．離散フィルトレー
ション {Fjh } を Fjh = Fjh によって定義する．高頻度データ解析において，大数の法則の役割を果たすのが次の補題である．
補題 1. （Genon-Catalot and Jacod, 1993）{Fjh } - 適合過程 Mjh , t ∈ [0, 1] に対し，確率収束[1/h]h
E[|Mjh |2 |Fj−1] → 0 (h → 0)

j=1が成立するとき，確率収束
[t/h]
[t/h]


h
h
h
sup Mj −E[Mj |Fj−1 ] → 0 (h → 0)t∈[0,1]  j=1

j=1が成立する．
これは離散マルチンゲールに Lenglart 不等式を適用すればただちに得られる．中心極限定理に相当するものを記述するために，まず安定収束の概念を導入する．
定義 1. （安定収束）部分 σ 代数 G ⊂ F に対して，完備可分距離空間 E に値をとる確率変数
列 {Zn } が，確率変数 Z に G - 安定収束するとは，任意の有界な G 可測確率変数 A に対して結合分布 (Zn , A) が (Z, A) に分布収束することである．
安定収束は分布収束より強く確率収束より弱い．分布収束するが安定収束しない例：X, Y を非退化独立同分布な確率変数としたとき，G = σ(X, Y ), Z2k+1 = X, Z2k = Y .

統計数理

78

第 65 巻

第 1 号 2017

定理 1. （Genon-Catalot and Jacod, 1993）M を {Ft } - 局所マルチンゲール，M⊥ を M と直交する有界マルチンゲール全体とする．{Fjh } - 適合過程 Mjh と連続過程 Q に対し，確率収束
[t/h]hh
1. sup E[Mj |Fj−1 ] → 0, (h → 0)t∈[0,1]  j=1

[t/h]

2.

h
E[|Mjh |2 |Fj−1] → Qt , (h → 0), ∀t ∈ [0, 1]

j=1
[t/h]

3.

h
E[|Mjh |4 |Fj−1] → 0, (h → 0), ∀t ∈ [0, 1]

j=1
[t/h]

4.

h
E[Mjh (Mjh − M(j−1)h )|Fj−1] → 0, (h → 0), ∀t ∈ [0, 1]

j=1
[t/h]

5.

⊥
⊥
h
E[Mjh (Mjh
− M(j−1)h)|Fj−1] → 0, (h → 0), ∀t ∈ [0, 1], ∀M ⊥ ∈ M⊥

j=1が成立するとき，D[0, 1] 値確率変数列としての F - 安定収束[·/h]

Mjh → ŴQ , (h → 0)

j=1が成立する．ここで Ŵ は確率空間 (Ω, F , P, {Ft }t≥0 ) の適当な拡張の上で定義された，F と独立な標準 Brown 運動．
ちなみにデータのサンプリングが等間隔でない場合には，より一般的な収束定理が必要であ
参照．
る．これについては Fukasawa（2010, 2011）
5.

差分二次形式に対する中心極限定理

フィルター付き確率空間 (Ω, F , P, {Ft }t≥0 ) は通常の仮定を満たすとし，W は {Ft } - Brown
運動，さらに {Ft } - 適合過程 U , V に対して X は（1.3）を満たすとする．ドリフト U に対して
（連続とは限らない）
伊藤過程であるこ
は U/V が局所有界であること，拡散項 V に対しては，とを要請する．つまり V = A + M と分解され，局所マルチンゲール M は可予測二次変分 Mを持ち，さらに適当な局所有界過程 a, q によって t
 t
At = A0 +as ds, M t =
qs ds
0

0

と表現できる．ドリフト U に対する仮定は V が 0 を取りうるモデルでは制約的だが，例えば
U も伊藤過程である，といった仮定で代用できる．データは

1 jh
Xjh =
Xt dt, j = 1, 2, . . . , [1/h]h (j−1)hと記述できる．
（1.4）で定義される V̄ に対す
さて，任意に与えた α = (α0 , . . . , αm−1 ) ∈ Rm , m ∈ N に対し，
る，以下の形の推定量 V̂h (α) を考える：

高頻度データに対する Whittle 推定[1/h]V̂h (α) =(5.1)

h(Xih − Xi−1
)

i=m+1

m−1

79

h
h
αj (Xi−j− Xi−j−1).

j=0

m = 1, α0 = 1 の場合が，数値微分 Xjh による離散二次変分の安直な近似である．m = 1,による推定量（2.4）である．第 3 節で見たように，数値微分を
α0 = 3/2 の場合が Gloter（2000）h, i = 1, 2, . . . の隣接項間に相関が生じている．その相関かとることで，その差分列 Xih − Xi−1
らも情報を取り出すことが，m ≥ 2 を考える動機である．次の定理は Gloter（2000）の中心極限を拡張するものである．
定理（m = 1 の場合）

定理 2. V̂h (α) が V̄ の一致推定量であるための必要十分条件は
2
1
α0 + α1 = 1,
3
6

(5.2)

であり，このとき F - 安定収束
h−1/2 (V̂h (α) − V̄ ) → MN

0, g(α) 1
0

Vt2 dt(h → 0)が成立する．ここでg(α) = α02 +

m−1

m−2

m−3

8
1
1 2 1
4
1
αj2 +
αj αj+1 +αj αj+2 .
α0 α1 + α0 α2 +
α1 +
9
9
36
2 j=1
9 j=1
18 j=1

証明. Girsanov - 丸山変換によって U = 0 と仮定してよい．また局所化によって q は有界と仮定してよい．n = [1/h],
 ih
Vih =
Vt dt
(i−1)hとおく．条件
（5.2）
の下で
n

V̄ =
i=m+1
n

Vih +


=
i=m+1
n

=
i=m+1



 mh
0

Vt dt +

 1
nh

Vt dt

2
1
α0 + α1 Vih + Op (h)
3
6
n

1
1
1
h
+
α0 + α1 Vi−1
α0 Vih + Op (h)
3
6
3
i=m+1

と書ける．そこで


1
1
h
h
h
h
h
Zi,j
= h−1/2 (Xih − Xi−1)(Xi−j− Xi−j−1) − δ0,j (Vih + Vi−1
) − δ1,j Vih
3
6

と置けば
h−1/2 (V̂n (α) − V̄ ) =

n

m−1

h
αj Zi,j+ Op (h1/2 )i=m+1 j=0となる．これが離散マルチンゲールの形をしていないのが技術的な問題で，これを定理 1 が使える状況に書き直すことがポイントである．


1 ih
1 h
h
h
−1/2
h
h
h
E[Zi,0 |Fi−1 ] = h] + h|Ui−1|2 −E[Vs |Fi−1]ds − Vi−1
hE[|Uih |2 |Fi−1
3 (i−1)h
3

統計数理

80

=h

−5/2



 ih
(i−1)h

=h

−1/2



第 65 巻

(ih − s)2 −h2
3

第 1 号 2017h]dt
E[Vt |Fi−11 hh+ h−1/2 h|Ui−1|2 − Vi−13h|2 −h|Ui−11 h
Vi−1
3+ O(h3/2 ),


1 h
h
h
h
h
h
E[Zi,1|Fi−1] = h−1/2 hUi−1(Ui−2+ Ui−1) − Vi−1
,
6
h
h
h
h
h
|Fi−1
] = h1/2 Ui−1
(Ui−j−1+ Ui−j), j ≥ 2
E[Zi,jという計算をして，
n

m−1

n

h
αj Zi,j=i=m+1 j=0m−1

i=m+1 j=0
n−1

m−1

=
i=m+1 j=0
n−1

n

h
h
h
αj (Zi,j− E[Zi,j
|Fi−1
]) +

m−1

h
h
αj E[Zi,j|Fi−1]

i=m+1 j=0
n−1

h
h
h
αj (Zi,j− E[Zi,j
|Fi−1
]) +





m−1

h
αj E[Zi+1,j|Fih ] + Op (h1/2 )i=m+1 j=0

1 h
1 h
h
=
α0
Vi + |Uih |2 −
Vi
3h
3h
i=m+1


1 h
h
h
h
+ α1 Uih (Ui−2+ Ui−1) + Uih (Ui−1+ Uih ) −
Vi
6h
m−1

+

1/2

h
Uih + |Uih |2 −
2Ui−1

h
h
h
h
αj (Uih (Ui−j−1+ Ui−j) + Uih (Ui−j+ Ui−j+1))

+ Op (h1/2 )

j=2
n−1

=:
i=m+1

h1/2

m−1

h
αj Fi,j+ Op (h1/2 )j=0を得る．これは近似的に離散マルチンゲールの形をしている：
hhE[Fi,j|Fi−1] = O(h),hhE[Fi,j(Xih − X(i−1)h )|Fi−1] = O(h),

4
5
h
h
h
|2 E[Vih |Fi−1] + 2 E[|Vih |2 |Fi−1] + O(h1/2 )
|Ui−1
3h
9h
1
h 2
h
h
h 2
h
h
h
h
h
E[|Fi,1| |Fi−1]=+ Ui−1) + Ui−1(Ui−2+ Ui−1) + |Ui−1|2 )E[Vih |Fi−1]((Ui−23h
5h+
E[|Vih |2 |Fi−1] + O(h1/2 )
36h2
1
h
h
h
h
h
h
h
h
E[Fi,0
Fi,1
|Fi−1
]=
(Ui−2
+ Ui−1) + |Ui−1|2 )E[Vih |Fi−1](2Ui−13h
2h+ 2 E[|Vih |2 |Fi−1] + O(h1/2 )
9h
1 h
h
h
h
h
h
h
h
h
E[Fi,0
Fi,j
|Fi−1
]=
+ Ui−j) + Ui−j+ Ui−j+1}E[Vih |Fi−1] + O(h1/2 )Ui−1 {2(Ui−j−1
3h
1
h
h
h
h
h
h
h
h
h
E[Fi,1
Fi,j
|Fi−1
]=
+ Ui−1){2(Ui−j−1+ Ui−j) + Ui−j+ Ui−j+1}{(Ui−2
6h
h
h
h
h
h
h
{Ui−j−1+ Ui−j+ 2(Ui−j+ Ui−j+1)}}E[Vih |Fi−1] + O(h1/2 )+ Ui−1h 2
h| |Fi−1]=
E[|Fi,0hhhE[Fi,j
Fi,k
|Fi−1
]=

1
h
h
h
h
h
h
+ Ui−j){2(Ui−k−1+ Ui−k) + Ui−k+ Ui−k+1}{(Ui−j−16h高頻度データに対する Whittle 推定

81

h
h
h
h
h
h
h
+ (Ui−j+ Ui−j+1){Ui−k−1+ Ui−k+ 2(Ui−k+ Ui−k+1)}}E[Vih |Fi−1]+ O(h1/2 ),ここで j, k ≥ 2. すると確率収束
n−1

h 2
h
E[|Fi,0| |Fi−1
]→

h

 1
0

i=m+1
n−1

h 2
h
E[|Fi,1| |Fi−1
]→

h
i=m+1
n−1

Vs2 ds,

19
36

 1

h
h
h
E[Fi,0
Fi,1
|Fi−1
]→

4
9

h
h
h
E[Fi,0
Fi,2
|Fi−1
]→

1
18

h
i=m+1

0

 1
0

Vs2 ds,Vs2 ds,さらに j ≥ 3 に対して
n−1

h
i=m+1
n−1

h

 1
0

Vs2 ds,hhhE[Fi,0
Fi,j
|Fi−1
]→0

i=m+1

が従う． 一方 j, k ≥ 2 に対しては，
n−1

h

h
h
h
E[Fi,j
Fi,k
|Fi−1
]

i=m+1


→

1
2
2
1
1
δj,k + δj+1,k + δj,k+1 +δj+1,k−1 +δj−1,k+1
2
9
9
36
36

 1
0

Vs2 dsで，この確率収束は j = k = 1 の場合を除けば j, k ≥ 1 において正しい．かくして
⎡
⎤
2
n−1
m−1



h
h
h
E ⎣
αj Fi,j  |Fi−1 ⎦


i=m+1
j=1

 1
m−1
m−2
m−2
19 2 14
1→
αj2 +αj αj+1 +αj αj+2
Vs2 ds
α1 +
36
2 j=2
9 j=1
18 j=1
0

 1
m−1
m−2
m−3
4
1
1 2 1
αj2 +
αj αj+1 +αj αj+2
Vs2 ds.
=
α1 +
36
2 j=1
9 j=1
18 j=10
あとは定理 1 を
Mih = h1/2

m−1hαj Fi,j, M =X

j=0

に対して適用すればよい．
定理 2 より，m = 1 のときは一致性のために α0 = 3/2 が必要で，そのとき g(α) = 9/4 = 2.25,
これは Gloter（2000）の結果である．m = 2 の場合，（5.2）の下

統計数理

82

第 65 巻

第 1 号 201719 2 8
α1 + α0 α1
36
9
19
8
2
= α0 +(6 − 4α0 )2 + α0 (6 − 4α0 )
36
9
53 2
= α0 − 20α0 + 19
9

2
90
107
53
+
=
α0 −
.
9
53
53

g(α) = α02 +最適化すると α0 = 90/53, α1 = −42/53 で，そのとき最小値 g(α) = 107/53 ≈ 2.018868 となり，より良いものが構成できたことになる．さらに m = 3 の
この時点で既に先行研究の推定量（2.4）
場合，（5.2）の下
19 2 81
1
4
α1 + α0 α1 + α0 α2 + α22 + α1 α2
36
9
9
2
9
53 2
1
1 2 4
= α0 − 20α0 + 19 + α0 α2 + α2 + (6 − 4α0 )α2
9
9
2
9


2
53
9
81 2
26
107
5
=
+
α0 −
10 + α2
α2 −
α2 +
9
53
6
212
159
53



2
2
9
81
52
77327
5
53
+
+
α0 −
10 + α2
α2 −
=
9
53
6
212
243
38637

g(α) = α02 +で，最適化すると α2 = 52/243, α0 = 7420/4293, α1 = −3922/4293 によって最小値 g(α) =77327/38637 ≈ 2.001372 をとる．m を上げると計算量は上がるが，それでも二次形式の計算なので，通常のデータ量なら計算コストを考慮する必要はない．ではどこまで漸近分散を下げられるだろうか．
補題 2.


inf



2
1
g(α)  α = (α0 , α1 , . . . , αm−1 ) ∈ Rm , α0 + α1 = 1, m ∈ N
3
6

= 2.

また
(5.3)

α0 =

√

√ √
3, αj = 2 3( 3 − 2)j , j = 1, 2, . . . , m − 1とすると g(α) → 2 (m → ∞).
証明. ある m ∈ N, α ∈ Rm が存在して g(α) < 2 であるとする．定理 2 の中心極限定理は
U = 0, V = θ (定数), h = 1/n の場合にも勿論成立するから，
h−1/2 (V̂ (α) − θ) → N (0, g(α)θ2 )となる．これは定常 Gauss 過程の未知パラメータ θ の推定量で，誤差分散が Cramér-Rao の下
界 2θ2 /n を下回るものが得られたことになり矛盾である．直接 g(α) の定義に代入すれば（5.3）
に対して g(α) → 2 が分かり，下限が 2 であることが従う．
下限を達成する α の列
（5.3）
は Whittle 推定量（3.3）の形をもとに発見したものである．我々のモデル（1.3）は一般に非定常で非 Gauss 過程であるが，それでも高頻度データに対しては定常
Gauss 過程に対する Whittle 推定のアイデアが応用できることを初めて示した次の定理が本研
究の主な貢献である．
（3.3）
によって定義する．このとき h → 0 で F - 安定収束定理 3. 推定量 V̂hw を高頻度データに対する Whittle 推定

h

−1/2

(V̂hw − V̄ ) → MN



 1
0, 2
0

Vt2 dt83



が成立する．
証明. Girsanov - 丸山変換によって U = 0 と仮定してよい．また局所化によって q, V は有界と仮定してよい．
[1/h]V̂hm =i=2√

h3(Xih − Xi−1
)2 + 2[1/h] (m−1)∧(i−2)
i=2

√ √
h
h
h
3( 3 − 2)k (Xih − Xi−1)(Xi−k− Xi−k−1
)

k=1

とおけば
lim sup h−1 E[|V̂hw − V̂hm |2 ] = 0(5.4)

m→∞ h>0である．また，（5.3）
で与えた α に対して (5.1) によって V̂h (α) を定義すると，任意の m に対して確率収束h−1/2 (V̂hm − V̂h (α)) → 0
が成立する．これと定理 2 を合わせると
h−1/2 (V̂hm − V̄ ) → MN

(5.5)


0, g(α) 1
0

Vt2 dtとなる．補題 2 より m → ∞ の時 g(α) → 2 だから（5.4）と（5.5）を合わせて主張を得る．
謝辞24684006 及び基盤研究（B）
24300107 の助成を受けた．
本研究は科学研究費補助金若手研究
（A ）

参

考

文

献

Bartlett, M.S. (1946). On the theoretical speciﬁcation and sampling properties of autocorrelated timeseries, Supplement to the Journal of the Royal Statistical Society, 8(1), 27–41.
Brockwell, P.J., Davis, R.A. and Yang, Yu (2007). Continuous-time Gaussian autoregression, Statistica
Sinica, 17, 63–80.
Clement, E., Delattre, S. and Gloter, A. (2013). An inﬁnite dimensional convolution theorem with applications to the eﬃcient estimation of the integrated volatility, Stochastic Processes and Their
Applications, 123, 2500–2521.
Ditlevsen, S. and Sørensen, M. (2004). Inference for observations of integrated diﬀusion processes,
Scandinavian Journal of Statistics, 31, 417–429.
Dzhaparidze, K. (1986). Parameter Estimation and Hypothesis Testing in Spectral Analysis of Stationary Time Series, Springer-Verlag, New York.
Fukasawa, M. (2010). Realized volatility with stochastic sampling, Stochastic Processes and Their Applications, 120, 829–852.
Fukasawa, M. (2011). Discretization error of stochastic integrals, Annals of Applied Probability, 21,
1436–1465.
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́, Probabilités et Statistiques,29, 119–151.

84

統計数理

第 65 巻

第 1 号 2017

Gittes, F. and Schmidt, C. (1998). Signals and noise in micromechanical measurements, Methods in
Cell Biology, 55, 129–156.
Gloter, A. (2000). Discrete sampling of an integrated diﬀusion process and parameter estimation of thediﬀusion coeﬃcient, ESAIM Probability and Statistics, 4, 205–224.
Gloter, A. (2001). Parameter estimation for a discrete sampling of an integrated Ornstein-Uhlenbeck
process, Statistics: A Journal of Theoretical and Applied Statistics, 35(3), 225–243.
Gloter, A. (2006). Parameter estimation for a discretely observed integrated diﬀusion process, Scandinavian Journal of Statistics, 33, 83–104.
Gloter, A. and Gobet, E. (2008). LAMN property for hidden processes: The case of integrated diﬀusions,Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 44(1), 104–128.
Gobet, E. (2001). Local asymptotic mixed normality property for elliptic diﬀusion: A Malliavin calculus
approach, Bernoulli, 7(6), 899–912.
Jacod, J. and Protter, P. (1998). Asymptotic error distributions for the Euler method for stochastic
diﬀerential equations, Annals of Probability, 26(1), 267–307.
Li, T., Kheifets, S., Medellin, D. and Raizen, M.G. (2010). Measurement of the instantaneous velocity
of a Brownian particle, Science, 328, 1673–1675.
Mattingly, J.C., Stuart, A.M. and Higham, D.J. (2002). Ergodicity for SDEs and approximations: Locally Lipschitz vector ﬁelds and degenerate noise, Stochastic Processes and Their Applications,101, 185–232.
Nerburgh, R., Peidle, J. and Reuckner, W. (2006). Einstein, Perrin, and the reality of atoms: 1905
revisited, American Journal of Physics, 74, 478.
Pandit, S.M. and Wu, S.M. (1975). Unique estimates of the parameters of a continuous stationary
stochastic process, Biometrika, 62(2), 497–501.
Papavasiliou, A., Pavliotis, G.A. and Stuart, A.M. (2009). Maximum likelihood drift estimation for
multiscale diﬀusions, Stochastic Processes and Their Applications, 119, 3172–3210.
Pavliotis, G.A. and Stuart, A.M. (2007). Parameter estimation for multiscale diﬀusions, Journal of
Statistical Physics, 127(4), 741–781.
Pokern, Y., Stuart, A.M. and Wiberg, P. (2009). Parameter estimation for partially observed hypoelliptic diﬀusions, Journal of the Royal Statistical Society: Series B, 71, Part 1, 49–73.
Rootzén, H. (1980). Limit distributions for the error in approximations of stochastic integrals, Annalsof Probability, 8(2), 241–251.
Tsai, H. and Chan, K.S. (2005). Quasi-maximum likelihood estimation for a class of continuous-timelong-memory processes, Journal of Time Series Analysis, 26(5), 691–713.
Whittle, P. (1953). Estimation and information in time series analysis, Arkiv för Matematik, 2(5),423–434.

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

Whittle Estimation for High-frequency Data
Masaaki Fukasawa1,2
1 Graduate School of Engineering Science, Osaka University

2 Center for Mathematical Modeling and Data Science, Osaka University

In this paper, we consider a statistical estimation problem of the diﬀusion term of acontinuous Itô process based on high-frequency data. This problem has been extensively
studied in the literature, where it has been typically assumed that the Itô process itself is observed discretely. In particular, it is well-known that a quasi-likelihood based on
the Euler-Maruyama approximation yields an asymptotically eﬃcient estimator. Here, we
study the case where the Itô process is hidden but its integrated process is observed at
high-frequency. It is known that a naive method that simply uses numerical derivativesof observed integrated processes results in an inconsistent estimation. We prove a central
limit theorem for quadratic forms of ﬁrst-order diﬀerences of the numerical derivatives.
Using Whittle’s approximation of the inverse of a covariance matrix, we construct a consistent estimator of the quadratic variation of the Itô process of which the asymptotic
variance is smaller than those of previously proposed estimators.

Key words: High-frequency data, Whittle estimation, central limit theorem, stable convergence,
Langevin model.