Econometrics101 - 復習

Table of Contents

• 1. Linear Regressionの解釈
  • 解釈 1: Linear Conditional Expectations
  • 解釈 2: The CEFの”Best”線形近似
  • 解釈 3: Causal Model
    • ATT and selection bias
• 2. Ordinary Least Squares Estimator
  • The residualとThe error termの違い
  • Demeaned Regressors
  • Projection Matrix
  • Residual Regression: FWL theorem
  • Leverage Values
    • the leverageの３つの性質
    • Leave-One-Out Regression
  • Goodness of fit measure
    • RegressorとR-squaredの関係
• 3. Property of OLS Estimator
  • Assumptions
  • Unbiasedness
  • Conditional Variance of OLS estimator
  • Gauss-Markov Theorem
• 4. Asymptotic Property of OLS Estimator

1. Linear Regressionの解釈

$(Y_i, X_i, U_i)$をrandom vectorで、$Y_i, U_i \in \mathbf R$, $X_i \in \mathbf R^{k+1}$とします. $X$の最初の要素は定数1とします、つまり、

このとき、以下の線型回帰問題を考えます.

この$\beta$は以下の推定量を解くことで得られます.

(1-1) のFOCをとると

この$\beta$の解釈をまず解説します.

解釈 1: Linear Conditional Expectations

The CEFが線形であると仮定します, i.e., $E[Y_i|X_i] = X_i’\beta^*$ and $U_i = Y_i - E[Y_i|X_i]$. また、$E[|Y|] < \infty$とします

Theorem 1-1: THE CEF Decomposition Property

このとき、$E[U_i|X_i] = 0$である. また、$U_i$はいかなる$X_i$の関数とは無相関である.

Proof

またこのことより、

定理の後者の命題については、

Theorem 1-3: THE CEF Prediction Property

$m(X_i)$をany function of $X_i$とします. The CEFは以下の最小化問題の解と一致する:

Proof

Theorem 1-3: THE Linear CEF theorem

the CEFがLinearのとき、\(X_i'\beta\)はthe CEFそのものである.

Proof

$E[Y_i|X_i] = X_i’\beta^*$と仮定します. このとき、定理1-1より

よって \(\beta^* = E[X_iX_i']^{-1}E[X_iY_i] = \beta\)

このとき、the CEFが推定できたことを以て、\(\beta_k\)を$X_{ik}$が\(Y_i\)に与える因果の意味での効果(the centeris paribus)と解釈したくなりますが、これは間違いです. あくまで\(\beta\)は \((Y_i, X_i)\) の同時分布の情報を集約した数値であり、必ずしもそれは因果の関係性を意味するわけではありません. 例えば、\(Y_i\)を個人の年間収入、\(X_i\)を年間所得税納付金額としたとき、$E[Y_i|X_i]$が上昇すればするほど、\(X_i\)は上がると考えられますが、その逆はthe centeris paribusの意味において考えることは一般的には難しいです.

解釈 2: The CEFの”Best”線形近似

Theorem 2-1: The BEST MMSE linear approximation

The function \(X_i'\beta\)はthe CEFのthe MMSE linear approximationである, i.e.,

Proof

従って、

なお$\beta$がThe BEST MMSE linear approximationであることは以下の式からもわかります

Approxiamtionはsupportに依存する

$\text{Supp}(X)$をrandom vector $X$のサポートと定義します. このとき

よって、

• the approximationは$X$の分布に依存する
• finite sampleにおいて、$X$のdataが観測されないエリアではthe approximationのweightはゼロになる
• $\text{Supp}(X)$が非常に狭い範囲の場合、finite sampleにおいて the approximationは非常に不安定になる

解釈 3: Causal Model

$(Y_i, D_i)$のrandom vectorからなるデータを考えます. $Y_i$はthe outcome, $D_i$はtreatment statusを示すbinary変数とします. このとき、$D_i$と$Y_i$のthe causal relationshipはthe potential outcomeを用いて記述することができます;

このときのthe treatment effectは

と記述することができます. この設定の下、$(Y_i, D_i)$の単回帰モデルを考えてみます.

このとき、

• $\alpha = E[Y_{0i}]$
• $\beta = E[\theta_i]$
• $\eta_i = Y_{0i} - E[Y_{0i}]$

と対応付けて解釈することができます. The treatment statusごとにそれぞれのthe CEFを評価すると

よって、

とATE + selection biasの関係性が見えてきます. 仮に $\eta_i \perp D_i$が成立している場合、(1-3)のselection biasの項が消えるので、単回帰モデルでATEが推定できることがわかります.

ATT and selection bias

• $y$: the outcome variable
• $\mathbf x$: a set of explanatrory variables
• $s$: a treatment status indicator
• $\mu_0(\mathbf x) = E[y|\mathbf x, s= 0]$
• $\mu_1(\mathbf x) = E[y|\mathbf x, s= 1]$

このとき、treatment statusに応じたグループ間の差分はATT + selection biasで表現できます.

次に、$\mu_s(\mathbf x) = \mathbf x\beta_s$とします.

2. Ordinary Least Squares Estimator

The CEFが$E[Y_i|X_i] = X_i’\beta$だと仮定します. このときの$\beta$を推定するためには

を解けばよいこととなります. しかし、the joint distribution $(Y_i, X_i)$は事前には知られていないので、我々はこのthe sample analogで$\beta$を推定します.

$\mathbf Y = (Y_1, \cdots, Y_N)’, \mathbf X = (X_1, \cdots, X_N)’$とすると、(2-2)は次のような最小化問題に書き換えられます

(2-3)について、$b$のFOCをとると

$\mathbf X$がfull rankであるならば$\mathbf X’\mathbf X$の逆行列がとれるので

残差を $e_i \equiv Y - X_i’ \hat\beta$と定義したとき、(2-5)より

(2-6)が意味することは

• $\mathbf X, \mathbf e$は直交
• $\mathbf X, \mathbf e$は無相関

The residualとThe error termの違い

上記(2-6)の内容を深堀したいと思います.

• $x_{1i}, x_{2i}, \epsilon_i$ are mutually independent standard normal random variables

という回帰モデルを考えます. Finite sampleのもとでは、OLS residualはregressorsと直交していますが、the true error termは必ずしも直交していません. Pythonコードで見ていきたいと思います.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import numpy as np
import pandas as pd
import seaborn as sns

## set seed
np.random.seed(42)

## Data Generation
const = np.array([1]*1000, ndmin=2).T
data = np.random.normal(size = (1000, 3))
data = np.concatenate([const, data], axis = 1)
coefficient_matrix = np.array([1, 1, 1, 1], ndmin=2).T

y = data @ coefficient_matrix

## OLS
X = data[:,:-1]
beta_ols = np.linalg.inv(X.T @ X) @ X.T @ y
residual = y - X @ beta_ols

## regressor and residual
observed_data = np.concatenate([data, residual], axis = 1)

## Covariance matrix
sns.heatmap(np.cov(observed_data, rowvar=False), annot=True, fmt='g');
sns.heatmap(np.cov(data, rowvar=False), annot=True, fmt='g');

Demeaned Regressors

というモデルをOLSで推定したいときの結果を以下のように表すとします:

残差とFOCの性質より、

(2-8)の式より

(2-10)を(2-9)の式に代入すると$\sum X_i(Y_i- \bar Y + (X_i’ - \bar X)\hat\beta)=0$を得る.

これを$\hat\beta$について解くと

• この推定量をthe demeaned formula for the least squares estimatorといいます
• OLS regressionで得られる推定量はdemeaned dataについてinterceptなしでOLS regressionを実行した場合に得られる推定量と一致します
• sample analogではstandard errorの数値が自由度の関係から微妙に異なってしまう

standard errorの比較: OLS estimator vs demeaned OLS estimator

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import statsmodels.api as sm
import numpy as np

## READ DATA
duncan_prestige = sm.datasets.get_rdataset("Duncan", "carData")

## preprocessing
#----- model 1
Y = duncan_prestige.data['income']
X = duncan_prestige.data['education']

X_naive = sm.add_constant(X)

#------ demeaned data
Y_demeaned = Y - Y.mean()
X_demeaned = X - X.mean()

## Fitting
model = sm.OLS(Y,X_naive)
model_demeaned = sm.OLS(Y_demeaned,X_demeaned)
results, results_demeaned = model.fit(), model_demeaned.fit()

print(results.summary())
print(results_demeaned.summary())

Then, The results are

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 income   R-squared:                       0.525
Model:                            OLS   Adj. R-squared:                  0.514
Method:                 Least Squares   F-statistic:                     47.51
Date:                xxx, dd mmm yyyy   Prob (F-statistic):           1.84e-08
Time:                        01:43:38   Log-Likelihood:                -190.42
No. Observations:                  45   AIC:                             384.8
Df Residuals:                      43   BIC:                             388.5
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         10.6035      5.198      2.040      0.048       0.120      21.087
education      0.5949      0.086      6.893      0.000       0.421       0.769
==============================================================================
Omnibus:                        9.841   Durbin-Watson:                   1.736
Prob(Omnibus):                  0.007   Jarque-Bera (JB):               10.609
Skew:                           0.776   Prob(JB):                      0.00497
Kurtosis:                       4.802   Cond. No.                         123.
==============================================================================



                              Demeaned OLS Regression Results                                
=======================================================================================
Dep. Variable:                 income   R-squared (uncentered):                   0.525
Model:                            OLS   Adj. R-squared (uncentered):              0.514
Method:                 Least Squares   F-statistic:                              48.62
Date:                xxx, dd mmm yyyy   Prob (F-statistic):                    1.25e-08
Time:                        01:43:17   Log-Likelihood:                         -190.42
No. Observations:                  45   AIC:                                      382.8
Df Residuals:                      44   BIC:                                      384.6
Df Model:                           1                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
education      0.5949      0.085      6.972      0.000       0.423       0.767
==============================================================================
Omnibus:                        9.841   Durbin-Watson:                   1.736
Prob(Omnibus):                  0.007   Jarque-Bera (JB):               10.609
Skew:                           0.776   Prob(JB):                      0.00497
Kurtosis:                       4.802   Cond. No.                         1.00
==============================================================================

Projection Matrix

次に線形代数の射影の観点からOLS estimatorの性質を考えます. $\mathbf Y = \mathbf X \beta + \epsilon$というモデルを考えます.

このとき、$\mathbf X$で張られるfieldへの射影を実行する射影行列(projection matrix) $\mathbf P$を以下のように定義します:

また残差行列 $M$を次のように定義します

射影行列の性質

1: $\mathbf X$を$\mathbf P$で射影すると自分自身になる

2: $\mathbf X = [X_1, X_2]$としたときをsubvector $X_1$を$\mathbf P$で射影すると自分自身になる

射影行列$P$について

表記の簡略化のため

このとき

よって,

(2-13)の結果を用いると

Residual Regression: FWL theorem

というOLS estimatorを考えます. ここで$X_1$を$(X_2, \cdots X_K)$に回帰します.

(2-14)の両辺に $\hat V’$を左から掛けます:

よって $\hat V’\hat V$の逆行列が存在する時

このような残差回帰で$\hat\beta$のsubvectorを推定できる性質のことをFWL定理といいます. なお、射影行列を用いると$\hat V$は以下のように表現できます

(2-17)に(2-18)を代入すると

(2-19)から$Y$を$(X_2, \cdots, X_K)$に回帰した残差を$X_1$を$(X_2, \cdots, X_K)$に回帰した残差すると$\hat\beta_1$が得られることがわかります. また、(25)の結果得られる残差と$\hat e$も一致します.

Coded example

というregression modelを考えます. このとき, $\beta$をフルモデルのOLS, 残差回帰それぞれで以下推定します.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
mport statsmodels.api as sm
import numpy as np

## set seed
np.random.seed(42)

## parameter
N = 10000
mean = [0, 0, 0, 0]
cov = [[1, -1, 2, 0], [-1, 4, -1, 1], [2, -1, 6, -2], [0, 1, -2, 4]]
beta = 1
gamma = np.array([[2, 3, 4]]).T

## data generation
data = np.random.multivariate_normal(mean, cov, N)
x, z = np.hsplit(data,[1,])
u = np.random.normal(0, 1, N).reshape(N, 1)
y = x * beta + z @ gamma + u

## residuals: regress (y, x) on z and get each of the residuals
x_resid = x - z @ np.linalg.inv(z.T @ z) @ z.T @ x
y_resid = y - z @ np.linalg.inv(z.T @ z) @ z.T @ y

## estimation
reg_ols = sm.OLS(y,data)
results = reg_ols.fit()

reg_resid_1 = sm.OLS(y,x_resid)
results_resid_1 = reg_resid_1.fit()

reg_resid_2 = sm.OLS(y_resid,x_resid)
results_resid_2 = reg_resid_2.fit()

Then,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                      y   R-squared (uncentered):                   0.990
Model:                            OLS   Adj. R-squared (uncentered):              0.990
Method:                 Least Squares   F-statistic:                          2.466e+05
Date:                xxx, dd mmm yyyy   Prob (F-statistic):                        0.00
Time:                        03:49:15   Log-Likelihood:                         -14191.
No. Observations:               10000   AIC:                                  2.839e+04
Df Residuals:                    9996   BIC:                                  2.842e+04
Df Model:                           4                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1             0.9261      0.061     15.187      0.000       0.807       1.046
x2             1.9852      0.014    139.290      0.000       1.957       2.013
x3             3.0339      0.023    129.403      0.000       2.988       3.080
x4             4.0271      0.016    253.591      0.000       3.996       4.058
==============================================================================

                                 Residual Regression Part 1 Results                                
=======================================================================================
Dep. Variable:                      y   R-squared (uncentered):                   0.000
Model:                            OLS   Adj. R-squared (uncentered):              0.000
Method:                 Least Squares   F-statistic:                              2.315
Date:                xxx, dd mmm yyyy   Prob (F-statistic):                       0.128
Time:                        03:49:21   Log-Likelihood:                         -37200.
No. Observations:               10000   AIC:                                  7.440e+04
Df Residuals:                    9999   BIC:                                  7.441e+04
Df Model:                           1                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1             0.9261      0.609      1.521      0.128      -0.267       2.119
==============================================================================


                              Residual Regression Part 2 Results                                
=======================================================================================
Dep. Variable:                      y   R-squared (uncentered):                   0.023
Model:                            OLS   Adj. R-squared (uncentered):              0.022
Method:                 Least Squares   F-statistic:                              230.7
Date:                xxx, dd mmm yyyy   Prob (F-statistic):                    1.57e-51
Time:                        03:49:29   Log-Likelihood:                         -14191.
No. Observations:               10000   AIC:                                  2.838e+04
Df Residuals:                    9999   BIC:                                  2.839e+04
Df Model:                           1                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1             0.9261      0.061     15.189      0.000       0.807       1.046
==============================================================================

• Degree of Freedomがモデルによって異なるので$\hat\beta_1$の値が一致してもstandard errorが異なる
• R-squaredも $1 - \text{RSS}/\text{TSS} = \text{ESS}/\text{TSS}$という定義より、値が異なってくる

練習問題: Partially Linear Model, Robinson Estomator 1988

• $g(\cdot)$: a general function of $\mathbf x$
• $\beta$: $ 1 \times M$ vector

というモデルを考えます. このとき、以下を証明せよ

• $\tilde y \equiv y - E[y|\mathbf x]$
• $\tilde y \equiv \mathbf z - E[\mathbf z|\mathbf x]$

証明

$E[y|\mathbf x]$についてまず考えます. (28)より$E[u|\mathbf x] = E[E[y|\mathbf x, \mathbf z]|\mathbf x] =0$. 従って、

(26) - (29)より

この推定量は$E[y|\mathbf x], E[\mathbf z|\mathbf x]$をnonparametricに推定し、それを用いて$y, \mathbf z$のresidualをそれぞれ計算. そのresidualsを用いてOLS regressionをするという形で計算します. この推定量から得られる$\hat\beta$は$\sqrt N$-consistent, asymptotically normalであることが知られています.

Python simulation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
mport numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.sandbox.regression.predstd import wls_prediction_std

import warnings
warnings.filterwarnings('ignore')

from matplotlib import style
style.use("fivethirtyeight")

## set seed
np.random.seed(42)

## parameter
N = 1000
mean = [0, 0]
cov = [[1, -1], [-1, 3]]
beta = 1.5
gamma = np.array([[0.25, 0.75]]).T
lottery = np.random.uniform(0, 1, N)

## data generation
Z = np.random.multivariate_normal(mean, cov, N)/2
z_1, z_2 = np.hsplit(Z, [1,])
D_hat = np.exp(Z @ gamma)/ (1 + np.exp(Z @ gamma))
D = np.int64(lottery < D_hat.flatten()).reshape(N, 1)
e = np.random.normal(0, 1, N).reshape(N, 1)
y = D * beta + np.exp(z_1)**2 + np.sin(z_1 * np.pi) - z_1 **3 + np.cos(z_2 * np.pi)**2 + e

## DataFrame
df = pd.DataFrame(np.concatenate([y, D_hat, D, z_1, z_2], axis = 1), columns = ['y', 'D_hat', 'D', 'z_1', 'z_2'])

次に２グループの平均単純比較をします

1
2
3
4
5
6
## Naive estimator
naive = df.groupby('D')['y'].mean().values
naive = np.insert(naive, 2, naive[1] - naive[0])
res = pd.DataFrame([naive])
res.columns = ['control mean', 'treated mean', 'mean diff']
print(res)

Then

   control mean  treated mean  mean diff
0      2.159487      3.543704   1.384217

OLS estimatesは

1
2
model_naive = smf.ols('y ~ D + z_1 + z_2', data=df).fit()
model_naive.summary()

Then

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.567
Model:                            OLS   Adj. R-squared:                  0.566
Method:                 Least Squares   F-statistic:                     434.8
Date:                xxx, dd mmm yyyy   Prob (F-statistic):          1.71e-180
Time:                        23:45:28   Log-Likelihood:                -1945.3
No. Observations:                1000   AIC:                             3899.
Df Residuals:                     996   BIC:                             3918.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.9697      0.078     25.165      0.000       1.816       2.123
D              1.6919      0.112     15.043      0.000       1.471       1.913
z_1            3.7839      0.133     28.550      0.000       3.524       4.044
z_2            0.0196      0.080      0.247      0.805      -0.137       0.176
==============================================================================

最後にRobinson estimatorを実施します（簡易版）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
## Robinson Estimator
#--- kernel density estimator object
kernel_D = stats.gaussian_kde(np.vstack([D.ravel(), z_1.ravel(),z_2.ravel()]), bw_method='silverman')
kernel_Y = stats.gaussian_kde(np.vstack([y.ravel(), z_1.ravel(),z_2.ravel()]), bw_method='silverman')

#--- input data for kernel desity estimation
input_array = df.loc[:, ['z_1', 'z_2']].values
D_input_array_0 = np.hstack([np.zeros(N).reshape(N, 1), input_array])
D_input_array_1 = np.hstack([np.ones(N).reshape(N, 1), input_array])
y_range = np.linspace(np.round(min(y), 2)-0.01, np.round(max(y)+0.01, 2), 3000)
y_sample = len(y_range)
y_kernel_hat = []

#--- density estimation
for row in input_array:
    covariate_row_reshaped = np.tile(row, y_sample).reshape(y_sample, 2)
    input_array_y = np.hstack([y_range.reshape(-1, 1), covariate_row_reshaped]).T
    
    prob_array = kernel_Y(np.array(input_array_y.tolist()))
    conditional_mean = np.sum(y_range.ravel() * prob_array)/np.sum(prob_array)
    
    y_kernel_hat.append(conditional_mean)

df_fix = df.copy()
df_fix['y_kernel_hat'] = np.array(y_kernel_hat)
df_fix['d_kernel_hat'] = kernel_D(D_input_array_1.T)/(kernel_D(D_input_array_0.T) + kernel_D(D_input_array_1.T))
y_residual = df_fix['y'] - df_fix['y_kernel_hat']
d_residual = df_fix['D'] - df_fix['d_kernel_hat']

reg_robinson = sm.OLS(y_residual,d_residual)
results = reg_robinson.fit(use_t = True, cov_type='HC2')
print(results.summary())

Then,

                                OLS Regression Results                                
=======================================================================================
Dep. Variable:                      y   R-squared (uncentered):                   0.328
Model:                            OLS   Adj. R-squared (uncentered):              0.327
Method:                 Least Squares   F-statistic:                              505.5
Date:                xxx, dd mmm yyyy   Prob (F-statistic):                    6.58e-91
Time:                        23:41:23   Log-Likelihood:                         -1477.6
No. Observations:                1000   AIC:                                      2957.
Df Residuals:                     999   BIC:                                      2962.
Df Model:                           1                                                  
Covariance Type:                  HC2                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1             1.5968      0.071     22.483      0.000       1.457       1.736
==============================================================================

• もともとのthe error termの分散によってkernel density estimationの精度が大きく変わってしまう点に留意

Leverage Values

• $X: n\times k$ matrix
• $E[\epsilon|X] = 0$
• $\text{rank}(X) = k$

というモデルを考えます. 射影行列の議論より、

(31)の射影行列はthe hat matrixともいい、

(31)と(33)を比較すると、

つまり、$h_{ii}$はthe observed outcome, $\hat y_i$ が予測値に与える影響を数値化したものといえるため、$h_{ii}$ のことをthe leverageといいます.

REMARKS

• $h_{ii}$はデータポイントのx値とn個のデータポイントすべてのx値の平均との間の距離を表す指標と考えることもできます
• Outlier detectionで活用されることもあります
• $h_{ii}$ について簡単な計算をすると以下のことが示せます:

Python: The leverageとcovariateの関係

「$h_{ij}$はデータポイントのx値とn個のデータポイントすべてのx値の平均との間の距離を表す指標」であることをPythonを用いて確認したいと思います.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf

## generating data
data = [
        [1, 0.1, -0.0716],
        [2, 0.45401, 4.1673],
        [3, 1.09765,	6.5703],
        [4, 1.27936, 13.815],
        [5, 2.20611, 11.4501],
        [6, 2.50064, 12.9554],
        [7, 3.0403, 20.1575],
        [8, 3.23583, 17.5633],
        [9, 4.45308, 26.0317],
        [10, 4.1699, 22.7573],
        [11, 5.28474, 26.303],
        [12, 5.59238, 30.6885],
        [13, 5.92091, 33.9402],
        [14, 6.66066, 30.9228],
        [15, 6.79953, 34.11],
        [16, 7.97943, 44.4536],
        [17, 8.41536, 46.5022],
        [18, 8.71607, 50.0568],
        [19, 8.70156, 46.5475],
        [20, 9.16463, 45.7762],
        [21, 4, 40]
        ]
df = pd.DataFrame(data, columns = ['row', 'x', 'y'])
df['const'] = 1

## compute the Hat matrix
X = df.loc[:, ['const', 'x']].values.reshape(-1, 2)
H_matrix = X @ np.linalg.inv(X.T @ X) @ X.T
df['h_leverage'] = np.diag(H_matrix)

## the leverageとx valueの関係
ns.scatterplot(x='x', y='h_leverage', data=df)
plt.axvline(np.mean(df['x']), color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = plt.ylim()
plt.text(np.mean(df['x'])*1.1, max_ylim*0.9, 'x Mean: {:.2f}'.format(np.mean(df['x'])));

Then,

the leverageの３つの性質

Theorem

1. $\sum_{i}^N h_{ii} = k$
2. $0 \leq h_{ii} \leq 1$
3. $h_{ii} \geq 1/n$ if an intercept is included

Proof (1)

Proof (2)

$H$は $H = HH$という性質を持っているので

従って、$h_{ii} \geq 0$. また、(1)より $\sum h_{ii} = k \leq n$なので$h_{ii}\leq 1$

Proof (3)

$X = [1, \hat X_i]$の場合を考えます.

Leave-One-Out Regression

• $X_i$: k vector for ith observation
• $E[\epsilon|X] = 0$

というモデルを考えます. これをLeave-One-Out Regression(LOO)という手法で推定します. LOOとは、自分自身のobservationを除外したsub-sampleを用いてパラメータを推定し、そのパラメータを用いて自分自身のoutcomeを予測するという手法です.

$\hat\beta$と$\tilde\beta_{-i}$の関係

• $\tilde e_i$: LOOを用いて得られたresidual
• $\hat e_i$: OLSで得られたresidual
• $h_{i}$: the hat matrixのi-th diagonal element

このとき、いかが成り立ちます:

(38)の証明

まず以下の定理を所与とします：

• A: non-singular matrix
• b: a vector
• $\lambda$: a scalar

If

Then,

次に、(36)で確認した$\tilde\beta_{-i}$計算の各項を見ていきます.

(39)の両辺に右から$X_{i}$を掛けると

$\hat\beta$の関係式より

より

(37)の証明

(40), (41)より

$\hat e_i$と$\tilde e_{i}$の関係

(38)より

これを踏まえるとLOOのprediction errorとOLS MSEの関係性は以下のように書けます.

(42)の解釈としては

• the prediction standard errorはすべてのサンプルを用いで学習するOLS standard error(ある種のtraining error)よりも大きい

REMARKS

• OLS estimatorはinfluencial observationの影響を強く受ける推定量なので、training dataのoutlier detectionを実施する際に用いられることがある
• ただし、outlierを除外することは常に良いとは限らない（例：outlierが真のモデルに従っている場合など）

Goodness of fit measure

Estimation of Error Variance

というモデルに従ってデータが生成されたとします. このとき仮に$e_i$が観測されるならば, $\sigma^2$の推定量として

次に、$e_i$が観測されない場合は OLS residual $\hat e_i$を用いて以下のように推定します

(43), (44)を比較すると

• (45)は射影行列はpositive semi definiteより
• the feasible estimator $\hat \sigma^2$は数値的に理想的なestimatorよりも小さいことがわかる

Analysis of Variance

というregression modelを考えます. $(x_i, u_i)$はuncorrelatedの場合

(46)を変形すると

• $y_i$の変動は$x_i’\beta$で説明できる割合と$u_i$で説明できる割合に分けて考えることができることを(47)は示しています
• 特に, $\text{Var}(x_i’\beta)/\text{Var}(y_i)$をGoodness of fit, $\mathbf R^2$ として一般的には用います

次にこれをOLS estimatorの観点から考えてみます.

(48)から $\mathbf 1_n \bar Y$を引くと

Then,

よって、

$R^2$の解釈の留意点

• Econometricsにおける分析の主眼点はtreament effectのパラメーターをconsistentに推定することであって、$R^2$の数値の高さ自体はあくまでa secondary importance
• $R^2$の数値はregressorを加えれば加えるほど高くなる性質があり、regressorsの選択がtrue modelに基づいているか否かの情報は提供してくれない
• regressorを加えれば加えるほど、the residualの推定の自由度が下がるため、$R^2$の数値の数値のぶれが大きくなる
• causal modelの評価には有用ではない

二点目の性質に対処したスコアとしてAdjusted R-squaredという評価スコアがあります

• K: the number of regressors including the constant term
• TSSとSSRをそれぞれの自由度で除したものを計算し、それらの比率をもってGoodness of fitを計算しようというモチベーション

RegressorとR-squaredの関係

$R^2$の数値はregressorを加えれば加えるほど高くなることを確認します.

という２つのOLS regression modelを考えます. (52), (53)を合わせると

OLS regresionのresidualの性質より

(54), (55)から次式を得ます

(56)における$(X(\beta - \beta_0) - X_1\beta_1)’(X(\beta - \beta_0) - X_1\beta_1) \geq 0$より

3. Property of OLS Estimator

Assumptions

Unbiasedness

Proof (1)

Proof (2)

$D’D$は$d = Da$とすると、任意の0でないベクトル$a$について、$d’d \geq 0$が成立するので、$D’D \geq 0$

4. Asymptotic Property of OLS Estimator