はじめに

最近、こちらの書籍を読みました。

計量経済学---ミクロデータ分析へのいざない

作者:末石直也
日本評論社

Amazon

わかりやすい理論の説明に合わせて、実際のデータ分析行う際の視点に立った解説が端々にあり、非常に良い本でした。

その中で、第8章のブートストラップの紹介の中で、Wild Bootstrapが扱われていました。

少し気になったので、簡単にまとめます。

残差Bootstrap

はじめに対比のために残差Bootstrap法について記しておく。

以下の線形回帰モデルを考える。

$\displaystyle{ Y = X_i \beta + u_i \\ E[u_i | X_i ] = 0 }$

残差Bootstrapでは、説明変数を固定して誤差項をランダムサンプルする。

手順は以下の通り。

モデル $Y_i ^\ast = X_i \hat \beta + \hat u_i (i=1,...,n)$ を考えます。ここで、 $\hat \beta$ は最小二乗推定量です。

残差 $\hat u_i$ から経験分布関数 $F_n$ を形成し、 $F_n$ からの無作為復元抽出により、ブートストラップ残差 $u_1 ^\ast,..., u^\ast_n$ を得る
1に基づき、ブートストラップ標本 $Y_i ^\ast = X_i \hat \beta +u_i ^\ast (i=1,...,n)$ を計算する
各ペア $(X_i, Y_i)$ を $(X_i,Y_i^\ast)$ で置き換え、ブートストラップ推定値 $\hat \beta^\ast$ 、および $( \hat \sigma ^\ast ) ^2$ を計算する

しかし、このBootstrap法は分散が不均一の場合にブートストラップ推定値の推定精度が低下することが知られている。

Wild Bootstrap

Wild Bootstrapはランダムサンプリングではなく選択分布(pick distribution)を利用することで、不均一分散の場合でも精度よく推定することができる。

Wild Bootstrapでは次のようなモデルを検討する。

$\displaystyle{ Y^*_i = X_i \hat \beta + f(\hat u_i)v_i^* }$

ここで $f()$ は残差 $\hat u_i$ の変換関数、 $v_i^*$ はある確率分布に従う変数である。

変換関数 $f()$ はいくつか提案されている。

$\displaystyle{ f(\hat u_i) = \sqrt{n/(n − k)} \hat u_i}\\ f(\hat u_i) = \frac{\hat u_i}{(1 − h_i)^{1/2}} \\ f(\hat u_i) = \frac{\hat u_i}{1 − h_i}$

ここで $h_i$ はハットマトリックス $P_X ≡ X ( X ^ T X ) ^ {−1} X ^ T$ の $i$ の対角要素である。

$v_i^*$ は、以下の性質を満たす確率分布であることが求められる。

$\displaystyle{ E[v_i^*] = 0\\ E[v_i^{*2} ] = 1\\ E[v_i^{*3} ] = 1 }$

MammenはLiuの条件を満たすような条件分布として以下のMammenの二項分布を提案している。

$\displaystyle{ P(v_i^* = \frac{1+1\sqrt{5}}{\sqrt{5}}) = \frac{\sqrt{5}-1}{2\sqrt{5}}\\ P(v_i^* = \frac{1-1\sqrt{5}}{\sqrt{5}}) = \frac{\sqrt{5}+1}{2\sqrt{5}} }$

二つの値しか取らないという状態であるが、今までのBootstrap法よりも精度良いサンプルリングが行えることがわかっている。

またDavidson and Flachaire (2008) は Rademacher分布を用いる方法を提案している。

$\displaystyle{ P(v_i^* = 1) = \frac{1}{2}\\ P(v_i^* = -1) = \frac{1}{2} }$

この他にもWild Bootstrap用の分布が提案されているようです。　

そして、なぜこのように二項分布で良いのか気になったので、Mammenの論文に軽く目を通したのですが、よくわらず...

もう少しちゃんと読んでみます。

projecteuclid.org

Rで実行

RでWild Bootstrap法を実行するにはlmboot packageで行うことができます。

www.rdocumentation.org

用いたデータは以下のモデルに従い生成されるデータとしました。

$\displaystyle{ Y = X_1 + X_2 +X_3 + \epsilon_i, \\ X_i \sim N(0, 1.2), \\ \epsilon_i \sim N(0, 0.6|X_2|) }$

誤差項の分散が不均一にならないような設定となっています。

基本的なWild Bootstrap法を実行するにはwild.boot関数を用いる。

wild_boot_model <-wild.boot(y ~ ., B = 1000, data = sample_df %>% select(-eps))

結果を確認する。

> quantile(wild_boot_model$bootEstParam[, "x1"], probs=c(.025, .975))
       2.5%       97.5% 
-0.08068078  0.11693998 
> quantile(wild_boot_model$bootEstParam[, "x2"], probs=c(.025, .975))
      2.5%      97.5% 
-0.1127211  0.2875900 
> quantile(wild_boot_model$bootEstParam[, "x3"], probs=c(.025, .975))
     2.5%     97.5% 
0.4301977 0.6467196

得られたBootstrap標本のヒストグラムを可視化してみます。

ggboot_histogram <- function(BD, X, bw = 1/50){
  tmp <- BD %>% as.data.frame() 
  
  wide <- max(tmp %>% pull(X))-min(tmp %>% pull(X))
  BN <- nrow(BD)
  coef <- quantile(tmp[,X], probs=c(.025, .975))
  ggplot(tmp,aes_(x = as.name(X)))+
    geom_histogram(binwidth = wide*bw)+
    geom_density(eval(bquote(aes(y = ..count..*.(wide*bw)))), fill = 1, alpha = 0.3)+
    geom_vline(xintercept=mean(tmp[,X]), linetype="dashed", colour = "red", size=.2) +
    ggtitle(paste0(X, " coeff by ", BN, " Bootstrap regression \n",
                   "Mean:", round(mean(tmp[,X]),4), "   2.5%:", round(coef[1],4), " - 97.5%:",round(coef[2],4)))+
    ylab("")
}

ggboot_histogram(BD=wild_boot_model$bootEstParam, "x3")