データの端に興味があるときの分位点回帰

東京工業大学/株式会社Nospare リサーチャーの栗栖です．
この記事は最近の研究成果「データの端に興味がある場合のノンパラメトリック分位点回帰」(Kurisu and Otsu(2020) の続編です．この記事では前回紹介したノンパラメトリック極値分位点回帰の数値実験と実データ分析の例について紹介します．

ノンパラメトリック極値分位点回帰

数値実験の結果を紹介する前に，局所線形分位点回帰について復習しておきます.
確率変数 $(Y,X) \in \mathbb{R} \times \mathbb{R}^{d}$ に対し, $X = x$ の時の $Y$ の条件付分布関数を $F_{Y}(y|x)$ とします．
分位点回帰では，データ $(Y_{1},X_{1}),\cdots, (Y_{n}, X_{n}) \in \mathbb{R} \times \mathbb{R}^{d}$ が与えられたときに，これらを使って $Y$ の条件付き分位点

\theta_{\alpha}(x) = \inf_{y \in \mathbb{R}}\{y: F_{Y}(y|x) > \alpha\},\ \alpha \in (0,1)

を推定する問題を考えます．これは以下の最適化問題を解くことに相当します：

\hat{\beta}_{n}^{(\alpha)} = \arg\min_{\beta \in \mathbb{R}^{d+1}}\sum_{j=1}^{n}K\left({X_{j} - x_{0} \over \delta_{n}}\right)\rho_{\alpha}\left(Y_{j} - X_{j}(x_{0},\delta_{n})'\beta\right)　\tag{2}

ここで，$K: \mathbb{R}^{d} \to \mathbb{R}$ はカーネル関数，$x(x_{0},\delta_{n}) = (1,(x_{1}-x_{0,1})/\delta_n,\dots,(x_{d}-x_{0,d})/\delta_n)'$，$\beta = (\beta_{0},\dots,\beta_{d})'$, $\hat{\beta}_{n}^{(\alpha)} = (\hat{\beta}_{n,0}^{(\alpha)},\dots,\hat{\beta}_{n,d}^{(\alpha)})'$ であり，$\delta_n \to 0$ ($n \to \infty$) はバンド幅と呼ばれるチューニングパラメータです．この値により$X_j$，$j=1,\dots,n$ のうち，どの程度$x_0$に近い値をとるデータを推定に利用するかが決まります．Kurisu and Otsu (2020) で提案しているノンパラメトリック極値分位点回帰では"extremal quantile" の設定，即ち，推定したい分位点のレベル $\alpha$ がサンプルサイズ $n$ に依存して $0$ に収束するような状況

n\delta_n^d \alpha_n \to k \in (0,\infty),\ n \to \infty

を考えています．このような設定を考えることで，$\alpha$ が $0$ に近い場合の分析に興味がある場合に通常の分位点回帰の推定精度が悪くなってしまうという問題に対処できるようになります．

上記の最適化問題 (2) において，$\hat{\beta}_{n,0}^{(\alpha_n)}$ が推定したい条件付分位点 $\theta_{\alpha_n}(x_0)$ の推定値に対応しています．以下ではこの値を $\hat{\theta}_{\alpha_n}(x_0) := \hat{\beta}_{n,0}^{(\alpha_n)}$ と書くことにします．このとき，$k(m-1)> d+1$ を満たす $m>1$ に対して$(\hat{\theta}_{\alpha_n}(x_0) - \theta_{\alpha_n}(x_0))/(\hat{\theta}_{m\alpha_n}(x_0)-\hat{\theta}_{\alpha_n}(x_0))$ がある分布に収束することを利用して，Kurisu and Otsu (2020) ではサブサンプリング法を用いて $\theta_{\alpha_n}(x_0)$ の $100(1-\tau)$% 信頼区間を構成することを提案しています．

サンプリング法

$\theta_{\alpha_{n}}(x_{0})$ の $100(1-\tau)$% 信頼区間を構成するアルゴリズムは以下の通りです
(詳細については Kurisu and Otsu(2020) を適宜参照してください)：

(Step1) データ $(Y_{1},X_{1}),\cdots, (Y_{n},X_{n})$ から $\hat{\beta}_{n,0}^{(\alpha_{n})} = \hat{\theta}_{\alpha_{n}}(x_{0})$, $\hat{r}_{n}:=1/(\hat{\theta}_{m\alpha_n}(x_0)-\hat{\theta}_{\alpha_n}(x_0))$ を計算．

(Step 2) データを長さが $b$ $(<n)$ の$B_{n} = n-b+1$ 個のブロック に分割:
$\{(Y_{1},X_{1}),\cdots,(Y_{b},X_{b})\}$, $\cdots$, $\{(Y_{n-b+1},X_{n-b+1}),\cdots,(Y_{n},X_{n})\}$．

(Step 3) $\alpha_{b}$ を実際に推定したい分位点レベル $\alpha_{n}$ より大きくとり ($\alpha_{n} < \alpha_{b}$), $\alpha = \alpha_{b}$ として $i$ 番目のブロックの $b$ 個のデータを用いて最適化問題 $(2)$ を解いて局所線形回帰の定数項の推定量 $\hat{\beta}_{i,b,0}^{(\alpha_{b})}$ を求める (バンド幅を $\delta_b=(k/b\alpha_b)^{1/d}$，$k=n\delta_n^{d}\alpha_n$ として計算)．さらに各ブロックごとに $\hat{r}_{b}:=1/(\hat{\beta}_{i,b,0}^{(m\alpha_{b})}-\hat{\beta}_{i,b,0}^{(\alpha_{b})})$ も計算．また $n$ 個のデータ全て利用して $\hat{\beta}_{n,0}^{(\alpha_{b})}$ を計算する．

(Step 4) それぞれの $\hat{\beta}^{(\alpha_{b})}_{i,b,0}$ に対して $T_{n,b} = \hat{r}_{b}(\hat{\beta}^{(\alpha_{b})}_{i,b,0} - \beta_{n,0}^{(\alpha_{b})})$ を計算し， $T_{n,1},\dots,T_{n,B_{n}}$ の経験分布の $\tau$-分位点 ($\tau \in (0,1)$) を $q_{n}(\tau)$ とすると，$\theta_{\alpha_{n}}(x_{0})$ の $100(1-\tau)$% 信頼区間は

C_{1-\tau}(\alpha_n) = [\hat{\theta}_{\alpha_{n}}(x_{0}) - q_{n}(\tau/2)\hat{r}_{n},\hat{\theta}_{\alpha_{n}}(x_{0}) - q_{n}(1-\tau/2)\hat{r}_{n}]

で計算できます．上記の信頼区間の構成方法は $B_{n}$，$b_{n}$ が大きく，$b_{n}/n$ が小さい場合に妥当性が示せます．

数値実験

ここでは $\theta_{\alpha_n}(x_0)$ の $100(1-\tau)$% 信頼区間

C_{1-\tau}(\alpha_n) = [\hat{\theta}_{\alpha_n}(x_r) - q_n(\tau/2)\hat{r}_n, \hat{\theta}_{\alpha_n}(x_r) - q_n(1-\tau/2)\hat{r}_n]

の有限標本におけるパフォーマンスについて紹介します．

データ生成過程

データ生成過程として以下のモデルを考えます．

Y_j = 0.5\sin(X_j) + \sqrt{2.5 +0.5X_j^2}U_{\ast,j},\ j=1,\dots,n.

ここで，$n \in \{2000, 5000\}$,

• $X_j$ は $[-1,0]$ 上の一様分布に従う i.i.d. データ，
• $U_{\ast,j}$ は (i) 自由度 $3$ or $30$ の $t$ 分布または (ii) シェイプパラメータ $3$ or $30$ の Weibull 分布に従う i.i.d. データ

とします．

チューニングパラメータの選択

$\hat{\theta}_{\alpha_n}(x_0)$や信頼区間の計算にはバンド幅 $\delta_n$ を決める必要があり，この値によって数値実験やデータ分析の結果が変化します．この意味で $\delta_n$ の選択は重要な問題です．ここでは leave-one-out cross validation (LOOCV) を用いて $\delta_n$ を選択することにします．さらに，信頼区間の構成においてサブサンプルサイズ $b$ の選択も結果に影響を与えるチューニングパラメータなのですが，この値については，予備的な数値実験により，$b \approx n/10$ ととれば比較的安定した結果が得られることがわかったので，数値実験と実データ分析では $b = [n/10]$ ($n/10$ の整数部分) としています．その他のチューニングパラメータは $k = n\delta_n^d \alpha_n$, $\alpha_b=n\alpha_n/b$, $\delta_b = (k/b\alpha_b)^{1/d}$ としています．

以上の設定のもとで，$x_0 = -0.5$ における $\hat{\theta}_{\alpha_n}(x_r)$ を計算し，$\theta_{\alpha_n}(x_0)$ の 90%, 95% 信頼区間を繰り返し計算し，その経験被覆確率 (empirical coverage probability) を計算したものが以下の表です (繰り返し計算の回数は250 回)．またこの計算ではカーネル関数として biweight カーネル

K(w) = 
\begin{cases}
{15 \over 16}(1-|w|^2) & |w| \leq 1 \\
0 & |w|>1
\end{cases}

を使用し，$k = n\delta_n\alpha_n$, $b = 200$, $\alpha_b = n\alpha_n/b$ としています．この結果から，極値分位点の情報を取り込んで構成した信頼区間 $C_{1-\tau}(\alpha_n)$ を用いることで，分位点レベル $\alpha_n$ が小さい場合でも nominal coverage probability $1-\tau$ に近い被覆確率が得られることがわかります．

正規近似との比較　

上記のデータ生成過程の下で(i)提案手法と(ii)正規分布に基づく信頼区間のパフォーマンスを比較しましょう．Fan et al.(1994) では $\alpha \in (0,1)$ を固定した場合の局所線形分位点回帰に基づく$\theta_{\alpha}(x_0)$ の推定量 $\hat{\theta}_{\alpha_n}(x_0)$ 漸近正規性 (中心極限定理) が示されています (Theorem 3)．この結果をもとに $\theta_{\alpha}(x_0)$ の $100(1-\tau)$% 信頼区間を計算することができます：

C_{1-\tau}^{N}(\alpha) = \left[\hat{\theta}_{\alpha}(x_0) - q(\tau/2)\sqrt{\tau^2(x_0) \over n\delta_n}, \hat{\theta}_{\alpha}(x_0) + q(1-\tau/2)\sqrt{\tau^2(x_0) \over n\delta_n}\right],

ここで，$q(\tau)$ は標準正規分布の $(1-\tau)$-分位点，

\tau^2(x_0) = {\alpha(1-\alpha)\int K^2(w)dw \over f_{X}(x_0)g_{Y}^{2}(\theta_{\alpha}(x_0)|x_0)}

であり，$f_{X}$は$X$の密度関数，$g_Y(\cdot|x)$ は $X=x$ を与えたときの $Y$ の条件付き密度関数です．$f_X(x_0)$, $g_Y(\theta_{\alpha}(x_0)|x_0)$ はそれぞれカーネル密度推定を用いることで推定することができます．またそれぞれの推定量の計算に必要なバンド幅は LOOCV で選択することができます．下の表は $n = 2000$ の場合における $C_{1-\tau}(\alpha_n)$ と $C_{1-\tau}^{N}(\alpha_n)$ の経験被覆確率を計算した結果です (繰り返し計算の回数は250回)．この結果から，$\hat{\theta}_{\alpha}(x_0)$ の中心極限定理に基づく信頼区間の構成方法では $\alpha$ の値が小さい場合，適切な被覆確率が得られないことがわかります．

実データ分析

ノンパラメトリック極値分位点回帰の方法をポンド(GBP)とオーストラリアドル(AUD)の為替レートのデータ分析に応用してみましょう．下の図は2006年3月22日から2008年8月30日の間に3時間間隔で観測されたGBP-AUD 為替レート$\{R_j\}_{j=1}^{n+1}$, $n=5053$ を

Y_j = 100(\log R_{j+1} - \log R_{j}),\ j=1,\dots,5053

と変換してプロットした結果です．このデータがノンパラメトリックな定常AR(1)型の確率過程

Y_j = m(Y_{j-1}) + \sigma(Y_{j-1})U_{\ast,j}

に従っているとして $Y_{j-1} = 0$ における条件付き分位点 $\theta_{\alpha_n}(x_0)$, $\alpha_n \in \{0.01, 0.005\}$ を推定し，その 90%, 95% 信頼区間を計算してみましょう．

チューニングパラメータの選択

時系列データでは，一般には独立なデータに対して用いられる LOOCV が適用できない(未来のデータを使って過去のデータを validate してしまう)ため，ここでは LOOCV の代わりに Yu and Jones (1998) で提案された rule-of-thumb を用いてバンド幅を選択します．サブサンプルサイズ $b$ については数値実験と同様 $b = [n/10]$ としています．またカーネル関数や $k$, $\alpha_b$ も数値実験と同様に設定しています．

下の表が $\theta_{\alpha_n}(x_0)$ の推定値とその 90%, 95% 信頼区間の計算結果です．この結果を見ると，特に 95% 信頼区間は $\theta_{\alpha_n}(x_0)$ の推定値よりも小さな値の方向に広くなっていることがわかります．これは $r_n (\hat{\theta}_{\alpha_n}(x_0) - \theta_{\alpha_n}(x_0))$ の漸近分布が正規分布のような原点について左右対称な分布ではなく，非対称な分布であることに起因する結果です．

まとめ

この記事では「データの端に興味があるときのノンパラメトリック分位点回帰」の続編として，Kurisu and Otsu (2020) で提案したノンパラメトリック極値分位点回帰について，数値実験と実データ分析について紹介しました．株式会社Nospareでは極値統計学に限らず，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．

参考文献
[1] Fan, J., Hu, T.-C. and Truong, Y. K. (1994) Robust non-parametric function estimation. Scandinavian Journal of Statistics 21, 433-446.
[2] Kurisu, D and Otsu, T. (2020) Nonparametric inference for extremal conditional quantiles. R&R at Econometric Theory.
[3] Yu, K. and Jones, M. C. (1998) Local linear quantile regression. Journal of the American Statistical Association, 93, 228-237.