データサイエンス
福島キャンプ 2019
(新しい時系列計量分析の理論と応用)1

国友直人2 (編)2019 年 8 月

1

学術振興会・科学研究プロジェクト「新しい時系列計量分析の理論と応用」(2017 年度〜2020 年度) が主催した（2019 年 8 月 7 日, 福島大学) 研究集会 (明治大学先端数理科学インスティテュート (MIMS) 協賛）
における講演集。
2
明治大学政治経済学部

概要
経済や金融 (ファイナンス) などで観察されるデータなどを分析対象とする統計的時系列分析では幾つかの検討するべき基本的な問題がなお存在している。例えば経済において「通常の常識では起こりにくいとされる事象」についてのリスク解析や対策の重要性についての認識が高まっているが、実際に 2008 年に起きたリーマンショックや 2011 年ごろに発生したヨーロッパ諸国の金融危機などが顕著な例として挙げることができる。国際的に連動している現代の経済・金融市場
では従来の計量分析ではほとんど考慮されてこなかった変動を経験しているのである。社会・経済を理解し、より良いものにしていくには、事前には予想が困難である自然災害や経済変動における稀ではあるが実際に起きると大きな影響のある不確実な事象を科学的に分析し、有効な対策を考察することが必要であり重要となっている。また近年に明らかになりつつあるように高頻度
金融データやマクロ経済データの分析においては、経済変数間の関係の統計的分析が重要である
が、確率過程の統計分析や非定常多次元時系列分析にはなお解決すべき課題が少なくない。実は考察の範囲を広げて自然や生物を研究対象とする分野でも似たような現象が観察され、共通する問題と固有の統計的問題がある。
科学研究プロジェクト「新しい時系列計量分析の理論と応用」では主に経済・金融・社会における新しい動向を背景として、近年の日本経済・社会の理解の方法として重要な新しい計量分析の方法を開発、応用を検討している。特に金融現象やマクロ経済の統計分析では、時々起きる大きな経済変動は重要であるにもかかわらず、なお研究の蓄積が不十分であり、様々な研究の可能性がある。また、ミクロ金融の分析ではジャンプや確率過程の一般理論を踏まえた金融時系列分
析はなお十分とは言えず、マイクロマーケット分析や確率過程の計量分析の方法を確立が望まれる。こうした経済・金融におけるマクロ分野、ミクロ分野と云う二つの時系列分析において新しい分析の枠組みを構築する必要がある。またマクロ経済データが象徴的であるように、経済・金融分野で観察されるデータは統計的には非定常性が見られるとともに本来的に多次元時系列であり、次元数も必ずしも小さくない場合も研究対象である。
本年度は研究プロジェクトの三年目であり、特にファイナンス市場の分析では重要な高次元デー
タ分析、マイクロマーケットの実証分析、から計量経済分析の問題やデータサイエンスの問題をなど含めて活発な議論を行った。ここに収録した研究報告が統計学・データサイエンス、計量経
済学、計量ファイナンス、経済時系列分析、などの理論と応用を展開する一助になることを期待する。

2019 年 9 月
編者

1

研究集会・プログラム

科学研究プロジェクト「新しい時系列計量分析の理論と応用」
日程：2019 年 8 月 7 日 (水)会場：福島大学経済経営棟 5F
オーガナイザー：国友直人
プロジェクト参加者：国友直人 (明治大学)・大屋幸輔 (大阪大学)・佐藤整尚 (東京大学)・栗栖大輔 (東京工業大学)
会場責任者：井上健 (福島大学)
協賛：明治大学先端数理科学インスティテュート (MIMS)< セッション I：高次元データ分析 >Chair: 大屋幸輔
13:00〜13:50「高頻度データにおける高次元共分散行列の統計推測」
小池裕太 (東京大学)

13:50〜14:20「ノイズを含む高頻度高次元ボラティリティ分析」
国友直人< 休憩 >< セッション II：Financial Economics and Econometrics>
Chair: 三崎広海 (筑波大学)
14:30〜15:10「Estimation of risk aversion for Japanese stock market using implied moments」
大屋幸輔

15:10〜15:50 「指値注文市場における気配設定」
太田亘 (大阪大学)< 休憩 >< セッション III：Econometrics>
Chair:国友直人
16:00〜16:40「A two-sample alternative to using instruments in regressions with omitted variables」
蛭川雅之 (龍谷大学)

16:40〜17:10「Nonparametric estimation of density functions from repeated measurements」
栗栖大輔

17:10〜17:40「周波数回帰と経済データへの応用」
佐藤整尚

2

高頻度データにおける高次元共分散行列の統計推測
小池祐太
東京大学 MI センター・大学院数理科学研究科
CREST JST2019 年 8 月 7 日
データサイエンス・福島キャンプ 2019@ 福島大学小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

1 / 49

1

背景2

Chernozhukov-Chetverikov-Kato (CCK) 理論

3

主結果

4

残差スパース性検定への応用

5

数値実験

6

実データ解析例

7

まとめと今後の課題

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

2 / 49背景

Yt = (Yt1 , . . . , Ytd )⊤ (t ∈ [0, 1]): d 個の資産の対数価格過程
Y = (Yt )t∈[0,1] は確率基底 (Ω, F, (Ft ), P) 上で定義された連続 Itô セ
ミマルチンゲールと仮定:∫ t∫ t
Yt = Y0 +µs ds +
σs dBs ,t ∈ [0, 1].
00µ = (µs )s∈[0,1] : d 次元 (Ft )-発展的可測過程
σ = (σs )s∈[0,1] : Rd×r 値 (Ft )-発展的可測過程▶ B = (B )s s∈[0,1] : r 次元標準 (Ft )-Brown 運動▶ 適当なモーメント条件を課す▶
▶

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

3 / 49

背景Y は等間隔な時点 th = thn = h/n, h = 0, 1, . . . , n で観測されているとする目的 観測データ (Yth )nh=0 に基づいて, Y の累積共分散行列∫ 1[Y , Y ]1 =Σt dt,

Σt := σt σt⊤0

を n → ∞ のときに推定する (高頻度観測の設定)
最も自然な [Y , Y ]1 の推定量として, 実現共分散行列n\[Y, Y ]1 :=n∑(Yth − Yth−1 )(Yth − Yth−1 )⊤h=1が考えられる小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

4 / 49背景
近年, 高頻度データ解析の分野でも, n → ∞ のとき d → ∞ なる高次元の設定で [Y , Y ]1 を推定する研究が活発に行われている▶
▶

応用: 最小分散ポートフォリオの重みの計算など
マーケットマイクロストラクチャー効果の影響で, サンプル数 n はせ
いぜい 400 程度しかとれないが, もし d = 100 ならば, 推定すべきパラメーター数は d(d + 1)/2 ≈ 5000 と n に比べて膨大となる!!

上の設定は特に d 変量正規確率変数を n 個独立に観測する状況を含n\み, その場合 [Y, Y ]1 は標本共分散行列に対応するから, 前述と同様の問題がやはり生じる
他方, 問題の解決方法も基本的には i.i.d. の設定のアナロジーで行われる: 累積共分散行列に何らかの意味でのスパース性を課す!!
次頁に代表的な先行研究を挙げる (他にもたくさんある)小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

5 / 49代表的な先行研究
Wang and Zou (2010)[Y , Y ]1 の成分の多くが 0, もしくは非常に小さいというスパース性を課して, バンド化/閾値法によってナイーブな推定量を正則化する▶ Bickel and Levina (2008a,b) の高頻度版▶

Fan, Furger, and Xiu (2016)
Y が連続時間のファクターモデルに従うと仮定し, その残差過程の共分散行列がスパースであると仮定して, ファクター部分は線形回帰, 残差部分は正則化で共分散行列を推定し, それを合算する▶ Fan et al. (2011) の高頻度版▶

Brownlees, Nualart, and Sun (2018)[Y , Y ]1 の逆行列にスパース性を課して, graphical Lasso によって[Y , Y ]1 の逆行列を推定する▶ Yuan and Lin (2007) の高頻度版 (高次元の理論は Ravikumar et al.
(2011) に依拠)▶

. . . しかし, そのようなスパース性の仮定の妥当性はどのように検証すればよいか?
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

6 / 49

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

図 1: S&P100 指数構成銘柄の実現相関行列 (2018 年 3 月)
AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1

15 分足リターンで計算. 長方形はセクターを表す

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

7 / 49

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

図 2: マーケットファクターで回帰後の残差過程の実現相関行列 (2018 年 3 月)
AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1

15 分足リターンで計算. 長方形はセクターを表す. マーケッ
トファクターは SPY で代理

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

8 / 49連続時間ファクターモデルの残差スパース性検定
ここでは Fan et al. (2016) の仮定に着目してみる (Wang and Zou(2010) の仮定は特殊ケースと考えられる)

d 番目の資産 Y d をファクター過程とみなし, 以下の連続時間 1 ファ
クターモデルを考える:
j = 1, . . . , d := d − 1.

Y j = βj Y d + R j ,(1)

β j : 定数 (因子負荷量)▶ R j : [R j , Y d ] ≡ 0 を満たすセミマルチンゲール (残差過程)▶

各ペア (i, j) について, 以下の統計的仮説検定を考える:
(i,j)

H0

: [R i , R j ]1 = 0 a.s.

小池祐太 (東大 MI センター)

vs

(i,j)

H1

高次元共分散行列の統計推測: [R i , R j ]1 ̸= 0

a.s.

2019 年 8 月 7 日(2)

9 / 49連続時間ファクターモデルの残差スパース性検定
ある特定のペア (i, j) に着目した場合の検定 (2) は Bibinger and
Mykland (2016) においてすでに考察されており, 次の検定統計量が提案されている:
√ ij
nT̂n
(i,j)
Tn := √
,
ij
V̂n
ここに,
n

n

n

n

\
i , Y d ] [Y\j , Y d ] − [Yd, Y d]i , Y j ] [Y\T̂ijn := [Y\
1
1
1
1
ij

ij

で V̂n は T̂n の漸近分散の推定量→L N(0, 1) となり,
(i,j)
(i,j)
H1 の下で |Tn | →p ∞ となることが証明されている(i,j)

適当な正則条件を仮定すると, H0

小池祐太 (東大 MI センター)(i,j)の下で Tn高次元共分散行列の統計推測

2019 年 8 月 7 日

10 / 49

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

図 3: S&P100 指数構成銘柄の実現相関行列 (2018 年 3 月): 5%水準で有意な成分=1, それ以外=0
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1

15 分足リターンで計算. 長方形はセクターを表す

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

11 / 49

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

図 4: マーケットファクターで回帰後の残差過程の実現相関行列 (2018 年 3 月):
5%水準で有意な成分=1, それ以外=0
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1

15 分足リターンで計算. 長方形はセクターを表す. マーケッ
トファクターは SPY で代理

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

12 / 49

連続時間ファクターモデルの残差スパース性検定
残差過程の共分散行列の正則化推定のために必要なスパース性を検証するには, 仮説検定 (2) を (i, j) ∈ Λn := {(i, j) : 1 ≤ i < j ≤ d} について同時に実行する必要がある▶

どの成分が有意かどうかも知りたい ⇒ 多重検定の問題

特に, この多重検定問題の family-wise error rate (FWER) をコントロールしたい▶

FWER = 正しい帰無仮説を 1 つでも誤って棄却してしまう確率

FWER のコントロールを達成するために, Romano and Wolf (2005)
で提案されたステップダウン法を利用したい▶

アルゴリズムの詳細については Romano and Wolf (2005) を参照⇒ L ⊂ Λn が与えられたとき, 統計量 max(i,j)∈L |Tn(i,j) | の分布を近似する方法が必要

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

13 / 49

この目的は以下のようなタイプの CLT を証明すると達成できる (“高次元漸近混合正規性”):()
b1/2
sup P (Ξn Sn ≤ y ) − P Ξn C
ζ
≤
y
→0
(3)
n
n
y ∈Rm

√

(
)
n
\
n vec [Y, Y ]1 − [Y , Y ]1

▶

Sn :=

▶

bn は Sn の漸近共分散行列の推定量:Cbn := n
C

n
∑
h=1

)
n ∑(
⊤
χh χ⊤
h+1 + χh+1 χh
2
n−1

χh χ⊤
h −

h=1

[
]
ただし, χh := vec (Yth − Yth−1 )(Yth − Yth−1 )⊤▶

ζn は F と独立な d 2 次元標準正規確率ベクトル▶

Ξn は m × d 2 ランダム行列で, 行ごとにスパースなもの (厳密な定式化は後述); m は n に依存しうる小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

14 / 49

Jacod’s Paradise
もし m, d ともに n に依存しないのであれば, (3) 式は Jacod の安定収束理論の帰結として得られる

定理 1 (Jacod and Protter (2012), Theorem 5.4.2)∫1
確率 1 で 0 (∥µs ∥2ℓ2 + ∥Σs ∥2ℓ2 )ds < ∞ が成り立つならば, Sn は n → ∞ のとき C̄1/2 ζ に安定収束する. ただし, ζ は F と独立な d 2 次元標準正規確率ベクトル, C̄ = (C̄pq )1≤p,q≤d 2 は d 2 次非負定値対称行列に値をとるランダム行列で, その成分は次式で与えられる:C̄d(p−1)+q,d(p ′ −1)+q ′

∫ T(
=

)
′ qq ′
pq ′ p ′ q
Σpp
Σ
+
Σ
Σ
dt,
t
t
t
t

0

p, q, p ′ , q ′ = 1, . . . , d.

定理 1 の仮定は µs ≡ 0 の場合は必要条件!!
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

15 / 49

Paradise Lost

しかし, n → ∞ のとき d → ∞ となる高次元の設定下では, Jacod の理論は完全に破綻する
Jacod の定理の証明は, Skorohod 位相に関する緊密性の議論とあるマ
ルチンゲール問題の解の一意性に強く依拠している▶ しかし, d → ∞ のとき, S は通常の位相に関してはもはや緊密でない!
n▶

高次元設定下での CLT の証明は, 実質的には分布間の距離に対する定量的評価の問題となり, 一歩先の評価が必要▶
▶

次元に対してリーズナブルな上界が欲しい「分布間の距離としてどれを採用するか」も本質的

本研究の貢献 d → ∞ のときに (3) 式のような収束を証明するための別の手段を構築

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

16 / 49

Chernozhukov-Chetverikov-Kato (CCK) 理論

最近, Chernozhukov, Chetverikov, and Kato (2013, 2015, 2017) の一連の研究で, (3) のような形の収束を証明するための新たな理論が開発された▶

Chernozhukov-Chetverikov-Kato 理論, 略して CCK 理論と呼ぶことにする

オリジナルの CCK 理論は, Ξn Sn が独立な確率ベクトルの和, もしく
b1/2
は Gauss 型で, かつ Ξn C
n が非ランダムでなければ適用できない
しかし, 証明のほとんどの部分では上述の数学的構造は使われない▶

証明の一部を修正することで我々の状況にも適用できる!

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

17 / 49CCK 理論: Gaussian comparison の場合
以下 m ≥ 3 と仮定する

定理 2 (Chernozhukov et al. (2015), Theorem 2)
F , Z を平均 0 の m 次元 Gauss 型確率変数とし, ある定数 σ > 0 が存在して minj E [Zj2 ] ≥ σ 2 が成り立つと仮定する. このとき, σ のみに依存する
定数 C > 0 が存在して,sup |P (F ≤ x) − P (Z ≤ x)| ≤ C ∆1/3 log2/3 mx∈Rmが成り立つ. ただし,∆=

小池祐太 (東大 MI センター)

max |E [Fi Fj ] − E [Zi Zj ]|1≤i,j≤m2

高次元共分散行列の統計推測

2019 年 8 月 7 日

18 / 49CCK 理論: 問題の分割
集合 A ⊂ R と ε > 0 に対して以下のように定義:

Aε := {u ∈ R : |u − v | ≤ ε for some v ∈ A}
次の補題は初等的議論で示せる:

補題 1
F , Z を m 次元確率変数とし, ある定数 ε, η > 0 が存在して, 任意のy ∈ Rm と Borel 集合 A ⊂ R に対して()()
P max (Fj − yj ) ∈ A ≤ P max (Zj − yj ) ∈ Aε + η
1≤j≤m

1≤j≤m

が成り立つと仮定する. このとき次が成り立つ:()
sup |P(F ≤ x) − P(Z ≤ x)| ≤ sup P 0 ≤ max (Zj − yj ) ≤ ε + η.
x∈Rm

小池祐太 (東大 MI センター)y ∈Rm

高次元共分散行列の統計推測

1≤j≤m

2019 年 8 月 7 日

19 / 49CCK 理論: 反集中不等式
Z が Gauss 型の場合, 補題 1 の最後の不等式の右辺第 1 項は次の形の反集中不等式を用いて次元に関してリーズナブルに抑えられる:

補題 2 (Nazarov の不等式 (Chernozhukov et al. (2017), Lemma A.1))
Z を m 次元 Gauss 型確率変数とし, ある定数 σ > 0 が存在してminj E [Zj2 ] ≥ σ 2 が成り立つと仮定する. このとき, σ のみに依存する定数
C > 0 が存在して, 任意の ε > 0 に対して()√
sup P 0 ≤ max (Zj − yj ) ≤ ε ≤ C ε log m
y ∈Rm

1≤j≤m

が成り立つ.

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

20 / 49CCK 理論: 最大値関数の近似
β > 0 に対して, 関数 Φβ : Rm → R をm∑Φβ (x) = β −1 log e βxj (x = (x1 , . . . , xm )⊤ ∈ Rm )
j=1

で定義する
任意の x ∈ Rm に対して

0 ≤ Φβ (x) − max xj ≤ β −1 log m

(4)

1≤j≤m

が成り立つ (cf. Chernozhukov et al. (2015), Eq.(1))

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

21 / 49CCK 理論: 指示関数の近似
任意に ε > 0, y ∈ Rm と Borel 集合 A ⊂ R をとる

β = ε−1 log m (従って ε = β −1 log m) とおけば, (4) より()
P max(Fj − yj ) ∈ A ≤ P(Φβ (F − y ) ∈ Aε ) = E [1Aε (Φβ (F − y ))]
j

が成り立つ

補題 3 (Chernozhukov et al. (2016), Lemma 5.1)
任意の ε > 0 と任意の Borel 集合 A ⊂ R に対して, 次の条件を満たす C ∞関数 g : R → R が存在する:(i) ある普遍定数 K > 0 が存在して ∥g ′ ∥∞ ≤ ε−1 , ∥g ′′ ∥∞ ≤ K ε−2 および ∥g ′′′ ∥∞ ≤ K ε−3 が成り立つ.
(ii) 任意の x ∈ R に対して 1A (x) ≤ g (x) ≤ 1A3ε (x) が成り立つ.
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

22 / 49CCK 理論: 期待値の差の評価への帰着
補題 3 を A = Aε として適用すれば,

E [1Aε (Φβ (F − y ))] ≤ E [g (Φβ (F − y ))],E [g (Φβ (Z − y ))] ≤ E [1A4ε (Φβ (Z − y ))]
[
)]
(
)
(
5ε
≤ E 1A5ε max(Zj − yj )= P max(Zj − yj ) ∈ A
j

j

を得る
以上より, σ のみに依存する定数 c1 > 0 が存在して,sup |P(F ≤ x) − P(Z ≤ x)|√≤ c1 ε log m + sup |E [g (Φβ (F − y ))] − E [g (Φβ (Z − y ))]|

x∈Rm

y ∈Rm

が成り立つ
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

23 / 49

以下に注意:▶
▶

以上の議論は任意の確率変数 F に対して成立!
g ◦ Φβ は C ∞ 級関数!

更に, g ◦ Φβ の高次偏導関数の ℓ1 -ノルムは「次元に関して効率的に」評価できる:

補題 4 (Chernozhukov et al. (2014), Lemma 4.3)
g : R → R を C 3 級関数とする. このとき, 任意の x ∈ Rm に対してm∑
∂ 2 (g ◦ Φβ )(x) ≤ ∥g ′′ ∥∞ + 2β∥g ′ ∥∞ ,∂xi ∂xj

i,j=1
m∑i,j,k=1∂ 3 (g ◦ Φβ )(x) ≤ ∥g ′′′ ∥∞ + 6β∥g ′′ ∥∞ + 6β 2 ∥g ′′ ∥∞∂xi ∂xj ∂xkが成り立つ
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

24 / 49CCK 理論: Stein の方法

Stein の方法によって |E [g (Φβ (F − y ))] − E [g (Φβ (Z − y ))]| を評価する

補題 5 (Stein の等式)
G を平均 0 の m 次元 Gauss 型確率変数とする. このとき, 有界な偏導関数をもつ任意の C 1 級関数 f : Rm → R に対して
[
]
m
∑
∂f
E [Gi f (G )] =E [Gi Gj ]E(G )(i = 1, . . . , m)∂xj
j=1

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

25 / 49CCK 理論: Slepian 補間 (スマート・パス法)
一般性を失わずに F と Z は互いに独立と仮定してよい
補題 4 より ψ := g ◦ Φβ の 2 階の偏導関数はすべて有界√√
関数 Ψ : (0, 1) → R を Ψ(t) = E [ψ( tF + 1 − tZ − y )] (t ∈ (0, 1))で定義すると, Ψ は微分可能であり, 各 t ∈ (0, 1) について1∑
Ψ (t) =
E
2
m

′

i=1

[(

F
Z
√i − √ i
t
1−t

)

√
∂ψ √
( tF + 1 − tZ − y )∂xi]

F が Gauss 型ならば, Stein の補題より
[ 2
]
m
√
1 ∑
∂ ψ √
′
Ψ (t) ={E [Fi Fj ] − E [Zi Zj ]} E
( tF + 1 − tZ − y )2∂xi ∂xji,j=1

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

26 / 49CCK 理論: 証明の完成
従って, 微分積分学の基本定理と補題 4 より, ある普遍定数 c2 > 0 が存在して,|E [g (Φβ (F − y ))] − E [g (Φβ (Z − y ))]| = |Ψ(1) − Ψ(0)|≤ c2 (ε−2 + ε−1 β)∆≤ 2c2 ε−2 (log m)∆

以上より√
sup |P (F ≤ x) − P (Z ≤ x)| ≤ c1 ε log m + 2c2 ε−2 (log m)∆

x∈Rm

を得る

ε = ∆1/3 log1/6 m として定理 2 の結論を得る
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

27 / 49

Malliavin 解析
上記の議論からわかるように, F が Gauss 型であるという仮定はStein の方法の適用にのみ必要
言い換えれば, F が Stein の方法を有効に適用できるクラスに属してさえいれば, 上の議論はそのまま成立する ⇒ Malliavin 解析!!
Malliavin の部分積分公式は Stein の等式の無限次元版 (Sakamoto andYoshida, 1999, 2004)▶ “Malliavin-Stein method” (Nourdin and Peccati, 2009)▶

∫1∫1
F = 0 us dBs と書ける場合を考える (u: E [ 0 us2 ds] < ∞ なる Rm×r
値発展的可測過程)

補題 6
G ∈ D1,2 (Rm ) ならば, 有界な偏導関数をもつ任意の C 1 級関数
f : Rm → R に対して
[(∫ 1
]
)
m
∑
∂f
i·
E [Fi f (G )] =E
us · Ds Gj ds(G )(i = 1, . . . , m)
∂xj
0
j=1

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

28 / 49

Nualart-Yoshida 補間
極限の混合正規性は問題となるか?
▶

Z = C1/2 ζ, ζ ∼ N(0, 1) ⊥⊥ F, C: m 次半正定値 F-可測ランダム行列▶

F と Z はもはや独立でない！
上の議論をそのまま成立させるには, C が対角行列でないとうまくい
かない (cf. Nourdin et al., 2016)⇒ 問題となる！
▶

この問題は, E [ψ(F − y )] と E [ψ(Z − y )] の補間を周波数領域で行えば解決できる (Nualart and Yoshida, 2017)▶

F と Z の特性関数どうしを補間して, 反転公式を使う

θ ∈ [0, 1] と z ∈ Rm に対して,√z⊤ Czλ(θ; z) := θ −1z · F − (1 − θ2 ),2φ(θ; z) := E [e λ(θ;z) ]とおくと, φ(1; z), φ(0; z) がそれぞれ F , Z の特性関数となる

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

29 / 49

Nualart-Yoshida 補間√∂φ(θ; z) = E [e λ(θ;z) ( −1z · F + θz⊤ Cz)]∂θであり, 右辺第 1 項に補題 6 を適用して整理すると,∂φ(θ; z) = θE e λ(θ;z)
∂θ
−

√

m
∑

(
zi zj

−1

2

E e λ(θ;z)

∫ 1


)
usi· · Ds Fj ds + Cij 

0

i,j=1



1 − θ2√( −1)2m∑(∫ 1
zi zj zk

i,j,k=1

)usi· · Ds Cjk ds  (5)0(Fj , Cjk ∈ D1,2 と仮定)

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

30 / 49

Nualart-Yoshida 補間ψ が十分滑らかならば, Fourier の反転公式と微分積分学の基本定理より
E [ψ(F − y )] − E [ψ(Z − y )]∫√1ψ̂(z)e −1z·y {φ(1; z) − φ(0; z)}dz
=
m
(2π)
Rm
)
∫ (∫ 1√1−1z·y ∂φ=(θ; z)dθ dz
ψ̂(z)e(2π)m Rm∂θ0(ψ̂ は ψ の Fourier 変換)

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

31 / 49

(5) 式を代入したあと再び反転公式を適用して以下を得る:
E [ψ(F − y )] − E [ψ(Z − y )]


(∫ 1
)
∫ 1
m
2
√
∑
∂
ψ
=
θE 
(θF + 1 − θ2 Z )usi· · Ds Fj ds − Cij  dθ
∂xi ∂xj
0
0
i,j=1


(∫ 1
)
∫ 1
m
√
1 − θ2  ∑∂3ψ+E(θF + 1 − θ2 Z )usi· · Ds Cjk ds  dθ2∂xi ∂xj ∂xk
0
0i,j,k=1(6)
従って, 補題 4 より, ある普遍定数 c3 > 0 が存在して,|E [ψ(F − y )] − E [ψ(Z − y )]|[]∫ 1≤ c3 ε−2 (log m)E max
usi· · Ds Fj ds − Cij1≤i,j≤m−3+ c3 ε小池祐太 (東大 MI センター)[

2(log m) E0usi· · Ds Cjk ds

max

1≤i,j≤m]∫ 10

高次元共分散行列の統計推測

2019 年 8 月 7 日

32 / 49CCK 理論: 混合正規極限の場合
まとめると以下を得る:

定理 3
F , C, m はパラメーター n ∈ N に依存して, 以下の条件を満たすとする:(i) F ∈ D1,2 (Rm ), C ∈ D1,2 (Rm×m )[]∫1(ii) (log m)2 E max1≤i,j≤m 0 usi· · Ds Fj ds − Cij → 0 (n → ∞)[]∫1(iii) (log m)7/2 E max1≤i,j≤m 0 usi· · Ds Cjk ds → 0 (n → ∞)(iv) limb↓0 lim supn→∞ P(min1≤j≤m Cjj < b) = 0
このとき, supx∈Rm |P (F ≤ x) − P (Z ≤ x)| → 0 (n → ∞) が成り立つ

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

33 / 49CCK 理論: 混合正規極限の場合

定理 3 に関するコメント

Sn は伊藤積分で書けるので定理 3 は適用可能だが, (3) を得るには不十分⇒ Sn は実際には二重伊藤積分なので, 補題 6 を 2 回適用して評価を精密にする必要がある▶
▶

詳細は K. (2019a) 参照 (Ξn は上と同様の議論で対処できる)
もしくは, (6) 式を Nualart and Yoshida (2017) に従って整理すると,(3) を得るのに十分な結果を得られるかもしれない

定理 3 において, 実際には F は Skorohod 積分であれば十分(anticipative でもよい)▶

特に F のマルチンゲール性は使われていない (これは非効率的?)

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

34 / 49

主結果: 記号
確率変数 ξ と p > 0 に対して ∥ξ∥p := (E [|ξ|p ])1/p
ランダム配列 F = (Fi1 ,...,iq )1≤i1 ,...,iq ≤r と p > 0 に対して

v
u ∑
u r
∥F ∥p,ℓ2 := t

Fi21 ,...,iq

i1 ,...,iq =1

m × l 行列 A に対して |||A|||∞ := max1≤i≤m

p

∑l

ij
j=1 |A |正方行列 A の対角成分の最小値を min diag(A) で表す
以下では B に関する (偏)Malliavin 解析を考える
ある非ランダムな m × d 2 行列 Υn と m × d 2 ランダム行列 X n が存在して Ξn = Υn ◦ X n と書けると仮定 (◦ は成分ごとの積の意) (Ξn の「スパース性」に関する仮定)

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

35 / 49

定理 4 (K. (2019a), Theorem 4.1)
すべての t ∈ [0, 1] について µt ∈ D1,∞ (Rd ), σt ∈ D2,∞ (Rd×r ) および2
X n ∈ D2,∞ (Rm×d ) であると仮定する. また, |||Υn |||∞ ≥ 1, |||Υn |||5∞ = O(nϖ ) と()
max max ∥X ijn ∥p + sup ∥Dt X ijn ∥p,ℓ2 + sup ∥Ds,t X ijn ∥p,ℓ2 = O(1),1≤i≤m 1≤j≤d 2

0≤t≤1

max sup∥µit ∥p + sup ∥Ds µit ∥p,ℓ2(

1≤i≤d 0≤t≤1

max sup1≤i≤d 0≤t≤1

sup
0≤t≤1− n10≤s,t≤1)= O(1),

0≤s,t≤1()∥Σiit ∥p + sup ∥Ds σti· ∥p,ℓ2 +
0≤s,t≤1max

1≤k,l≤d∥Ds,t σui· ∥p,ℓ2= O(1),

0≤s,t,u≤1= O(n−γ )kl
Σkl
t+ 1 − Σt
n

sup

2

がすべての p ∈ [1, ∞) とある ϖ ∈ (0, 12 ), γ ∈ (0, 12 ] について成り立つと仮定する. さらに, 以下を仮定する:

bn Ξ⊤
lim lim sup P(min diag(Ξn Cn ) < b) = 0.
b↓0

n→∞

このとき, ある c > 0 について d + m = O(nc ) ならば, (3) 式が成立する
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

36 / 49

残差スパース性検定への応用
定理 4 より, L ⊂ Λn が与えられたとき, 適当な条件下で統計量(i,j)
max(i,j)∈L |Tn | の分布はある混合 Gauss 型確率ベクトルの最大値の分布で近似できる⇒ 後者の分布はシミュレーション (ブートストラップ法) で評価できる
標準的なブートストラップ法は非エルゴード性の問題によって機能しない (cf. Dovonon et al., 2013; Hounyo, 2017)⇒ Gauss 型 MA(1) ワイルドブートストラップを用いる:
n,∗i, Y j]
[Y\
1 :=

n
∑
∗
(ηh∗ − ηh−1)(Ytih − Ytih−1 )(Ytjh − Ytjh−1 ).
h=1i.i.d.

ここに, ηh∗ ∼ N(0, 0.5) ⊥⊥F▶

本質的に Hounyo (2017) の WBB ブートストラップの特殊ケース

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

37 / 49

残差スパース性検定への応用

n,∗

n

n

n,∗

i , Y d ] [Y\j , Y d ] + [Y\i , Y d ] [Y\j, Y d]T̂ijn,∗ := [Y\
1
1
1
1
n,∗

n

n

n,∗

\
i , Y j ] [Y\d , Y d ] − [Yi , Y j ] [Y\d, Y d]
− [Y\
1
1
1
1
とし,
(i,j)
Tn,∗ :=(i,j)

とおくと, max(i,j)∈L |Tn
ることが証明できる

√ ij
nTn,∗
√
V̂ijn
(i,j)

| の分布を max(i,j)∈L |Tn,∗ | の分布で近似でき

K. (2019a, Proposition 4.2) 参照

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

38 / 49

数値実験
モデル (1) を以下の設定でシミュレーション (cf. Fan et al., 2016):▶

ファクター過程 Y d は Heston モデル:

dYtd = µdt +√vt dBtd ,)√√ (dvt = κ(θ − vt )dt + η vt ρdBtd + 1 − ρ2 dBtd+1 .

パラメーターは µ = 0.05, κ = 3, θ = 0.09, η = 0.3, ρ = −0.6 とし, 初
期値 v0 は定常分布から生成j▶ 残差過程: dR = γ ⊤ dB (j = 1, . . . , d)ttjdB t = (Bt1 , . . . , Bt )⋆ γ1 , . . . , γd は B と独立な d 次元確率ベクトルで, Γ := (γ ⊤ γj )1≤i,j≤d はi
同サイズのブロック 10 個をもつブロック行列. 各ブロックの対角成分
は [0.2, 0.5] 上の一様乱数で生成し, ブロック内相関は定数 ργ とする⋆▶

β 1 , . . . , β d は [0.25, 2.25] 上の一様乱数で生成

d = 100 として, ργ ∈ {0.25, 0.5, 0.75} と変動させ FWER と平均検出力を比較
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

39 / 49

表 1: 有意水準 5%での FWER

ργ
0.25
0.50
0.75

n = 26
0.022
0.023
0.026

n = 39
0.007
0.008
0.011

n = 78
0.003
0.004
0.006

n = 130
0.004
0.005
0.008

n = 195
0.009
0.009
0.014

n = 390
0.018
0.019
0.023

Romano-Wolf 法による多重検定の有意水準 5%での FWER(10,000 回のモ
ンテカルロ実験の結果から計算). ブートストラップ複製は 999 回.

表 2: 有意水準 5%での平均検出力

ργ
0.25
0.50
0.75

n = 26
0.000
0.000
0.004

n = 39
0.000
0.001
0.017

n = 78
0.000
0.037
0.393

n = 130
0.006
0.458
0.977

n = 195
0.048
0.956
1.000

n = 390
0.567
1.000
1.000

Romano-Wolf 法による多重検定の有意水準 5%での平均検出力 (10,000 回のモンテカルロ実験の結果から計算). ブートストラップ複製は 999 回.
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

40 / 49

実データ解析例
期間: 2018 年 3 月▶

この期間を区間 [0, 1] とみなす (オーバーナイトリターンは無視)

15 分足リターン
分析対象: S&P100 指数構成銘柄
ファクター過程: SPDR S&P 500 ETF (SPY)
データソース: ブルームバーグ
図 1–2 に対応する共分散行列のスパース性を検定するために,
Romano-Wolf のステップダウン法を実行する▶

棄却値の計算は 999 回のブートストラップ複製に基づく

Fan et al. (2016) で報告されたブロック対角構造を確認するために,セクターによって銘柄をソートしておく

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

41 / 49

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

図 5: S&P100 指数構成銘柄の実現相関行列 (2018 年 3 月): FWER が 5%水準で有意な成分=1, それ以外=0
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1

15 分足リターンで計算. 長方形はセクターを表す

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

42 / 49

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

AMZN
BKNG
CHTR
CMCSA
DIS
F
FOX
FOXA
GM
HD
LOW
MCD
NKE
SBUX
TGT
TWX
CL
COST
CVS
KHC
KO
MDLZ
MO
PEP
PG
PM
WBA
WMT
COP
CVX
HAL
KMI
OXY
SLB
XOM
AIG
ALL
AXP
BAC
BK
BLK
BRK.B
C
COF
GS
JPM
MET
MS
USB
WFC
ABBV
ABT
AGN
AMGN
BIIB
BMY
CELG
DHR
GILD
JNJ
LLY
MDT
MRK
PFE
UNH
BA
CAT
EMR
FDX
GD
GE
HON
LMT
MMM
RTN
UNP
UPS
UTX
AAPL
ACN
CSCO
FB
GOOG
GOOGL
IBM
INTC
MA
MSFT
NFLX
ORCL
PYPL
QCOM
TXN
V
DWDP
SPG
T
VZ
DUK
EXC
NEE
SO

図 6: マーケットファクターで回帰後の残差過程の実現相関行列 (2018 年 3 月):
FWER が 5%水準で有意な成分=1, それ以外=0
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1

15 分足リターンで計算. 長方形はセクターを表す.
マーケットファクターは SPY で代理

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

43 / 49

まとめ

サンプル数とともに次元が上昇するような高次元の設定下で実現共分散行列の漸近混合正規性を証明した▶

証明は Malliavin 解析と Nualart-Yoshida 補間による
Chernozhukov-Chetverikov-Kato 理論の変形

理論の応用例として, 連続時間 1 ファクターモデルの残差スパース性検定を説明
数値実験と実データ解析例によって理論の適用可能性を示した
論文▶

“Mixed-normal limit theorems for multiple Skorohod integrals inhigh-dimensions, with application to realized covariance”, Electronic
Journal of Statistics, 13 (2019), no.1, 1443–1522.

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

44 / 49

今後の課題より複雑な状況への適用
非同期観測の場合 (Hayashi-Yoshida 推定量は二重伊藤積分で書ける!)
マイクロストラクチャーノイズ
▶ ジャンプ▶
▶

ファクターを PCA 等で構成した場合の残差スパース性検定

FWER ではなく False Discovery Rate (FDR) や False Discovery
Proportion (FDP) のコントロールを考える
高次元精度行列の統計推測への応用
De-biased graphical Lasso を使うと, インプットの共分散行列推定量の高次元漸近混合正規性からほぼ自動的に精度行列推定量の高次元漸近
混合正規性が出てくる (Koike (2019b) 参照)▶ 最適収束レートを達成する精度行列推定法と組み合わせることは可能か？ (例えば Sun and Zhang (2013) の scaled Lasso)▶

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

45 / 49

References I
M. Bibinger and P. A. Mykland. Inference for multi-dimensional high-frequency
data with an application to conditional independence testing. Scand. J. Stat.,
43:1078–1102, 2016.
P. J. Bickel and E. Levina. Regularized estimation of large covariance matrices.
Ann. Statist., 36(1):199–227, 2008a.
P. J. Bickel and E. Levina. Covariance regularization by thresholding. Ann.
Statist., 36(6):2577–2604, 2008b.
C. Brownlees, E. Nualart, and Y. Sun. Realized networks. J. Appl. Econometrics,
33(7):986–1006, 2018.
V. Chernozhukov, D. Chetverikov, and K. Kato. Gaussian approximations andmultiplier bootstrap for maxima of sums of high-dimensional random vectors.
Ann. Statist., 41(6):2786–2819, 2013.
V. Chernozhukov, D. Chetverikov, and K. Kato. Gaussian approximation ofsuprema of empirical processes. Ann. Statist., 42(4):1564–1597, 2014.

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

46 / 49

References II
V. Chernozhukov, D. Chetverikov, and K. Kato. Comparison andanti-concentration bounds for maxima of Gaussian random vectors. Probab.
Theory Related Fields, 162:47–70, 2015.
V. Chernozhukov, D. Chetverikov, and K. Kato. Empirical and multiplierbootstraps for suprema of empirical processes of increasing complexity, andrelated Gaussian couplings. Stochastic Process. Appl., 126:3632–3651, 2016.
V. Chernozhukov, D. Chetverikov, and K. Kato. Central limit theorems andbootstrap in high dimensions. Ann. Probab., 45(4):2309–2353, 2017.
P. Dovonon, S. Gonçalves, and N. Meddahi. Bootstrapping realized multivariate
volatility measures. J. Econometrics, 172:49–65, 2013.
J. Fan, Y. Liao, and M. Mincheva. High-dimensional covariance matrix estimationin approximate factor models. Ann. Statist., 39(6):3320–3356, 2011.
J. Fan, A. Furger, and D. Xiu. Incorporating global industrial classification
standard into portfolio allocation: A simple factor-based large covariancematrix estimator with high-frequency data. J. Bus. Econom. Statist., 34(4):489–503, 2016.
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

47 / 49

References III
U. Hounyo. Bootstrapping integrated covariance matrix estimators in noisy
jump-diﬀusion models with non-synchronous trading. J. Econometrics, 197:
130–152, 2017.
J. Jacod and P. Protter. Discretization of Processes. Springer, 2012.
Y. Koike. Mixed-normal limit theorems for multiple Skorohod integrals inhigh-dimensions, with application to realized covariance. Electron. J. Stat., 13(1):1443–1522, 2019a.
Y. Koike. De-biased graphical Lasso for high-frequency data. Working paper.
arXiv:1905.01494, 2019b.
I. Nourdin and G. Peccati. Stein’s method on Wiener chaos. Probab. Theory
Related Fields, 145:75–118, 2009.
I. Nourdin, D. Nualart, and G. Peccati. Quantitative stable limit theorems on the
Wiener space. Ann. Probab., 44(1):1–41, 2016.
D. Nualart and N. Yoshida. Asymptotic expansion of Skorohod integrals. Working
paper. Available at arXiv: https://arxiv.org/abs/1801.00120, 2017.

小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

48 / 49

References IV
P. Ravikumar, M. J. Wainwright, G. Raskutti, and B. Yu. High-dimensional
covariance estimation by minimizing ℓ1 -penalized log-determinant divergence.
Electron. J. Stat., 5:935–980, 2011.
J. P. Romano and M. Wolf. Exact and approximate stepdown methods for
multiple hypothesis testing. J. Amer. Statist. Assoc., 100(469):94–108, 2005.
Y. Sakamoto and N. Yoshida. Malliavin calculus and statistical asymptotic theory(in Japanese). Proc. Inst. Statist. Math., 47(1):175–199, 1999.
Y. Sakamoto and N. Yoshida. Asymptotic expansion formulas for functionals of
ϵ-Markov processes with a mixing property. Ann. Inst. Statist. Math., 56(3):545–597, 2004.
T. Sun and C.-H. Zhang. Sparse matrix inversion with scaled Lasso. J. Mach.
Learn. Res., 14:3385–3418, 2013.
Y. Wang and J. Zou. Vast volatility matrix estimation for high-frequency financial
data. Ann. Statist., 38(2):943–978, 2010.
M. Yuan and Y. Lin. Model selection and estimation in the Gaussian graphicalmodel. Biometrika, 94(1):19–35, 2007.
小池祐太 (東大 MI センター)

高次元共分散行列の統計推測

2019 年 8 月 7 日

49 / 49

Detecting the number of factors of quadratic variationin the presence of microstructure noise :High-Dimension Case
Naoto Kunitomo
Meiji University

August 7, 2019
This is a joint work with Daisuke Kurisu.

Kunitomo

Rank test for QV

2019/8/7

1 / 40Outline

Introduction
Estimation of QV
Representing the problem as a Random Field model
Characteristic Roots in High-Dimension
A Mathematical Interpretation and Discussions
Simulation Results
Further Problems

Kunitomo

Rank test for QV

2019/8/7

2 / 40Introduction
One of important observations on the asset price movements has
been that many financial asset prices move in similar way in theirtrends, volatilities and jumps.
There is a question how to cope with many asset prices when thenumber of factors of volatilities or quadratic variation of asset prices
is less than the dimension of observed prices. Kunitomo and Kurisu(2019) is developping a statistical method for detecting the numberof factors.
The problem is related to the reduced rank regression problem, which
has been well-known in statistical multivariate analysis (Anderson(1984, AS) and Anderson (2003)).
There is an important question when the dimension of financial asset
prices is large (high-dimension) because we have a large number ofassets traded in financial markets. This is an extension of Kunitomo
and Kurisu (2019).
Kunitomo

Rank test for QV

2019/8/7

3 / 40Introduction

Figure: Intra-day movements of 5 Stock Prices at TSE(afternoon session of January 25, 2016).

We set the same values for the starting process since we focus on thevolatility structure (of quadratic variation) of asset prices.
Kunitomo

Rank test for QV

2019/8/7

4 / 40Introduction
The new feature in our formulation is the fact that we are dealing
with the continuous stochastic process as the hidden process with
diﬀusions and jump processes while we have discrete observations
with measurement errors.
We develop a new way to determine the number of factors ofquadratic covariation or the integrated volatility of asset prices basedon the separating information maximum likelihood (SIML) method
when the dimension of observed vectors is large.
The SIML method is originally developed by Kunitomo, Sato andKurisu (2018).
Kunitomo, N. and Kurisu, D. (2019). Detecting factors of quadratic
variation in the presence of microstructure noise. MIMS-RBP SDS-10.
Kunitomo

Rank test for QV

2019/8/7

5 / 40Estimation of QV
We adopt the the construction of the whole filtered probability space(Ω, F, (Ft )t∈[0,1] , P) where both the hidden process X and the noise V.
Let Y(tin ) = (Yj (tin )) be the (p-dimensional) observed (log-)prices atti ∈ [0, 1] and i = 1, · · · , n, which satisfyies
Y(tin ) = X(tin ) + v(tin ) (i = 1, · · · , n) ,X(tin ) (= (Xj (tin )): the p × 1 hidden stochastic vector process.
v(tin ) (= (vj (tin ))): a sequence of (mutually) independent market′
microstructure noises with E[v(tin )] = 0 and E[v(tin )v(tin ) ] = Σv
(a positive definite matrix).
We also assume that v is independent of X.

Kunitomo

Rank test for QV

2019/8/7

6 / 40Estimation of QV
We also assume that X is a p-dimensional continuous-time stochastic
process which is given by∫ t∫ t
σ(s)dW(s)X(t) = X(0) +b(s)ds +
0
0
∫ t∫
∫ t∫
+
∆(s, x)(µ − ν)(ds, dx) +∆(s, x)µ(ds, dx),
0

∥x∥<1

0

∥x∥≥1

b(s): the p-dimensional adapted drift process, σ(s): the p × q1
instantaneous predictable volatility process, c(s) = σ(s)σ ′ (s),W(s) = (Wj (s)): the q1 × 1 standard Brownian motions.
∆(ω, s, x) is a Rp -valued predictable function on Ω × [0, ∞) × Rq2 ,µ(·) is a Poisson random measure on [0, ∞) × Rq2 andν(ds, dx) = ds ⊗ λ(dx) is the predictable compensator of µ with a
σ-finite measure λ on (Rp , B p ).
Kunitomo

Rank test for QV

2019/8/7

7 / 40Estimation of QV
The jump terms are denoted as ∆Xs = (∆Xj (s))(∆Xj (s) = Xj (s) − Xj (s−) and Xj (s−) = limu↑s Xj (u) at any s ∈ [0, 1].
Here ∥ · ∥ is the Euclidean norm on Rp .
Assumption 2.1(a) The path t 7→ b(t, ω) is locally bounded.
∫ t+u(b) The process σ is continuous and t ||σ(s)||ds > 0 a.s. for allt, u > 0.
(c) For some positive K , the measurable functions b(s, x), σ(s, x), and∆(s, x) (0 ≤ s ≤ 1) satisfy∫∥b(s, x)∥2 + ∥σ(s, x)∥2 +∥∆(s, u, x)∥2 λ(du) ≤ K [ 1 + ∥x∥2 ] , x ∈ Rp ,U0

where U0 = {∥u∥ < 1}.

(d) The noise terms v(tin )(= (vj (tin ))′ , 1 ≤ i ≤ n, 1 ≤ j ≤ p) are asequence of i.i.d. random variables with E[vj (tin )] = 0,E[v(tin )v′ (tin )] = Σv (p.d.), E[vj4 (tin )] < ∞, v ⊥⊥ X.
Kunitomo

Rank test for QV

2019/8/7

8 / 40Estimation of QV
The fundamental quantity for the continuous-time Itô semimartingale withp ≥ 1 is the quadratic variation (QV) matrix, which is given by
∫ 1
∑
′
(x)
Σx =
c(s)ds +(∆Xs )(∆Xs ) = (σgh ) .
0

0≤s≤1

The SIML estimator of Σ̂x (Kunitomo, Sato and Kurisu (2018)) for the
integrated volatility is defined byn
1 ∑′zk zk ,
Σ̂x =
mn

m

k=1

where zk = (zgk ), j = 1, · · · , p; k = 1, · · · , mn , mn = O(nα ), 0 < α < 1,′which are constructed by the transformation from Yn = (yk ) (n × p) to′
Zn (= (zk )) by()
Zn = Kn Yn − Ȳ0
Kunitomo

Rank test for QV

2019/8/7

9 / 40Estimation of QV
where Kn =√nPn C−1n ,
10 ··· 00
 −1 10 ··· 0 
−1
0 ··· Cn = 
, 0 −1 1 0
0 −1 10000 −1 1
n×n

and
√
(n)
(n)
Pn = (pjk ) , pjk =

[
]
2
2π
1
1
cos
(k − )(j − ) .
2n + 1
2
2
n + 12Intuitively, the SIML estimator Σ̂x is an estimator of the spectral density
matrix of X at the origin, i.e. Σx = fX (0). (Indeed, this is true if we apply
the SIML method to state space models.)
Kunitomo

Rank test for QV

2019/8/7

10 / 40

Representing the problem as a Random Field
We consider the inference probelm when the dimension p is large. For
j = 1, · · · , p and at s = j/p, we set Xs∗ (t) = Xj (t) (0 ≤ t ≤ 1) , where
Xj (t) is the j-th component of X(t). We construct the random field Xs (t)on s ∈ [0, 1], t ∈ [0, 1] byj −1Xs∗ (t) = Xj−1 (t) + p[s −][Xj (t) − Xj−1 (t)] ,pj
for j−1p ≤ s ≤ p and j = 1, · · · , p. Then we have the continuous timeprocess∫ t∫ t
Xs∗ (t) = X0∗ (t) +bs∗ (u)du +σ ∗s (u)dBu
0
0
∫ ∫
∗
+
∆s (x, u)(µ(u) − ν(u))(du, dx, s)u ∥x∥<1∫ ∫+∆∗s (x, u)µ(du, dx, s).
uKunitomo∥x∥≥1

Rank test for QV

2019/8/7

11 / 40

We write the quadratic variation as∫ 1∑′
σx (s, t) =σ ∗s (u)σ ∗t (u)du +∆Xs (u)∆Xt (u) ,
0

0≤u≤1

which is a non-negative definite (random or non-random) matrix. Then, byusing Mercer Theorem,∑
σx (s, t) =µi βi (s)βi (t) ,i≥1where 0 ≤ µ1 ≤ µ2 ≤ · · · are the eigenvalues and βi (u) are the associatedeigenfunction.
In the same way, we construct the noise process
Vs∗ (t) (0 ≤ s ≤ 1, 0 ≤ t ≤ 1) and Ys∗ (t) (0 ≤ s ≤ 1, 0 ≤ t ≤ 1) byj −1Vs∗ (t) = Vj−1 (t) + p[s −][Vj (t) − Vj−1 (t)]
p
and
j −1
Ys∗ (t) = Yj−1 (t) + p[s −][Yj (t) − Yj−1 (t)]pKunitomo
Rank test for QV
for s ∈ [0,1] and t ∈ [0, 1].
2019/8/7
12 / 40

Then we haveYs∗ (t) = Xs∗ (t) + Vs∗ (t) (0 ≤ s ≤ 1, 0 ≤ t ≤ 1).

Kunitomo

Rank test for QV

2019/8/7

13 / 40

Characteristic Roots in High-Dimension
Let the random field gm (s, t) (0 ≤ s, t ≤ 1) from gm ( pi , pj ) = [Gm ]ij ,σx ( pi , pj ) = [Σx ]ij and σv ( pi , pj ) = [Σv ]ij where
Gm = ((Gm )ij )) (i, j = 1, · · · , p), Σx = ((Σx )ij )) (i, j = 1, · · · , p) andΣv = ((Σv )ij )) (i, j = 1, · · · , p). We denote
i j
i j
i j
σm ( , ) = σx ( , ) + am σv ( , ) .
p p
p p
p p
Let

Dm (λ) =gm ( p1 , p1 ) − λσv ( p1 , p1 )gm ( p2 , p1 ) − λσv ( p2 , p1 )..
.
gm ( pp , p1 ) − λσv ( pp , p1 )Kunitomo

gm ( p1 , p2 ) − λσv ( p1 , p2 )gm ( p2 , p2 ) − λσv ( p2 , p2 )..
.
gm ( pp , p2 ) − λσv ( pp , p2 )Rank test for QV

···
···
..
.
···

gm ( p1 , pp ) − λσv ( p1 , pp )gm ( p2 , pp ) − λσv ( p2 , pp )..
.
gm ( pp , pp ) − λσv ( pp , pp )2019/8/7.

14 / 40

Then the characteristic roots when we take H = Σv = σv2 Ip as the leadingcase. The underlying characteristic equation can be written as thesolutions of|Dm (λ)| = 0 .
We transform the present problem into the stochastic integral equation andj−1jirandom field. For i−1
p ≤ t ≤ p and p ≤ s ≤ p (i, j = 1, · · · , p), definegm (t, s)=i −1 j −1,)ppi j −1i −1 j −1+cp (t)[gm ( ,) − gm (
,
)]
p
p
p
p
i −1 ji −1 j −1, ) − gm (,)]+cp (s)[gm (
p
p
p
p
i j
i j −1i −1 ji −1 j −1+cp (t)cp (s)[gm ( , ) − gm ( ,) − gm (, ) + gm (
,
)] ,
p p
p
p
p
p
p
p
gm (

where cp (t) = p[t − (i − 1)/p] and cp (s) = p[s − (j − 1)/p].
Kunitomo

Rank test for QV

2019/8/7

15 / 40

Under some conditions such as the boundedness of fourth cordermoments. we have
[(∫ ∫
)2 ]
1
1
gm (t, s)dtds
E< +∞ .
00Hence we can use the standard arguments on L2 . For a fixed p, we write
i j|(gm ( , ) − am σv2 δ(i, j)) − (λ − am )σv2 Ip | = 0 .
p p
If we take n, m → ∞, theni j|plimm→∞ (gm ( , ) − am σv2 δ(i, j)) − (λ − am )σv2 Ip | = 0
p p
impliespλi − am −→ 0 (i = 1, · · · , rx ) .
For the case when p is large, we may set p = pm . In order to develop theasymptotic theory, we impose a condition that the number of eigenvalues
is less than m. Also we assume that the volatility and co-volatility
functions are well-defined.
Kunitomo

Rank test for QV

2019/8/7

16 / 40

Condition II : pmm −→ 0 as n, m −→ +∞.
It may be possible to extend the following results when pm /m → c andc > 0 as m → ∞.
Miyahara (1982) has pointed out the fundamental diﬃculty on L2
arguments in the present situation (See Section 4). The rank condition
should be replaced by the following condition of continuous integral
equation :Condition III : There exist eigenfunctions βk (t) (0 ≤ t ≤ 1, k ≥ 1) such
that
∫ 1
k
σx (t, s)βk (t)dt = 0 (at s = , k = 1, · · · , rx ) .
p0Because gm (t, s) is a symmetric non-negative definite operator on[0, 1] × [0, 1], the set of eigenvalues is at most denumerable andeigenfunctions corresponding to the diﬀerent eigenvalues are orthogonal by
the (standard) integral equation theory.

Kunitomo

Rank test for QV

2019/8/7

17 / 40

We construct σm (t, s), σx (t, s), σv (t, s) and bk (t) from Σm , Σx , Σv andbij for gm (t, s). We take bjk (j, k = 1, · · · , p) such that and construct
bk (t) from bjk for 0 ≤ t ≤ 1.

Kunitomo

Rank test for QV

2019/8/7

18 / 40

We consider the stochastic integral equation as
∫ 1
∫ 1
k
gm (s, t)bk (t)dt = λk
σv (t, s)bk (t)dt (at s = , k = 1, · · · , p) ,p00where λk (0 ≤ λ1 ≤ λ2 ≤ · · · , 1 ≤ k ≤ ∞). This equation corresponds tothe continuous version. By using σm (t, s) = σx (t, s) + am σv (t, s), we
re-write∫ 1∫ 1[gm (t, s) − am σv (t, s)] bk (t)dt = [λk − am ]σv (t, s)bk (t)dt .
00

For a fixed p, we take the probability limit as m −→ ∞.
Since we have the rank condition (R), E[gm (t, s)] = σx (t, s) + am σv (t, s),∫ 1∫ 1
plimm [gm (t, s)−am σv (t, s)]βk (t)dt = plimm [λk −am ]σv (t, s)βk (t)dt00

we find that as m −→ ∞,pλj − am −→ 0 (j = 1, · · · , rx )and 0 < rx = p − qx .
Kunitomo

Rank test for QV

2019/8/7

19 / 40

We summarize our heuristic arguments.
Theorem 3.1 : Assume some regularity conditions and takem = mn = [nα ] (0 < α < 1), ln = [nβ ] with 0 < α < β < 1. Also assume
that the market micro-structure noises are mutually independent Gaussian
random variables. For a fiexd p, as m → ∞,
[
]
rx
m 1 ∑
w
[ (λi − am )]2 −→ χ2 (1)am (2) 2rx
i=1

where rx = p − qx and 1 ≤ rx < ∞.
Since Theorem 3.1 does not depend on p and rx , we obtain the following
result, which has some important meanings in applications.
Corollary 3.2 : In Theorem 3.1 as p → ∞,
[
]
rx
m 1 ∑
w
2
[ (λi − am )] −→ χ2 (1)am (2) 2rx
i=1

where rx = p − qx and 1 ≤ rx ≤ ∞.
Kunitomo

Rank test for QV

2019/8/7

20 / 40

A Mathematical Interpretation and Discussions
There is a fundamental problem of our formulation and derivations in theprevious sections. We have treated the high-dimension case when(1)(2)
p −→ ∞ and qx = qx + qx −→ ∞. However, when p = +∞ there canbe an infinite number of Brownian motions as the limit, and it is not clearthat we have the proper mathematical meaning of our formulation.
Fortunately, however, there is an earlier work on the cylindrical Brownian
motion on a Banach space by Miyahara (1982) and we can use hisarguments. We are reviewing the first part of Miyahara (1982). For the
probability space (Ω, F, P) and an increasing family of σ−field Ft (t ≥ 0),
let H be a real separable Hilbert space with norm denoted by ∥ · ∥. A′
mapping Bt (h, ω) : [0, ∞) × H × Ω → R is called a cylindrical Brownian
motion (c.s.m.) on H if if satisfies the following conditions. (i)B0 (h, ·) = 0 and Bt (h, ·) is F−adapted, (ii) For any h ∈ H, h ̸= 0,
Bt (h.·)/∥h∥ is a one dimensional Brownian motion, (iii) For any t ∈ [0, ∞)and α, β ∈ R1 and h, k ∈ H, Bt (αh + βk) = αBt (h) + βB(k) (P − a.s.).
Kunitomo

Rank test for QV

2019/8/7

21 / 40

′

Let V be a Banach space which is a dense set of H, and let V be the dual′
Banach space of V. Then we have V ⊂ H ⊂ V and we denote the′canonical bilinar norm on V × V by <, >V ×V ′ or simply by <, >. A′′
V − valued process B̃t is called an H-Brownian motion on V (H-B.M.) if′
it satisfies (i) B̃· (ω) : R+ → V is almost surely continuous and B̃0 (ω) = 0,(ii) For each y ∈ V, y ̸= 0, < y , B̃t > /∥y ∥ is a one-dimensional Brownian
motion.
Then Miyahara (1982) gives the following fundamental result for us.
′

Theorem M (Miyahara 1982) : Let B̃t be a H−B.M. on V . Then there
exists a unique cylindrical Brownian motion Bt on H which satisfies thefollowing equality for any y ∈ V and t ∈ R+
< y , B̃t >= Bt (y ) (P − a.s..)

Kunitomo

Rank test for QV

2019/8/7

22 / 40

We can interpret any real-valued Brownian motions as the bilinear norm of
the cylindrical Brownian motions on H. (We have used such Brownian
motions in (2.2).) Miyahara (1982) has also shown that even when we
have p = ∞, we have the corresponding Banach-valued Brownian motions
and the solutions of stochastic diﬀerential equations.
We denote Xs∗ (t) = (Xj (t)) at s = j/p (j = 1, · · · , p) with
E[∥X (t)∥2 ] < ∞.
In a similar framework a Corollary of Araujo (1978) gives the following
CLT for a sequence of Banach-valued random variables in H1 .
Let Wu (s) = (Wi/n (j/p)) at u = i/n, s = j/p , i = 1, · · · , n; j = 1, · · · , p
be i.i.d. random vectors with E[Wu (j/p)] = 0 andW[Uu (j/p)2 ] = σu2 (> 0).
Proposition 4.1 : Let Wu (s) are i.i,d. random variables at
u = i/n (i = 1, · · · , n) andj −1
Wu (s) = Wu ((j − 1)/p) + [s −][Wu (i/p) − Wu ((i − 1)/p)] .
pAssume2+δ
Kunitomo
Rank (i/p)|
2019/8/7
23 / 40
sup E[|Wu test for QV ] < ∞We construct the random field Wu (s, t) (0 ≤ u ≤ 1, 0 ≤ s, t ≤ 1) from thediscrete process Uu (i/p.j/p) for 0 ≤ u ≤ 1; i, j = 1, ·, p. It can be doneik−1ksuch as that for j−1
p ≤ r ≤ p and p ≤ s ≤ p (j, k = 1, · · · , p). Define
Wu (s, t)=j −1 k −1,)pp
j k −1j −1 k −1+cp (s)[W ( ,) − W(,)]
ppnp
j −1 kj −1 k −1+cp (t)[W (, ) − W(
,
)]
n
p
p
p
j k
j k −1j −1 kj −1 k −1+cp (s)cp (t)[W ( , ) − W ( ,) − W(, ) + W(
,
)] ,
p p
p
p
p
p
p
p
W(

where cp (r ) = p[s − (i − 1)/p] and cp (t) = p[s − (j − 1)/p].

Kunitomo

Rank test for QV

2019/8/7

24 / 40

Then a Corollary of Araujo (1978) implies the following CLT for asequence of Banach-valued random variables in H2 , which would be usefulfor our analysis.
Proposition 4.2 : Define Wu (s, t) be i.i,d. random variables at
u = i/n (i = 1, · · · , n) as above. Assume
sup E[|Wu (s, t)|2+δ ] < ∞
1≤i≤n

√ ∑
for some δ > 0. Then (1/ n) ni=1 Wi/n (s, t) converges weakly to a
Gaussian random variable.

Kunitomo

Rank test for QV

2019/8/7

25 / 40

Simulation Results′

Let X = (X(t))t≥0 (= (X (1) (t), · · · , X (p) (t)) )t≥0 be the vector of Itôsemimartingale satisfying
(j)

(j)

(j)

(j)

dX (j) (t) = σt dBt , j = 1, · · · , q1,xdX (3) (t) = Zt dNt

j = 1, · · · , q2,x ,where B = (B (1) , · · · , B (q1,x ) ) is the two dimensional (standard) Brownian
motion vector, N = (N (1) , · · · , N (q2,x ) is the Poisson process with intensity
λj (j = 1, · · · , q2,x ) as 10. We assume that N is independent of W
(qx = q1,x + q2,x ) and Z = (Zt )t≥0 is the jump sizes with Zt ∼ N(0, 5−2 ).
(See Cont and Tankov (2004) for the generation of jump processes.) For
the volatility process σ of the diﬀusion part, we set
(j)

(j)

(j)

(j)

d(σt )2 = aj (µj − (σt )2 )dt + κj σt dWtσ ,j = 1, · · · , q1,x ,where σ (1) and σ (2) are independent.
Kunitomo

Rank test for QV

2019/8/7

26 / 40

When q1,x = 2 and q1,x = 1, we set a1 = 2, a2 = 3, µ1 = 0.8, µ2 = 0.7,(j)(j)
κ1 = κ2 = 0.5, E[dBtσ dBt ] = ρj dt, ρ1 = ρ2 = −0.5. (We also haveq1,x = 5 and q1,x = 5,)0 In our simulation, we consider the following twomodels :′Model 1 : Y(t) = Γ1 (X (1) (t), X (2) (t)) + v(t)and′Model 2 : Y(t) = Γ1 (X (1) (t), X (2) (t)) + Γ2 X (3) (t) + v(t).
Here we denote the coeﬃcients matrices (p × q1,x and p × q2,x ,
(1)
(p) ′
(1)
(p) ′
respectively) as Γ1 = (γ 1 , . . . , γ 1 ) , Γ2 = (γ 2 , . . . , γ 2 ) and they are(j)(j)
sampled as γ 1 ∼ [0.25, 1.75] and γ 2 ∼ [0.25, 1.75], j = 1, . . . , p.

Kunitomo

Rank test for QV

2019/8/7

27 / 40

The observation vectors are′Y(i/n) = (Y (1) (i/n), . . . , Y (p) (i/n)) , i = 1, . . . , n
and we set ∆ = ∆n = 1/n. As the market microstructure noise vectors,′
we set vi/n = (v (1) (i/n), . . . , v (p) (i/n)) , and use independent Gaussian
noises for each component, that is,′(v (1) (i/n), . . . , v (p) (i/n)) ∼ i.i.d.Np (0, cIp ) (i = 1, . . . , n) ,with a pre-specified value c. In all simulations, we set p = 100.

Kunitomo

Rank test for QV

2019/8/7

28 / 40

Simulation Results
Let N be the number of Monte Carlo iterations. We plotted the meanvalue of the eigenvalues of the SIML estimator for the quadratic variationin Figures . To compute the SIML estimators, we set mn /n = 0.06 andln /n = 0.09.
n
1
1 ∑
′
zj zj , Σ̂v =
mn
ln

m

Σ̂x =

j=1

where akn = 4n sin2

[ (
π
2

2k−1
2n+1

n
∑

′

−1
ajn
zj zj ,

j=n−ln +1)].

In the following Figures, we set 0 ≤ λ̂1 ≤ . . . ≤ λ̂p are eigenvalues of−1Σ̂v Σ̂x .

Kunitomo

Rank test for QV

2019/8/7

29 / 40

In Model-1 and Molde-2 we have p=100 dimensions observation vectors(p = 100). Model-1 has two factor of diﬀusion type (qx = 2, rx = 98)
while Model-2 has two diﬀusion type factors and one jumps factor(qx = 3, rx = 97). Figures show that the estimated characteristic roots
reflect the true rank of hidden stochastic process. They show the
distributions of the test statistic we are developed in this paper. It seems
that our method of evaluating the rank condition of hidden volatility
factors based on the characteristic roots and the SIML estimation detects
the number of factors properly in these two numerical simulations.

Kunitomo

Rank test for QV

2019/8/7

30 / 40

Figure 5.1 : Mean of estimated log characteristic roots (log eigen values)of Model 1 when c = 10−6 (left) and c = 10−8 (right).
We set ∆ = 1/20000, mn /n = 0.06 and ln /n = 0.09. The number ofMonte Carlo iteration is 300.
Kunitomo

Rank test for QV

2019/8/7

31 / 40

Figure 5.2 : Mean of estimated log characteristic roots (log eigen values)of Model 2 when c = 10−6 (left) and c = 10−8 (right).
We set ∆ = 1/20000, mn /n = 0.06 and ln /n = 0.09. The number ofMonte Carlo iteration is 300.
Kunitomo

Rank test for QV

2019/8/7

32 / 40

Figure 5.3 : Empirical distributions of test statistic Tn of Model 1 whenrx = 97(left), rx = 98(center) and rx = 99(right).
We set ∆ = 1/20000, c = 10−8 , mn /n = 0.06, and ln /n = 0.09.
The number of Monte Carlo iteration is 300. The red line is the density of
the chi square distribution with 1 degree of freedom.
Kunitomo

Rank test for QV

2019/8/7

33 / 40

Figure 5.4 : Empirical distributions of test statistic Tn of Model 2 whenrx = 96(left), rx = 97(center) and rx = 98(right).
We set ∆ = 1/20000, c = 10−8 , mn /n = 0.06, and ln /n = 0.09.
The number of Monte Carlo iteration is 300. The red line is the density of
the chi square distribution with 1 degree of freedom.
Kunitomo

Rank test for QV

2019/8/7

34 / 40

Figure 5.5 : Mean of estimated log characteristic roots (log eigen values)of Model 1 when c = 10−6 (left) and c = 10−8 (right).
We set ∆ = 1/20000, mn /n = 0.06 and ln /n = 0.09.
The number of Monte Carlo iteration is 300.
Kunitomo

Rank test for QV

2019/8/7

35 / 40

Figure 5.6 : Mean of estimated log characteristic roots (log eigen values)of Model 2 when c = 10−6 (left) and c = 10−8 (right).
We set ∆ = 1/20000, mn /n = 0.06 and ln /n = 0.09.
The number of Monte Carlo iteration is 300.
Kunitomo

Rank test for QV

2019/8/7

36 / 40

Figure 5.7 : Empirical distributions of test statistic Tn of Model 1 whenrx = 89(left), rx = 90(center) and rx = 91(right). We set ∆ = 1/20000,
c = 10−8 , mn /n = 0.06, and ln /n = 0.09.
The number of Monte Carlo iteration is 300.
The red line is the density of the chi square distribution with 1 degree ofKunitomo
Rank
test for QV
2019/8/7
37 / 40
freedom.

Further Problems
We can develo a detecting procedure for the number of hidden factors
by using the SIML method when the true hidden stochastic process isa class of Itô semimartingales, there can be market microstructurenoises and the dimension is large at the same time.
Our procedure is partly based on the Banach-valued random variables
and CLT(central limit theorems).
We have derived the asymptotic distributions of characteristic roots.
From our limited simulations and an empirical application, ourapproach works well in practical situations.
The results in Kunitomo and Kurisu (2019) can be generalized when
the dimension is high. (Kunitomo, N. and Kurisu, D. (2019).
Detecting factors of quadratic variation in the presence of
microstructure noise. MIMS-RBP SDS-10.)
There are several related problems of Ito semi-martingales withmarket microstructure noises when the dimension is large.
Kunitomo

Rank test for QV

2019/8/7

38 / 40References
Anderson, T.W. (1984), Estimating Linear Statistical Relationships.
Ann. Statist. 12 1-45.
Christensen, K., M. Podolskij and M. Vetter (2013), On Covariation
Estimation for Multivariate Continuous Itô Semimartingales with
Noise in Non-Synchronous Observation Schemes. J. Multivariate
Anal. 120 59-84.
Kunitomo, N., Sato and D. Kurisu (2018), Separating Information
Maximum Likelihood Estimation for High Frequency Financial Data.
Springer.
Araujo, A. P. (1978), ”On the Central Limit Theorem in BanachSpaces,” Journal of Multivariate Analysis, 8, 598-613.
Ito, K. and Nishio, M. (1968), ”On the Convergence of Sum of
Independent Banach Space Valued Random Variables,” Osaka Journalof Mathematics, 35-48.
Miyahara, Yoshio (1982), Stochastic Evolution Equations and White
Noise Analysis, Carleton Mathematical Lecture Notes No.42.
Kunitomo

Rank test for QV

2019/8/7

39 / 40

Mandrekar and Büdinger (2015), Stochastic Integration in BanachSpace, Springer.
Malliavin, P. and M. E. Mancino (2009), ”A Fourier Transform
Method for Nonparametric Estimation of Multivariate Volatility,” The
Annals of Statistics, 37-4, 1993-2010.

Kunitomo

Rank test for QV

2019/8/7

40 / 40

Estimation of risk aversion for Japanese stock market
using implied moments

Kosuke Oya
Graduate School of Economics & MMDS, Osaka Univsersity「データサイエンス・福島キャンプ 2019」2019/08/07, 福島大学

1

OutlineWhat is risk neutral probablity density?
Two probability measures P and Q
Volatility index, variance risk premium, volatility spread
Moment under P
Sample moment & realized moment
Estimation for relative risk aversion
Summary and future studies

2

What is risk neutral probablity
density?

Risk neutral valuation• S: a finite set of state• πP (s)：probability of state s• pc(s): price of contingent claim that pays !1 in state s• x(s): asset whose payoff for state sp(x)：price of asset x#pc(s)p(x) =
pc(s) x(s) =πP (s)
x(s) = EP [ m(s) x(s) ]πP (s)s∈S
s∈S$!!%&!!'
!

"

!

f

f

R : risk free rate, 1/R =

!
s∈S

!"

m(s)

pc(s) =!
s∈S

πP (s)m(s) = EP [m(s)]#m(s)πP (s)11 !
1p(x) = fx(s) = fπQ (s)x(s) = f EQ [ x(s) ]R s∈S EP [m(s)]R s∈SR$!!!!!!!!%&!!!!!!!!'πQ (s)

3

What is risk neutral pdf?
• S(t): stock index price at time t,• R(t) = ln S(t + 1)/S(t): one period return,• P: physical probability measure, Q: risk neutral prob. measure,• p(r): (physical) prob. density func. of R(t), q(r): risk neutral pdf,
Bakshi et al. (2003)• power utility with relative risk aversion γ

suppl.

• pricing kernel m(r) = e−γrq(r) = (m(r) p(r)m(r)p(r)dr=(e−γr p(r)e−γr p(r)dr(1)• risk neutral γ = 0 ⇔ q(r) = p(r)• risk averse γ > 0 ( q(r) = m(r)p(r) for simplicity ),q(r) > p(r) for r < 0, and vice versa
4

Italian Government Bond Futures: Fig.6 in Fornari and Mele (2001)
This is not example of q(r) = e−γ r p(r)5Statistical devices
Estimation for risk neutral probablity density• under the assumption of ( specific parametric model )specific distribution p(r) and functional form m(r)• Aı̈t-Sahalia and Lo (1998)：Nonparametric estimation
moment (provides useful & enough information)• 1st order moment: risk neutral valuation• 2nd order moment: (variance, volatility) risk reference•)

Model free implied variance : volatility index ( VIX, VXJ, VI )• “implied” means ... impled from observed option priceHereafter, k th order moments under P and Q
(k)

(k)

µP = EP [R(t)k ] and µQ = EQ [R(t)k ], respectively
6

Implied moment：Volatility index• VIX は Cboe Global Markets が開発した index で，S&P500 の将来の収益率の変動の大きさの期待値として算出されている• 我が国では大阪大学 CSFI（現：MMDS）が 2008 年 7 月に学術研究目的として VXJ (Volatility Index Japan) を公表，2010 年 11 月に日
本経済新聞社が日経平均 VI の公表開始し，2012 年 2 月に大阪取引所において日経平均 VI 先物の取引が開始• vol. index は株価指数との間に負の相関をもつことから恐怖指数（fear gauge）としても知られる

7

Nikkei average

25000

20000

15000

10000

Volatility Index Japan

2010/01

2012/01

2014/01

2016/012018/01

2016/012018/01

60

40

20

2010/01

2012/01

2014/01

出所： 大阪大学 mmds

8

Implied moment：Volatility index• MFIV は，put と call option の価格 P(K) と C(K) に関する項を権利行
使価格 K に関して積分する以下の式から(1)

EQ [(R(t) − µQ )2 ] = 2

* F
0

P(K)
dK + 2
K2

* ∞
F

C(K)
dK
K2

• 実際には K は離散的で，P(K) と C(K) も対応する K のすべての領域
で観測できるわけではないため近似計算が必要• この積分を離散近似によって算出しているのが VIX• 大阪大学が公開している VXJ はこの積分を Gatheral (2006) の方法
で適切に計算したもの (Fukasawa et al., 2011)• Carr and Madan (2001) や Fukasawa (2012) は 2 階連続微分可能なS(t) の関数の Q のもとでの期待値を導出(4)• µ(3)，µQ ，そして歪度と尖度も option 価格を使って同様にQ

9

Implied moment ：分散リスク・プレミアム（VRP）+ (1) ,2(2)
それぞれの測度のもとでの分散：σ2P = µP − µP ,+ (1) ,2(2)
σ2Q = µQ − µQ

分散リスク・プレミアム (VRP)： VRP = σ2Q − σ2P

VRP を金融市場に関連する予測変数とする研究• Bollerslev et al. (2009)：市場収益率に対する予測力
long-run risk model の拡張（追加的な不確実性の源泉：vol. の変動の risk）• Bollerslev et al. (2014)：
US を含む主要 8 カ国の VRP と共通する global VRP• 大屋 (2011)：景気指標に対する予測力• Ubukata and Watanabe (2014)：
予測力あり（credit spread，景気指標），予測力なし（市場収益率）• 渡部 (2016)：実現分散に Heterogenous AR モデルとその拡張モデル
10

Implied moment：P のもとでの moment との関係
Pricing kernel m(R) = e−γR の仮定のもとsuppl.

• Bakshi and Madan (2006)：µP , (k = 2, 3, 4) と γ で µQ ，µQ を表現
(k)

(1)

(2)

γ2 (3)(1)(2)µQ = −γµP +µP + o(γ2 )2

γ2 (4)
(2)
(2)
(3)
(2)
µQ = µP − γµP +(µP − (µP )2 ) + o(γ2 )2• Bakshi et al. (2003)：P のもとでの歪度 θP ，尖度 κP で θQ を表現θQ ≈ θP − γ (κP − 3) ( σ2P )1/2
• risk neutral pdf q(r) は，リスク回避度 γ > 0，κP > 3 のとき，現実の確率密度関数 p(r) よりも負の方向に歪む• R の（P で）正規分布に従っているなら，q(r) も歪度=011Implied moment：ボラティリティ・スプレッド
Theorem ( Bakshi and Madan (2006): Volatility spread )σ2Q − σ2P
σ2P

≈ −γ

( σ2P )1/2$!!%&!!'

2× θP + γ ×$%&'

2nd order

3rd

σ2P

2
$%&'
2nd

× ( κP −
$%&'

3

)

4th

vol. スプレッドが拡大するのは• 分布の裾が厚い ( κP > 3 ) ，負の歪度 ( θP < 0 )• リスク回避度 ( γ > 0 )スプレッドがゼロとなるのは• リスク回避度 γ = 0(1)• R ∼ N(µP , σ2P ) ⇒ θP = 0, κP = 3
12

Moment under P確率測度 P のもとでの moment
低頻度観測データ• 月次の分散：日次収益率の利用（月内の営業日数は小標本）• 日次の分散：日次収益率系列に対する rolling window での標本分散
高頻度観測データ• 実現（realized）moment：tick data によるモーメント推定• 実現分散（Realized Variance: RV）：σ2P の推定量• 高頻度観測の収益率の期待値はゼロと仮定 µP = 0(1)

実現 moment（Realized moment: RM）
(k)

RM

=

n
!

r(ti )ki=1

13

実現 moment の漸近的性質• 期間 [0, T] において，価格 S(ti )，0 = t0 < t1 < · · · < tn = T が観測• W(s) は std Wiener 過程，σ(s) は実数値適合過程，J(t) は jump 過程• 収益率 r(ti ) = p(ti ) − p(ti−1 )，時点 s での jump ∆p(s) = p(s) − p(s−)
対数価格 p(t) = ln S(t) が以下のように表されていると考える* tln S(t) = p(t) = p(0) +n → ∞ のとき
実現分散：RM

(2)

p

−→
p

RM(k) −→0* T
0σ(s)dW(s) + J(t)2σ(s) ds +!

!

(∆p(s))20<s≤T(∆p(s))k , k = 3, 4

0<s≤T

14

実現 moment の漸近的性質• RM(k) , (k = 3, 4) は R(t) が従っている確率過程の連続部分に起因する• 3 次の場合は収益率と分散変動の共分散であるレバレッジ効果• 4 次の場合は分散変動の分散を捉えることができていない• ただし，実現 moment の漸近的な性質は ti − ti−1 → 0, (n → ∞) だが，実際には microstructure noise の影響を考慮して，5 分間隔といった時間間隔が選ばれることが多い (生方・渡部, 2011)• microstructure noise に関しては
Ubukata and Oya (2009); Aı̈t-Sahalia and Jacod (2014)• 5 分間隔の RM(k) , (k = 3, 4) では，R(t) の jump 部分に起因する項が支配的になっていると考えられる• 次節での実証では，RM(k) を使った相対的リスク回避度 γ の推定と長期の低頻度観測データによる推定の両方を行い，推定結果を総合的に検討する
15

Estimation for relative risk
aversion

Previous studies• Aı̈t-Sahalia and Lo (1998) and Aıt-Sahalia and Lo (2000)• オプション価格を利用した状態価格密度の nonparametric estimatorsとリスク回避度の推定• リスク回避度で調整された nonparametric VaR• Bliss and Panigirtzoglou (2004)• オプション価格を利用した q(r) の推定とリスク回避度の推定• Bakshi and Madan (2006)• ボラティリティ・スプレッドとリスク回避度 γ の推定• Yoon and Byun (2012)• BM06 による S&P 500, 日経 225 and KOSPI 200 での γ の推定• γ̂： S&P 500 > 日経 225 > KOSPI 200• 大屋 (2017)• BM06 の vol. スプレッドの紹介と実現 moment を使っての γ の推定
16Volatility spread, skewness, kurtosis
Monthly volatility spread• σ2Q には VXJ の 2 乗値の月中平均，σ2P には月次実現分散
実現分散，歪度，尖度 (Amaya et al., 2015)• t 月の第 i 日の第 j 番目の収益率を ri (j), j = 1, . . . , ni
-i
√ -ni
ni
!
ni j=1 ri (j)3
ni nj=1
ri (j)4
2• DRVi =ri (j) , DRSi =, DRSi =
3/2
DRVi2
DRV
j=1
i
• RVt =

Ni
!
i=1

Ni

DRVi ,1 !
RSi =DRSi ,
Ni
i=1

Ni

1 !
RKi =
DRKi
Ni
i=1

月次標本歪度・標本尖度• rolling window による日次収益率を使った標本歪度・標本尖度
17日経平均株価：月次収益率，RV，RS，RK
実現分散

0

-.4

.01

-.2

.02

0

.03

.04

.2

収益率

2000m1

2005m1

2010m1

2015m1

2020m1

2000m1

2005m1

2015m1

2020m1

2015m1

2020m1

実現尖度

5

-1

10

0

15

1

20

2

25

3

実現歪度

2010m1

2000m1

2005m1

2010m1

2015m1

2020m1

2000m1

2005m1

2010m118

日経平均株価：標本歪度，標本尖度
標本尖度(60日window)

0

-3

5

-2

-1

10

0

15

1

20

標本歪度(60日window)

2000m1

2005m1

2010m1

2015m1

2020m1

2000m1

2010m1

2015m1

2020m1

標本尖度(90日window)

0

-3

5

-2

-1

10

0

15

1

20

標本歪度(90日window)

2005m1

2000m1

2005m1

2010m1

2015m1

2020m1

2000m1

2005m1

2010m1

2015m1

2020m119

表 1：歪度・尖度の推定値（期間別）

期間 1
98/Mar – 04/Oct

期間 2
04/Nov – 11/Jun

期間 3
11/Jul – 18/Feb

RS

SS60

SS90

RK

SK60

SK90

mean
se

0.077
0.046

-0.0320.087-0.057
0.096

7.022
0.388

3.502
0.197

3.728
0.225

mean
se

0.254
0.066

-0.2480.084-0.316
0.102

11.679
0.723

3.707
0.372

4.336
0.597

mean
se

0.452
0.081

-0.102
0.090

-0.14
0.097

17.474
0.496

4.235
0.412

4.732
0.577

RS，RK は実現歪度，実現尖度，SS は標本歪度，SK は標本尖度で，添字はウィンドウの大きさを表す。

• 実現歪度は期間 2，3 でその標本平均が有意に正• 標本歪度はいずれの期間でも負だが，有意なのは期間 2 のみ

• 尖度はどちらも，期間 1 から 3 になるにつれて，その標本平均の値を増加させている• 実現尖度は標本尖度よりも大きな値をとっている
検討結果へ

20

分散リスク・プレミアム，ボラティリティ・スプレッド

.04
.03
.02
.01
0

0

.01

.02

.03

.04

.05

インプライド分散

.05

実現分散

2000m1

2005m1

2010m1

2015m1

2020m1

2000m1

2010m1

2015m1

2020m1

ボラティリティ・スプレッド

-.5

-.01

0

0

.5

.01

1

.02

1.5

.03

分散リスク・プレミアム

2005m1

2000m1

2005m1

2010m1

2015m1

2020m1

2000m1

2005m1

2010m1

2015m1

2020m121

リスク回避度の推定

GMM 推定ε(t) =σ2Q (t) − σ2P (t)σ2P (t)

+γ

+

,1/2

σ2P (t)γ2 2
θP (t) −σ (t) (κP (t) − 3)2 P

• θP (t) と κP (t)• 推定 HF：実現歪度，実現尖度• 推定 LF60：60 日 rolling window による標本歪度，標本尖度• 推定 LF90：90 日 rolling window による標本歪度，標本尖度• 操作変数ベクトル Z(t)set 1： Z(t) = (1, σ2Q (t))′ ， set 2： Z(t) = (1, σ2Q (t), σ2Q (t − 1))′• JT (γ) = T gT (γ)′ WT gT (γ) を最小化する γ• WT は Newey-West の HAC ウエイト行列
T
1!
• gT (γ) ≡Z(t − 1)ε(t)
T t=1

22

表 2a：推定 HF

期間

推定

操作
変数

γ

HAC
標準誤差

Hansen’s
J 統計量

p値

set1

5.498
4.352
2.825

0.348
0.358
0.312

1.293
1.435
1.330

0.256
0.231
0.249

set2

5.311
4.439
2.823

0.288
0.332
0.312

2.099
1.650
1.336

0.350
0.438
0.513

1: 98/Mar–04/Oct
2: 04/Nov–11/Jun

HF

3: 11/Jul–18/Feb
1: 98/Mar–04/Oct
2: 04/Nov–11/Jun
3: 11/Jul–18/Feb

HF• 操作変数 set 1，set 2 どちらの推定結果もほぼ同じ• 期間 1 から 3 に向けて，低下傾向23表 2b：推定 LF

期間

推定

1: 98/Mar–04/Oct
2: 04/Nov–11/Jun
3: 11/Jul–18/Feb

LF60

操作
変数

γ

HAC
標準誤差

Hansen’s
J 統計量

p値

set1

7.243
9.153
7.008

3.129
2.868
1.856

2.467
5.081
1.469

0.116
0.024
0.226

2.407
1.473
1.857

5.747
1.499
1.528

0.056
0.473
0.466

1: 98/Mar–04/Oct
2: 04/Nov–11/Jun
3: 11/Jul–18/Feb

LF60

set2

6.079
6.730
7.067

1: 98/Mar–04/Oct
2: 04/Nov–11/Jun
3: 11/Jul–18/Feb

LF90

set1

7.604
6.830
6.373

2.564
1.736
1.601

2.504
2.132
1.435

0.114
0.144
0.231

1: 98/Mar–04/Oct
2: 04/Nov–11/Jun
3: 11/Jul–18/Feb

LF90

set2

5.792
5.872
6.370

2.335
1.084
1.597

2.762
2.389
1.444

0.251
0.303
0.486

• 標準誤差と比べると各期間の推定値の差は顕著ではない• 操作変数 set 2，LF90 の方が標準誤差が小さい

24

表 3：制約付き推定（操作変数 set 2）

推定

γ

HAC
標準誤差

Hansen’s
J 統計量

p値

θP = 0

1
2
3

HF

4.817
3.972
2.508

0.367
0.344
0.331

3.616
1.841
0.859

0.164
0.398
0.651

θP = 0

1
2
3

LF90

10.042
7.711
7.919

2.182
0.771
1.565

5.707
1.705
1.515

0.058
0.426
0.469

κP = 31
2
3

HF-85.169-13.970-10.383

64.583
6.526
3.179

0.711
1.427
0.556

0.701
0.490
0.757

κP = 3

1
2
3

LF90

4.111
11.392
9.968

4.589
4.111
4.903

3.869
1.679
2.054

0.145
0.432
0.358

制約

期間

• 制約 θP = 0 については，推定 HF では制約なしの結果とほぼ同じで，推定 LF でも差は小さい• 制約 κP = 3 は明らかに間違っている

25

検討結果

期間 1

表 1

より歪度は RS でも SS60 ，SS90 でもゼロからの乖離は有意ではないので，推定結果は推定 HF(IV set2) を採用

期間 2 表 1 より歪度は，RS は有意に正，SS60 ，SS90 ともに有意に負を示しているが，表 3 での歪度=0 の制約を課した推定結果と制約を課さない推定結果に大きな差はないため，s.e. の小さい
推定 HF(IV set2) を採用
期間 3 表 1 より RS は有意に正となっているが，表 3 の結果より，推定 HF(IV set2) を採用

26

ローリング・ウィンドウによる推定結果: 推定 HF

γ の推定値
95%信頼区間

0

1

2

3

4

5

6

7

8

9

10

推定HF (IV set 2) による γ の推定値

2001m1 2003m1 2005m1 2007m1 2009m1 2011m1 2013m1 2015m1 2017m1

27

ローリング・ウィンドウによる推定結果: 推定 LF

30

推定LF90 (IV set 2) による γ の推定値

-10

0

10

20

γ の推定値
95%信頼区間2001m1 2003m1 2005m1 2007m1 2009m1 2011m1 2013m1 2015m1 2017m1

28

θQ
θP

-9 -8 -7 -6 -5 -4 -3 -2 -1

0

1

2

3

4

ローリング・ウィンドウによる推定結果: P と Q のもとでの歪度2001m1 2003m1 2005m1 2007m1 2009m1 2011m1 2013m1 2015m1 2017m1

Bakshi et al. (2003)：θQ ≈ θP − γ (κP − 3) ( σ2P )1/2
29

Summary and future studies

まとめと今後の検討課題

まとめ

1. リスク中立確率測度から市場参加者のリスクに対するスタンスに関して，リスク回避度をひきだす Bakshi and Madan (2006) の方法
を日本市場の日経平均株価に適用
2. 金融危機以前はリスク回避度は低下傾向を示していたが，金融危機を機に一旦，リスク回避度の上昇がみられ，その後，再び低下傾向を示していることが明らかになった
今後の検討課題• θP ，κP の推定• pricing kernel(utility function) の特定化
30

参考文献

ReferencesAı̈t-Sahalia, Yacine and Jean Jacod (2014) High-frequency financial econometrics :Princeton University Press.
Aı̈t-Sahalia, Yacine and Andrew W Lo (1998) “Nonparametric estimation of state-pricedensities implicit in financial asset prices,” Journal of Finance, Vol. 53, No. 2, pp.
499–547.
Aıt-Sahalia, Yacine and Andrew W Lo (2000) “Nonparametric risk management andimplied risk aversion,” Journal of econometrics, Vol. 94, No. 1-2, pp. 9–51.
Amaya, Diego, Peter Christoffersen, Kris Jacobs, and Aurelio Vasquez (2015) “Doesrealized skewness predict the cross-section of equity returns?” Journal of Financial
Economics, Vol. 118, No. 1, pp. 135–167.
Bakshi, Gurdip and Dilip Madan (2006) “A theory of volatility spreads,” Management
Science, Vol. 52, No. 12, pp. 1945–1956.

31

Bakshi, Gurdip, Nikunj Kapadia, and Dilip Madan (2003) “Stock return characteristics,
skew laws, and the differential pricing of individual equity options,” Review of
Financial Studies, Vol. 16, No. 1, pp. 101–143.
Bliss, Robert R and Nikolaos Panigirtzoglou (2004) “Option-implied risk aversion
estimates,” Journal of Finance, Vol. 59, No. 1, pp. 407–446.
Bollerslev, Tim, George Tauchen, and Hao Zhou (2009) “Expected stock returns andvariance risk premia,” Review of Financial Studies, Vol. 22, No. 11, pp. 4463–4492.
Bollerslev, Tim, James Marrone, Lai Xu, and Hao Zhou (2014) “Stock return
predictability and variance risk premia: statistical inference and international
evidence,” Journal of Financial and Quantitative Analysis, Vol. 49, No. 3, pp.
633–661.
Carr, Peter and Dilip Madan (2001) “Optimal positioning in derivative securities.”
Fornari, Fabio and Antonio Mele (2001) “Recovering the probability density function ofasset prices using GARCH as diffusion approximations,” Journal of Empirical
Finance, Vol. 8, No. 1, pp. 83–110.

32

Fukasawa, Masaaki (2012) “The normalizing transformation of the implied volatility
smile,” Mathematical Finance, Vol. 22, No. 4, pp. 753–762.
Fukasawa, Masaaki, Isao Ishida, Nabil Maghrebi, Kosuke Oya, Masato Ubukata, andKazutoshi Yamazaki (2011) “Model-free implied volatility: from surface to index,”
International Journal of Theoretical and Applied Finance, Vol. 14, No. 4, pp.
433–463.
Gatheral, Jim (2006) The volatility surface: a practitioner’s guide : John Wiley & Sons.
Ubukata, Masato and Kosuke Oya (2009) “Estimation and testing for dependence inmarket microstructure noise,” Journal of Financial Econometrics, Vol. 7, No. 2, pp.
106–151.
Ubukata, Masato and Toshiaki Watanabe (2014) “Market variance risk premiums inJapan for asset predictability,” Empirical Economics, Vol. 47, No. 1, pp. 169–198.
Yoon, Sun-Joong and Suk Joon Byun (2012) “Implied risk aversion and volatility risk
premiums,” Applied Financial Economics, Vol. 22, No. 1, pp. 59–70.

33

生方雅人・渡部敏明 (2011) 「実現ボラティリティ–ボラティリティの計測方法の発展とリ
スクマネジメントへの応用可能性 (特集ボラティリティと証券投資)」，『証券アナリスト
ジャーナル』，第 49 巻，第 8 号，16–26 頁．
大屋幸輔 (2011) 「ボラティリティの景気予測力 -バリアンス・リスクプレミアムの検証から-」，『浅子和美・飯塚信夫・宮川努編『世界同時不況と景気循環分析』, 第 7 章』，
141–157 頁．

(2017) 「ボラティリティ・スプレッド」，『先物・オプションレポート』，第 29 巻，第 12 号．
渡部敏明 (2016) 「日経 225 分散リスク・プレミアムの予測力」，『先物・オプションレ
ポート』，第 28 巻，第 11 号．

34

補足：効用関数と pricing kernel

Find pt ( the value at time t ) of a payoff xt+1
• Assume that the investor can freely buy or sell as much of thepayoff xt+1 (which is r.v.) as he wishes, at the price pt .
• Let et be the original consumption level at time t and ξ be the amount of
the asset he choose to buy.
• the consumptions at t and t + 1 are ct = et − pt ξ and ct+1 = et+1 + xt+1 ξ• Investor’s utility function U(ct , ct+1 ) = u(ct ) + β Et [ u(ct+1 ) ]• u(·) is increasing, reflecting a desire for more consumption, and concave,reflecting the declining marginal value of additional consumption• β captures impatience, called the subjective discount factor

max u(ct ) + β E[ u(ct+1 ) ],ξ

s.t. ct = et − pt ξ and ct+1 = et+1 + xt+1 ξ

The first order condition.
/u′ (ct+1 )u′ (ct+1 (xt+1 ))pt = Et β ′xt+1 = Et [ mt+1 (xt+1 ) xt+1 ], mt+1 (xt+1 ) = βu (ct )u′ (ct )35

Under the assumption of• time-additive utility func. u(St ) + β u(St+1 ), β: subjective discount factor• power utility u(S) =S1−γ − 1, γ is the relative risk-aversion
1−γ

Pricing kernel (stochastic discount factor) is given asu′ (St+1 )mt+1 = β ′= β e−γ (log St+1 −log St ) = β e−γ Rt+1
u (St )
where Rt+1 = log St+1 − log St . Then the risk-neutral density q(R) is
q(R) = (m(R) p(R)m(R) p(R)dR= p(R) × (e−γ Re−γ R p(R)dR• risk-averse investors pay more attention to unpleasant states incomparison to the physical counterpart.
slide 436

Representaion µ(1) in terms of µ(k)PQ

• Denote m.g.f. of R under p(R) and q(R) as MP (λ) and MQ (λ).
• Assume µP = 0*
MP (λ) =eλ r p(r)dr* 0
1
1 2 2 1 3 3
1 4 4=
1 + λ r + λ r + λ r + λ r p(r)dr + o(λ4 )2624
λ2 (2) λ3 (3) λ4 (4)= 1 + µP + µP + µP + o(λ4 )
2
6
24
(
*
eλ r e−γ r p(r) dr MP (λ − γ)λr
MQ (λ) =e q(r)dr = (=MP (−γ)e−γ r p(r)dr
(1)

*

µQ =

r q(r)r =22
MP (λ − γ)′ 2λ=0
MP (−γ)slide 11

37

Theorem 1 in BM(2006): Variance spread
For m(R) satisfying m(R) = 1 − A1 R + (1/2) A2 R2 + O(R3 ) where2222
A1 = − ∂m/∂R2R=0 and A2 = ∂2 m/∂R2 2R=0 ,σ2Q − σ2P
σ2P

≈ − A1 ( σ2P )1/2 × θP +• Suppose that m(R) =* ∞0

A2 σ2P
2⎛⎞
2 A21 ⎟⎟⎜⎜⎟⎟ .
× ⎜⎜⎝ κP − 1 −A ⎠2e−zR ν(dz) for some measure ν on R+ .

• This m(R) encompasses HARA marginal utility class and loss-aversion
utilities.
• For example z ∼ G(γ, 1/b), then E[z] = γb, E[z2 ] = γ(γ + 1)b2
m.g.f. of z is M(λ) = (1 − bλ)−γ .

• Then we have m(R) = M(−R) = (1 + bR)−γ which represents HARA
marginal utility.
38

Quotes competition in limit order markets
Wataru Ohta
Graduate School of Economics
and Center for Mathematical Modeling and Data Science
Osaka University

Outline

1. The purpose of this study
2. Extant empirical studies
3. Extant theoretical studies
4. Model5. Equilibrium
6. Empirical study in the Tokyo Stock Exchange
7. Summary2

1. Purpose• Order submission strategies in limit order markets• no designated market makers• matching orders submitted by public traders• How transactions take place in limit order markets• Is there a systematic pattern to the ask quote and bid quote?
• What is the optimal order submission strategy under a certain market
condition?
• a limit order undercutting the best quote by one tick
• a limit order undercutting the best quote by more than one tick
• a limit order at the best quote• a limit order submitted behind the best quote• What affects the size of holes?
• How does the tick size affect quote competition of limit orders?

32-1. Empirical studies• Higher volatility can relate with lower liquidity.
• Egginton, Van Ness, Van Ness (2016)• Quote stuffing relates with decreased liquidity and increased
short-term volatility.

• Hasbrouck (2018)• Higher short-term volatility relates with recurrent cycles of undercutting
similar to the Edgeworth cycles.

42-2. June 24, 2002, Code 7269 in Tokyo

1475

1470

1465

1460

14 500

15 00015 500

16 00016 50017 00017 500• Horizontal axis denotes time in seconds from the opening of the day,and vertical axis denotes price.
• Red line is the ask, blue line is the bid, black dots indicate
transactions.

52-3. June 4, 2002, Code 8273 in Tokyo

1983.01982.5

1982.01981.5

1981.01980.516 00018 000

20 000• Traders queue at the same quote.
• There is no quote competition.
• Tick size is the binding constraint in choosing quotes of limit orders.

63-1. Theoretical studies• There can be too many states of the limit order book.
• Foucault (1999, Journal of Financial Markets)• assumes that limit orders are automatically canceled in one period
after their submission in a discrete time setting.

• This paper assumes that limit orders expire in two periods after theirsubmission.

• Foucault, Kadan, and Kandel (2005, Review of Financial Studies)• Sellers and buyers arrive at the market alternately.
• Traders cannot submit limit orders at or behind the best quotes.
• This paper assume sellers and buyers arrive at the market randomly
and traders can choose any price under the tick size restriction.

• Parlour (1998, Review of Financial Studies)• assumes that the bid-ask spread is one tick

• Goettler, Parlour, and Rajan (2005, Journal of Finance; 2009,
Journal of Financial Economics)• numerical solutions• Rosu (2009, Review of Financial Studies; 2009, Mimeo)• continuous time setting with zero tick size

7

quoted price3-2. Theoretical results of this paper
16
14
12
10
8
6

16
14
12
10
8
6
4
0

10

20
period3040• An example of quote dynamics on the equilibrium path when limit
orders expire in two periods• Red lines denote ask quotes, blue lines denote bid quotes, dots
indicate transactions.
• Cyclic quote dynamics• There is possibility of queuing under the large tick size.

84-1. Model• Discrete and infinite periods• In each period, one potential trader arrives at the exchange andsubmits an order.
• The probability of a trader arriving is α ∈ (0, 1], and the probability ofa seller conditional on a trader arriving is β ∈ (0, 1).
• α and β are exogenous.
• α is the trader arrival rate.
• β represents the order imbalance.
• A trader is a seller with the probability of π s = αβ, and a buyer with theprobability of π b = α(1 − β).

• No asymmetric information

9

4-2. Model• A trader can submit an order with a volume of one share only when
he arrives at the market.
• The order is either a limit order or a market order.
• A limit order: price contingent• chooses the limit price under the restriction of the tick size k.
• cannot cancel or modify his order once he submits it.
• a limit order automatically expires in two periods after their submission.

• A market order: not price contingent• immediately executed at the best available price on the market

104-3. Model• A seller holds one share of the asset and evaluates it as v L ≥ 0.
• A payoff for a seller is P − v L if he sells a share at a price P.

• A buyer holds no share, and evaluates a share as v H (v H > v L ).
• The payoff for a buyer is v H − P if he buys a share at a price P.

• A trader receives zero payoff if he does not trade.
• The discount rate is assumed to be zero.
• Traders maximize expected payoffs.
• The market is transparent: traders observe the book before
submitting orders.

114-4. Trading rule
• The price-time precedence rule
• Among limit buy (sell) orders, the highest (lowest)-priced orders havepriority.

• Among limit orders at the same price, priority is given to orders
submitted earlier.

• The discriminatory pricing rule
• a transaction price is the price of the matched limit order.

• The order choice depends on the state of the book and the futureorder flow.
• Trade-off between price and execution probability in submitting anorder.

• a market order: executed immediately at a low price• a limit order: execution is uncertain though the price is favorable oncethe order is executed.

12

4-5. Pure-strategy Markov Perfect Equilibrium• Markov strategy which depends only on the type of trader and on thebook.
• Pure strategy.

13

quoted price5-1. An example of equilibrium quote dynamics
16
14
12
10
8
6

16
14
12
10
8
6
4
0

10

20
period30

40• v L = 0, v H = 21, k = 1, the trader arrival rate α = 1, and the orderimbalance β = 1/2.
• At first, one-tick quote-cutting• At a certain level, the best quote jumps more than one tick.
• A limit sell order behind the best ask is submitted.

14

5-2. Edgeworth cycle• An equilibrium with cyclic quote dynamics exists when the tick size k
is small enough.
• Edgeworth cycle• Edgeworth (1925)• Maskin and Tirole (1988, Econometrica)• Noel (2007, Journal of Industrial Economics) in Canadian retailgasoline markets

15

5-3. Edgeworth cycle• A seller who arrives at an empty book submits a limit sell order at a
relatively high ask, and allows the next seller to undercut his ask.
• Hindering quote-cutting by the aggressive ask is costly.

• The next seller undercuts the best ask by only one tick.
• expects further quote-cutting.
• one-tick quote-cutting minimizes the cost in price to get higher priority
against the limit orders in the book.

• When the best ask reaches a certain level after one-tick
quote-cutting, a seller submits a very aggressive ask (quote jump)• Jumping quote is reasonable because hindering further quote cutting
raises the execution probability.

• The next seller submits a limit sell order far behind the best ask.
• gives up further quote-cutting because the price has already become
too low.

• The execution probability of a limit order behind the best ask is small.
• The high price compensates for the large loss in the execution
probability.

16

5-4. The trader arrival rate α and the critical asks
1.0

0.8

Al
Af

0.6

0.4

Au

0.2

Am
Α

0.0

0.2

0.4

0.6

0.8

1.0

• Horizontal axis is the trader arrival rate α and vertical axis is price.
• A f is the first ask submitted to an empty book.
• Au is the end of the range of one-tick quote-cutting.
• A m is the most aggressive ask.
• A l is the least aggressive ask.

17

5-5. The order imbalance β and the critical asks
1.0

Al

0.8

Af

0.6

Au

0.4

Am
0.2

Β
0.2

0.4

0.6

0.8

1.0

• When sellers outnumber buyers (when β is greater),• Sellers submit more aggressive limit sell orders• narrow the size of holes on the ask side.
• shorten the steps of one-tick quote-cutting• Holes created by the jumping quote on the ask side are smaller.
• Cyclic quote dynamics on the ask side are more frequently observed

185-6. Multiple equilibria: Queuing• When the tick size is large, traders queue at the same quote.
• The cost of quote-cutting to obtain price priority is large.
• Queuing does not occur when the tick size is small.

• The higher the arrival rate of sellers (buyers) is, the more likely
sellers (buyers) queue at the same quote.
• The benefit of quote-cutting is small when the number of traders on the
opposite side of the market is small.

• Order imbalance caused by liquidity needs increases the depth.

196-1. Empirical study in the Tokyo Stock Exchange• Data: Nikkei Needs tick data
• Period: from August 2014 to December 2015
• Stocks in the TOPIX100• Tick size reduction after July 22, 2014
• Price range from 500 yen to 1,000 yen (the tick size is 0.1 yen) andfrom 5000 yen and 10,000 yen (the tick size is 1 yen)

• 26 stocks with more than 200 trading days satisfying the above pricerange condition• Divide the trading period into 5-minute intervals• Proportion of the number of specific types of order submission to thetotal number of quote changes for each interval• Cycle: two consecutive quote cutting followed by quote rebounding• OIB: (number of market and limit buy order - number of market andlimit sell order)/(number of market and limit buy order + number ofmarket and limit sell order).
• Winsorize variables at the 2nd and 98th percentiles

206-2. Frequency of cycles
Frequency of cycles, Relative spread, and RV
25
20
15
10
0

0

0

5

.05

.1

.05
.1
.15
RV (15-second interval)

.2

.15

Frequency of cutting, jumping, rebounding and queueing

900

1130
5-minute interval

1500• Black: proportion of cycles to thetotal number of quote changes• Blue: relative bid-ask spread• Green: realized variance

900

1130
5-minute interval

1500• Black: proportion of quotecutting• Blue: proportion of quote
jumping• Green: proportion of quoterebounding• Red: proportion of queuing21

6-3. Quote cycles in the 5-min interval from 9:05 to 9:09
• Total: 8268 stock*days• 4913 stock*days have no cycle, 3355 stock*days have cycles
N of cycles

Mean
SD
Min
Max
Mean
SD
Min
Max

Tick spread

Relative Realized
spread varianceIntervals without cycles (4913 stock*days)
0
4.20
0.06
0.09
0
1.98
0.03
0.08
0
1.41
0.19
0.006
0
17.46
0.26
0.82
Intervals with cycles (3355 stock*days)
3.87
5.81
0.08
0.18
6.69
2.79
0.04
0.13
1
1.46
0.02
0.006
138
17.46
0.26
0.82

Realized
Kurtosis
4.21
1.87
1.86
15.63
4.13
1.86
1.86
15.63

22

6-4. Bid-ask spread and moments of quote-midpoint returns• Bid-ask spread and moments of quote-midpoint returns classified byquintiles of OIB for the 5-minute intervals from 9:05 to 9:09

Quintile 1 (low OIB)
Quintile 2
Quintile 3
Quintile 4
Quintile 5 (high OIB)
Total
F-valuep-valueTickspread
4.90
4.81
4.77
4.90
4.89
4.85
1.00
0.41

Relativespread
0.069
0.068
0.067
0.069
0.068
0.068
0.66
0.62

Realizedvariance
0.129
0.124
0.127
0.129
0.135
0.129
2.29
0.06

Realized
skewness-1.031-0.449
0.024
0.498
1.039
0.019
1091.1
0.00

Realized
kurtosis
4.266
4.141
4.032
4.175
4.286
4.179
4.87
0.00

• F-value and p-value are for F-test of equality of quintile means23

6-5. Order submission• Proportion (percent) of specific sell orders to the total number ofquote changes classified by quintiles of OIB for the 5-minuteintervals from 9:05 to 9:09

Quintile 1 (low OIB)
Quintile 2
Quintile 3
Quintile 4
Quintile 5 (high OIB)
Total
F-valuep-value

Cycle onthe ask side
0.12
0.07
0.04
0.04
0.02
0.06
90.82
0.00

One-tick
cutting
11.06
9.19
7.96
6.89
5.36
8.09
831.6
0.00

Quote
jumping
2.49
2.19
2.04
1.97
1.71
2.07
89.21
0.00

Quote
rebounding
1.54
1.17
0.96
0.85
0.67
1.04
158.6
0.00

Queuing
23.41
19.75
17.18
14.70
11.17
17.23
2645
0.00

• F-value and p-value are for F-test of equality of quintile means24

7. Summary• In limit order markets, cyclic quotes dynamics can be observed, andholes emerge in the book during cycles.
• When the trader arrival rate is high, cyclic quote dynamics are more
frequently observed and the holes are larger.
• When sellers outnumber buyers, cycles on the ask side are more
frequently observed and the holes on the ask side are smaller.
However, sellers queue at the same ask if the tick size is relatively
large.

25

References
Egginton, J. F., B. F. Van Ness, and R. A. Van Ness, 2016, Quote Stuffing, Financial
Management 45(3): 583-608.
Foucault, Thierry, 1999, Order flow composition and trading costs in a dynamic limit ordermarket, Journal of Financial Markets 2, 99-134.
Foucault, Thierry, Ohad Kadan, and Eugene Kandel, 2005, Limit order book as a market for
liquidity, Review of Financial Studies 18, 1171-1217.
Goettler, Ronald L., Christine A. Parlour, and Uday Rajan, 2005, Equilibrium in a dynamic
limit order market, Journal of Finance 60, 2149-92.
Goettler, Ronald L., Christine A. Parlour, and Uday Rajan, 2009, Informed traders and limit
order markets, Journal of Financial Economics 93, 67-87.
Hasbrouck, J., 2018, High-Frequency Quoting: Short-Term Volatility in Bids and Offers,
Journal of Financial and Quantitative Analysis 53(2): 613-641.
Maskin, Eric, and Jean Tirole, 1988, A theory of dynamic oligopoly, ii: Price competition,kinked demand curves, and Edgeworth cycles, Econometrica 56, 571-99.
Noel, Michael D., 2007, Edgeworth price cycles: Evidence from the Toronto retail gasoline
market, Journal of Industrial Economics 55, 69-92.
Parlour, Christine A., 1998, Price dynamics in limit order markets, Review of Financial
Studies 11, 789-816.
Rosu, Ioanid, 2008, Liquidity and information in order driven markets, Mimeo.
Rosu, Ioanid, 2009, A dynamic model of the limit order book, Review of Financial Studies
22, 4601-41.

26

A Two-Sample Alternative
to Using Instruments
in Regressions with Omitted Variables
Masayuki Hirukawa Irina Murtazashvili Artem Prokhorov
Ryukoku University
Drexel University
University of Sydney
7 August 2019
Fukushima University

Fukushima University

1

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Motivation
Relevant variables that are omitted lead to a challenging problem.

Examples of the missing regressor:
1. An ability measure in Mincer’s (1974) wage regression.
– Card (1995) argues that the estimation result su¤ers from the “abilitybias” unless the regression includes a variable representing ability asa regressor.
– A test score is typically considered as a variable representing ability.
– Test scores are unavailable in the Current Population Survey (CPS)and the Panel Study of Income Dynamics (PSID).

7 August 20191

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

2. Work experience in gender wage gap (e.g., Zabalza and Arrufat, 1985;
Black, Trainor, and Spencer, 1999).
– No actual work experience is included in the General Household Survey (GHS) and CPS.
– Actual work experience is available in the National Longitudinal Survey (NLS) and PSID.

In the absence of a proxy, the most common solution is to look for an
instrument, but...
1. Weak instruments?
2. Finite-sample properties of IV estimates?

7 August 20192

Fukushima University

1.1

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Our Proposal
Why not turn to a two-sample approach using a proxy?

Our two-sample estimation procedure is to run ordinary least squares (OLS)
after replacing the missing regressor with a nonparametric estimate of itsconditional mean.
– This estimator does not rely on the combined sample constructed via
the nearest-neighbor matching (NNM).
– Hirukawa and Prokhorov (2018) have already proposed a NNM-based
two-sample estimation (to be discussed shortly).

7 August 20193

Fukushima University

1.2

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Contributions

1. Nonparametric imputation of the missing regressor:
Imposing a parametric form on the conditional mean of the missing regressor (e.g., Fang, Keane and Silverman, 2008; Flavin and Nakagawa,2008) has at least two problems:(a) Bias due to misspeci…cation of the conditional mean.
(b) Lack of adjustments in standard errors under the assumption that
imputed predictors are free of estimation errors.
2. Nonparametrically generated regressors:
Our approach is analogous to Pagan (1984).
The regressor is not imputed within the same sample in our estimation
procedure.
7 August 20194

Fukushima University

A Two-Sample Alternative to IV

Plan of Talk
1. Motivation
2. Setup
3. Strategy 1: IV Estimation
4. Strategy 2: Two-Sample Estimation
5. Finite-Sample Performance
6. Application to Estimation of Return to Schooling
7. Concluding Remarks
7 August 20195

Hirukawa, Murtazashvili & Prokhorov

Fukushima University

2

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Setup

True Model.
0
Y = 0 + X10 1 + X20 2 + X3I
3 + u;

where X1 2 Rd1 , X2 2 Rd2 and X3I 2 Rd3I .
Either 1 or 3 is the parameter of interest.

Standard Situation.

A single complete sample (Y; X1; X2; X3I ) is available.

If E ( uj X1; X2; X3I ) = 0, then OLS for this regression is consistent.

7 August 20196

Fukushima University

Our Case.

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

X2 is missing (in the primary sample).

Two Datasets.
Y
X

X1
X

X2

X3
S1
X
S2
X
X
Note: X3 = (X3I ; X3E )Estimation Strategies.
1. IV estimation with X2 left omitted, where only S1 is used.
2. Two-sample estimation with (a proxy of) X2 imputed from a di¤erent
data source, where both S1 and S2 are used.
7 August 20197

Fukushima University

3

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Strategy 1: IV Estimation
The true regression0
Y = 0 + X10 1 + X20 2 + X3I
3+ucan be rewritten as a “short” regression000
Y = 0 + X10 1 + X3I
3 + u + X2 2 = XS S + v;
00where XS := 1; X10 ; X3I.

– X2 is treated as the vector of omitted variables.
– (X1; X3I ) and X2 may be correlated.

7 August 20198

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

It is still possible to estimate S (and thusof instruments Z 2 RdZ is available.

1 ) consistently if the vector

The estimation is based on the moment restriction
E (ZS v ) := E

("

1
Z

#

u + X20 2)= 0(dZ +1) 1:– Even if dZ
d1 + d3I is satis…ed, consistency of this estimation
strategy is built on somewhat restrictive assumptions such as
Z q X2 and E (X2)0 2 = 0:

– A simple IV estimator in the just-identi…ed case (dZ = d1 + d3I ) is
1 10
1
n
1X
0
^ IV;S := Q
^ 1R
^ ZY := @A@
ZS;iXS;iZS;iYiA :
ZX
n i=1
n i=1
0

7 August 2019

n
1X

9

Fukushima University

4

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Strategy 2: Two-Sample Estimation
What if instruments may be unavailable or there may be other reasons foravoiding Strategy 1 (e.g., weak instruments)?

Suppose that econometricians have access to the auxiliary data set
S2 = S2m =nom
X2j ; X3Ij ; X3Ej:
j=1

S1 and S2 contain common variables X3 = (X3I ; X3E ).

– X3I and X3E are included in and excluded from the regression0
Y = 0 + X10 1 + X20 2 + X3I
3 + u:

7 August 201910

Fukushima University

4.1

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Additional Notations
The vector of matching variables X3 consists of continuous (“C ”) anddiscrete (“D”) variables so that
d3 := dim (X3) = dim (X3C ) + dim (X3D ) =: d3C + d3D :– Assume that d3C > 0, i.e., S1 and S2 have at least one continuous
variable in common.
Let X3 := X3C X3D , where X3C := supp (X3C ) and X3D := supp (X3D ).
De…neg (X3) ==7 August 2019""

g1 (X3)
g2 (X3)
1
2

#

:=

#

"
11

:=
X1
X2

"

E ( X1j X3)
E ( X2j X3)
g1 (X3)
g2 (X3)

#

:

#

and

Fukushima University

4.2

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

The Estimator
Reformulate the true “long” regression0
Y = 0 + X10 1 + X20 2 + X3I
3+uas
Y = X0 + ;

where
X :=

0
0
0
0
1; X1; g2 (X3) ; X3I and:= u + 02 2:

The estimator fundamentally takes the form of an OLS-type regression of
Y on X .
– Actually g2 ( ) is unknown and thus must be estimated nonparametrically.

7 August 201912

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Given g^2 ( ), some consistent nonparametric estimator of g2 ( ), the entire
estimation procedure takes the following two steps:
1. Regardnom
X2j ; X3jin S2 and fX3ign
i=1 in S1 as m data and
j=1

n design points, respectively, and obtain n nonparametric estimatesfg^2 (X3i)gn
i=1 .

0
0
0
0
^
2. Run OLS for the regression of Yi on Xi := 1; X1i; g^2 (X3i) ; X3Ii .

The estimator ofin the …nal form is

1 1
n
X
1
^ P ILS := Q
^ iX
^ 0A
^ iYi:^ 1 R^ ^ := @
X
X
i
^X
^ XY
X
n i=1
n i=1
0

n
1X

– Because this estimator is the OLS with g^2 (X3i) imputed in place ofthe missing regressor X2i, we call it the plug-in least squares (PILS)
estimator hereinafter.
7 August 201913

Fukushima University

4.3

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Kernel Regression Smoother for Mixed Data
From among many nonparametric methods, we adopt the kernel regression
smoother for mixed continuous and categorical data by Racine and Li(2004) to estimate g2 ( ).
The kernel is W (t3; x3; h; ) := K (t3C ; x3C ; h) L (t3D ; x3D ; ).
K:
L:
t3:
x3 :
h; :

product kernel for the continuous component of X3.
product kernel for the discrete component of X3.
data point.
design point.
bandwidths.

The nonparametric regression estimator of g2 ( ) is de…ned as

7 August 2019Pm
j=1 X2j W X3j ; X3i; h;g^2 (X3i) := P
m W X ; X ; h;
3j
3i
j=1
14

; i = 1; : : : ; n:Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Which univariate kernel to choose for the continuous component?
– A standard symmetric kernel such as the Epanechnikov kernel; or– The beta kernel by Chen (1999), which takes the form oftx=b (1 t)(1 x)=b
KB(x;b) (t) =1 f0
B fx=b + 1; (1 x) =b + 1gt

1g

for the design point x 2 [0; 1] and the smoothing parameter b.

Figure: Shapes of the Beta Kernel for b = 0:2 andx = 0:00; 0:25; 0:50; 0:75; 1:00
7 August 201915

Fukushima University

4.4

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Regularity Conditions

Assumption 1. The two random samples (S1; S2) = (S1n; S2m) are drawn
independently from the joint distribution of (Y; X1; X2; X3) with …nite fourthorder moments.

Assumption 2. X3C is continuously distributed with a convex and compact
support X3C , and its density is bounded and bounded away from zero on X3C .

Assumption 3. (i) E ( uj X1; X3) = 0 and 2u (X1; X3) := E u2 X1; X3 2
(0; 1). (ii) E 1 02 = 0. (iii) g2 ( ) is non-constant on X3C if X3E contains at least one continuous variable, and g2 ( ) is strictly nonlinear on X3Cotherwise.

7 August 201916

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Assumption 4. (i) Let f ( ) be the marginal pdf of X3C . Then, the secondorder derivatives of f ( ) and f ( ) g2 ( ) with respect to X3C are continuousand bounded uniformly on X3C . (ii) There exist some constants 2 (0; 1)
and C 2 [1; 1) so that supx32X3 E jY j2+ X3 = x3
C.
Assumption 5. The univariate continuous kernel is either (a) a symmetricand bounded pdf that satis…es the …rst-order Lipschitz condition or (b) thebeta kernel.

Assumption 6. Sequences of the smoothing parameters hp (= hp (m) > 0),bp (= bp (m) > 0) and q (= q (m) 2 (0; 1)) and the boundary parameterp = p (m) > 0 for p = 1; : : : ; d3C and q = 1; : : : ; d3D satisfy one of
the following conditions as m ! 1: (a) for a symmetric kernel, hp; q ! 0
Qd3C
and log m= m p=1
hp ! 0; and (b) for the beta kernel, bp; q ; p ! 0,!
r
Qd3C
bp= p ! 0 and log m= m
p=1 bp p ! 0.

7 August 201917

Fukushima University

4.5

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

General Asymptotic Results of PILS

Theorem 1.

Theorem 2.

p

If Assumptions 1-6 hold, then ^ P ILS !

If Assumptions 1-6 hold, then1
N 0(d+1) 1; VX := N 0(d+1) 1; QXX^ 1 1
Bg2 = Q
^X
^n
X

n h
X

i=1

^ i fg2 (X3i)
X^ ) fg^2 (X3i)
E (Xas n; m ! 1.

p

dn ^ P ILS
Bg2 !
1
QXX
as n; m ! 1, whereg^2 (X3i)g0

E (^g2 (X3i)j X3i)g

n
X
1
^ iX
^ 0;
^^^ =
X
QXX = E XX 0 ; Q
i
XX
n i=1

0

i

2;

n
1 X
0 ; and= lim piin;m!1 ni=1^^ ) [^E (X
g2 (X3i) E f g^2 (X3i)j X3ig]0 2:
i = Xi i

7 August 201918

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

In addition, the asymptotic variance VX can be consistently estimated by^ 1 ^ ^ ^Q
^ 1 ;V^X := Q
^ ^
^ ^
XX

XX

where
n
X
1
^^^ =
^ i ^ 0i;
n i=1

0

n
1X

^i = X^ i^i + @^i^i = Yi

^i =

0
@

n i=1

^ 0 ^ P ILS ;
X
i

n
1 X

n k=1

1

^iA ;

1

^ k A g^2 (X3i)0 ^ P ILS;2;X

and ^ P ILS;2 is the PILS estimate of2.

Remark. PILS is a partial mean estimator of Newey (1994), and thus ^ ^ ^can be obtained straightforwardly by his approach.
7 August 201919

Fukushima University

4.6p

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorovn-Asymptotic Normality vs Curse of Dimensionality

Assumption 6’. Sequences of the smoothing and boundary parameters satisfy
one of the followings as m ! 1: (a) for a symmetric kernel, there is a bandwidth h (= h (m) > 0) so that h1; : : : ; hd3C / h, 1; : : : ; d3D / h2, h !

0, and log m= mhd3C ! 0; and (b) for the beta kernel, there are a smoothing parameter b (= b (m) > 0) and a boundary parameter (= (m) > 0)so that b1; : : : ; bd3C / b, 1; : : : ; d3D / b, 1; : : : ; d3C / , b; ! 0,nb= ! 0, and log m= m (b )d3C =2o! 0.

Corollary 1. Let h / (log m=m) and b / (log m=m)2 for some constant> 0. Also suppose that one of the following divergence patterns in (n; m)is true: (i) n=m !
2 (0; 1); (ii) n=m ! 0; or (iii) n=m ! 1 and
p ^
d
4
n=m
! 0. Then, under Assumptions 1-5 and 6’, n P ILS!

N 0(d+1) 1; VX holds for d3C
7 August 20193 as n; m ! 1.
20

Fukushima University

4.7

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

A Comparison with the MSII Estimator

NNM.

The matched-sample indirect inference (MSII) estimation by Hirukawa
and Prokhorov (2018) starts from obtaining a combined sample by meansof NNM on the observed variables.

For each observation of X3 in S1, K -NNM picks out K closest matchesof X2 from S2 through …nding …rst to K th closest matches of X3 in S2
with respect to some distance function such as the Mahalanobis distance.
– The resulting combined sample can be written as
S=

7 August 2019non
Yi; X1i; X2j1(i); : : : ; X2jK (i); X3Ii; X3Ei
:
i=1

21

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Matched-Sample OLS (MSOLS).

The OLS estimation (MSOLS) for the regression of Yi on
Xi;j(i) :=00001; X1i; X2j(i); X3Ii ;

where
K
1 X
X
;
X2j(i) :=K k=1 2jk (i)generates a non-vanishing, classical measurement error bias.
– The bias is attributed to using X2j(i) as a proxy for X2i.
– The source of attenuation (= bias toward zero) is

7 August 2019

22

2=E

0
2 2 .

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Bias-Corrected Estimation.
As is well known in the literature on errors-in-variables models, the bias
cannot be corrected in general without imposing additional identi…cation
conditions.
– It can be corrected analytically with no such extra conditions because
S2 serves as repeated measurements.
The MSII is a bias-corrected estimator de…ned as
0

where

^ II := @

n
1X

n i=1

0
Xi;j(i)Xi;j(i)

n

1 1
n
1 ^A 1 X
Xi;j(i)Yi;
K
n i=1
o

^ := diag 0^(d1 +1) (d1 +1); 2 ; 0d3I d3I ;and ^ 2 is the di¤erence-based variance estimator that can be obtainedby reordering S2 with respect to X3.

7 August 201923

Fukushima University

A Two-Sample Alternative to IV

Two-Step Bias Correction forp

Hirukawa, Murtazashvili & Prokhorovn-Asymptotic Normality.

p

MSII can attain n-asymptotic normality when the number of matchingvariables is only one.

Hirukawa and Prokhorov (2018) propose a two-step estimator called thefully-modi…ed MSII (MSII-FM) estimator.
– In its second step, MSII-FM eliminates the second-order bias due tothe so-called matching discrepancy (Abadie and Imbens, 2006) asymptotically by means of a polynomial approximation.
– The estimator can achieve parametric convergence when the numberof matching variables is four or less.

7 August 201924

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Novel Features of PILS.
1. PILS does not rely on NNM.
It requires no bias correction that is a key ingredient in MSII(-FM).
The asymptotic variance of PILS does not much exceed that of OLS.
2. The asymptotic analysis for PILS explicitly incorporates discrete matchingvariables.
3. The role of excluded continuous matching variables for identi…cation isclari…ed.
Hirukawa and Prokhorov (2018) maintain the assumption that all common variables enter the regression and are used for both estimation andmatching.
7 August 201925

Fukushima University

5

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Finite-Sample Performance

Regressions.
[L] : Y = 0 + 1X1 + 2X2 + u; and[S ] : Y = 0 + 1X1 + (u + 2X2) := 0 + 1X1 + v:
where
0=

1=

2 = 1:

1=2 =(Data Generation.

1. Z; X3EC

2. X3ED
7 August 2019

iid

iid

U [ 2; 2].

Bernoulli (1=2)

26

1=2 wp 1=2.
1=2 wp 1=2

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov2
3. X1 = 1 + X3E+ X3ED + ( Z + ), whereC
f1:00; 0:70; 0:15g.

corr (X1; Z ) =

n

s

240 2
240 2 + 481

o8>< 0:58 for

0:44 for>: 0:11 for

4. X2 = h X3EC + X3ED + 2, where 2h ( x) =iid

N (0; 1) and

3 =2
2 =2
>
x
x
2x
>
>
>
>
>
jxj 1
>
>
: x + (5= ) (x= )N (0; 1) and

2= 3(5=2) f (2= )iid

27

N (0; 1).

2

= 1:00= 0:70 := 0:15

8
>
x
>
>
>
>
fexp (px) exp ( x)g = fexp (x) + exp ( x)g>>>< (3=4) x + 2 1

5. Y = 1 + X1 + X2 + u, where u
7 August 2019iid

1=2g ;[Model A][Model B][Model C]:[Model D][Model E]= 3=4 [Model F]

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Samples.
1. An unobserved complete sample:
S = f(Yi; X1i; X2i; X3i; Zi)gni=1 :2. Two observable samples:
S1 = f(Yi; X1i; X3i; Zi)gn
i=1 and S2 =nom
X2j ; X3j:
j=1

3. The matched sample via NNM with respect to X3:
S=

n

on
Yi; X1i; X2j1(i); : : : ; X2jK (i); X3i; X3j1(i); : : : ; X3jK (i); Zi;
i=1

where the number of matches K 2 f1; 2; 4; 8; 16; 32; 64; 128g.

Remark. Sample sizes are n = m 2 f500; 1000g, and the number of replications is 1000.
7 August 201928

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Estimators.

1. OLS*:
2. OLS-S:
3. IV-S:
4. MSOLS:
5. MSII:
6. PILS-E:
7. PILS-B:
8. PARA

Infeasible OLS estimator for [L] using S .
(Inconsistent) OLS estimator for [S ] using S1.
IV estimator for [S ] using S1.
(Inconsistent) MSOLS estimator for [L] using S .
MSII estimator for [L] using S .
PILS estimator for [L] using S1, S2 and the Epanechnikov kernel.
PILS estimator for [L] using S1, S2 and the beta kernel.
A parametric analogue to PILS for [L] using S1 and S2 by Fang,
Keane and Silverman (2008).

Bandwidth Selection.
^ = ^ X (log m=m)0:3 for the Epanechnikov kernel.
h^b = ^ U (log m=m)0:6 for the beta kernel, where U := (X + 2) =4 2 [0; 1].
^ = (log m=m)0:6 for the discrete kernel.
7 August 201929

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

A Few Words on PARA.
PARA takes the following two steps:
1. Run OLS for the regression of X2j on 1; X3EC ;j ; X3ED ;j using S2and obtain n parametric predictorsnon^ 2i = ^ 0 + ^ 1X3E ;i + ^ 2X3E ;i
X
C
D
i=1

using the OLS estimates (^ 0; ^ 1; ^ 2).

^ 2i .
2. Run OLS for the regression of Yi on 1; X1i; X
Under our Monte Carlo design, PARA is consistent if and only if h( ) is
odd.
– Consistency holds forh ( x) =7 August 2019(x
fexp (x)[Model A]:
exp ( x)g = fexp (x) + exp ( x)g [Model B]
30

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Model A: h (x) = x (n = m = 500).

1.00

0.70

0.15

Estimator
M ean
SD
RM SE
SE
CR
M ean
SD
RM SE
SE
CR
M ean
SD
RM SE
SE
CR

OLS*
0.99840.02240.02250.022395%
0.99850.02450.02460.024595%
0.99920.02700.02710.027295%

OLS-S
1.0611
0.04250.07440.04341.07400.04700.08760.04791.09200.05290.1061
0.0536IV-S
0.9956
0.07370.07380.073895%
0.99340.10630.10650.106295%
0.81162.72552.73205.591799%

MSOLS
1.02270.03730.04360.0391
1.02780.04130.04980.04291.03530.04680.05860.0476MSII
0.99750.04270.04280.039593%
0.99750.04760.04770.043592%
0.99790.05420.05430.048392%

PILS-E
1.00690.03470.03540.031893%
1.00880.03890.03990.035092%
1.01180.04450.04610.038890%

PILS-B
1.00500.03450.03490.031894%
1.00650.03860.03920.034992%
1.0091
0.04430.04520.038891%

PARA
0.99890.03180.03180.031795%
0.99920.03510.03510.034895%
1.00000.03940.03940.038694%

Notes: M ean = simulation average of the parameter estimate; SD = simulation standard
deviation of the parameter estimate; RM SE = root mean-squared error of the parameter
estimate; SE = simulation average of the standard error; and CR = coverage rate for thenominal 95% con…dence interval.
7 August 201931

Fukushima University

A Two-Sample Alternative to IV

Model E: h (x) = jxj

1.00

0.70

0.15

Estimator
M ean
SD
RM SE
SE
CR
M ean
SD
RM SE
SE
CR
M ean
SD
RM SE
SE
CR

OLS*
0.99810.02370.02380.023996%0.99810.02640.02640.026696%
0.99870.02980.02990.0302
96%

Hirukawa, Murtazashvili & Prokhorov

1 (n = m = 500).

OLS-S
1.22820.03320.2306
0.03371.27530.03650.2777
0.03681.33970.04080.34210.0405IV-S
0.99390.06200.06230.062895%0.99050.08960.09010.090696%
0.90131.32401.32761.915498%

MSOLS
1.1647
0.03620.16860.04771.20550.04050.20940.05791.26580.04630.26980.0751MSII
0.96660.13970.1436
0.117397%0.93080.46230.46750.378696%
0.97281.17351.17381.912495%

PILS-E
1.01160.03980.04140.039594%
1.01640.04670.04950.045993%
1.02620.0582
0.06380.055891%

PILS-B
1.00610.03970.04020.039694%
1.0091
0.04670.0476
0.046194%
1.01530.05860.0606
0.056294%

PARA
1.17700.03340.18010.0337
1.2166
0.03700.21980.0371
1.2726
0.04160.27580.0413

Notes: M ean = simulation average of the parameter estimate; SD = simulation standard
deviation of the parameter estimate; RM SE = root mean-squared error of the parameter
estimate; SE = simulation average of the standard error; and CR = coverage rate for thenominal 95% con…dence interval.
7 August 201932

Fukushima University

6

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Application to Estimation of Return to Schooling
We estimate the earnings equationlog (earnings) = 0 + 1education + 2ability+ 3experience + 4experience2+ 5married + 6black + 7south + 8urban + u
using the following two samples:
1. S1 from the 1972 wave of PSID (n = 2430).
– An ability measure (IQ score) is included (but this may not be agood ability proxy).
2. S2 from the CARD dataset provided by Wooldridge (2013) (m = 1102).
– Another ability measure (“Knowledge of the World of Work” or
KWW score) is included.

7 August 201933

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

IV Estimation.
Father’s education is used as the IV for individual’s education in the shortregression with ability omitted.
The sample correlation between these variables is 0.434.

Two-Sample Estimation.
List of common variables:(

included: education (C), married (D), black (D), south (D), urban (D)excluded: age (C), south while growing up (D)A single match (K = 1) and the Mahalanobis metric are adopted for
NNM.
7 August 201934

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

(1)

(2)

(3)

(4)

(5)

(6)

(7)

OLS*

OLS-S

IV-S

MSOLS

MSII-FM

PILS

PARA

0:0635(0:0056)
0:0809(0:0043)
0:0017(0:0001)
0:0313(0:0078)
0:3717(0:0539)
0:1302(0:0323)
0:0921(0:0288)
0:1363(0:0282)

0:0718(0:0053)
0:0818(0:0043)
0:0017(0:0001)

0:0953(0:0049)
0:0830(0:0044)
0:0017(0:0001)( )0:3793(0:0536)
0:1741(0:0316)
0:0983(0:0287)
0:1499(0:0284)( )0:3844(0:0535)
0:1249(0:0304)
0:0814(0:0293)
0:1278(0:0288)

0:0727(0:0059)
0:0826(0:0049)
0:0017(0:0001)
0:0012(0:0027)
0:3799(0:0535)
0:1849(0:0393)
0:0979(0:0286)
0:1538(0:0297)

0:0685(0:0080)
0:0766(0:0065)
0:0016(0:0001)
0:0029(0:0081)
0:3777(0:0535)
0:1504(0:0806)
0:0989(0:0289)
0:1404(0:0390)

0:0598(0:0070)
0:0825(0:0043)
0:0017(0:0001)
0:0075(0:0029)
0:3954(0:0541)
0:1749(0:0316)
0:1036(0:0287)
0:1541(0:0281)

0:2850(0:0891)
0:1170(0:0094)
0:0017(0:0001)
0:2444(0:0609)
0:2226(0:0605)
0:1819(0:0309)
0:1372(0:0299)
0:0138(0:0421)

No

No

No

Yes

Yes

Yes

Yes

2430

2430

2430

2430
1102

2430
1102

2430
1102

2430
1102

Estimation Results.

Education
Experience
Experience2
Ability
M arried
Black
South
U rbanData Combination?
Sample Size:nm

7 August 201935

Fukushima University

7

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

Concluding Remarks
The PILS estimation procedure is developed for models where endogenous
regressors enter the model due to an omitted variable.
– Convergence properties of PILS are explored.
Monte Carlo study con…rms attractive …nite-sample properties of PILS.
– It is numerically demonstrated that PILS is generally more e¢ cient thanIV and dominates other competing estimators (e.g., MSII) in terms ofmean squared error.
The curse of dimensionality in continuous matching variables is still a concern.
– We may adopt propensity score matching as a means of dimension
reduction using multiple matching variables.

7 August 201936

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov

References[1] Abadie, A., and G. W. Imbens (2006): “Large Sample Properties of Matching Estimators for Average Treatment E¤ects,” Econometrica, 74, 235 267.

[2] Black, B., M. Trainor, and J. E. Spencer (1999): “Wage Protection Systems, Segregation and Gender Pay Inequalities: West Germany, theNetherlands and Great Britain,” Cambridge Journal of Economics, 23,449 - 464.

[3] Card, D. (1995): “Using Geographic Variation in College Proximity to
Estimate the Return to Schooling,”in L. N. Christophides, E. K. Grant, and
R. Swidinsky (eds.), Aspects of Labour Market Behavior: Essays in Honour
of John Vanderkamp. Toronto: University of Toronto Press, 201 - 222.

7 August 201937

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov[4] Chen, S. X. (1999): “Beta Kernel Estimators for Density Functions,”Computational Statistics and Data Analysis, 31, 131 - 145.
[5] Fang, H., M. P. Keane, and D. Silverman (2008): “Sources of Advantageous Selection: Evidence from the Medigap Insurance Market,” Journalof Political Economy, 116, 303 - 350.
[6] Flavin, M., and S. Nakagawa (2008): “A Model of Housing in the Presenceof Adjustment Costs: A Structural Interpretation of Habit Persistence,”
American Economic Review, 98, 474 - 495.
[7] Hirukawa, M., and A. Prokhorov (2018): “Consistent Estimation of Linear
Regression Models Using Matched Data,” Journal of Econometrics, 203,
344 - 358.
[8] Mincer, J. A. (1974): Schooling, Experience and Earnings. New York: National Bureau of Economic Research.
7 August 201938

Fukushima University

A Two-Sample Alternative to IV

Hirukawa, Murtazashvili & Prokhorov[9] Newey, W. K. (1994): “Kernel Estimation of Partial Means and a General
Variance Estimator,” Econometric Theory, 10, 233 - 253.
[10] Pagan, A. (1984): “Econometric Issues in the Analysis of Regressions with
Generated Regressors,” International Economic Review, 25, 221 - 247.
[11] Racine, J. S., and Q. Li (2004): “Nonparametric Estimation of Regression
Functions with Both Categorical and Continuous Data,”Journal of Econometrics, 119, 99 - 130.
[12] Wooldridge, J. M. (2013): Introductory Econometrics: A Modern Approach, 5th Edition. Mason, OH: South-Western Cengage Learning.
[13] Zabalza, A., and J. L. Arrufat (1985): “The Extent of Sex Discriminationin Great Britain,” in A. Zabalza and Z. Tzannatos (eds.), Women andEqual Pay: The E¤ects of Legislation on Female Employment and Wages
in Britain, Cambridge, U.K.: Cambridge University Press, 70 - 96.
7 August 201939

Nonparametric estimation of density functions fromrepeated measurements
Daisuke Kurisu
School of Engineering,Tokyo Institute of Technology
Joint works with
Karun Adusumilli (University of Pennsylvania)
Taisuke Otsu (London School of Economics)
Yoon-Jae Whang (Seoul National University)

August 7, 2019 @福島大学

D. Kurisu

Inference from repeated measurents

August 7, 20191 / 38

Outline

Introduction
Estimation of density functions
Estimation of distribution functions under symmetric error
Estimation of distribution functions under asymmetric error
Conclusion

D. Kurisu

Inference from repeated measurents

August 7, 20192 / 38

Introduction

Kurisu, D and Otsu, T. (2019). On the uniform convergence of
deconvolution estimators from repeated measurements. submitted.
STICERD Econometrics Paper Series EM604.
Adusumilli, K., Kurisu, D., Otsu, T. and Whang, Y.-J. (2019).
Inference on distribution functions under measurement error.
R&R for J. Econometrics.
STICERD Econometrics Paper Series EM594.

D. Kurisu

Inference from repeated measurents

August 7, 20193 / 38

Estimation of density functions
We consider the following repeated measurements model.
Consider a bivariate i.i.d. sample {Yj,1 , Yj,2 }nj=1 of (Y1 , Y2 ), which is
generated by
Y1 = X + ϵ1 , Y2 = X + ϵ2 ,where (X , ϵ1 , ϵ2 ) are unobservables.
X: an error-free variable of interest, (ϵ1 , ϵ2 ) are measurement error forX.
We are interested in estimating the densities of X , ϵ1 , and ϵ2 .
Applications:(1) Household income and expenditure survey data (AKOW(2019))→ reinterview studies (to assess the credibility of the survey data).
(2) Medical data
→ measurements of patient’s status by diﬀerent methods (to assess
the compatibility of measurements).
D. Kurisu

Inference from repeated measurents

August 7, 20194 / 38

Estimation of density functions

Assumption M(ϵ1 , ϵ2 ) are independent copies of a random variable ϵ.
X is independent of (ϵ1 , ϵ2 ).
X and ϵ have square integrable Lebesgue densities fX and fϵ ,respectively.
The characteristic functions φX (·) = E [e i·X ] and φϵ (·) = E [e i·ϵ ]vanish nowhere, and E [ϵ] = 0.

D. Kurisu

Inference from repeated measurents

August 7, 20195 / 38

Estimation of density functions

Defineψ(u1 , u2 ) = E [e i(u1 Y1 +u2 Y2 ) ] = φX (u1 + u2 )φϵ (u1 )φϵ (u2 ).
Under the condition E [|Y1 |] < ∞, Kotlarski’s indentity gives us a explicit
identification formula of φX , that is(∫ u)∂ψ(0, u2 )/∂u1
φX (u) = exp
du2 .
ψ(0, u2 )0
The Li and Vuong’s (LV) estimator (Li and Vuong(’98, JMA)) of φX isdefined by taking its sample counterpart.

D. Kurisu

Inference from repeated measurents

August 7, 20196 / 38

Estimation of density functions(Definition of LV estimator):∫ uφ̂X (u) = exp0∂ ψ̂(0, u2 )/∂u1ψ̂(0, u2 )

du2 ,

∑
where ψ̂(u1 , u2 ) = n1 nj=1 e i(u1 Yj,1 +u2 Yj,2 ) and∑∂ ψ̂(u1 ,u2 )= n1 nj=1 iYj,1 e i(u1 Yj,1 +u2 Yj,2 ) .
∂u1
Based on this estimator, the density fX of X can be estimated by∫
1fˆX (x) =e −iux φ̂X (u)φK (hu)du,
2π R∫
where φK (u) = R e iux K (x)dx is the Fourier transform of a kernel
function K and h = hn is a sequence of positive numbers(bandwidths) such that hn → 0 as n → ∞.
D. Kurisu

Inference from repeated measurents

August 7, 20197 / 38

Estimation of density functions

Based on the expression φϵ (u) = ψ(0, u)/φX (u), the characteristic
function φϵ of ϵ can also be estimated byφ̂ϵ (u) =ψ̂(0, u).
φ̂X (u)The estimator fˆϵ of the density fϵ is given by replacing φ̂X in thedefinition of fˆX with φ̂ϵ .

D. Kurisu

Inference from repeated measurents

August 7, 20198 / 38

Estimation of density functions (CK)We also consider a regularized version of the LV estimator proposed in
Comte and Kappus(’15, JMA).
(Definition of CK estimator): Their main idea is to regularize φ̂X as
φ̃X (u) =

φ̃mod
X (u)
,
max{1, |φ̃mod
X (u)|}

where

∫ u
φ̃mod
X (u) = exp0∂ ψ̂(0, u2 )/∂u1du2
ψ̃(0, u2 )withψ̃(0, u2 ) =D. Kurisuψ̂(0, u2 ).
√min{1, n|ψ̂(0, u2 )|}Inference from repeated measurents

August 7, 20199 / 38

Estimation of density functions (CK)√The additional term, min{1, n|ψ̂(0, u2 )|} circumvents unfavorable
eﬀects caused by small values of the denominator.
The denominator of φ̃X , max{1, |φ̃mod
X (u)|} is introduced to improve
the quality of the estimator by imposing that the estimand is acharacteristic function, which should not take values lager than one.
The CK estimator f˜X of the density fX is defined by replacing φ̂Xwith φ̃X .
Also, φϵ can be estimated byφ̃ϵ (u) =φ̃X (u)ψ̂(0, u)√, where φ̄X (u) =.
φ̄X (u)min{1, n|φ̃X (u)|}The CK estimator f˜ϵ of fϵ is given by replacing φ̂X in the definition offˆX with φ̃ϵ .
D. Kurisu

Inference from repeated measurents

August 7, 201910 / 38

Estimation of density functions

To estimate densities, we need to choose the kernel function K , andimpose the following conditions.
Assumption K∫∫
The kernel function K∫satisfies R K (x)dx = 1, R x ℓ K (x)dx = 0 forℓ = 1, . . . , p − 1, and R |x|p K (x)dx < ∞ with a positive even integerp.
Also, φK (u) = 0 for any |u| > 1.
This assumption says that K is a p-th order kernel function.

D. Kurisu

Inference from repeated measurents

August 7, 201911 / 38

Estimation of density functions(Examples of kernel functions)
Kernel functions satisfying Assumption K is typically constructed by
specifying its Fourier transform φK . Examples of φK include(Finite order kernel):φK (u) = (1 − u 2 )k I{u ∈ [−1, 1]}, k ≥ p + 2.
(Infinite order kernel 1):



1 {
} if |u| ≤ c0 ,−b exp(−b/(|u|−c0 )2 )if c0 < |u| < 1,
φK (u) = exp(|u|−1)2

0
if 1 ≤ |u|,for 0 < c0 < 1 and b > 0.
(Infinite order kernel 2): The sinc kernel function.
φK (u) = I{u ∈ [−1, 1]}, K (x) = sin(x)/x.
D. Kurisu

Inference from repeated measurents

August 7, 201912 / 38

Estimation of density functions
We consider the following two scenarios for the characteristic functions φXand φϵ .
Assumption OS (ordinary smooth): For some positive constants
βx > 1, Cx ≥ cx , ωx , βϵ > 1, Cϵ ≥ cϵ , and ωϵ , it holds
cx |u|−βx ≤ |φX (u)| ≤ Cx |u|−βx , for all |u| ≥ ωx ,cϵ |u|−βϵ ≤ |φϵ (u)| ≤ Cϵ |u|−βϵ , for all |u| ≥ ωϵ .
The conditions βx , βϵ > 1 are introduced to guarantee the consistencyof the density functions.
We need to use the lower and upper bounds of the characteristic
functions to obtain suitable bounds of the stochastic anddeterministic bias terms of the estimators.
Examples: Laplace, Gamma.
D. Kurisu

Inference from repeated measurents

August 7, 201913 / 38

Estimation of density functions

Assumption SS (super smooth): For some positive constants ρx ,Cx ≥ cx , ωx , ρϵ , Cϵ ≥ cϵ , ωϵ , and some constants βx , βϵ ∈ R, it holds
cx |u|βx exp(−|u|ρx /µx ) ≤ |φX (u)| ≤ Cx |u|βx exp(−|u|ρx /µx ),for all |u| ≥ ωx ,cϵ |u| exp(−|u| /µϵ ) ≤ |φϵ (u)| ≤ Cϵ |u|βϵ exp(−|u|ρϵ /µϵ ),βϵ

ρϵfor all |u| ≥ ωϵ .
We use lower and upper bounds to control estimation errors of
estimators.
Examples: Normal, Mixture Normal, Cauchy.

D. Kurisu

Inference from repeated measurents

August 7, 201914 / 38

Estimation of density functions

Theorem (KO(2019))
Suppose that Assumption M holds true, and E [|Y1 |2+η ] < ∞ for someη > 0.
(i) : Under Assumption OS and n−1/2 Tn3βx +2βϵ +1 log Tn → 0 as n → ∞,it holds()−1/2 2βx +2βϵ +1
sup |φ̂X (u) − φX (u)| = Op n
Tn
log Tn ,u∈[−Tn ,Tn ]

sup

(
)
|φ̂ϵ (u) − φϵ (u)| = Op n−1/2 Tn3βx +2βϵ +1 log Tn .

u∈[−Tn ,Tn ]D. Kurisu

Inference from repeated measurents

August 7, 201915 / 38

Estimation of density functions

Theorem (conti.)(i) Additionally, suppose that Assumption K holds for p ≥ max{βx , βϵ }.
Set Tn = h−1 . Then)(sup |fˆX (x) − fX (u)| = Op n−1/2 h−2βx −2βϵ −2 log h−1 + hβx −1 ,|u|≤h−1()
sup |fˆϵ (x) − fϵ (u)| = Op n−1/2 h−3βx −2βϵ −2 log h−1 + hβϵ −1 .

|u|≤h−1

D. Kurisu

Inference from repeated measurents

August 7, 201916 / 38

Estimation of density functions
Theorem (conti.)(ii) Under Assumption SS andn−1/2 Tn1−3βx −2βϵ (log Tn ) exp (3Tnρx /µx + 2Tnρϵ /µϵ ) → 0 as n → ∞,it holds
sup|φ̂X (u) − φX (u)|u∈[−Tn ,Tn ]

(
( ρx
))
2Tn
2Tnρϵ
−1/2 1−2βx −2βϵ
= Op nTn(log Tn ) exp
+
,
µx
µϵ
sup |φ̂ϵ (u) − φϵ (u)|u∈[−Tn ,Tn ]

(
))
( ρx
2Tnρϵ
3Tn
−1/2 1−3βx −2βϵ
= Op nTn+.
(log Tn ) expµxµϵ

D. Kurisu

Inference from repeated measurents

August 7, 201917 / 38

Estimation of density functions
Theorem (conti.)(ii) Additionally, suppose that Assumption K holds and that there exists
0 < c ≤ 1 such that φK (x) = 1 for |x| ≤ c. Set Tn = h−1 . Thensup |fˆX (u) − fX (u)||u|≤h−1()−ρϵρx −ρx2h−ρx+ 2hµ− c µh−1/2 2βx +2βϵ −2
−1
ρ
/q−β
−1
x
x
µ
ϵ +h
x
= Op nh(log h )e x
e,sup |fˆϵ (u) − fϵ (u)||u|≤h−1

(
)
−ρϵ
2h−ρx
c ρx h−ρx+ 2hµρ /q−βx −1 − µx
ϵ +h x= Op n−1/2 h3βx +2βϵ −2 (log h−1 )e µx
e,where q = 1 when βx , βϵ > 0, and q > 1 when βx , βϵ ≤ 0.
D. Kurisu

Inference from repeated measurents

August 7, 201918 / 38

Estimation of density functions
The same uniform convergence results in (i) and (ii) hold true even if
we replace the LV estimator (φ̂X , φ̂ϵ , fˆX , fˆϵ ) with the CK estimator(φ̃X , φ̃ϵ , f˜X , f˜ϵ ).
To obtain the uniform consistency of φ̂ (or φ̃) on [−Tn , Tn ], it issuﬃcient to assumen−1/2 Tn2βx +2βϵ +1 log Tn → 0 under Assumption OS,
ρ2Tn x

n−1/2 Tn1−2βx −2βϵ (log Tn )e µx

ρ
2T ϵ
ϵ

+ µn

→ 0 under Assumption SS,

as n → ∞.
We can take c = 1 in Theorem (ii) if we use the sinc kernel function.
We can also derive uniform convergence rates when X is ordinarysmooth and ϵ is super smooth and vise versa since the decay rates of
the characteristic functions are not essential in our proofs.
D. Kurisu

Inference from repeated measurents

August 7, 201919 / 38

Estimation of density functions(Comparison with Li and Vuong (’98, JMA))
We note that LV established the uniform convergence rates of theirestimators under the assumption that both X and ϵ have bounded
support. Our theorem does not require such boundedness.
Also, the convergence rates obtained in our theorem are typicallyfaster thanin LV. For( those obtained) example, if we set
α/2(1+β+β)xϵTn = O (n/ log log n)with 0 < α < 1/2 as in Lemma
3.1 of LV, our Theorem (i) implies that(()α)− 1 +α−22(1+βx +βϵ )nsup |φ̂X (u) − φX (u)| = Op,
log log nu∈[−Tn ,Tn ]and this convergence rate is faster than that given in LV, i.e.,()− 1 +α2n. Similar comments apply to other cases.
log log n
D. Kurisu

Inference from repeated measurents

August 7, 201920 / 38

Estimation of density functions(Comparison with Bonhomme and Robin(’10, Rev. Econ. Stud.))
The convergence rates in our theorem are also faster than those given
in Bonhomme and Robin (’10). For example, under Assumption OS,
Bonhomme and Robin (’10, Theorem 1) implies that()
sup |φ̂X (u) − φX (u)| = Op n−1/2 Tn3βx +3βϵ +2 log Tn ,u∈[−Tn ,Tn ]

sup

)
(
|φ̂ϵ (u) − φϵ (u)| = Op n−1/2 Tn3βx +3βϵ +2 log Tn .

u∈[−Tn ,Tn ]It also should be noted that our assumption on (Y1 , Y2 ) is weaker
than an assumption in Bonhomme and Robin (’10) since we do notneed the existence of moment generating functions of Y12 and Y1 Y2 .
More precisely, the same convergence rate given in Lemma 1 of theirpaper can be obtained under weaker condition by proving a
multivariate version of Neumann and Reiss (’09, Theorem 4.1).
D. Kurisu

Inference from repeated measurents

August 7, 201921 / 38

Estimation of density functions

We obtained the same uniform convergence rates for the LV estimatorand the CK estimator.
It is open whether the LV estimator can achieve the L2 convergence
rate in Comte and Kappus(’15, JMA).
To control the L2 risk of the LV-type estimators which are defined by
the ratio of the (regularized) empirical averages, it seems crucial tointroduce some regularization as in Comte and Kappus(’15).

D. Kurisu

Inference from repeated measurents

August 7, 201922 / 38

Estimation of density functions

Multivariate version of Theorem 4.1 in Neumann and Reiß(’09, Bernoulli)
The relaxation of assumptions in the previous works can be achievedby a multivariate version of theorem 4.1 in Neumann and Reiß(’09)
on the normalized empirical characteristic function process.
Let {Yj = (Yj,1 , . . . , Yj,d )′ }nj=1 be Rd -valued i.i.d. random variables.
For t = (t1 , . . . , td )′ ∈ Rd , define
1 ∑ it·Yje,nnψ(t) = E [e it·Y ], ψ̂(t) =j=1

D. Kurisu

Inference from repeated measurents

August 7, 201923 / 38

Estimation of density functions
and for k = (k1 , . . . , kd )′ ∈ Nd , |k| =∑dj=1 kj ,√1 ∑ it·Yj
Cn (t) = √(e− E [e it·Y1 ]) = n(ψ̂(t) − ψ(t)),
n
n

j=1

n
1 ∑
∂ |k|
(k)
Cn (t) = √(e it·Yj − E [e it·Y ]),
kd
k1
n
·
·
·
∂t
∂t
1
d
j=1
[
]
(
)
(k)
(k)
E [∥Cn ∥L∞ (w ) ] = E sup w (∥t∥)|Cn (t)| ,t∈Rdwhere w (t) = (log(e + |t|))−1/2−δ for some δ > 0 and ∥ · ∥ is the
Euclidean norm on Rd .

D. Kurisu

Inference from repeated measurents

August 7, 201924 / 38

Estimation of density functions
For example,n
1 ∑ ∂ it·Yj
Cnk (t) = √(e− E [e it·Y ]), k = (0, 1)′ ,
∂t1
n
j=1


n

√ 1 ∑(g (Yj ) − E [g (Y1 )]) , g ∈ G1,k ,
= n
n

j=1
{
}
∂
d
G1,k := y 7→(cos(t · y ) + i sin(t · y )) : t ∈ R .
∂t1
From Corollary 19.35 in van der Vaart (1998), we have that
(√
)
(k)
2k
sup E [∥Cn ∥L∞ (w ) ] ≤ CJ[·]E [Y1 ], G1,kn≥1

for a universal constant C and where J[·] (δ, G ) =∫δ√log(N[·] (ϵ, G ))dϵ is0

the L2 -bracketing number of G1,k .
D. Kurisu

Inference from repeated measurents

August 7, 201925 / 38

Estimation of density functions

Lemma 3(KO(2019))[()2+η ]
∏d
(k
∨1/2)
j
Assume E< ∞ for some η > 0. Thenj=1 |Yj |(k)sup E [∥Cn ∥L∞ (w ) ] < ∞.
n≥1

We can also apply this result to examples in Bonhomme and Robin(’10) to
weaken their assumptions on the existence of m.g.f.s.

D. Kurisu

Inference from repeated measurents

August 7, 201926 / 38

Estimation of density functions
Since
(k)

∥Cn ∥L∞ (w ) ≥

√

∂ |k|

n sup
∥u∥≤Tn∂u1k1 · · · ∂udkd()
ψ̂(u) − ψ(u)inf w (∥u∥),∥u∥≤TnLemma 3 implies that
[
]
(k)
(
)
supn≥1 E [∥Cn ∥L∞ (w ) ]
∂ |k|
E
sup
ψ̂(u) − ψ(u)
≤ √
kd
k1
n inf ∥u∥≤Tn w (∥u∥)∥u∥≤Tn ∂u1 · · · ∂ud()= O n−1/2 log Tn .
Then we can show
sup∥u∥≤Tn(∂ |k|∂u1k1 · · · ∂udkd())
ψ̂(u) − ψ(u) = Op n−1/2 log Tnfrom Markov’s inequality.
D. Kurisu

Inference from repeated measurents

August 7, 201927 / 38

Estimator of distribution function under symmetric error(We use the notations in AKOW(2019))
Suppose we observe a random sample {Xj }nj=1 generated fromX = X ∗ + ϵ,
where X ∗ is an unobservable variable of interest and ϵ is its measurement
error. We assume that ϵ is independent of X ∗ .
More specifically, suppose that we observe
Xj,k = Xj∗ + ϵj,k for k = 1, 2 and j = 1, . . . , n.
where {Xj∗ }nj=1 and {ϵj,k }j,k are two independent sets of i.i.d. randomvariables.

D. Kurisu

Inference from repeated measurents

August 7, 201928 / 38

Estimator of distribution function under symmetric error
If the pdf fϵ of ϵ is “symmetric”, its Fourier transform fϵft can be estimatedby (Delaigle, Hall and Meister (’08, AoS))
1∑fˆϵft (ω) =cos{ω(Xj,1 − Xj,2 )}
n
n

1/2

j=1

We can estimate the cdf of FX ∗ by
F̂

X∗

(
)
∫
n
2
t − Xj,k
1 1 ∑∑
1 1 sin(ωu) K ft (ω)(t) = −L̂
, L̂(u) =
dω,
2 2n
h
2π −1 ω fˆϵft (ω/h)
j=1 k=1

where K ft is the Fourier transform of a kernel function K .

D. Kurisu

Inference from repeated measurents

August 7, 201929 / 38

Estimator of distribution function under symmetric error

We are concerned with approximation for the distribution of the
maximal deviation
Tn = sup |F̂X ∗ (t) − FX ∗ (t)|,t∈Twhere T is a compact interval specified by the researcher.
A direct use of such approximation is construction of the confidenceband for FX ∗ over T .

D. Kurisu

Inference from repeated measurents

August 7, 201930 / 38

Estimator of distribution function under symmetric error
Why distribution function?
1
2
Let {Xj }nj=1
and {Yj }nj=1
be two independent samples of X and Y . X
and Y are generated as

X = X ∗ + ϵ, Y = Y ∗ + δ,
To test the hypothesis of the (first-order) stochastic dominance
H0 : FX ∗ (t) ≤ FY ∗ (t) for all t ∈ T ,against the negation of H0 , density-based methods are cumbersome
to be handled.
In AKOW, we applied the stochastic order test to the Korea
Household Income and Expenditure Survey data to investigate welfarechanges of diﬀerent population subgroups between 2006 and 2012.
D. Kurisu

Inference from repeated measurents

August 7, 201931 / 38

Estimator of distribution function under symmetric error
Bootstrap approximation: the bootstrap version of F̂X ∗ is given by
(
)
♯
n
2
t − Xj,k
1
1 ∑∑♯F̂X ∗ (t) = −
L̂
,
2 2n
h
j=1 k=1♯
is randomly drawn from the pooled observations {Xj,k }where Xj,kwith equal weights. The bootstrap counterpart of Tn is obtained asTn♯ = sup |F̂X♯ ∗ (t) − F̂X ∗ (t)|.
t∈T

Let ĉα denote the (1 − α)-th quantile of the bootstrap statistic Tn♯ .
Theorem 1(AKOW(2019))
Under some regularity conditions and suppose either Assumption OS or SS
holds. Then,
P(Tn ≤ ĉα ) ≥ 1 − α + o(1).
D. Kurisu

Inference from repeated measurents

August 7, 201932 / 38

Estimator of distribution function under symmetric error
Bootstrap based implementation of stochastic dominance test: to testthe (first-order) stochastic dominance
H0 : FX ∗ (t) ≤ FY ∗ (t) for all t ∈ T ,against the negation of H0 , we can use a Kolmogorov-type test.
The test statistics and its bootstrap counterpart are given by
Dn1 ,n2 = sup{F̂X ∗ (t) − F̂Y ∗ (t)},
t∈T
{
}
♯
Dn1 ,n2 = sup F̂X♯ ∗ (t) − F̂Y♯ ∗ (t) − {F̂X ∗ (t) − F̂Y ∗ (t)} ,t∈T1
where F̂X♯ ∗ and F̂Y♯ ∗ are computed using bootstrap samples {Xi♯ }ni=1
♯ n2
n1
n2
and {Yi }i=1 from {Xi }i=1 and {Yi }i=1 , respectively.

D. Kurisu

Inference from repeated measurents

August 7, 201933 / 38

Estimator of distribution function under symmetric error

Let ĉαD denotes the (1 − α)-th quantile of the bootstrap statistics Dn♯1 ,n2 .
Theorem 6(AKOW(2019))
Suppose either Assumption OS or SS holds. Then, under some regularity
conditions, we have the following results.
Under H0 ,P(Dn1 ,n2 > ĉαD ) ≤ α + o(1).
Under the alternative H1 (i.e. H0 is false),
P(Dn1 ,n2 > ĉαD ) → 1, as n → ∞.

D. Kurisu

Inference from repeated measurents

August 7, 201934 / 38

Estimator of distribution function under asymmetric error
Let m < n. We can estimate FX ∗ under possibly asymmetric error
distribution by()
m2t − Xj,k
1
1 ∑∑AF̌X ∗ (t) = −
Ľ
,
2 2m
h
j=1 k=1
[
]
∫
e iωu
1 1 1
Ľ(u) =Im ft
K ft (ω)dω,π 0 ω
fˇϵ (ω/h)where fˇϵft is the Comte and Kappus’ estimator (computed from “n”
observations).
Let TnA = supt∈T |F̌XA∗ (t) − FX ∗ (t)|, and ĉαA denote the (1 − α)-thquantile of the bootstrap statistic TnA♯ = supt∈T |F̌XA♯∗ (t) − F̌XA∗ (t)|, where
(
)
♯
m
2
t − Xj,k
1 ∑∑1A♯F̌X ∗ (t) = −
Ľ
.
2 2m
h
j=1 k=1

D. Kurisu

Inference from repeated measurents

August 7, 201935 / 38

Estimator of distribution function under asymmetric error
Theorem 2(AKOW(2019))
Under some regularity conditions and suppose(i) Assumption OS holds true, and mnζ0 −1 h2βϵ −2βx −1 → 0 as n → ∞ forsome 0 < ζ0 < 1, or(ii) Assumption SS holds true, and mnζ1 −1 (log n)2βx +1 (log log n)2 → 0 as
n → ∞ for some 0 < ζ1 < 1.
Then
P(TnA ≤ ĉαA ) ≥ 1 − α − o(1).
We can construct an asymptotic confidence band for FX ∗ over T withlevel α as [F̂X ∗ (t) ± ĉα ] (Theorem 1) or [F̌XA∗ (t) ± ĉαA ] (Theorem 2).
Other applications and real data analysis are given in AKOW(2019).
D. Kurisu

Inference from repeated measurents

August 7, 201936 / 38

Conclusion
We investigated the uniform convergence of deconvolution estimatorsfrom repeated measurements.
We achieved faster uniform convergence rates than those obtained inexisting results under weaker assumptions.
Our results can be applied to uniform inference on distributionfunctions under possibly asymmetric measurement error.
(i) Kurisu, D and Otsu, T. (2019). On the uniform convergence of
deconvolution estimators from repeated measurements.
STICERD DP EM604. submitted.
(ii) Adusumilli, K., Kurisu, D., Otsu, T. and Whang, Y.-J. (2019).
Inference on distribution functions under measurement error.
R&R for J. Econometrics. STICERD DP EM594.
D. Kurisu

Inference from repeated measurents

August 7, 201937 / 38

References
Bonhomme, S. and Robin, J.-M. (2010). Generalized non-parametricdeconvolution with an application to earnings dynamics. Rev. Econ.
Stud. 77, 491-533.
Comte, F. and Kappus, J. (2015). Density deconvolution fromrepeated measurement without symmetry assumption on the errors.
J. Multivariate Anal. 140, 31-46.
Delaigle, A., Hall, P. and A. Meister (2008). On deconvolution withrepeated measurements. Ann. Statist. 36, 665-685.
Li, T. and Vuong, Q. (1998). Nonparametric estimation of themeasurement error model using multiple indicators. J. Multivariate
Anal. 65, 139-165.
Neumann, M. and Reiß, M. (2009). Nonparametric estimation forLévy processes from low-frequency observations. Bernoulli, 15,
223-248.
D. Kurisu

Inference from repeated measurents

August 7, 201938 / 38

Fukushima41

1 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

データサイエンス・サマーキャンプ福島佐藤整尚（東京⼤学・経済）(データサイエンスの業界でよく使われているJupyter Notebook を使って発表したいと思います。）

where with
is an
that

,
,
-matrix representing the first-order difference, andwith a diagonal matrix. That is,

To detect a trend component in , we use
only with
with smaller contains information about a longer cyclic component.

is an

,
-orthogonal matrix suchfor an estimation, because

2019/08/29 16:33

Fukushima41

2 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [33]: n <- length(shouhi)-1
Pmat <- sqrt(2/(n+0.5))*cos(outer(seq(n)-0.5,2*pi/(2*(n)+1)*(seq(n)-0.5)))z <- Pmat %*% diff(log(shouhi))parm(2,1)
tsplot(shouhi)#abline(v=82)
tsplot(z)abline(h=0,col=2)

trend series:

for :2019/08/29 16:33

Fukushima41

3 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [35]: parm(2,1)
m <- 2zT <- zzT[(m+1):n] <- 0
plot.ts(z,type="n")lines(zT)#abline(h=0,col=2)
plot.ts(shouhi)
Trend <- exp(cumsum(c(log(shouhi[1]),t(Pmat) %*% zT)))lines(Trend,col=2)
title(paste("Trend (m=",m,")"))

2019/08/29 16:33

Fukushima41

4 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

季節周波数付近のｚのみを逆変換すればよい

2019/08/29 16:33

Fukushima41

5 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [36]: parm(2,1)
m2 <- 2n2 <- as.integer(n/2)
zS <- z*0
zS[(n2-m2):(n2+m2)] <- z[(n2-m2):(n2+m2)]
zS[(n-m2):(n)] <- z[(n-m2):(n)]plot.ts(z,type="n")lines(zS)#abline(h=0,col=2)
plot.ts(log(shouhi),lty=2,ylab="",yaxt="n")par(new=T)
Seasonal <- cumsum(c(0,t(Pmat) %*% zS))
plot.ts(Seasonal,col=2)
title(paste("Seasonal (m2=",2,")"))

2019/08/29 16:33

Fukushima41

6 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

ｚにおいて季節性に対応する部分を０にすれば良い。

In [37]: parm(2,1)
m2 <- 2n2 <- as.integer(n/2)
zAD <- zzAD[c((n2-m2):(n2+m2),(n-m2):(n))] <- 0
plot.ts(z,type="n")lines(zAD)plot.ts(shouhi,lty=2)
SA <- cumsum(c(log(shouhi)[1],t(Pmat) %*% zAD))lines(exp(SA),col=2)
title(paste("S.A. (m2=",m2,")"))

2019/08/29 16:33

Fukushima41

7 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

ダミー変数をPで変換して回帰するが、トレンド成分か季節成分以外の周波数において推定する。
In [7]: #Additive Outlier
D.ao <- rep(0,n+1)
D.ao[82] <- 1parm(2,1)
plot.ts(D.ao)
z.ao <- Pmat %*% diff(D.ao)
plot.ts(z.ao)
abline(h=0,col=2)

2019/08/29 16:33

Fukushima41

8 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [10]: plot.ts(z,type="n")abline(h=0,col=3)
lines(zAD)lines(z.ao/10,col=2)
abline(v=10,lty=2)

2019/08/29 16:33

Fukushima41

9 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [14]: #Level Shift
D.ls <- rep(0,n+1)
D.ls[83:(n+1)] <- 1parm(2,1)
plot.ts(D.ls)
z.ls <- Pmat %*% diff(D.ls)
plot.ts(z.ls)
abline(h=0,col=2)

2019/08/29 16:33

Fukushima41

10 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [15]: plot.ts(z,type="n")abline(h=0,col=3)
lines(zAD)lines(z.ls/10,col=2)
abline(v=10,lty=2)

2019/08/29 16:33

Fukushima41

11 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [20]: #ver.2012.3.21
DDAIC <- 4
TCRATE <- 0.8 #^(12/frequency)
len <- length#dyn.load("decomp6r.dll")len <- lengthparm <- function(m,n) {par(mfrow=c(m,n))}tsplot <- function(x,y=NULL,...){
if(is.null(y)) {plot.ts(x,...)}
else {
plot.ts(cbind(x,y),plot.type="single",...,col=1:4)}}"regsiml"<function(data, reg = NULL, m1 = 2, period = 4, sorder = 1, log = 0,pb=2,pa=2,hh=-20){
pb <- pb*(sorder > 0)
pa <- pa*(sorder > 0)
odata <- data
if(log > 0) data <- log(data)n0 <- len(data)if(sorder > 0) {dx <- diff(data)dx <- c(0,rep(dx[1:period],pb),dx,rep(dx[(n0-period):(n0-1)],pa))
x.tmp <- cumsum(dx)x.tmp <- x.tmp + (data[1]-x.tmp[period*pb+1])
data <- x.tmp}
n <- len(data)n1 <- n-1
if(!is.null(reg)) {reg <- cbind(reg)k <- ncol(reg)}dx <- diff(data)if(hh < 0) {dy <- dx
iddd1 <- NULL
if(!is.null(reg)) {if(any(dimnames(reg)[[2]]=="ao")) {iddd <- (cbind(reg[,dimnames(reg)[[2]]=="ao"]) == 1)
if(pb > 0) {for(ii in seq(pb)) {iddd1 <- rbind(iddd1,iddd[(1:period)*(ii > 0),,drop=F])}}
iddd1 <- rbind(iddd1,iddd)

2019/08/29 16:33

Fukushima41

12 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [39]: #length(shouhi)#ts(shouhi,start=c(1994,1),frequency=4)zz <- x12siml(shouhi, reg = NULL, trend = 10, ar=0, ilog = 1, frequency= 4, start =c(1994, 1), iplot = T, sorder = 1,hh=0,ao=c(2014,1))[1] 0.3541929 0.3930816
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "0.0103""0.1115" "1.1298" "1" "9" "0.3155"
Mean Sum Sq R Squared F-value
Df 1 Df 2 Pr(>F)[1,] "0.0099""0.1657" "21.4431" "1" "108" "0"

2019/08/29 16:33

Fukushima41

13 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [40]: shouhi.o <- shouhi
shouhi.o[75] <- shouhi[75]*100
zz <- outlier(shouhi.o,c(2008,1,2015,4),type="ao",trend=17,ilog=0)
plot.ts(shouhi,type="n")lines(shouhi.o)lines(zz$trend,col=2)

2019/08/29 16:33

Fukushima41

14 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

[1] "make mat"[1] "make mat"[1] "make mat"[1] 5243278 6336719
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "166397.1" "0.1445" "2.7024" "1" "16" "0.1197"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "914164.6" "1e-04""0.0083" "1" "109" "0.9276"[1] "411.721884555715" "2008""1"[1] 5243278 4842078[1] "410.156568702371" "2008""2"[1] 5243278 4764030[1] 5243278 5171799[1] 5243278 5765884
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "175171.4" "0.0519" "0.8757" "1" "16" "0.3633"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "914098.2" "2e-04""0.0241" "1" "109" "0.8768"[1] 5243278 5245184
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "179900.2" "0""2e-04" "1" "16" "0.9884"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "913666.2" "0.0012" "0.1273" "1" "109" "0.722"[1] 5243278 5875371
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "172549.8" "0.0801" "1.3924" "1" "16" "0.2552"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "912168.8" "0.0044" "0.4858" "1" "109" "0.4873"[1] 5243278 4967688[1] "408.454530909937" "2009""4"[1] 5243278 3945295[1] "402.189922619381" "2010""1"[1] 5243278 5266835
Mean Sum Sq R Squared F-value
Df 1 Df 2 Pr(>F)[1,] "111472.9" "0.6161" "25.6726" "1" "16" "1e-04"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "913306.9" "0.002""0.2131" "1" "109" "0.6452"[1] "398.101638787292" "2010""2"[1] 5243278 7607983
Mean Sum Sq R Squared F-value
Df 1 Df 2 Pr(>F)[1,] "116852.2" "0.5781" "21.9241" "1" "16" "2e-04"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "911475.3" "0.006""0.6525" "1" "109" "0.421"[1] 5243278 4504819[1] 5243278 4714534[1] 5243278 5106308[1] 5243278 4000143[1] 5243278 6285875
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "145489.3" "0.346""8.464" "1" "16" "0.0102"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)[1,] "912557.3" "0.0036" "0.3926" "1" "109" "0.5322"[1] 5243278 5772231
Mean Sum Sq R Squared F-value
Df 1 Df 2 Pr(>F)[1,] "105169.8" "0.6582" "30.8174" "1" "16" "0"
Mean Sum Sq R Squared F-value Df 1 Df 2 Pr(>F)

2019/08/29 16:33

Fukushima41

15 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

2019/08/29 16:33

Fukushima41

16 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

2019/08/29 16:33

Fukushima41

17 / 17ﬁle:///C:/Users/Seisho Sato/AppData/Local/Packages/Microsoft...

In [ ]:

2019/08/29 16:33