作者：chen_h
微訊號 & QQ：862251340
微信公眾號：coderpai

簡介

線性迴歸是分析兩個或者多個變數之間關係的標準工具

在本文中，我們將使用 Python 包 statsmodels 來估計，解釋和視覺化線性迴歸模型。在此過程中，我們將討論各種主題，包括：

• 簡單和多元線性迴歸
• 視覺化
• 內生性和遺漏變數偏差
• 二階最小二乘法

簡單線性迴歸

我們希望確定機構中的差異是否有助於解釋觀察到的經濟結果。我們如何衡量制度差異和經濟結果？

在文字中，經濟結果用 1995 年人均國內生產總值代表，並根據匯率進行調整。制度差異由政治風險服務小組建立的 1985 - 1995 年平均徵收保護指數代表。本文中使用的這些變數和其他資料可在 Daron Acemoglu 的網頁上下載。

我們將使用 pandas 的 .read_stata() 函式將 .dta 檔案中包含的資料讀入資料幀。

import pandas as pd

df1 = pd.read_stata('https://github.com/QuantEcon/QuantEcon.lectures.code/raw/master/ols/maketable1.dta')
df1.head()

讓我們用一個散點圖來看一下人均 GDP 與徵收指數保護之間是否存在明顯的關係。

import matplotlib.pyplot as plt
plt.style.use('seaborn')

df1.plot(x='avexpr', y='logpgp95', kind='scatter')
plt.show()

該圖顯示了徵收與人均 GDP 之間存在相當強的正相關關係。

具體而言，如果對徵收的更高保護是衡量制度質量的指標，那麼更好的制度似乎與更好的經濟結果正相關（人均 GDP 更高）。

鑑於這種情節，選擇線性模型來描述這種關係似乎是一個合理的假設。

我們可以把我們的模型寫成：

$logpgp95_{i} = \beta_{0} + \beta_{1}avexpr_{i} + \mu_{i}$

其中，

• $\beta_{0}$ 是 y 軸上線性趨勢線的截距；
• $\beta_{1}$ 是線性趨勢線的斜率，表示保護風險對人均 GDP 的邊際效應
• $\mu_{i}$ 是隨機誤差項（由於未包含在模型中的因素，觀測值與線性趨勢的偏差）

在視覺上，這個線性模型涉及選擇最合適資料的直線，如下：

import numpy as np

# Dropping NA's is required to use numpy's polyfit
df1_subset = df1.dropna(subset=['logpgp95', 'avexpr'])

# Use only 'base sample' for plotting purposes
df1_subset = df1_subset[df1_subset['baseco'] == 1]

X = df1_subset['avexpr']
y = df1_subset['logpgp95']
labels = df1_subset['shortnam']

# Replace markers with country labels
plt.scatter(X, y, marker='')

for i, label in enumerate(labels):
    plt.annotate(label, (X.iloc[i], y.iloc[i]))

# Fit a linear trend line
plt.plot(np.unique(X),
         np.poly1d(np.polyfit(X, y, 1))(np.unique(X)),
         color='black')

plt.xlim([3.3,10.5])
plt.ylim([4,10.5])
plt.xlabel('Average Expropriation Risk 1985-95')
plt.ylabel('Log GDP per capita, PPP, 1995')
plt.title('Figure 2: OLS relationship between expropriation risk and income')
plt.show()

估計線性模型引數（$\beta$）的最常用技術是普通最小二乘法（OLS）。

顧名思義，OLS 模型通過找到最小化殘差平方和的引數來解決。即：

$\underset{\hat{\beta}}{min} \sum_{i=1}^{N}\hat{\mu_{i}}^{2}$

其中，$\hat{\mu}$ 觀察值和因變數預測值之間的差異。

為了估計常數項 $\beta_{0}$ ，我們需要在資料集中新增一列 1（如果用 $\beta_{0}$ 和 $x_{i}=1$ 替代 $\beta_{0}$，則考慮方程式 ）

df1['const'] = 1

現在我們可以使用 OLS 函式在 statsmodels 中構建我們的模型。

我們將 pandas 函式幀與 statsmodels 一起使用，但是標準陣列也可以用作參考。

import statsmodels.api as sm

reg1 = sm.OLS(endog=df1['logpgp95'], exog=df1[['const', 'avexpr']], missing='drop')
type(reg1)

statsmodels.regression.linear_model.OLS

到目前為止，我們已經構建了我們的模型。我們需要使用 .fit() 來獲得引數估計 $\hat{\beta_{0}}$ 和 $\hat{\beta_{1}}$ 。

results = reg1.fit()
type(results)

statsmodels.regression.linear_model.RegressionResultsWrapper

我們現在將擬合的迴歸模型儲存在結果中。要檢視 OLS 迴歸結果，我們可以呼叫 .summary() 方法。

print(results.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:               logpgp95   R-squared:                       0.611
Model:                            OLS   Adj. R-squared:                  0.608
Method:                 Least Squares   F-statistic:                     171.4
Date:                Mon, 19 Nov 2018   Prob (F-statistic):           4.16e-24
Time:                        17:03:03   Log-Likelihood:                -119.71
No. Observations:                 111   AIC:                             243.4
Df Residuals:                     109   BIC:                             248.8
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.6261      0.301     15.391      0.000       4.030       5.222
avexpr         0.5319      0.041     13.093      0.000       0.451       0.612
==============================================================================
Omnibus:                        9.251   Durbin-Watson:                   1.689
Prob(Omnibus):                  0.010   Jarque-Bera (JB):                9.170
Skew:                          -0.680   Prob(JB):                       0.0102
Kurtosis:                       3.362   Cond. No.                         33.2
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

從結果中，我們可以看到：

• 截距 $\hat{\beta_{0}} = 4.63$
• 斜率 $\hat{\beta_{1}} = 0.53$
• 正的 $\hat{\beta_{1}}$ 引數估計意味著制度質量對經濟結果有積極影響
• $\hat{\beta_{1}}$ 的 p 值為 0.000 意味著制度對 GDP 的影響具有統計顯著性（使用 p< 0.05 作為拒絕規則）
• R平方值 0.611

使用我們的引數估計，我們現在可以將我們的估計關係寫為：

$\hat{logpgp95_{i}} = 4.63+0.53avexpr_{i}$

上述的方程是對我們資料的最好的擬合，正如圖二描述的。

我們能使用這個等式去對我們的資料進行預測，比如當我們的 avexpr 為 7.07 時，我們就可以得到對數 GDP 為 8.38。

mean_expr = np.mean(df1_subset['avexpr'])
mean_expr

6.515625

predicted_logpdp95 = 4.63 + 0.53 * 7.07
predicted_logpdp95

8.3771

獲得此結果的一種更加簡單（也更加準確）的方法是使用 .predict() 並設定 constat=1 和 avexpri = mean_expr 。

results.predict(exog=[1, mean_expr])

array([8.09156367])

我們可以通過在結果上呼叫 .predict() 來獲取資料集中 avexpri 的每個值的預測 logpgp95i 陣列。

根據 avexpri 繪製預測值，表明預測值位於我們上面擬合的線性線上。

為了比較目的，還繪製了 logpgp95i 的觀察值。

# Drop missing observations from whole sample

df1_plot = df1.dropna(subset=['logpgp95', 'avexpr'])

# Plot predicted values

plt.scatter(df1_plot['avexpr'], results.predict(), alpha=0.5, label='predicted')

# Plot observed values

plt.scatter(df1_plot['avexpr'], df1_plot['logpgp95'], alpha=0.5, label='observed')

plt.legend()
plt.title('OLS predicted values')
plt.xlabel('avexpr')
plt.ylabel('logpgp95')
plt.show()