FWL定理(Frisch–Waugh–Lovell theorem) をシミュレーションで確認します。
以降の数式、導出は https://qiita.com/ykawakubo/items/327dbdb0379f87aa9248 を参照引用しています。
\(\mathbf{y}(n\times1)\) と \(\mathbf{X}(n\times k)\) の線形回帰モデルを考えます。 \[\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\epsilon},\quad \boldsymbol{\epsilon}\sim N\left(0,\sigma^2\mathbf{I}_n\right)\tag{1}\] ここで \(\mathbf{X}\) を \(\mathbf{X}_1(n\times k_1)\) と \(\mathbf{X}_2(n\times k_2)\) に分割します( \(k=k_1+k_2\) )。 \[\begin{eqnarray}\mathbf{y}&=&\mathbf{X}_1\boldsymbol{\beta}_1+\mathbf{X}_2\boldsymbol{\beta}_2+\boldsymbol{\epsilon}\\\boldsymbol{\beta}&=&\left(\boldsymbol{\beta}_1^{\mathsf{T}},\boldsymbol{\beta}_2^{\mathsf{T}}\right)^{\mathsf{T}}\end{eqnarray}\] よって、OLS推定量 \(\hat{\boldsymbol{\beta}}\) は \[\hat{\boldsymbol{\beta}}=\left(\hat{\boldsymbol{\beta}}_1^{\mathsf{T}},\hat{\boldsymbol{\beta}}_2^{\mathsf{T}}\right)^{\mathsf{T}}\] と表され、また、\(\mathbf{y}\) を\(\mathbf{X}_1\)に回帰させたときのOLS推定量 \(\hat{\boldsymbol{\beta}}_1\) は \[\hat{\boldsymbol{\beta}}_1=\left(\mathbf{X}_1^{\mathsf{T}}\mathbf{X}_1\right)^{-1}\mathbf{X}_1^{\mathsf{T}}\mathbf{y}\] であるため、残差ベクトル \(\tilde{\mathbf{y}}\) は \[\tilde{\mathbf{y}}=\mathbf{y}-\mathbf{X}_1\tilde{\boldsymbol{\beta}}_1=\mathbf{y}-\mathbf{X}_1\left(\mathbf{X}_1^{\mathsf{T}}\mathbf{X}_1\right)^{-1}\mathbf{X}_1^{\mathsf{T}}\mathbf{y}\] となり、ここで \[\mathbf{P}_1=\mathbf{X}_1\left(\mathbf{X}_1^{\mathsf{T}}\mathbf{X}_1\right)^{-1}\mathbf{X}_1^{\mathsf{T}}\] として、 \[\tilde{\mathbf{y}}=\left(\mathbf{I}_n-\mathbf{P}_1\right)\mathbf{y}\] と表す。
続いて、\(\mathbf{X}_2\) の各列毎に \(\mathbf{X}_1\)に回帰させたときの回帰係数ベクトル \(\left(k_1\times 1\right)\) を \(k_2\) 列並べた行列 \(\left(k_1\times k_2\right)\) は
\[\mathbf{X}_1\left(\mathbf{X}_1^{\mathsf{T}}\mathbf{X}_1\right)^{-1}\mathbf{X}_1^{\mathsf{T}}\] となるため、残差行列 \(\tilde{\mathbf{X}}_2\) は \[\tilde{\mathbf{X}}_2 = \mathbf{X}_2 - \mathbf{X}_1\left(\mathbf{X}_1^\mathsf{T} X_1\right)^{-1} \mathbf{X}_1^\mathsf{T} \mathbf{X}_2 = \left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\] となる。
さらに、\(\tilde{\mathbf{y}}\) を \(\tilde{\mathbf{X}}_2\) に回帰させたときのOLS推定量を \(\tilde{\boldsymbol{\beta}}_2\) とすると \[\begin{eqnarray}\tilde{\boldsymbol{\beta}}_2&=&\left(\tilde{\mathbf{X}}_2^\mathsf{T}\tilde{\mathbf{X}}_2\right)^{-1}\tilde{\mathbf{X}}_2^\mathsf{T}\tilde{\mathbf{y}}\\&=&\left(\left(\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)^\mathsf{T}\left(\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)\right)^{-1}\left(\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)^\mathsf{T}\tilde{\mathbf{y}}\end{eqnarray}\] となるが、ここで \[\mathbf{M}=\mathbf{I}_n-\mathbf{P}_1\] とすると \(\mathbf{M}\) は対称冪等行列であるため \[\mathbf{M}^\mathsf{T}\mathbf{M}=\mathbf{M}\mathbf{M}=\mathbf{M}\]こちらを参照してください。
となり、\(\mathbf{I}_n-\mathbf{P}_1\) が対称冪等行列であることを利用すると \(\tilde{\boldsymbol{\beta}}_2\) は \[\begin{eqnarray}\tilde{\boldsymbol{\beta}}_2&=&\left(\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n - \mathbf{P}_1\right)^\mathsf{T}\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)^{-1}\left(\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)^\mathsf{T}\tilde{\mathbf{y}}\\&=&\left(\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)^{-1}\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n - \mathbf{P}_1\right)^\mathsf{T}\tilde{\mathbf{y}}\\&=&\left(\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n - \mathbf{P}_1\right)\mathbf{X}_2\right)^{-1}\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n - \mathbf{P}_1\right)\tilde{\mathbf{y}}\end{eqnarray}\tag{2}\] となる。
モデル(1)のOLS推定量 \(\hat{\boldsymbol{\beta}}\) は \[\begin{eqnarray}\hat{\boldsymbol{\beta}}&=&\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)^{-1}\mathbf{X}^{\mathsf{T}}\mathbf{y}\\\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)\hat{\boldsymbol{\beta}}&=&\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)^{-1}\mathbf{X}^{\mathsf{T}}\mathbf{y}\\\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)\hat{\boldsymbol{\beta}}&=&\mathbf{I}\mathbf{X}^{\mathsf{T}}\mathbf{y}\\\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)\hat{\boldsymbol{\beta}}&=&\mathbf{X}^{\mathsf{T}}\mathbf{y}\\\end{eqnarray}\] であることから \[\begin{bmatrix}\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1&\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_2\\\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_1&\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_2\end{bmatrix}\begin{bmatrix}\hat{\boldsymbol{\beta}}_1\\\hat{\boldsymbol{\beta}}_2\end{bmatrix}=\begin{bmatrix}\mathbf{X}_1^\mathsf{T}\mathbf{y}\\\mathbf{X}_2^\mathsf{T}\mathbf{y}\end{bmatrix}\tag{3}\] と展開され、(3)の1行目は \[\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\hat{\boldsymbol{\beta}}_1+\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_2\hat{\boldsymbol{\beta}}_2=\mathbf{X}_1^\mathsf{\mathsf{T}}\mathbf{y}\] であるので \[\begin{eqnarray}\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\hat{\boldsymbol{\beta}}_1&=&\mathbf{X}_1^\mathsf{T}\mathbf{y}-\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_2\hat{\boldsymbol{\beta}}_2\\\left(\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\right)^{-1}\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\hat{\boldsymbol{\beta}}_1&=&\left(\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\right)^{-1}\left(\mathbf{X}_1^\mathsf{T}\mathbf{y}-\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_2\hat{\boldsymbol{\beta}}_2\right)\\\mathbf{I}_n\hat{\boldsymbol{\beta}}_1&=&\left(\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\right)^{-1}\mathbf{X}_1^\mathsf{T}\mathbf{y}-\left(\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\right)^{-1}\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_2\hat{\boldsymbol{\beta}}_2\\\hat{\boldsymbol{\beta}}_1&=&\left(\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\right)^{-1}\mathbf{X}_1^\mathsf{T}\left(\mathbf{y}-\mathbf{X}_2\hat{\boldsymbol{\beta}}_2\right)\end{eqnarray}\] となり、求めた \(\hat{\boldsymbol{\beta}}_1\) を(3)の2行目である \[\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_1\hat{\boldsymbol{\beta}}_1+\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_2\hat{\boldsymbol{\beta}}_2=\mathbf{X}_2^\mathsf{T}\mathbf{y}\] に代入すると \[\begin{eqnarray}\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_1\left(\mathbf{X}^{\mathsf{T}}_1\mathbf{X}_1\right)^{-1}\mathbf{X}_1^\mathsf{T}\left(\mathbf{y}-\mathbf{X}_2\hat{\boldsymbol{\beta}}_2\right)+\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_2\hat{\boldsymbol{\beta}}_2&=&\mathbf{X}_2^\mathsf{T}\mathbf{y}\\\mathbf{X}^{\mathsf{T}}_2\mathbf{P}_1\left(\mathbf{y}-\mathbf{X}_2\hat{\boldsymbol{\beta}}_2\right)+\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_2\hat{\boldsymbol{\beta}}_2&=&\mathbf{X}_2^\mathsf{T}\mathbf{y}\\\mathbf{X}^{\mathsf{T}}_2\mathbf{P}_1\mathbf{y}-\mathbf{X}^{\mathsf{T}}_2\mathbf{P}_1\mathbf{X}_2\hat{\boldsymbol{\beta}}_2+\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_2\hat{\boldsymbol{\beta}}_2&=&\mathbf{X}_2^\mathsf{T}\mathbf{y}\\-\mathbf{X}^{T}_2\mathbf{P}_1\mathbf{X}_2\hat{\boldsymbol{\beta}}_2+\mathbf{X}^{\mathsf{T}}_2\mathbf{X}_2\hat{\boldsymbol{\beta}}_2&=&\mathbf{X}_2^\mathsf{T}\mathbf{y}-\mathbf{X}^{\mathsf{T}}_2\mathbf{P}_1\mathbf{y}\\\left(\mathbf{X}^{\mathsf{T}}_2\left(\mathbf{I}_n-\mathbf{P}_1\right)\mathbf{X}_2\right)\hat{\boldsymbol{\beta}}_2&=&\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n-\mathbf{P}_1\right)\mathbf{y}\\\hat{\boldsymbol{\beta}}_2&=&\left(\mathbf{X}^{\mathsf{T}}_2\left(\mathbf{I}_n-\mathbf{P}_1\right)\mathbf{X}_2\right)^{-1}\mathbf{X}_2^\mathsf{T}\left(\mathbf{I}_n-\mathbf{P}_1\right)\mathbf{y}\end{eqnarray}\tag{4}\] と \(\tilde{\boldsymbol{\beta}}_2\) が求められ、式{2}の \(\tilde{\boldsymbol{\beta}}_2\) と式(4)の \(\hat{\boldsymbol{\beta}}_2\) は一致していることが確認でき、FWL定理 が証明されます。
続いて実際の数値を利用して FWL定理 を確認します。
始めに線形回帰モデル(1)のサンプルデータを作成します。
library(dplyr)
set.seed(20240713)
k <- 10
n <- 20
X <- rnorm(n = k * n, mean = 0, sd = 1) %>% matrix(nrow = n)
X[, 1] <- 1
beta <- rnorm(n = k) %>% matrix(ncol = 1)
e <- rnorm(n = k, mean = 0, sd = 2)
y <- X %*% beta + ey %>% tibble::tibble()# A tibble: 20 × 1
.[,1]
<dbl>
1 -0.337
2 4.24
3 7.74
4 1.79
5 0.510
6 0.632
7 2.68
8 0.198
9 0.915
10 3.76
11 3.47
12 1.03
13 3.38
14 0.248
15 -1.88
16 2.15
17 6.38
18 2.34
19 -1.60
20 -3.64 X %>% tibble::tibble()# A tibble: 20 × 1
.[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.524 0.374 -0.546 -0.420 1.59 -0.433 -0.385 -0.0414
2 1 1.48 0.162 1.28 0.0728 -0.363 2.13 -0.384 -0.826
3 1 -0.0723 1.81 1.43 1.05 -0.905 3.05 -0.720 -1.14
4 1 -0.762 0.0781 1.77 -0.531 -0.591 1.67 1.13 1.45
5 1 -0.156 -1.49 1.22 -0.302 0.0744 -1.08 -0.00726 0.0734
6 1 0.947 -0.242 -1.70 0.917 1.10 -0.690 0.826 -1.20
7 1 1.20 -0.324 0.0497 -0.936 0.568 0.0349 0.224 0.478
8 1 0.00752 0.586 -0.311 -0.633 1.29 1.24 -1.29 -0.932
9 1 0.0761 -0.378 -1.25 0.567 0.147 0.387 -0.932 -0.485
10 1 -0.636 -1.23 0.109 -0.631 -0.968 0.243 -1.35 0.956
11 1 -0.408 0.873 -2.43 -0.727 -1.19 0.500 -1.01 -0.908
12 1 -0.381 0.355 0.951 1.44 0.536 -1.50 0.976 0.753
13 1 -1.21 -0.698 -1.05 0.324 -0.178 -0.689 -1.18 -0.745
14 1 -1.78 0.370 -1.02 1.41 -0.403 -1.61 -0.172 -1.15
15 1 -0.590 1.21 1.27 -0.560 0.352 -1.03 -0.473 0.898
16 1 -1.74 -0.463 -1.28 -2.17 -1.14 1.82 0.113 -0.0156
17 1 1.42 0.0629 0.159 1.26 -1.00 -0.941 0.616 -0.185
18 1 0.205 0.0283 0.759 -0.794 -0.936 1.68 1.29 -1.01
19 1 0.0924 0.946 -1.44 0.968 -0.140 -0.763 -0.328 -1.01
20 1 -0.260 -0.790 -0.599 -0.323 0.961 -2.53 1.86 0.708
# ℹ 1 more variable: .[10] <dbl>線形回帰モデル(1)の \(\hat{\boldsymbol{\beta}}\) とその部分ベクトル \(\hat{\boldsymbol{\beta}}_2\) を求めます。
始めに \(k_1=(1,\cdots,7),\quad k_2=(8,9,10)\) として分割します。
split.col <- 7
hat.beta <- solve(t(X) %*% X) %*% t(X) %*% y
hat.beta2 <- hat.beta[(split.col + 1):k, , drop = F]hat.beta %>% tibble::tibble()# A tibble: 10 × 1
.[,1]
<dbl>
1 1.76
2 1.00
3 -0.624
4 -0.322
5 1.66
6 -1.68
7 1.19
8 -0.0301
9 0.437
10 0.914 hat.beta2 %>% tibble::tibble()# A tibble: 3 × 1
.[,1]
<dbl>
1 -0.0301
2 0.437
3 0.914 続いて、残差回帰推定量 \(\tilde{\boldsymbol{\beta}}_2\) を求め、\(\hat{\boldsymbol{\beta}}_2\) と比較します。
X1 <- X[, 1:split.col]
X2 <- X[, (split.col + 1):k]
tilde.y <- y - X1 %*% solve(t(X1) %*% X1) %*% t(X1) %*% y
tilde.X2 <- X2 - X1 %*% solve(t(X1) %*% X1) %*% t(X1) %*% X2
tilde.beta2 <- solve(t(tilde.X2) %*% tilde.X2) %*% t(tilde.X2) %*% tilde.y
data.frame(hat.beta2, tilde.beta2) hat.beta2 tilde.beta2
1 -0.03011793 -0.03011793
2 0.43684544 0.43684544
3 0.91385313 0.91385313\(\hat{\boldsymbol{\beta}}_2\) と \(\tilde{\boldsymbol{\beta}}_2\) は一致しています。
もう一つサンプルデータを作成します。
今度は \(k_1=(1,\cdots5),\quad k_2=(6,\cdots10)\) として分割します。
split.col <- 5
hat.beta2 <- hat.beta[(split.col + 1):k, , drop = F]
X1 <- X[, 1:split.col]
X2 <- X[, (split.col + 1):k]
tilde.y <- y - X1 %*% solve(t(X1) %*% X1) %*% t(X1) %*% y
tilde.X2 <- X2 - X1 %*% solve(t(X1) %*% X1) %*% t(X1) %*% X2
tilde.beta2 <- solve(t(tilde.X2) %*% tilde.X2) %*% t(tilde.X2) %*% tilde.y
data.frame(hat.beta2, tilde.beta2) hat.beta2 tilde.beta2
1 -1.67833743 -1.67833743
2 1.19270432 1.19270432
3 -0.03011793 -0.03011793
4 0.43684544 0.43684544
5 0.91385313 0.91385313こちらも一致しています。
以上です。
