Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

FWL定理

証明

残差への回帰Y~=X~1β1\tilde{Y} = \tilde{X}_1 \beta_1は

Y−X2γ^=(X1−X2δ^)β1Y - X_2 \hat{\gamma} = (X_1 - X_2\hat{\delta}) \beta_1
(3)

であり、

γ^=(X2TX2)−1X2TYδ^=(X2TX2)−1X2TX1\hat{\gamma} = (X_2^T X_2)^{-1} X_2^T Y\\ \hat{\delta} = (X_2^T X_2)^{-1} X_2^T X_1
(4)

を代入するとそれぞれ

Y−X2γ^=Y−X2(X2TX2)−1X2TY=[I−X2(X2TX2)−1X2T]YX1−X2δ^=X1−X2(X2TX2)−1X2TX1=[I−X2(X2TX2)−1X2T]X1\begin{align} Y - X_2 \hat{\gamma} &= Y - X_2 (X_2^T X_2)^{-1} X_2^T Y\\ &= [I - X_2 (X_2^T X_2)^{-1} X_2^T] Y\\ X_1 - X_2\hat{\delta} &= X_1 - X_2 (X_2^T X_2)^{-1} X_2^T X_1\\ &= [I - X_2 (X_2^T X_2)^{-1} X_2^T] X_1 \end{align}
(5)

となる。

これらをそれぞれ被説明変数、説明変数として最小二乗推定すると

β^1=[(X1−X2δ^)T(X1−X2δ^)]−1(X1−X2δ^)T(Y−X2)={X1T[I−X2(X2TX2)−1X2T]X1}−1X1[I−X2(X2TX2)−1X2T]Y\begin{align} \hat{\beta}_1 &= [(X_1 - X_2\hat{\delta})^T (X_1 - X_2\hat{\delta})]^{-1} (X_1 - X_2\hat{\delta})^T (Y - X_2) \\ &= \{ X_1^T [I - X_2 (X_2^T X_2)^{-1} X_2^T] X_1 \}^{-1} X_1 [I - X_2 (X_2^T X_2)^{-1} X_2^T] Y \\ \end{align}
(6)

（I−X2(X2TX2)−1X2TI - X_2 (X_2^T X_2)^{-1} X_2^Tは冪等で対称な行列のため2行目で計算が簡略化できている）

次に、すべての説明変数XXをYYに回帰したとき

Y=X1β1+X2β2Y = X_1 \beta_1 + X_2 \beta_2
(7)

のX1X_1の係数β1\beta_1の推定量を求める。二乗誤差の最小化の解

(β^1,β^2)=arg⁡min⁡ ⁡β1,β2(Y−X1β1+X2β2)T(Y−X1β1+X2β2)(\hat{\beta}_1, \hat{\beta}_2) = \operatorname*{\arg \min\ }_{\beta_1, \beta_2} (Y - X_1 \beta_1 + X_2 \beta_2)^T (Y - X_1 \beta_1 + X_2 \beta_2)
(8)

の1階の条件は

X1TY−X1TX1β^1−X1TX2β^2=0X2TY−X2TX1β^1−X2TX2β^2=0X_1^T Y - X_1^T X_1 \hat{\beta}_1 - X_1^T X_2 \hat{\beta}_2 = 0\\ X_2^T Y - X_2^T X_1 \hat{\beta}_1 - X_2^T X_2 \hat{\beta}_2 = 0
(9)

となる。β^2\hat{\beta}_2を消すためには、2つ目の式に左から−X1TX2(X2TX2)−1-X_1^T X_2(X_2^T X_2)^{-1}を掛けて

−X1TX2(X2TX2)−1(X2TY−X2TX1β^1−X2TX2β^2)=−X1TX2(X2TX2)−1X2TY+X1TX2(X2TX2)−1X2TX1β^1+X1TX2(X2TX2)−1X2TX2β^2=−X1TX2(X2TX2)−1X2TY+X1TX2(X2TX2)−1X2TX1β^1+X1TX2β^2\begin{align} &-X_1^T X_2(X_2^T X_2)^{-1} (X_2^T Y - X_2^T X_1 \hat{\beta}_1 - X_2^T X_2 \hat{\beta}_2) \\ &= - X_1^T X_2(X_2^T X_2)^{-1} X_2^T Y + X_1^T X_2(X_2^T X_2)^{-1} X_2^T X_1 \hat{\beta}_1 + X_1^T X_2(X_2^T X_2)^{-1} X_2^T X_2 \hat{\beta}_2 \\ &= - X_1^T X_2(X_2^T X_2)^{-1} X_2^T Y + X_1^T X_2(X_2^T X_2)^{-1} X_2^T X_1 \hat{\beta}_1 + X_1^T X_2 \hat{\beta}_2 \end{align}
(10)

これを第1式に足すと

X1TY−X1TX1β^1−X1TX2β^2−X1TX2(X2TX2)−1X2TY+X1TX2(X2TX2)−1X2TX1β^1+X1TX2β^2=X1TY−X1TX2(X2TX2)−1X2TY−[X1TX1−X1TX2(X2TX2)−1X2TX1]β^1=0X_1^T Y - X_1^T X_1 \hat{\beta}_1 - X_1^T X_2 \hat{\beta}_2 - X_1^T X_2(X_2^T X_2)^{-1} X_2^T Y + X_1^T X_2(X_2^T X_2)^{-1} X_2^T X_1 \hat{\beta}_1 + X_1^T X_2 \hat{\beta}_2 \\ = X_1^T Y - X_1^T X_2(X_2^T X_2)^{-1} X_2^T Y - [X_1^T X_1 - X_1^T X_2(X_2^T X_2)^{-1} X_2^T X_1] \hat{\beta}_1 = 0
(11)

となる。

これを変形すると

β^1=[X1TX1−X1TX2(X2TX2)−1X2TX1]−1X1TY−X1TX2(X2TX2)−1X2TY={X1T[I−X1TX2(X2TX2)−1X2T]X1}−1X1T[I−X2(X2TX2)−1X2T]Y\begin{align} \hat{\beta}_1 &= [X_1^T X_1 - X_1^T X_2(X_2^T X_2)^{-1} X_2^T X_1]^{-1} X_1^T Y - X_1^T X_2(X_2^T X_2)^{-1} X_2^T Y \\ &= \{X_1^T [I - X_1^T X_2(X_2^T X_2)^{-1} X_2^T] X_1\}^{-1} X_1^T [I - X_2(X_2^T X_2)^{-1} X_2^T] Y \end{align}
(12)

となり、残差回帰の解と一致する。

数値例

import numpy as np

n = 100
d = 4
k = 1
np.random.seed(0)
beta = np.array(list(range(1, d + 1)))
print(f"{beta=}")
X = np.random.uniform(size=(n, d))
e = np.random.normal(size=n) * 0.1
Y = X @ beta + e

X1 = X[:, :k]
X2 = X[:, k:]
beta1 = beta[:k]
beta2 = beta[k:]
assert np.allclose(Y, X1 @ beta1 + X2 @ beta2 + e)
beta=array([1, 2, 3, 4])

# I - X_2 (X_2^T X_2)^{-1} X_2^T は冪等で対称な行列
I = np.eye(n)
M = (I - X2 @ np.linalg.inv(X2.T @ X2) @ X2.T)
assert np.allclose(M.T, M)
assert np.allclose(M, M @ M)
delta = np.linalg.inv(X2.T @ X2) @ X2.T @ X1
gamma = np.linalg.inv(X2.T @ X2) @ X2.T @ Y

Y_tilde = (Y - X2 @ gamma)
X_tilde = (X1 - X2 @ delta)

# 残差回帰
beta_hat2 = np.linalg.inv(X_tilde.T @ X_tilde) @ X_tilde.T @ Y_tilde
beta_hat2.round(1)
array([1.])
# 通常の線形回帰
beta_hat1 = np.linalg.inv(X.T @ X) @ X.T @ Y
beta_hat1.round(1)
array([1. , 2. , 2.9, 4. ])

残差の作図

Source
import matplotlib.pyplot as plt
import japanize_matplotlib

fig, axes = plt.subplots(figsize=[12, 3], ncols=3)
fig.subplots_adjust(wspace=0.3)
axes[0].scatter(X1, Y)
axes[0].set(ylabel=r"$Y$", xlabel=r"$X_1$", title="元の散布図")

axes[1].scatter(X_tilde, Y_tilde)
axes[1].set(ylabel=r"$\tilde{Y}$", xlabel=r"$X_1$", title=r"$X_1$以外の変動を抜いた$\tilde{Y}$と$\tilde{X}$の散布図")

beta_ = np.linalg.inv(X_tilde.T @ X_tilde) @ X_tilde.T @ Y_tilde
xrange = np.linspace(X_tilde.min(), X_tilde.max(), 100)
axes[1].plot(xrange, xrange.reshape(-1, 1) @ beta_)

axes[2].scatter(X1, Y_tilde)
axes[2].set(ylabel=r"$\tilde{Y}$", xlabel=r"$X_1$", title=r"$Y$から$X_1$以外の変動を抜いた$\tilde{Y}$との散布図")
beta_ = np.linalg.inv(X1.T @ X1) @ X1.T @ Y_tilde
xrange = np.linspace(X1.min(), X1.max(), 100)
axes[2].plot(xrange, xrange.reshape(-1, 1) @ beta_)

fig.show()
/tmp/ipykernel_88307/1189295102.py:22: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
<Figure size 1200x300 with 3 Axes>

外生性・直交条件からの説明

浅野・中村では誤差が直交することから証明するスタイルで、有斐閣の本はより直接的に証明している

OLSのPartialling out解釈

y=β0+β1x1+β2x2+εy = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon
(13)

というモデルのβ1\beta_1は

  1. yyをx1x_1に回帰する：y=β0+β1x1+β2z2+εy = \beta_0 + \beta_1 x_1 + \beta_2 z_2 + \varepsilon

  2. yyをx~1\tilde{x}_1に回帰する（ここでx~1\tilde{x}_1はx1x_1をx2x_2に回帰した残差）

  3. y~\tilde{y}をx~1\tilde{x}_1に回帰する（ここでy~\tilde{y}はyyをx2x_2に回帰した残差）

の3つの方法のいずれかで推定することができる

数値例

import numpy as np

n = 1000
np.random.seed(0)
x0 = np.ones(shape=(n, ))
x2 = np.random.uniform(0, 1, size=n)
x1 = 3 * x2 + np.random.uniform(0, 1, size=n) + np.random.normal(0, 1, size=n)
X = np.array([ x0, x1, x2 ]).T
e = np.random.normal(0, 1, size=n)

beta = np.array([10, 5, 7])  # 真のbeta
y = X @ beta + e
np.corrcoef(X, rowvar=False)
/home/mitama/notes/.venv/lib/python3.10/site-packages/numpy/lib/function_base.py:2897: RuntimeWarning: invalid value encountered in divide
  c /= stddev[:, None]
/home/mitama/notes/.venv/lib/python3.10/site-packages/numpy/lib/function_base.py:2898: RuntimeWarning: invalid value encountered in divide
  c /= stddev[None, :]
array([[ nan, nan, nan], [ nan, 1. , 0.64601675], [ nan, 0.64601675, 1. ]])
class OLS:
    def fit(self, X, y):
        self.beta_ = np.linalg.inv(X.T @ X) @ X.T @ y
        return self

    def predict(self, X):
        return X @ self.beta_

1. yyをx1x_1に回帰する

y=β0+β1x1+β2z2+εy = \beta_0 + \beta_1 x_1 + \beta_2 z_2 + \varepsilon

# 1. 𝑦 を 𝑥1 に回帰する
OLS().fit(X, y).beta_
array([9.98226817, 5.03640436, 6.80111723])

2. yyをx~1\tilde{x}_1に回帰する

ここでx~1\tilde{x}_1はx1x_1をx2x_2に回帰した残差

# 2. 𝑦 を 𝑥̃ 1 に回帰する
X_ = X[:, [0, 2]] # x0, x2だけのX
x1_res = x1 - OLS().fit(X_, x1).predict(X_)

# 切片を付けた場合
X_ = np.array([ x0, x1_res ]).T
OLS().fit(X_, y).beta_
array([23.29878836, 5.03640436])
# 切片を付けない場合
OLS().fit(x1_res.reshape(-1, 1), y).beta_
array([5.03640436])

3. y~\tilde{y}をx~1\tilde{x}_1に回帰する

ここでy~\tilde{y}はyyをx2x_2に回帰した残差

# 3. 𝑦̃ を 𝑥̃ 1 に回帰する
X_ = X[:, [0, 2]] # x0, x2だけのX
y_res = y - OLS().fit(X_, y).predict(X_)

# 切片を付けない場合
OLS().fit(x1_res.reshape(-1, 1), y_res).beta_
array([5.03640436])
# 切片を付けた場合
X_ = np.array([ x0, x1_res ]).T
OLS().fit(X_, y_res).beta_
array([4.27435864e-15, 5.03640436e+00])

partialling out

OLS推定では説明変数と残差の共分散はゼロになる。

y=β0+β1x1+⋯+βjxj+⋯+βdxd+εy = \beta_0 + \beta_1 x_1 + \cdots + \beta_j x_j + \cdots + \beta_d x_d + \varepsilon
(14)

というモデルがあったとき、説明変数xjx_jを他のすべての説明変数に回帰すると、その残差x~j\tilde{x}_jは説明変数xjx_jの分散の情報を残す一方で他の説明変数とは無相関になる。

したがってyyをこの残差x~j\tilde{x}_jに回帰すると、その回帰係数βj\beta_jはyyに対するxjx_jの影響のみを示す。

→他の変数の影響を排除（partialling out）できる

FWL定理の応用

データの可視化

元のモデルにdd次元の説明変数があったとしても、x~j\tilde{x}_jとyyの関係へと次元を削減することができるため、グラフに表示しやすい。

statsmodelsで簡単に実行できる

statsmodelsのplot_regress_exog関数は残差の分析に関する4つの図をまとめて出力でき、そのうちのひとつ「Partial regression plot」がpartialling outした残差同士の散布図になっている。

import statsmodels.api as sm
results = sm.OLS(y, X).fit()
fig = sm.graphics.plot_regress_exog(results, 'x2')
Source
import matplotlib.lines as mlines
import matplotlib.collections as mcoll

def change_alpha(fig, alpha = 0.2):
    """plot_regress_exogの透明度を変える"""
    for ax in fig.axes:
        # scatter() の点（0番だけ）
        for coll in ax.collections:
            if isinstance(coll, mcoll.PathCollection):
                coll.set_alpha(alpha)

        # plot(..., marker='o', linestyle='None') の点（1〜3番）
        for line in ax.lines:
            if isinstance(line, mlines.Line2D):
                has_marker = line.get_marker() not in (None, "", "None")
                no_line = line.get_linestyle() in ("None", " ", "", None)
                if has_marker and no_line:
                    line.set_alpha(alpha)

import statsmodels.api as sm
results = sm.OLS(y, X).fit()
fig = sm.graphics.plot_regress_exog(results, 'x2')
change_alpha(fig, alpha = 0.2)
<Figure size 640x480 with 4 Axes>

計算の高速化

PyHDFEパッケージのような高次元データの分析において活用されているらしい

統計的因果推論

などで用いられる。DMLはモデルの学習アルゴリズムに機械学習を許容するので、過学習や正則化によるバイアスに対処するためのcross-fittingという推定方法を提案している。

歴史：Yule-Frisch-Waugh-Lovell Theorem

[2307.00369] The Yule-Frisch-Waugh-Lovell Theorem

FWLの前にYuleがいたらしい

通常、Frisch and Waugh (1933) と Lovell (1963) の名前をとって Frisch-Waugh-Lovell Theoremと呼ぶ。しかしこの論文によれば Yule (1907) も重要な貢献をしており、計量経済学では知名度がないものの統計学分野では注目されているため、FWL定理ではなくYFWL定理と呼ぶことを提案している。

参考文献

FWL

（考察）再帰的なpartialling outで任意のパラメータを任意の次元で推定できないか？

→ できなかった

うまくいってAdditive modelとのつながりが見えれば面白かったんだが

import numpy as np

# データを生成
n = 1000
np.random.seed(0)
x0 = np.ones(shape=(n, ))
x2 = np.random.uniform(0, 1, size=n)
x3 = np.random.uniform(0, 1, size=n)
x1 = x2 + x3 + np.random.normal(0, 1, size=n)
X = np.array([ x0, x1, x2, x3 ]).T
e = np.random.normal(0, 1, size=n)

beta = np.array([10, 5, 7, 3])  # 真のbeta
y = X @ beta + e


class OLS:
    def fit(self, X, y):
        self.beta_ = np.linalg.inv(X.T @ X) @ X.T @ y
        return self

    def predict(self, X):
        return X @ self.beta_

1. yyをx1x_1に回帰する

y=β0+β1x1+β2z2+εy = \beta_0 + \beta_1 x_1 + \beta_2 z_2 + \varepsilon

# 1. 𝑦 を 𝑥1 に回帰する
OLS().fit(X, y).beta_
array([9.95538609, 5.03200117, 6.87766976, 3.05740035])

2. yyをx~1\tilde{x}_1に回帰する

ここでx~1\tilde{x}_1はx1x_1を残りの説明変数に回帰した残差

x~1=x1−(β^0+β^2x2+β^3x3)\tilde{x}_1 = x_1 - (\hat{\beta}_0 + \hat{\beta}_2 x_2 + \hat{\beta}_3 x_3)
(15)
# 2. 𝑦 を 𝑥̃ 1 に回帰する
X_ = X[:, [0, 2, 3]]
x1_res = x1 - OLS().fit(X_, x1).predict(X_)


# 切片を付けない場合
OLS().fit(x1_res.reshape(-1, 1), y).beta_
array([5.03200117])
import matplotlib.pyplot as plt
plt.scatter(x1_res, x1)
<Figure size 640x480 with 1 Axes>
from scipy.stats import pearsonr
pearsonr(x1_res, x1).statistic.round(3)
0.933

説明変数1つずつでpartialling outできないか？

  1. x~1=x1−β^0x0\tilde{x}_1 = x_1 - \hat{\beta}_0 x_0

  2. x~1=x~1−β^2x2\tilde{x}_1 = \tilde{x}_1 - \hat{\beta}_2 x_2

  3. x~1=x~1−β^3x3\tilde{x}_1 = \tilde{x}_1 - \hat{\beta}_3 x_3

for all nuisance parameter indices j∈Jj \in J

  1. β^j=arg⁡min⁡βjE[(y−βjxj)2]\renewcommand{\argmin}{\mathop{\arg\min}} \hat{\beta}_j = \argmin_{\beta_j} E[(y - \beta_j x_j)^2]

  2. x~j=y−β^jxj\tilde{x}_j = y - \hat{\beta}_j x_j

# 2. 𝑦 を 𝑥̃ 1 に回帰する
x1_res = x1.copy()
nuisance_indices = [0, 2, 3]
for i in nuisance_indices:
    X_ = X[:, [i]]
    x1_res = x1_res - OLS().fit(X_, x1_res).predict(X_)
    print(f"corr(x1_res, x1): {pearsonr(x1_res, x1).statistic:.3f}")
    print(f"beta1           : {OLS().fit(x1_res.reshape(-1, 1), y).beta_[0]:.3f}")

OLS().fit(x1_res.reshape(-1, 1), y).beta_
corr(x1_res, x1): 1.000
beta1           : 5.750
corr(x1_res, x1): 0.998
beta1           : 3.435
corr(x1_res, x1): 0.998
beta1           : 2.925
array([2.92542108])
# 2. 𝑦 を 𝑥̃ 1 に回帰する
x1_res = x1.copy()
nuisance_indices = [0, 2, 3]
for i in nuisance_indices[::-1]:
    X_ = X[:, [i]]
    x1_res = x1_res - OLS().fit(X_, x1_res).predict(X_)
    print(f"corr(x1_res, x1): {pearsonr(x1_res, x1).statistic:.3f}")
    print(f"beta1           : {OLS().fit(x1_res.reshape(-1, 1), y).beta_[0]:.3f}")

OLS().fit(x1_res.reshape(-1, 1), y).beta_
corr(x1_res, x1): 0.887
beta1           : 7.193
corr(x1_res, x1): 0.876
beta1           : 2.998
corr(x1_res, x1): 0.876
beta1           : 4.613
array([4.61310126])
References
  1. Lovell, M. C. (1963). Seasonal Adjustment of Economic Time Series and Multiple Regression Analysis. Journal of the American Statistical Association, 58(304), 993–1010. 10.1080/01621459.1963.10480682
計量経済学・因果推論
対数変換
計量経済学・因果推論
パネルデータ分析