2021-09-21

計量経済学メモ:行列を用いた最小二乗推定量の導出

・本稿の内容
最少二乗推定量の行列を用いた導出過程をメモします。重回帰モデルの説明変数の数は定数項を含めて $k$ 個、データのサイズは $n$ 個( $n\ge k$ )として話を進めます。

Ⅰ:準備
Ⅱ:最少二乗推定量の導出
Ⅲ:具体例
Ⅳ:参考文献

Ⅰ:準備

$k$ 変数関数 $z=f(\beta_{1},\beta_{1},\cdots ,\beta_{k})$ を $\boldsymbol{\beta}=(\beta_{1},\beta_{1},\cdots ,\beta_{k})'$ の関数と考える。 $f(\boldsymbol{\beta})$ の $\boldsymbol{\beta}$ による偏微分は以下のように定義される。

$\dfrac{\partial f(\boldsymbol{\beta}) }{\partial \boldsymbol{\beta}} = \begin{pmatrix} \dfrac{\partial f(\boldsymbol{\beta}) }{\partial \beta_{1}}\\ \dfrac{\partial f(\boldsymbol{\beta}) }{\partial \beta_{2}}\\ \vdots \\ \dfrac{\partial f(\boldsymbol{\beta}) }{\partial \beta_{k}} \end{pmatrix}$

$\boldsymbol{a}=(a_1,a_2,\cdots,a_k)'$ 、 $\boldsymbol{\beta}=(\beta_1,\beta_2,\cdots,\beta_k)'$ としたときに、以下の公式が成り立つ。

公式1

$\begin{eqnarray} \dfrac{\partial \boldsymbol{a}'\boldsymbol{\beta} }{\partial \boldsymbol{\beta}} &=& \boldsymbol{a}\\ &=& (\boldsymbol{a}')' \end{eqnarray}$

公式2

$\begin{eqnarray} \dfrac{\partial \boldsymbol{\beta}' \boldsymbol{a} }{\partial \boldsymbol{\beta}} &=& \boldsymbol{a} \end{eqnarray}$

公式1と公式2が成り立つことを確認する。単純化のために、 $k=2$ (つまり、 $\boldsymbol{a}=(a_1,a_2)'$ 、 $\boldsymbol{\beta}=(\beta_1,\beta_2)'$ )とする。

公式1の証明

$\begin{eqnarray} \boldsymbol{a}'\boldsymbol{\beta} &=& \begin{pmatrix} a_1&a_2 \end{pmatrix} \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}\\ &=& a_1\beta_1+a_2\beta_2 \end{eqnarray}$

$\begin{eqnarray} \dfrac{\partial \boldsymbol{a}'\boldsymbol{\beta} }{\partial \boldsymbol{\beta}} &=& \begin{pmatrix} \dfrac{\partial (a_1\beta_1+a_2\beta_2)}{\partial \beta_{1}}\\ \dfrac{\partial (a_1\beta_1+a_2\beta_2) }{\partial \beta_{2}}\\ \end{pmatrix}\\ &=& \begin{pmatrix} a_1\\ a_2\\ \end{pmatrix}\\ &=& \boldsymbol{a}\\ &=& (\boldsymbol{a}')' \end{eqnarray}$

よって、公式1が成り立つ。

公式2の証明

$\begin{eqnarray} \boldsymbol{\beta}'\boldsymbol{a} &=& \begin{pmatrix} \beta_1&\beta_2 \end{pmatrix} \begin{pmatrix} a_1\\ a_2 \end{pmatrix}\\ &=& a_1\beta_1+a_2\beta_2 \end{eqnarray}$

$\begin{eqnarray} \dfrac{\partial \boldsymbol{\beta}'\boldsymbol{a} }{\partial \boldsymbol{\beta}} &=& \begin{pmatrix} \dfrac{\partial (a_1\beta_1+a_2\beta_2)}{\partial \beta_{1}}\\ \dfrac{\partial (a_1\beta_1+a_2\beta_2)}{\partial \beta_{2}}\\ \end{pmatrix}\\ &=& \begin{pmatrix} a_1\\ a_2\\ \end{pmatrix}\\ &=& \boldsymbol{a} \end{eqnarray}$

よって、公式2が成り立つ。

$\boldsymbol{\beta}=(\beta_1,\beta_2,\cdots,\beta_k)'$ 、 $C=\begin{pmatrix}c_{11}& c_{12} & \cdots & c_{1k}\\c_{21}& c_{22} &\cdots & c_{2k}\\\vdots &\vdots&\ddots& \vdots \\c_{k1}& c_{k2} & \cdots & c_{kk}\end{pmatrix}$ としたときに、以下の公式が成り立つ。

公式3

$\begin{eqnarray} \dfrac{\partial \boldsymbol{\beta}'C\boldsymbol{\beta} }{\partial \boldsymbol{\beta}} &=& (C+C')\boldsymbol{\beta}\\ \end{eqnarray}$

公式4
$C$ が対称行列であれば

$\begin{eqnarray} \dfrac{\partial \boldsymbol{\beta}'C\boldsymbol{\beta} }{\partial \boldsymbol{\beta}} &=& 2C\boldsymbol{\beta}\\ \end{eqnarray}$

公式3と公式4が成り立つことを確認する。単純化のために、 $k=2$ (つまり、 $\boldsymbol{\beta}=(\beta_1,\beta_2)'$ 、 $C=\begin{pmatrix}c_{11}& c_{12}\\c_{21}& c_{22}\end{pmatrix}$ )とする。

公式3の証明

$\begin{eqnarray} \boldsymbol{\beta}'C\boldsymbol{\beta} &=& \begin{pmatrix} \beta_1&\beta_2 \end{pmatrix} \begin{pmatrix} c_{11}& c_{12}\\ c_{21}& c_{22} \end{pmatrix} \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}\\ &=& c_{11}\beta_{1}^{2}+c_{12}\beta_{1}\beta_{2}+c_{21}\beta_{1}\beta_{2}+c_{22}\beta_{2}^{2} \end{eqnarray}$

$\begin{eqnarray} \dfrac{\partial \boldsymbol{\beta}'C\boldsymbol{\beta} }{\partial \boldsymbol{\beta}} &=& \begin{pmatrix} \dfrac{\partial (c_{11}\beta_{1}^{2}+c_{12}\beta_{1}\beta_{2}+c_{21}\beta_{1}\beta_{2}+c_{22}\beta_{2}^{2})}{\partial \beta_{1}}\\ \dfrac{\partial (c_{11}\beta_{1}^{2}+c_{12}\beta_{1}\beta_{2}+c_{21}\beta_{1}\beta_{2}+c_{22}\beta_{2}^{2})}{\partial \beta_{2}}\\ \end{pmatrix}\\ &=& \begin{pmatrix} 2c_{11}\beta_{1}+(c_{12}+c_{21})\beta_{2}\\ (c_{12}+c_{21})\beta_{1}+2c_{22}\beta_{2}\\ \end{pmatrix}\\ &=& \begin{pmatrix} 2c_{11}&(c_{12}+c_{21})\\ (c_{12}+c_{21})&2c_{22}\\ \end{pmatrix} \begin{pmatrix} \beta_{1}\\ \beta_{2} \end{pmatrix}\\ &=& \left\{ \begin{pmatrix} c_{11}& c_{12}\\ c_{21}& c_{22}\\ \end{pmatrix} + \begin{pmatrix} c_{11}& c_{21}\\ c_{12}& c_{22}\\ \end{pmatrix} \right\} \begin{pmatrix} \beta_{1}\\ \beta_{2} \end{pmatrix}\\ &=& (C+C')\boldsymbol{\beta} \end{eqnarray}$

よって、公式3が成り立つ。

公式4の証明

$C=C'$ であるから、公式3より、

$\begin{eqnarray} \dfrac{\partial \boldsymbol{\beta}'C\boldsymbol{\beta} }{\partial \boldsymbol{\beta}} &=& 2C\boldsymbol{\beta}\\ \end{eqnarray}$

となる。よって、公式4が成り立つ。

Ⅱ:最少二乗推定量の導出

$k$ 個の説明変数を持ち、データサイズが $n$ 個である重回帰モデルを考える。

$y_1= \beta_1x_{11}+\cdot\cdot\cdot+\beta_kx_{1k}+e_1\\ y_2= \beta_1x_{21}+\cdot\cdot\cdot+\beta_kx_{2k}+e_2\\ \cdot\\ \cdot\\ \cdot\\ y_n= \beta_1x_{n1}+\cdot\cdot\cdot+\beta_kx_{nk}+e_n$

行列表記を用いると、

$\left(\begin{array}{c} y_1\\ y_2\\ \vdots \\ y_n \end{array}\right) = \begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1k}\\ x_{21}& x_{22} & \cdots & x_{2k}\\ \vdots &\vdots&\ddots& \vdots \\ x_{n1} & x_{n2} & \cdots & x_{nk} \end{pmatrix} \begin{pmatrix} \beta_1\\ \beta_2\\ \vdots \\ \beta_k \end{pmatrix} + \begin{pmatrix} e_1\\ e_2\\ \vdots \\ e_n \end{pmatrix}$

と書ける。本稿では定数項を持つ重回帰モデルを考えるため、以下では $x_{i1}(i=1,2\cdots n)$ は $1$ とする。定数項を持つ場合の重回帰モデルは以下のようになる。

$y_1= \beta_1+\beta_2x_{12}+\cdot\cdot\cdot+\beta_kx_{1k}+e_1\\ y_2= \beta_1+\beta_2x_{22}+\cdot\cdot\cdot+\beta_kx_{2k}+e_2\\ \cdot\\ \cdot\\ \cdot\\ y_n= \beta_1+\beta_2x_{n2}+\cdot\cdot\cdot+\beta_kx_{nk}+e_n$

行列表記を用いると、

$\underbrace{ \left(\begin{array}{c} y_1\\ y_2\\ \vdots \\ y_n \end{array}\right) }_{n×1} = \underbrace{ \begin{pmatrix} 1& x_{12} & \cdots & x_{1k}\\ 1& x_{22} & \cdots & x_{2k}\\ \vdots &\vdots&\ddots& \vdots \\ 1& x_{n2} & \cdots & x_{nk} \end{pmatrix} }_{n×k} \underbrace{ \begin{pmatrix} \beta_1\\ \beta_2\\ \vdots \\ \beta_k \end{pmatrix} }_{k×1} + \underbrace{ \begin{pmatrix} e_1\\ e_2\\ \vdots \\ e_n \end{pmatrix} }_{n×1}$

と書ける。

$\boldsymbol{y} =\left(\begin{array}{c} y_1\\ y_2\\ \vdots \\ y_n \end{array}\right), X=\begin{pmatrix} 1 & x_{12} & \cdots & x_{1k}\\ 1 & x_{22} & \cdots & x_{2k}\\ \vdots &\vdots&\ddots& \vdots \\ 1 & x_{n2} & \cdots & x_{nk} \end{pmatrix}, \boldsymbol{\beta} =\begin{pmatrix} \beta_1\\ \beta_2\\ \vdots \\ \beta_k \end{pmatrix}, \boldsymbol{e}=\begin{pmatrix} e_1\\ e_2\\ \vdots \\ e_n \end{pmatrix}$

とおいて、表記を簡略化する。

$\boldsymbol{y}=X\boldsymbol{\beta}+\boldsymbol{e}・・・①$

①式を整理すると、

$\boldsymbol{e}=\boldsymbol{y}-X\boldsymbol{\beta}・・・②$

と書ける。残差二乗和 $\boldsymbol{e}'\boldsymbol{e}$ が最小になるような $\boldsymbol{\beta}$ を求めていく。*1 $\boldsymbol{e}'\boldsymbol{e}$ は

$\begin{eqnarray} \boldsymbol{e}'\boldsymbol{e} &=& (\boldsymbol{y}-X\boldsymbol{\beta})'(\boldsymbol{y}-X\boldsymbol{\beta})\\ &=&(\boldsymbol{y}'-\boldsymbol{\beta}'X')(\boldsymbol{y}-X\boldsymbol{\beta})\\ &=&(\boldsymbol{y}'\boldsymbol{y}-\boldsymbol{y}'X\boldsymbol{\beta}-\beta'X'\boldsymbol{y}+\boldsymbol{\beta}'X'X\boldsymbol{\beta})・・・③ \end{eqnarray}$

と書ける。③式を $\boldsymbol{\beta}$ で微分する。

$\dfrac{\partial \boldsymbol{e}'\boldsymbol{e}}{\partial \boldsymbol{\beta}}=\dfrac{\partial (\boldsymbol{y}'\boldsymbol{y}-\boldsymbol{y}'X\boldsymbol{\beta}-\beta'X'\boldsymbol{y}+\boldsymbol{\beta}'X'X\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}・・・④$

④式の右辺を項別に見ていく。
1. $\boldsymbol{y}'\boldsymbol{y}$ について
$\boldsymbol{y}'\boldsymbol{y}$ は $\boldsymbol{\beta}$ と関係がないため、

$\dfrac{\partial \boldsymbol{y}'\boldsymbol{y}}{\partial \boldsymbol{\beta}}=\boldsymbol{0}$

となる。

2. $\\-\boldsymbol{y}'X\boldsymbol{\beta}$ について
Ⅰ:準備で確認した公式1を用いて、

$\\-\dfrac{\partial \boldsymbol{y}'X\boldsymbol{\beta}}{\partial \boldsymbol{\beta}}= -(\boldsymbol{y}'X)'=-X'\boldsymbol{y}$

となる。

3. $\\-\boldsymbol{\beta}'X'\boldsymbol{y}$ について
Ⅰ:準備で確認した公式2を用いて、

$\\-\dfrac{\partial \boldsymbol{\beta}'X'\boldsymbol{y}}{\partial \boldsymbol{\beta}}=-X'\boldsymbol{y}$

となる。

4. $\boldsymbol{\beta}'X'X\boldsymbol{\beta}$ について
$X'X$ は対称行列であることに注意し、Ⅰ:準備で確認した公式4を用いて、

$\dfrac{\partial \boldsymbol{\beta}'X'X\boldsymbol{\beta}}{\partial \boldsymbol{\beta}}=2X'X\boldsymbol{\beta}$

となる。

1～4の結果をまとめると

$\dfrac{\partial \boldsymbol{e}'\boldsymbol{e}}{\partial \boldsymbol{\beta}}=-X'\boldsymbol{y}-X'\boldsymbol{y}+2X'X\boldsymbol{\beta}・・・⑤$

となる。

$\boldsymbol{\beta}$ が残差二乗和を最小化するとき、⑤式は $\boldsymbol{0}$ である。そのときの $\boldsymbol{\beta}$ を $\boldsymbol{\hat{\beta}}$ と書くことにすると⑤式は

$\\-X'\boldsymbol{y}-X'\boldsymbol{y}+2X'X\boldsymbol{\hat{\beta}} = \boldsymbol{0}・・・⑥$

と書ける。

⑥式を整理すると

$\\-2X'\boldsymbol{y}+2X'X\boldsymbol{\hat{\beta}} = \boldsymbol{0}\\ \\-X'\boldsymbol{y}+X'X\boldsymbol{\hat{\beta}} = \boldsymbol{0}\\ X'X\boldsymbol{\hat{\beta}}=X'\boldsymbol{y}・・・⑦\\ \boldsymbol{\hat{\beta}}=(X'X)^{-1}X'\boldsymbol{y}・・・⑧$

と書ける。⑦式は正規方程式という。⑧式が最少二乗推定量である。*2

Ⅲ:具体例

データサイズ $n$ を5、定数項を含めた説明変数の数 $k$ を3とし、具体的な数値例で⑦式を使用してみる。数値例として用いるデータセットは以下の表の通りである。*3

表1 データセット

重回帰モデルは以下のようになる。

$y_1= \beta_1+\beta_2 x_{12}+\beta_3 x_{13}+e_1\\ y_2= \beta_1+\beta_2 x_{22}+\beta_3 x_{23}+e_2\\ y_3= \beta_1+\beta_2 x_{32}+\beta_3 x_{33}+e_3\\ y_4= \beta_1+\beta_2 x_{42}+\beta_3 x_{43}+e_4\\ y_5= \beta_1+\beta_2 x_{52}+\beta_3 x_{53}+e_5\\$

行列表記を用いると、

$\underbrace{ \left(\begin{array}{c} y_1\\ y_2\\ y_3\\ y_4\\ y_5 \end{array}\right) }_{5×1} = \underbrace{ \begin{pmatrix} 1& x_{12} & x_{13}\\ 1& x_{22} & x_{23}\\ 1& x_{32} & x_{33}\\ 1& x_{42} & x_{43}\\ 1& x_{52} & x_{53}\\ \end{pmatrix} }_{5×3} \underbrace{ \begin{pmatrix} \beta_1\\ \beta_2\\ \beta_3 \end{pmatrix} }_{3×1} + \underbrace{ \begin{pmatrix} e_1\\ e_2\\ e_3\\ e_4\\ e_5 \end{pmatrix} }_{5×1}$

と書ける。⑧式を用いて、最少二乗推定量を求める。

$\begin{eqnarray} \begin{pmatrix} \hat{\beta_1}\\ \hat{\beta_2}\\ \hat{\beta_3} \end{pmatrix} &=& \left\{ \begin{pmatrix} 1&1&1&1&1\\ x_{12}&x_{22}&x_{32}&x_{42}&x_{52}\\ x_{13}&x_{23}&x_{33}&x_{43}&x_{53} \end{pmatrix} \begin{pmatrix} 1& x_{12} & x_{13}\\ 1& x_{22} & x_{23}\\ 1& x_{32} & x_{33}\\ 1& x_{42} & x_{43}\\ 1& x_{52} & x_{53}\\ \end{pmatrix} \right\}^{-1} \begin{pmatrix} 1&1&1&1&1\\ x_{12}&x_{22}&x_{32}&x_{42}&x_{52}\\ x_{13}&x_{23}&x_{33}&x_{43}&x_{53} \end{pmatrix} \left(\begin{array}{c} y_1\\ y_2\\ y_3\\ y_4\\ y_5 \end{array}\right)\\ &=& \left\{ \begin{pmatrix} 1&1&1&1&1\\ 5&7&10&11&12\\ 4&8&6&2&5 \end{pmatrix} \begin{pmatrix} 1& 5 & 4\\ 1& 7 & 8\\ 1& 10 & 6\\ 1& 11 & 2\\ 1& 12 & 5\\ \end{pmatrix} \right\}^{-1} \begin{pmatrix} 1&1&1&1&1\\ 5&7&10&11&12\\ 4&8&6&2&5 \end{pmatrix} \left(\begin{array}{c} 24\\ 20\\ 42\\ 60\\ 54 \end{array}\right) \end{eqnarray}$

上式の計算を続ければ最少二乗推定量は求まるが、手計算だと時間がかかる。本稿では $R$ 言語を用いて、残りの計算を実施する。以下のような $R$ コードを記述する。

#被説明変数Y_iを5×1行列にする
y<-matrix(data=c(24,20,42,60,54),nrow=5)
y

#説明変数X_i1,X_i2,X_i3を5×3行列にする
#※定数項ありのモデルだから、X_i1は要素がすべて1の5×1行列にする
x1<-matrix(data=c(1,1,1,1,1),nrow=5)
x2<-matrix(data=c(5,7,10,11,12),nrow=5)
x3<-matrix(data=c(4,8,6,2,5),nrow=5)
X<-cbind(X1,X2,X3)
X

#Xの転置行列
X_t<-t(X)
X_t

#X_tとXの積の逆行列
X_inv<-solve(X_t%*%X)
X_inv

#推定された係数
betahat<-X_inv%*%X_t%*%Y
betahat