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.

OLS推定量の性質

重回帰モデル

NN個のサンプルがあり、ii番目のサンプルについての回帰式を次のように表記する

Yi=XiTβ+ui,i=1,,NY_i = X_i^T \beta + u_i, \hspace{1em} i=1,\dots, N

ここでXiRDX_i \in \mathbb{R}^Dは説明変数の行列XRN×DX\in \mathbb{R}^{N\times D}を1行取り出したもので、Yi,uiRY_i, u_i \in \mathbb{R}も1つのサンプルの被説明変数と誤差項である。

行列表記にすると

Y=Xβ+uY = X \beta + u

と表すことができる。

OLS推定量

目的関数は残差の二乗和であるため、

L(β)=uTu=(YXβ)T(YXβ)=YTYYTXβ(Xβ)TY+(Xβ)TXβ=YTY2βTXTY+βTXTXβ\begin{align} L(\beta) &= u^T u\\ &= (Y - X \beta)^T (Y - X \beta)\\ &= Y^T Y - Y^T X \beta - (X\beta)^T Y + (X \beta)^T X \beta\\ &= Y^T Y - 2 \beta^T X^T Y + \beta^T X^T X \beta\\ \end{align}

である。これを微分してゼロとおくと

L(β)β=2XTY+2XTXβ=0\frac{\partial L(\beta)} {\partial \beta} = - 2 X^T Y + 2 X^T X \beta = 0

となり、β\betaについて解くと

β=(XTX)1XTY\beta = (X^T X)^{-1} X^T Y
import numpy as np

x0 = np.array([1, 1, 1])
x1 = np.array([1, 2, 3])
x2 = np.array([2, 8, 9])
X = np.array([x0, x1, x2]).T
beta = np.array([3, 5, 7]) # 真のbeta
y = X @ beta

# OLS推定量
beta_hat = np.linalg.inv(X.T @ X) @ (X.T @ y)
beta_hat.round(1)
array([3., 5., 7.])

OLS推定量の別表記

一致性や不偏性の議論のための準備として、OLS推定量を変形する。

重回帰モデルY=Xβ+uY = X\beta + uをOLS推定量β^=(XTX)1XTY\hat{\beta} = (X^T X)^{-1} X^T Yに代入して変形すると

β^=(XTX)1XTY=(XTX)1XT(Xβ+u)=(XTX)1XTXβ+(XTX)1XTu=β+(XTX)1XTu=β+(1NXTX)11NXTu\begin{align} \hat{\beta} &= (X^T X)^{-1} X^T Y\\ &= (X^T X)^{-1} X^T (X\beta + u)\\ &= (X^T X)^{-1} X^T X\beta + (X^T X)^{-1} X^T u\\ &= \beta + (X^T X)^{-1} X^T u\\ &= \beta + \left(\frac{1}{N} X^T X\right)^{-1} \frac{1}{N} X^T u \end{align}

となる。

単回帰でいうと

β^=β+Cov(X,u)Var(X)\hat{\beta} = \beta + \frac{Cov(X, u)}{Var(X)}

である。

単回帰の導出

単回帰モデル Y=α+βX+uY=\alpha+\beta X+u は両辺をXXと共分散をとると

Cov(X,Y)=Cov(X,α+βX+u)=Cov(X,α)定数との共分散0+βCov(X,X)=Var(X)+Cov(X,u)=βVar(X)+Cov(X,u)\begin{aligned} \operatorname{Cov}(X, Y) &= \operatorname{Cov}(X, \alpha + \beta X + u)\\ &= \underbrace{ \operatorname{Cov}(X, \alpha) }_{定数との共分散 → 0} + \beta \underbrace{ \operatorname{Cov}(X, X) }_{=\operatorname{Var}(X) } +\operatorname{Cov}(X, u) \\ &= \beta \operatorname{Var}(X) + \operatorname{Cov}(X, u)\\ \end{aligned}

傾き係数のOLS推定量に代入すると

β^=Cov(X,Y)Var(X)=βVar(X)+Cov(X,u)Var(X)=β+Cov(X,u)Var(X)\begin{aligned} \hat{\beta} &= \frac{\operatorname{Cov}(X, Y)}{\operatorname{Var}(X)}\\ &= \frac{ \beta \operatorname{Var}(X) + \operatorname{Cov}(X, u) }{\operatorname{Var}(X)}\\ &= \beta + \frac{\operatorname{Cov}(X, u)}{\operatorname{Var}(X)} \end{aligned}

XXに内生性がある、すなわちCov(X,u)0Cov(X, u) \neq 0であるとβ^β\hat{\beta} \neq \betaとなる

OLS推定量のバリアンス

β^=β+(XTX)1XTu\hat{\beta} = \beta + (X^T X)^{-1} X^T uuN(0,σ2I)u\sim N(0, \sigma^2 I)の仮定より、β^N(β,σ2(XX)1)\hat{\beta} \sim N(\beta, \sigma^2 (X^\top X)^{-1})となる。

よって Var[β^]=σ2(XX)1\operatorname{Var}[\hat{\beta}] = \sigma^2 (X^\top X)^{-1} となる

証明

以下の定理を使う

定理

uuを確率変数ベクトルとし、μRn\mu\in\mathbb{R}^nbRpb\in\mathbb{R}^pCRp×nC\in\mathbb{R}^{p\times n}rankC=p\operatorname{rank} C = pとする。

uN(μ,Σ)u \sim N(\mu, \Sigma) のとき、Cu+bN(Cμ+b,CΣC)Cu + b \sim N(C\mu + b, C\Sigma C^\top)

C:=(XTX)1XTC := (X^T X)^{-1} X^Tとおけば

CT=[(XTX)1XT]T=X[(XTX)1]T((AB)T=BTAT)=X[(XTX)T]1((A1)T=(AT)1)=X(XTX)1((XTX)T=XTX)\begin{aligned} C^T = [(X^T X)^{-1} X^T]^T &= X [(X^T X)^{-1}]^T \quad (\because (AB)^T = B^T A^T) \\ &= X [(X^T X)^T]^{-1} \quad (\because (A^{-1})^T = (A^T)^{-1} ) \\ &= X (X^T X)^{-1} \quad (\because (X^T X)^T=X^T X ) \\ \end{aligned}

であるため、uN(0,σ2I)u\sim N(0, \sigma^2 I) の仮定が満たされるとき、

β+CuN(β,σ2CCT)=β+(XTX)1XTuN(β,σ2(XTX)1XTX(XTX)1)=β+(XTX)1XTuN(β,σ2(XTX)1)\begin{aligned} &\beta + C u \sim N(\beta, \sigma^2 CC^T)\\ &= \beta + (X^T X)^{-1} X^T u \sim N(\beta, \sigma^2 (X^T X)^{-1} X^T X (X^T X)^{-1} )\\ &= \beta + (X^T X)^{-1} X^T u \sim N(\beta, \sigma^2 (X^T X)^{-1}) \end{aligned}

不偏性

外生性

単回帰モデルY=α+βX+uY = \alpha + \beta X + uを例にとる。

また、外生性の条件は別の表現もできる

証明
E(Xu)=EX[E(XuX)]=EX[XE(uX)=0]=EX(X0)=0E(u)=EX[E(uX)=0]=0E(X u) = E_{X}[ E(X u | X) ] = E_{X}[ X \underbrace{ E(u|X) }_{ =0 } ] = E_{X}(X\cdot 0) = 0 \\ E(u) = E_{X}[ \underbrace{ E(u|X) }_{ =0 } ] = 0

さらに、共分散との関係も導出できる

証明
Cov(X,u)=E(Xu)=0E(X)E(u)=0=0Cov(X, u) = \underbrace{ E(X u) }_{ =0 } - E(X) \underbrace{ E(u) }_{ =0 } = 0

OLS推定量の不偏性

単回帰モデルY=α+βX+uY = \alpha + \beta X + uのOLS推定量β\betaの確率極限は

Undefined control sequence: \DeclareMathOperator at position 1: \̲D̲e̲c̲l̲a̲r̲e̲M̲a̲t̲h̲O̲p̲e̲r̲a̲t̲o̲r̲{\cov}{\text{Co…

\DeclareMathOperator{\cov}{\text{Cov}}
\DeclareMathOperator{\var}{\text{Var}}
\DeclareMathOperator{\plim}{\text{plim}}
\plim \beta = \beta + \frac{\cov(X, u)}{\var(X)}

となる。外生性が満たされるとき\cov(X, u) = 0であるため、plimβ=β\plim \beta = \betaとなり、OLS推定量は母回帰係数の不偏推定量となる。

証明
Undefined control sequence: \DeclareMathOperator at position 1: \̲D̲e̲c̲l̲a̲r̲e̲M̲a̲t̲h̲O̲p̲e̲r̲a̲t̲o̲r̲{\cov}{\text{Co…

\DeclareMathOperator{\cov}{\text{Cov}}
\DeclareMathOperator{\var}{\text{Var}}
\DeclareMathOperator{\plim}{\text{plim}}
\begin{align}
\plim \beta
&= \frac{\cov(X, Y)}{\var(X)}\\
&= \frac{\cov(X, \alpha + \beta X + u)}{\var(X)}\\
&= \frac{\cov(X, \alpha) + \cov(X, \beta X) + \cov(X, u) }{\var(X)}\\
&= \frac{\beta \var(X) + \cov(X, u) }{\var(X)}\\
&= \beta + \frac{\cov(X, u)}{\var(X)}
\end{align}

一致性

異常値がない(X,uX, uは4次までのモーメントを持つ)という仮定と大数の法則により以下が成立する

1NXTX=1Ni=1NXiXiTpE(XiXiT)1NXTu=1Ni=1NXiuipE(Xiui)=0\begin{align} \frac{1}{N} X^T X &= \frac{1}{N} \sum^N_{i=1} X_i X_i^T \overset{p}{\longrightarrow} E(X_i X_i^T) \\ \frac{1}{N} X^T u &= \frac{1}{N} \sum^N_{i=1} X_i u_i \overset{p}{\longrightarrow} E(X_i u_i) = 0 \end{align}

多重共線性がないという仮定により(E(XiXiT))1(E(X_i X_i^T))^{-1}が存在する

(1NXTX)11NXTup0\left( \frac{1}{N} X^T X \right)^{-1} \frac{1}{N} X^T u \overset{p}{\longrightarrow} 0

よって

β^pβ\hat{\beta} \overset{p}{\longrightarrow} \beta

漸近正規性

OLS推定量

β^=β+(1NXTX)11NXTu\hat{\beta} = \beta + \left(\frac{1}{N} X^T X\right)^{-1} \frac{1}{N} X^T u

を整理して以下の形にする

N(β^β)=(1NXTX)11NXTu\sqrt{N} (\hat{\beta} - \beta) = \left( \frac{1}{N} X^T X \right)^{-1} \frac{1}{\sqrt{N}} X^T u

1NXTu\frac{1}{\sqrt{N}} X^T u1Ni=1NXiui\frac{1}{\sqrt{N}} \sum^N_{i=1} X_i u_iと書くことができる。OLSの仮定より

E(Xiui)=0Var(Xiui)=E(u2XiXiT)\begin{align} E(X_i u_i) &= 0\\ Var(X_i u_i) &= E(u^2 X_i X^T_i) \end{align}

なので、中心極限定理により

1NXTu=1Ni=1NXiuidN(0,E(ui2XiXiT))\frac{1}{\sqrt{N}} X^T u = \frac{1}{\sqrt{N}} \sum^N_{i=1} X_i u_i \overset{d}{\longrightarrow} N\left( 0, E(u_i^2 X_i X_i^T) \right)

となる。

一致性のときに導出した

1NXTXpE(XiXiT)\frac{1}{N} X^T X \overset{p}{\longrightarrow} E(X_i X_i^T)

を使うと、スルツキーの定理を用いて

N(β^β)=(1NXTX)11NXTud(E(XiXiT))1×N(0,E(ui2XiXiT))=N(0,(E(XiXiT))1E(ui2XiXiT)(E(XiXiT))1)=N(0,V)\begin{align} \sqrt{N} (\hat{\beta} - \beta) &= \left( \frac{1}{N} X^T X \right)^{-1} \frac{1}{\sqrt{N}} X^T u \overset{d}{\longrightarrow} \left( E(X_i X_i^T) \right)^{-1} \times N\left( 0, E(u_i^2 X_i X_i^T) \right)\\ &= N\left( 0, \left( E(X_i X_i^T) \right)^{-1} E(u_i^2 X_i X_i^T) \left( E(X_i X_i^T) \right)^{-1} \right)\\ &= N(0, V)\\ \end{align}

となる。

スルツキーの定理

確率変数の行列YN,Y,XN,XRN×NY_N, Y, X_N, X \in \mathbb{R}^{N\times N}、正則行列CRN×NC \in \mathbb{R}^{N\times N}があるとする。

NN\to \inftyのとき

XNdXYNdCX_N \xrightarrow{d} X\\ Y_N \xrightarrow{d} C\\

とする。

このとき、以下の結果が成り立ち、これを スルツキーの定理 という

  1. XN+YNdX+CX_N + Y_N \xrightarrow{d} X + C

  2. YNXNdCXY_N X_N \xrightarrow{d} C X

  3. YN1XNdC1XY_N^{-1}X_N \xrightarrow{d} C^{-1} X

VVは以下のように一致推定できる

V^=[1Ni=1NXiXiT]11Ni=1Nu^i2XiXiT[1Ni=1NXiXiT]1=(1NXTX)11NXTU^X(1NXTX)1\begin{aligned} \hat{V} & =\left[\frac{1}{N} \sum_{i=1}^N X_i X_i^T\right]^{-1} \frac{1}{N} \sum_{i=1}^N \hat{u}_i^2 X_i X_i^T\left[\frac{1}{N} \sum_{i=1}^N X_i X_i^T\right]^{-1} \\ & =\left(\frac{1}{N} X^T X\right)^{-1} \frac{1}{N} X^T \hat{U} X\left(\frac{1}{N} X^T X\right)^{-1} \end{aligned}

ただし、U^\hat{U}は対角要素にu^12,,u^N2\hat{u}_1^2, \dots, \hat{u}_N^2を並べた対角行列である。

バイアスとバリアンス

最小二乗推定量はすべての線形不偏推定量の中で最もバリアンスが小さい(最良である)ことを示すガウス・マルコフの定理というものがある。

不偏性

任意のパラメータの線形結合θ=αβ\theta=\boldsymbol{\alpha}^\top {\boldsymbol{\beta}}を考える。例えばf(x0)=x0βf(x_0)=x_0^\top \betaがこの形である。

この最小二乗推定値は

θ^=αβ^=α(XX)1Xy\hat{\theta} = \boldsymbol{\alpha}^\top \hat{\boldsymbol{\beta}} = \boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}

で、期待値をとると

E[θ^]=E[αβ^]=E[α(XX)1Xy]=E[α(XX)1X(Xβ+ε)]=α(XX)1XXβ+α(XX)1XE[ε]=α(XX)1XXβ=αβ\begin{align} E[\hat{\theta}] &= E[\boldsymbol{\alpha}^\top \hat{\boldsymbol{\beta}}]\\ &= E[\boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}]\\ % &= \boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top E[\boldsymbol{y}]\\ &= E[\boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top (X\beta + \varepsilon) }] \\ &= \boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top X \beta + \alpha^\top (X^\top X)^{-1} X^\top} E[ \varepsilon ] \\ &= \boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top X \beta}\\ &= \boldsymbol{\alpha}^\top \boldsymbol{\beta} \end{align}

となり(※)、θ^\hat{\theta}が不偏推定量である(E[θ^]=θE[\hat{\theta}] = \theta)ことがわかる。

バリアンス

αβ\boldsymbol{\alpha}^\top \boldsymbol{\beta}に対して不偏のまた別の線形推定量cy\boldsymbol{c}^\top \boldsymbol{y}があるとする。

両者の差を

αβ^cy=[α(XX)1Xc]y=:dy\begin{align} \boldsymbol{\alpha}^\top \boldsymbol{\hat{\beta}} - \boldsymbol{c}^\top \boldsymbol{y} &= [ \boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X} - \boldsymbol{c}^\top ] \boldsymbol{y}\\ &=: \boldsymbol{d}^\top \boldsymbol{y} \end{align}

とおく。このとき、不偏性E[cy]=αβE[\boldsymbol{c}^\top \boldsymbol{y}] = \boldsymbol{\alpha}^\top \boldsymbol{\beta}から両者の差の期待値はゼロになるべきであり、

E[dy]=dXβ=0E[\boldsymbol{d}^\top \boldsymbol{y}] = \boldsymbol{d}^\top \boldsymbol{X} \boldsymbol{\beta} = 0

が任意のβ\boldsymbol{\beta}について成り立たなければならないため、

dX=0\boldsymbol{d}^\top \boldsymbol{X} = \boldsymbol{0}

が成り立つ。

次に、2つの確率変数X,YX, Yに対し

V[X+Y]=V[X]+2Cov[X,Y]+V[Y]V[X+Y] = V[X] + 2 \text{Cov}[X, Y] + V[Y]

が成り立つから、cy\boldsymbol{c}^\top \boldsymbol{y}の分散は

V[cy]=V[αβ^dy]=V[αβ^]2Cov[α(XX)1Xy,dy]+V[dy]\begin{align} V[\boldsymbol{c^\top y}] &= V[\boldsymbol{ \alpha^\top \hat{\beta} - d^\top y }]\\ &= V[\boldsymbol{ \alpha^\top \hat{\beta} }] - 2 \text{Cov} [\boldsymbol{ \alpha^\top (X^\top X)^{-1} X^\top y }, \boldsymbol{ d^\top y} ] + V[\boldsymbol{ d^\top y}] \end{align}

と表すことができる。

Cov(A,B)=E[(AE[A])(BE[B])]\text{Cov}(A, B) = E[(A - E[A])(B - E[B])^\top]
Cov(ay,by)=E[(ayE[ay])(byE[by])]\text{Cov}(a^\top y, b^\top y) = E[(a^\top y - E[a^\top y])(b^\top y - E[b^\top y])^\top]
Cov(ay,by)=E[(ayE[ay])(byE[by])]=E[(ayαβ)(byE[by])]\text{Cov}(a^\top y, b^\top y) = E[(a^\top y - E[a^\top y])(b^\top y - E[b^\top y])^\top]\\ = E[(a^\top y - \alpha^\top \beta)(b^\top y - E[b^\top y])^\top]

次に、2つの確率変数

ay=aiyi,by=biyi\boldsymbol{a^\top y} = \sum a_i y_i, \hspace{2em} \boldsymbol{b^\top y} = \sum b_i y_i

の共分散は、誤差項ε\boldsymbol{\varepsilon}が無相関・等分散の仮定V[ε]=σ2IV[\boldsymbol{\varepsilon}] = \sigma^2 \boldsymbol{I}を満たすとき、

Cov[ay,by]=Cov[aε,bε]=aibiσ2=(ab)σ2\text{Cov}[\boldsymbol{a^\top y, b^\top y}] = \text{Cov}[\boldsymbol{a^\top \varepsilon, b^\top \varepsilon}] = \sum a_i b_i \sigma^2 = (\boldsymbol{a^\top b}) \sigma^2

となることから

Cov[α(XX)1Xy,dy]=α(XX)1Xdσ2\text{Cov}[\boldsymbol{ \alpha^\top (X^\top X)^{-1} X^\top y, d^\top y }] = \boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top d } \cdot \sigma^2

となり、dX=0\boldsymbol{d}^\top \boldsymbol{X} = \boldsymbol{0}よりこれは0となる。

よって

V[cy]=V[αβ^]+V[dy]V[\boldsymbol{ c^\top y }] = V[\boldsymbol{ \alpha^\top \hat{\beta} }] + V[\boldsymbol{ d^\top y }]

が成り立ち、分散は非負なので

V[cy]V[αβ^]V[\boldsymbol{ c^\top y }] \geq V[\boldsymbol{ \alpha^\top \hat{\beta} }]

を意味する。

よってαβ^\boldsymbol{ \alpha^\top \hat{\beta} }は最良線形不偏推定量BLUEである。

OLS推定の幾何学的意味

OLS推定量

β^=(XX)1Xy\hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{y}

y^=Xβ^\hat{\boldsymbol{y}} = X \hat{\boldsymbol{\beta}}に代入すると

y^=X(XX)1XPy=Py\hat{\boldsymbol{y}} = \underbrace{ X (X^\top X)^{-1} X^\top }_{P} \boldsymbol{y} = P \boldsymbol{y}

つまり、ベクトルy\boldsymbol{y}を行列P=X(XX)1XP = X (X^\top X)^{-1} X^\topで射影したものとみなすことができる。

この行列PPは対称行列で、P2=PP^2=Pとなる。この2つの性質を満たす行列を射影行列という。

P2=PP=(X(XX)1X)(X(XX)1X)=X(XX)1(XX)(XX)1X=X(XX)1X=P\begin{aligned} P^2 & = PP \\ & =(X(X^{\top} X)^{-1} X^{\top})(X(X^{\top} X)^{-1} X^{\top}) \\ & =X(X^{\top} X)^{-1}(X^{\top} X)(X^{\top} X)^{-1} X^{\top} \\ & =X(X^{\top} X)^{-1} X^{\top} =P \end{aligned}

射影行列は、XXの列空間X\Im Xにベクトルを正射影するという性質がある。y\boldsymbol{y}X\Im Xへの射影がy^\hat{\boldsymbol{y}}で、垂線の足が誤差u\boldsymbol{u}となる。

よって、最小二乗法はy\boldsymbol{y}からX\Im Xへの射影を求める操作であると捉えることができる。

OLSとBLUE / BUE

OLSはBLUEかBUEか?

行列AAによってb=Ayb = Ayのような線形結合で表現される推定量bbは線形推定量という。OLSもA=(XX)1XA=(X^\top X)^{-1} X^\topとおけば同じ形になっていることがわかる。

OLSはガウス・マルコフの定理でBLUE(線形不偏推定量のなかで最良)だと示された。

Hansen (2022) は線形制約は不要で、線形と非線形の両方のすべての不偏推定量の中で最良(BUE)だと主張した。 一方で Pötscher & Preinerstorfer (2022) はHansen (2022)に対する批判を展開した。 Portnoy (2022) は「一般線形モデルの場合、十分に広い分布族内のすべての分布に対して、不偏な推定量は線形でなければならない」と述べている。

Hansen(2022) の「OLSが最良の不偏推定量(BUE)」が正しいとしても、Portnoy (2022)によれば「一般線形モデルの不偏推定量は線形推定量」なので結局BLUE=BUEになり、いずれの主張も間違っていないことになる。

参考:

参考文献

References
  1. Hansen, B. E. (2022). A Modern Gauss–Markov Theorem. Econometrica, 90(3), 1283–1294. 10.3982/ecta19255
  2. Pötscher, B. M., & Preinerstorfer, D. (2022). A Modern Gauss-Markov Theorem? Really? 10.48550/ARXIV.2203.01425
  3. Portnoy, S. (2022). Linearity of Unbiased Linear Model Estimators. The American Statistician, 76(4), 372–375. 10.1080/00031305.2022.2076743