【译】使用python做计量之线性回归

这是用python做计量经济学实证分析系列第一篇。

原文地址：https://python.quantecon.org/ols.html

[scode type="green"]本系列文章对原文进行了一定的修改，主要是作图工具从matplotlib变为plotly.express，因为后者更为简洁明了，如果你不习惯plotly，可以看原文的代码[/scode]

[scode type="yellow"]本系列文章需要你对计量经济学有一定的了解，并且熟悉python、pandas、numpy[/scode]

引入

文献：Acemoglu D, Johnson S, Robinson J A. The colonial origins of comparative development: An empirical investigation[J]. American economic review, 2001, 91(5): 1369-1401.

地址：https://economics.mit.edu/files/4123

本文强调了制度在经济发展中的重要性。其主要贡献是利用殖民者死亡率作为制度差异外生变化的来源。这样的变化是必要的，以确定制度是否导致了更快的经济增长，而不是相反。

安装相应库：

pip install linearmodels

导入相应库

import numpy as np
import pandas as pd
import plotly.express as px
import plotly
import statsmodels.api as sm
from statsmodels.iolib.summary2 import summary_col
from linearmodels.iv import IV2SLS
plotly.offline.init_notebook_mode(connected=True)

简单的线性回归

Acemoglu研究的是制度差异是否导致了经济产出的差异。但是，我们如何衡量制度和产出呢？

在论文中，

经济产出以经汇率调整后的1995年人均GDP的对数代表。
制度差异以Political Risk Services Group构建的1985-95年期间的平均征收保护指数（protection against expropriation on average）代表。

数据可以在Daron Acemoglu的个人网站上下载。https://economics.mit.edu/faculty/acemoglu/data/ajr2001

我们使用read_stata() 读取.dta格式文件。

df1 = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/ols/maketable1.dta')
print(df1.head())

对logGDP和平均征收保护指数（avexpr）做散点图

fig = px.scatter(df1,x='avexpr',y='logpgp95')
plotly.offline.plot(fig)

上图显示了征收保护指数与人均GDP的对数之间相当强的正相关关系。具体来说，如果更高的征收保护是衡量制度质量的一个指标，那么更好的制度似乎与更好的经济成果（人均国内生产总值）呈正相关。选择一个线性模型来描述这种关系似乎是合理的假设。

我们可以把模型写成

$$\operatorname{logpgp} 95_{i}=\beta_{0}+\beta_{1} \text { avexpr }_{i}+u_{i}$$

其中，

$\beta_{0}$是趋势线在y轴上的截距
$\beta_{1}$是线性趋势线的斜率，表示征收保护指数对log人均GDP的边际影响
u是一个随机误差项（由于模型中未包含的因素，观测值与线性趋势有偏差）

绘制趋势线

# 删除空值
df1_subset = df1.dropna(subset=['logpgp95', 'avexpr'])
# Use only 'base sample' for plotting purposes
df1_subset = df1_subset[df1_subset['baseco'] == 1]
fig = px.scatter(df1_subset, x="avexpr", y="logpgp95", trendline="ols",
                 title={'text': "Figure 2: OLS relationship between expropriation risk and income",'y':0.9,'x':0.5},
                 text="shortnam")
fig.update_traces(textposition='top center')
fig.update_xaxes(title_text='Average Expropriation Risk 1985-95')
fig.update_yaxes(title_text='Log GDP per capita, PPP, 1995')

plotly.offline.plot(fig)

最常用的估计$\beta$的方法是最小二乘法（OLS），顾名思义，OLS模型是通过寻找使残差平方和最小化的参数来求解的。即：

$$\min _{\hat{\beta}} \sum_{i=1}^{N} \widehat{u}_{i}^{2}$$

为了求解常数项$\beta_0$，我们需要在数据集中添加一列1（即用$\beta_0x_i,x_i=1$来替换$\beta_0$）

df1['const'] = 1

现在我们可以使用statsmodels中的OLS函数构建我们的模型。

reg1 = sm.OLS(endog=df1['logpgp95'], exog=df1[['const', 'avexpr']], \
    missing='drop')
results = reg1.fit() # 获取参数的估计值
print(results.summary()) # 打印回归结果

从上图中，我们可以看到

截距$\hat \beta_0=4.63$
斜率$\hat\beta_1 = 0.53$
$\hat\beta_1$为正，意味着制度质量对经济成果有正向作用。
$\hat\beta_1$的 p-value为0，意味着上述结论在统计学上是显著的（使用p<0.05为拒绝条件）
R平方值为0.611表明，人均国内生产总值的对数61%左右的变化可以由征收保护指数解释。

使用我们的参数估计，我们现在可以将我们的估计关系写为

$$\hat{logpgp} 95_{i}=4.63+0.53 \text { avexpr }_{i}$$

我们可以用这个方程来预测征收保护指数指数对的人均GDP水平的影响。

例如，对于一个指数值为7.07（数据集的平均值）的国家，我们发现其1995年人均对数国内生产总值的预测水平为8.38。

mean_expr = np.mean(df1_subset['avexpr'])
predicted_logpdp95 = 4.63 + 0.53 * 7.07

获得这个结果的一个更简单（更准确）的方法是使用predict()方法，设置常数项constant = 1 and $avexpr_i= mean\_expr$

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

同样地，我们可以使用predict()预测所有值作图并与实际值比较。

# Drop missing observations from whole sample
df1_plot = df1.dropna(subset=['logpgp95', 'avexpr'])
df1_plot['predict'] = results.predict()
df_long=pd.melt(df1_plot, id_vars=['avexpr'], value_vars=['logpgp95', 'predict'])
fig = px.scatter(df_long, x='avexpr', y='value', color='variable',
                 title={'text': "OLS predict",'y':0.9,'x':0.5},
                 labels={'logpgp95':'observed'},symbol="variable")

fig.update_yaxes(title_text='logpgp95')
plotly.offline.plot(fig)

线性回归的扩展-多元回归

到目前为止，我们只考虑了影响经济表现的制度因素——几乎可以肯定的是，还有许多其他影响GDP的因素没有包括在我们的模型中。剔除影响logpgp95的其他变量将导致遗漏变量偏差（omitted variable bias），产生有偏差和不一致的参数估计。我们可以通过加入其他可能影响logpgp95i的因素，将我们的二元回归模型扩展到多元回归模型。

考虑其他因素例如：

气候对经济表现的影响，用纬度来代表气候
经济表现和制度（如文化、历史等）的差异；用大陆虚拟变量（dummies）来控制

让我们使用maketable2.dta中的数据估计本文中考虑的一些扩展模型

df2 = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/ols/maketable2.dta')
df2['const'] = 1
# 创建参数列表
X1 = ['const', 'avexpr']
X2 = ['const', 'avexpr', 'lat_abst']
X3 = ['const', 'avexpr', 'lat_abst', 'asia', 'africa', 'other']
# OLS回归
reg1 = sm.OLS(df2['logpgp95'], df2[X1], missing='drop').fit()
reg2 = sm.OLS(df2['logpgp95'], df2[X2], missing='drop').fit()
reg3 = sm.OLS(df2['logpgp95'], df2[X3], missing='drop').fit()

现在我们已经拟合了模型，我们将使用summary_col在一个表中显示结果（模型编号与本文中的编号相对应）

info_dict={'R-squared': lambda x:f"{x.rsquared:.2f}", 
           'No.observations': lambda x: f"{int(x.nobs):d}"}
# f是去掉转义字符，nobs就是总观测数量
results_table = summary_col(results=[reg1,reg2,reg3],float_format='%0.2f',
                            stars = True,model_names=['Model 1','Model 3','Model 4'],
                            info_dict=info_dict,
                            regressor_order=['const','avexpr','lat_abst','asia','africa'])

results_table.add_title('Table 2 - OLS Regressions')
print(results_table)

内生性-两阶段最小二乘

OLS模型可能存在内生性问题，从而导致模型估计的偏差和不一致。也就是说，制度与经济成果之间可能存在双向关系：

较富裕的国家可能负担得起或更喜欢更好的制度
影响收入的变量也可能与制度差异相关
指数的构建可能有失偏颇：研究者可能倾向于研究收入较高、拥有更好制度的国家

为了处理内生性，我们可以使用两阶段最小二乘（2SLS）回归，这是OLS回归的扩展。此方法需要用以下变量替换内生变量avexpr：

与avexpr相关
与误差项无关（即不应直接影响因变量，否则由于omitted variable bias，他会与$u_i$相关；）

新的回归变量被称为工具变量，其目的是消除我们的制度差异代理中的内生性。原文献的主要贡献是利用殖民者的死亡率来衡量制度差异。他们假设，殖民者的高死亡率导致建立了更具extractive in nature的制度（对征用的保护更少，这些制度至今仍然存在）。

使用散点图，我们可以看到征收保护与殖民者死亡率呈负相关，符合作者的假设并满足有效工具变量的第一个条件。

df1_subset2 = df1.dropna(subset=['logpgp95', 'avexpr'])
# Use only 'base sample' for plotting purposes
df1_subset2 = df1.dropna(subset=['logem4', 'avexpr'])
fig = px.scatter(df1_subset2, x="logem4", y="avexpr", trendline="ols",
                 title={'text': "Figure 3: First-stage relationship between settler mortality and expropriation risk",'y':0.9,'x':0.5},
                 text="shortnam")
fig.update_traces(textposition='top center')
fig.update_xaxes(title_text='Log of Settler Mortality')
fig.update_yaxes(title_text='Average Expropriation Risk 1985-95')

plotly.offline.plot(fig)

如果17至19世纪的殖民者死亡率对当前的国内生产总值有直接影响（除了通过制度产生的间接影响外），第二个条件可能无法满足例如，定居者死亡率可能与一个国家当前的疾病环境有关，这可能影响当前的经济表现。原文献认为这不太可能，因为：

大多数定居者死亡是由于疟疾和黄热病，对当地居民的影响有限。
例如，非洲或印度当地人民的疾病负担似乎没有高于平均水平，这是殖民前这些地区相对较高的人口密度造成的。

这似乎是一个有效的工具，我们可以使用2SLS回归获得一致和无偏的参数估计。

第一阶段

第一阶段的因变量是内生变量（$avexpr_i$）。工具是我们模型中所有外生变量的集合（而不仅仅是我们替换的变量）。以模型为例，我们的工具只是一个常数和殖民者死亡率对数$\log e m 4_{i}$，因此，我们将第一阶段回归估计为

$$\text { avexpr }_{i}=\delta_{0}+\delta_{1} \log e m 4_{i}+v_{i}$$

我们需要估计这个方程的数据位于maketable4.dta中（仅用baseco = 1的数据进行估算）

df4 = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/ols/maketable4.dta')
df4 = df4[df4['baseco'] == 1]
# Add a constant variable
df4['const'] = 1
# Fit the first stage regression and print summary
results_fs = sm.OLS(df4['avexpr'],df4[['const', 'logem4']],missing='drop').fit()
print(results_fs.summary())

二阶段

我们需要使用predict()检索$\text {avexpr}_{i}$的预测值。然后，我们将原始线性模型中的内生变量$\text {avexpr}_{i}$替换为预测值$\widehat {avexpr}_{i}$。

因此，我们的第二阶段回归为：

$$\operatorname{logpgp} 95_{i}=\beta_{0}+\beta_{1} a \widehat{v e x p} r_{i}+u_{i}$$

df4['predicted_avexpr'] = results_fs.predict()
results_ss = sm.OLS(df4['logpgp95'],df4[['const', 'predicted_avexpr']]).fit()
print(results_ss.summary())

第二阶段的回归结果为我们提供了一个制度对经济结果的影响的无偏且一致的估计。结果表明，与OLS结果相比，两者之间存在更强的正相关关系。

请注意，虽然我们的参数估计是正确的，但标准误不是，因此，不建议“手动”（使用OLS分阶段）计算2SLS。我们可以使用linearmodels包一步正确地进行2SLS回归。

注意，当使用IV2SLS时，外生变量和工具变量在函数参数中是分开的

iv = IV2SLS(dependent=df4['logpgp95'],exog=df4['const'],endog=df4['avexpr'],
            instruments=df4['logem4']).fit(cov_type='unadjusted')
print(iv.summary)

鉴于我们现在有一致和无偏见的估计，我们可以从我们估计的模型中推断，制度差异（源自殖民时期建立的制度）有助于解释当今各国收入水平的差异。

原文献使用0.94的边际效应来计算智利和尼日利亚之间的指数差异（即制度质量）意味着收入差异高达7倍，强调了制度在经济发展中的重要性。

版权属于：作者名称
本文链接：https://www.sitstars.com/archives/88/
转载时须注明出处及本声明

页面

分类

【译】使用python做计量之线性回归

引入

简单的线性回归

线性回归的扩展-多元回归

内生性-两阶段最小二乘

添加新评论取消回复

【译】使用python做计量之线性回归