Python-因果推断-全-
Python 因果推断（全）
原文：
causal-methods.github.io/Book/Introduction.html译者：飞龙
引言
原文：
causal-methods.github.io/Book/Introduction.html译者：飞龙
作者：Vitor Kamada
电子邮件：econometrics.methods@gmail.com
最后更新日期：2020 年 8 月 15 日
这本书是使用 Python 进行因果推断的实用指南。我解释了出现在经济学最负盛名的期刊，如《美国经济评论》和《计量经济学》中的方法和技术。
我不假设任何技术背景，但我建议您熟悉我之前的书中的概念：Probability and Statistics with Python。
计量经济学课程材料
我的本科、硕士和博士课程中的教学大纲、幻灯片/笔记、Python 和 R 代码可以在我的另一个Github上找到。
一、国际象棋选手比其他人更理性吗？
译者：飞龙
电子邮件：econometrics.methods@gmail.com
最后更新日期：2020 年 10 月 28 日
看一下下面来自 Palacios-Huerta & Volij（2009）的蜈蚣游戏图。每个玩家只有两种策略：“停止”或“继续”。如果玩家 1 在第一轮停止，他会得到 4 美元，玩家 2 会得到 1 美元。如果玩家 1 继续，玩家 2 可以选择停止或继续。如果玩家 2 停止，他将获得 8 美元，玩家 1 将获得 2 美元。如果玩家 2 继续，玩家 1 再次可以决定“停止”或“继续”。从社会角度来看，最好是两名玩家在六轮中都选择“继续”，然后玩家 1 可以获得 256 美元，玩家 2 可以获得 64 美元。然而，如果玩家 2 是理性的，他永远不会在第 6 轮继续玩，“因为如果他停止，他可以获得 128 美元。知道通过反向归纳，玩家 1 在第 1 轮继续是不理性的。
来源：Palacios-Huerta & Volij（2009）
几项实验研究表明，几乎没有人在第一次机会停下来。根据 Palacios-Huerta & Volij（2009）的下表，只有 7.5%的学生在第一轮停下来。大多数学生（35%）在第 3 轮停下来。在 McKelvey & Palfrey（1992）中可以找到大样本量的类似结果。
来源：Palacios-Huerta & Volij（2009）
让我们打开 Palacios-Huerta & Volij（2009）的数据集，其中包含有关国际象棋选手如何玩蜈蚣游戏的信息。
import pandas as pd
path = "https://github.com/causal-methods/Data/raw/master/"
sheet_name = "Data Chess Tournaments"
chess = pd.read_excel(path + "Chess.xls", sheet_name)
chess
| 比赛 | 标题 1 | ELORating1 | 标题 2 | ELORating2 | 终节点 | |
|---|---|---|---|---|---|---|
| 0 | 贝纳斯克 | GM | 2634 | 0 | 2079 | 1 |
| 1 | 贝纳斯克 | GM | 2621 | 0 | 2169 | 1 |
| 2 | 贝纳斯克 | GM | 2562 | FM | 2307 | 1 |
| 3 | 贝纳斯克 | GM | 2558 | 0 | 2029 | 1 |
| 4 | 贝纳斯克 | GM | 2521 | 0 | 2090 | 1 |
| ... | ... | ... | ... | ... | ... | ... |
| 206 | 塞斯塔奥 | GM | 2501 | 0 | 2019 | 1 |
| 207 | 塞斯塔奥 | 0 | 2050 | 0 | 2240 | 2 |
| 208 | 塞斯塔奥 | 0 | 2020 | 0 | 2004 | 1 |
| 209 | 塞斯塔奥 | 0 | 2133 | GM | 2596 | 1 |
| 210 | 塞斯塔奥 | 0 | 2070 | FM | 2175 | 1 |
211 行×6 列
让我们统计每个类别中的国际象棋选手人数。在样本中，我们可以看到 26 名国际大师（GM），29 名国际大师（IM）和 15 名联邦大师（FM）。
chess['Title1'].value_counts()
0 130
IM 29
GM 26
FM 15
WGM 6
WIM 3
WFM 1
CM 1
Name: Title1, dtype: int64
让我们将所有其他国际象棋选手标记为“其他”。
dictionary = { 'Title1':
{ 0: "Other",
'WGM': "Other",
'WIM': "Other",
'WFM': "Other",
'CM': "Other" }}
chess.replace(dictionary, inplace=True)
我们在“其他”类别中有 141 名国际象棋选手。
chess['Title1'].value_counts()
Other 141
IM 29
GM 26
FM 15
Name: Title1, dtype: int64
让我们关注标有 1 的国际象棋选手的头衔和评分。让我们忽略他们对手的头衔和评分，也就是标有 2 的国际象棋选手。
title_rating1 = chess.filter(["Title1", "ELORating1"])
title_rating1
| 标题 1 | ELORating1 | |
|---|---|---|
| 0 | GM | 2634 |
| 1 | GM | 2621 |
| 2 | GM | 2562 |
| 3 | GM | 2558 |
| 4 | GM | 2521 |
| ... | ... | ... |
| 206 | GM | 2501 |
| 207 | 其他 | 2050 |
| 208 | 其他 | 2020 |
| 209 | 其他 | 2133 |
| 210 | 其他 | 2070 |
211 行×2 列
我们可以看到国际大师的平均评分为 2513，而国际大师的平均评分为 2411。
请注意，评分是预测谁将赢得国际象棋比赛的很好的指标。
极不可能在“其他”类别中的国际象棋选手击败国际大师。
评分可以被理解为理性和反向归纳的代理。
# Round 2 decimals
pd.set_option('precision', 2)
import numpy as np
title_rating1.groupby('Title1').agg(np.mean)
| ELORating1 | |
|---|---|
| 标题 1 | |
| --- | --- |
| FM | 2324.13 |
| GM | 2512.96 |
| IM | 2411.69 |
| 其他 | 2144.60 |
让我们将分析限制在只有国际大师。
grandmasters = chess[chess['Title1'] == "GM"]
grandmasters
| 比赛 | 标题 1 | ELORating1 | 标题 2 | ELORating2 | 终节点 | |
|---|---|---|---|---|---|---|
| 0 | 贝纳斯克 | GM | 2634 | 0 | 2079 | 1 |
| 1 | 贝纳斯克 | GM | 2621 | 0 | 2169 | 1 |
| 2 | Benasque | GM | 2562 | FM | 2307 | 1 |
| 3 | Benasque | GM | 2558 | 0 | 2029 | 1 |
| 4 | Benasque | GM | 2521 | 0 | 2090 | 1 |
| 5 | Benasque | GM | 2521 | FM | 2342 | 1 |
| 7 | Benasque | GM | 2510 | 0 | 2018 | 1 |
| 8 | Benasque | GM | 2507 | 0 | 2043 | 1 |
| 9 | Benasque | GM | 2500 | 0 | 2097 | 1 |
| 10 | Benasque | GM | 2495 | 0 | 2043 | 1 |
| 12 | Benasque | GM | 2488 | 0 | 2045 | 1 |
| 13 | Benasque | GM | 2488 | GM | 2514 | 1 |
| 15 | Benasque | GM | 2485 | 0 | 2092 | 1 |
| 16 | Benasque | GM | 2482 | 0 | 2040 | 1 |
| 18 | Benasque | GM | 2475 | 0 | 2229 | 1 |
| 19 | Benasque | GM | 2473 | 0 | 2045 | 1 |
| 35 | Benasque | GM | 2378 | 0 | 2248 | 1 |
| 156 | Leon | GM | 2527 | 0 | 2200 | 1 |
| 168 | Sestao | GM | 2671 | 0 | 2086 | 1 |
| 170 | Sestao | GM | 2495 | 0 | 2094 | 1 |
| 172 | Sestao | GM | 2532 | 0 | 2075 | 1 |
| 174 | Sestao | GM | 2501 | 0 | 2005 | 1 |
| 176 | Sestao | GM | 2516 | GM | 2501 | 1 |
| 177 | Sestao | GM | 2444 | 0 | 2136 | 1 |
| 193 | Sestao | GM | 2452 | 0 | 2027 | 1 |
| 206 | Sestao | GM | 2501 | 0 | 2019 | 1 |
所有 26 位特级大师都在第一轮停下了蜈蚣游戏。他们就是标准经济学教科书中描述的理性经济人。
不幸的是，在 78 亿人口的人口中，根据国际象棋联合会(FIDE)的说法，我们只有 1721 位特级大师。
grandmasters.groupby('EndNode').size()
EndNode
1 26
dtype: int64
让我们来检查一下国际大师(IM)。
international_master = chess[chess['Title1'] == "IM"]
# Return only 4 observations
international_master.head(4)
| Tournament | Title1 | ELORating1 | Title2 | ELORating2 | EndNode | |
|---|---|---|---|---|---|---|
| 6 | Benasque | IM | 2521 | 0 | 2179 | 1 |
| 11 | Benasque | IM | 2492 | 0 | 2093 | 1 |
| 14 | Benasque | IM | 2487 | IM | 2474 | 2 |
| 17 | Benasque | IM | 2479 | 0 | 2085 | 1 |
并非所有国际大师都在第一轮停下。
5 位国际大师在第二轮停下，还有 2 位在第三轮停下。
node_IM = international_master.groupby('EndNode').size()
node_IM
EndNode
1 22
2 5
3 2
dtype: int64
这 5 位在第二轮停下的国际大师代表了总国际大师人数的 17%。
只有 76%的国际大师按照新古典经济理论的预测行事。
length_IM = len(international_master)
prop_IM = node_IM / length_IM
prop_IM
EndNode
1 0.76
2 0.17
3 0.07
dtype: float64
让我们对联邦大师应用相同的程序。
在每个节点停下的联邦大师的比例与国际大师相似。
federation_master = chess[chess['Title1'] == "FM"]
node_FM = federation_master.groupby('EndNode').size()
length_FM = len(federation_master)
prop_FM = node_FM / length_FM
prop_FM
EndNode
1 0.73
2 0.20
3 0.07
dtype: float64
让我们将先前的描述性统计数据放在一个条形图中。
更容易地可视化每个结束蜈蚣游戏的节点上国际象棋选手的比例。
条形图表明，国际大师和联邦大师在玩蜈蚣游戏时的方式不同于国际大师和联邦大师。然而，国际大师和联邦大师在玩蜈蚣游戏时的方式看起来是相似的。
import plotly.graph_objects as go
node = ['Node 1', 'Node 2', 'Node 3']
fig = go.Figure(data=[
go.Bar(name='Grandmasters', x=node, y=[1,0,0]),
go.Bar(name='International Masters', x=node, y=prop_IM),
go.Bar(name='Federation Masters', x=node, y=prop_FM) ])
fig.update_layout(barmode='group',
title_text = 'Share of Chess Players that Ended Centipede Game at Each Node',
font=dict(size=21) )
fig.show()
让我们正式地测试一下特级大师(\(p_{g}\))的比例是否等于国际大师(\(p_{i}\))的比例。零假设(\(H_0\))是：
z 统计量为：
其中\(se(\hat{p}_{g}-\hat{p}_{i})\)是样本比例之间的标准误差：
其中\(n_g\)是特级大师的样本量，\(n_i\)是国际大师的样本量。
对于节点 1，我们知道：
然后：
z 统计量的 p 值为 0.007。因此，在显著性水平\(\alpha=1\%\)上拒绝零假设(\(H_0\))。
from statsmodels.stats.proportion import proportions_ztest
# I inserted manually the data from Grandmasters to
# ilustrate the input format
count = np.array([ 26, node_IM[1] ]) # number of stops
nobs = np.array([ 26, length_IM ]) # sample size
proportions_ztest(count, nobs)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
(2.68162114289528, 0.00732663816543149)
让我们还在节点 1 测试一下，国际大师(\(p_{i}\))的比例是否等于联邦大师(\(p_{f}\))的比例。零假设(\(H_0\))是：
z 统计量为 0.18，相应的 p 值为 0.85。因此，我们无法拒绝国际大师比例等于联邦大师比例的零假设。
count = np.array([ node_IM[1], node_FM[1] ])
nobs = np.array([ length_IM, length_FM ])
proportions_ztest(count, nobs)
(0.1836204648065827, 0.8543112075010346)
练习
1）使用 Palacios-Huerta & Volij（2009）的电子表格“Data Chess Tournaments”来计算停留在节点 1、节点 2、节点 3、节点 4 和节点 5 的其他国际象棋选手的比例。比例将总和为 100%。其他国际象棋选手是一个排除所有国际象棋选手（包括特级大师、国际大师和联邦大师）的类别。
2）这个问题涉及 Palacios-Huerta & Volij（2009）的电子表格“Data UPV Students-One shot”。
a) 在 Google Colab 中打开电子表格“Data UPV Students-One shot”。
b) 巴斯克大学（UPV）的学生玩了多少对蜈蚣游戏？
c) 计算停留在节点 1、节点 2、节点 3、节点 4、节点 5 和节点 6 的学生比例。比例将总和为 100%。
3）比较练习 1（国际象棋选手）和练习 2（c）（学生）的结果。为什么这两个子群体玩蜈蚣游戏的方式不同？推测并证明你的推理。
4）使用 Palacios-Huerta & Volij（2009）的电子表格“Data Chess Tournaments”来测试在蜈蚣游戏的节点 3 停止的国际大师比例是否等于在节点 3 停止的其他国际象棋选手的比例。
5）对于这个问题，使用 Palacios-Huerta & Volij（2009）的电子表格“Data Chess Tournaments”。创建一个条形图，比较国际大师和其他国际象棋选手在每个节点停止蜈蚣游戏的比例。
6）假设你是一位新古典经济学家。你如何证明在更强的理性和自利的假设下建立的标准经济理论与 Palacios-Huerta & Volij（2009）的论文中呈现的实证证据相矛盾？详细说明你的理由。
参考文献
Fey, Mark, Richard D. McKelvey, and Thomas R. Palfrey. (1996). An Experimental Study of Constant-Sum Centipede Games. International Journal of Game Theory, 25(3): 269–87.
Palacios-Huerta, Ignacio, and Oscar Volij. (2009). Field Centipedes. American Economic Review, 99 (4): 1619-35.
二、白人名字有助于获得工作面试吗？
原文：
causal-methods.github.io/Book/2%29_Does_a_White_Sounding_Name_Help_to_Get_Job_Interview.html译者：飞龙
电子邮件：econometrics.methods@gmail.com
最近更新：2020 年 9 月 15 日
让我们加载 Bertrand & Mullainathan（2004）的数据集。
import pandas as pd
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "lakisha_aer.dta")
df.head(4)
| id | ad | education | ofjobs | yearsexp | honors | volunteer | military | empholes | occupspecific | ... | compreq | orgreq | manuf | transcom | bankreal | trade | busservice | othservice | missind | ownership | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | b | 1 | 4 | 2 | 6 | 0 | 0 | 0 | 1 | 17 | ... | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |
| 1 | b | 1 | 3 | 3 | 6 | 0 | 1 | 1 | 0 | 316 | ... | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |
| 2 | b | 1 | 4 | 1 | 6 | 0 | 0 | 0 | 0 | 19 | ... | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |
| 3 | b | 1 | 3 | 4 | 6 | 0 | 1 | 0 | 1 | 313 | ... | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
4 行×65 列
让我们将分析限制在变量'call'和'race'上。
回拨：1 = 应聘者被召回面试；否则为 0。
种族：w = 白人，b = 黑人。
callback = df.loc[:, ('call', 'race')]
callback
| 回拨 | 种族 | |
|---|---|---|
| 0 | 0.0 | w |
| 1 | 0.0 | w |
| 2 | 0.0 | b |
| 3 | 0.0 | b |
| 4 | 0.0 | w |
| ... | ... | ... |
| 4865 | 0.0 | b |
| 4866 | 0.0 | b |
| 4867 | 0.0 | w |
| 4868 | 0.0 | b |
| 4869 | 0.0 | w |
4870 行×2 列
让我们计算按种族分组的观察数量（大小）和变量'call'的平均值。
我们有相同数量（2435）的黑人和白人申请者的简历。
只有 6.4%的黑人收到了回拨；而 9.7%的白人收到了回拨。
因此，白人申请者获得面试回拨的可能性约高出 50%。
换句话说，每 10 份白人申请者发送的简历中，才能得到 1 次面试机会，黑人申请者需要发送 15 份简历才能得到相同的结果。
# Round 2 decimals
pd.set_option('precision', 4)
import numpy as np
callback.groupby('race').agg([np.size, np.mean])
| 回拨 | |
|---|---|
| 大小 | |
| --- | --- |
| 种族 | |
| --- | --- |
| b | 2435.0 |
| w | 2435.0 |
有人可能会争辩说，这 3.2%的差异（9.65 - 6.45）不一定意味着对黑人的歧视。
你可以争辩说，白人申请者获得更多回拨，是因为他们拥有更多的教育、经验、技能，而不是因为肤色。
具体来说，你可以争辩说，白人申请者更有可能获得回拨，是因为他们更有可能拥有大学学位（资格的信号），而不是因为他们是黑人。
如果你从美国人口中随机抽取样本或检查美国人口普查数据，你会发现黑人获得大学学位的可能性比白人小。这是一个不容置疑的事实。
让我们检查 Bertrand & Mullainathan（2004）的数据集中拥有大学学位的黑人和白人的比例。
最初，大学毕业生在变量'education'中被编码为 4，3 = 一些大学，2 = 高中毕业，1 = 一些高中，0 = 未报告。
让我们创建变量'college' = 1，如果一个人完成了大学学位；否则为 0。
df['college'] = np.where(df['education'] == 4, 1, 0)
我们可以看到 72.3%的黑人申请者拥有大学学位。白人拥有大学学位的比例非常相似，为 71.6%。
为什么这些数字不代表美国人口，而且这些值彼此更接近？
因为数据不是从现实中随机抽取的样本。
college = df.loc[:, ('college', 'race')]
college.groupby('race').agg([np.size, np.mean])
| 大学 | |
|---|---|
| 大小 | |
| --- | --- |
| 种族 | |
| --- | --- |
| b | 2435 |
| w | 2435 |
Bertrand＆Mullainathan（2004）产生了实验数据。他们创建了简历。他们将黑人听起来的名字（例如：Lakish 或 Jamal）随机分配给一半简历，将白人听起来的名字（例如：Emily 或 Greg）分配给另一半。
通过姓名对种族的随机化使得两个类别白人和黑人在所有可观察和不可观察的因素上都相等（非常相似）。
让我们检查简历中其他因素的陈述。变量的名称是不言自明的，更多信息可以通过阅读 Bertrand＆Mullainathan（2004）的论文获得。
resume = ['college', 'yearsexp', 'volunteer', 'military',
'email', 'workinschool', 'honors',
'computerskills', 'specialskills']
both = df.loc[:, resume]
both.head()
| college | yearsexp | volunteer | military | workinschool | honors | computerskills | specialskills | ||
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 6 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
| 1 | 0 | 6 | 1 | 1 | 1 | 1 | 0 | 1 | 0 |
| 2 | 1 | 6 | 0 | 0 | 0 | 1 | 0 | 1 | 0 |
| 3 | 0 | 6 | 1 | 0 | 1 | 0 | 0 | 1 | 1 |
| 4 | 0 | 22 | 0 | 0 | 1 | 1 | 0 | 1 | 0 |
让我们使用不同的代码来计算整个样本（白人和黑人）以及黑人和白人之间的分样本变量的平均值。
请注意，整个样本的平均工作经验（yearsexp）为 7.84，黑人为 7.83，白人为 7.86。
如果您检查所有变量，您会发现黑人的平均值非常接近白人的平均值。这是随机化的结果。
我们还计算标准偏差（std），这是平均值周围变化的度量。请注意，整个样本和分样本之间的标准偏差几乎相同。就像平均值一样，在实验数据的情况下，你不应该看到标准偏差之间有太大的差异。
年龄经验变量的标准偏差为 5 年。我们可以粗略地说，大部分观察（约 68%）在平均值的 1 个标准差以下和 1 个标准差以上之间，即在[2.84，12.84]之间。
black = both[df['race']=='b']
white = both[df['race']=='w']
summary = {'mean_both': both.mean(), 'std_both': both.std(),
'mean_black': black.mean(), 'std_black': black.std(),
'mean_white': white.mean(), 'std_white': white.std()}
pd.DataFrame(summary)
| mean_both | std_both | mean_black | std_black | mean_white | std_white | |
|---|---|---|---|---|---|---|
| college | 0.7195 | 0.4493 | 0.7228 | 0.4477 | 0.7162 | 0.4509 |
| yearsexp | 7.8429 | 5.0446 | 7.8296 | 5.0108 | 7.8563 | 5.0792 |
| volunteer | 0.4115 | 0.4922 | 0.4144 | 0.4927 | 0.4086 | 0.4917 |
| military | 0.0971 | 0.2962 | 0.1018 | 0.3025 | 0.0924 | 0.2897 |
| 0.4793 | 0.4996 | 0.4797 | 0.4997 | 0.4789 | 0.4997 | |
| workinschool | 0.5595 | 0.4965 | 0.5610 | 0.4964 | 0.5581 | 0.4967 |
| honors | 0.0528 | 0.2236 | 0.0513 | 0.2207 | 0.0542 | 0.2265 |
| computerskills | 0.8205 | 0.3838 | 0.8324 | 0.3735 | 0.8086 | 0.3935 |
| specialskills | 0.3287 | 0.4698 | 0.3273 | 0.4693 | 0.3302 | 0.4704 |
为什么我们如此关心上表，显示平均白人和平均黑人申请者几乎相同？
因为白人申请者更有可能因为他们更高的教育水平而获得回电的论点是站不住脚的，如果白人和黑人申请者两组都具有类似的教育水平。
无论是未观察到的因素还是无法测量的因素，比如动机、心理特征等，都不能用来证明回电率的不同。
在一个实验中，只有处理变量（种族）是外生操纵的。其他一切都保持不变，因此结果变量（回电）的变化只能归因于处理变量（种族）的变化。
因此，实验研究消除了观察研究中出现的所有混杂因素。
实验是硬科学中的金标准。宣称因果关系的最严格方式。
调查、人口普查、直接从现实中提取的观察数据不能用来建立因果关系。它可能有助于捕捉关联，但不能证明因果关系。
形式上，我们可以写出下面的因果模型。这 3 行是等价的。我们只有在处理变量（\(T\)）被随机化时，才能声称\(\beta\)具有“因果”解释。在没有随机化的情况下，\(\beta\)只捕捉到相关性。
让我们使用普通最小二乘（OLS）方法估计上述模型。
在 Python 的 stasmodels 库中，截距是一个值为 1 的常数。
让我们创建变量“处理”：1 = 黑人申请者，0 = 白人申请者。
变量‘call’是结果变量（Y）。
df['Intercept'] = 1
df['Treatment'] = np.where(df['race'] =='b', 1, 0)
import statsmodels.api as sm
ols = sm.OLS(df['call'], df[['Intercept', 'Treatment']],
missing='drop').fit()
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
import pandas.util.testing as tm
让我们打印结果。
print(ols.summary().tables[1])
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.0965 0.006 17.532 0.000 0.086 0.107
Treatment -0.0320 0.008 -4.115 0.000 -0.047 -0.017
==============================================================================
现在我们可以将拟合模型写为：
我们已经在本节开头的代码中得到了上述系数，我将再次重现：
callback.groupby('race').agg([np.size, np.mean])
| 回电 | |
|---|---|
| 大小 | |
| --- | --- |
| 种族 | |
| --- | --- |
| b | 2435.0 |
| w | 2435.0 |
请注意，截距的值为 9.65%。这是接到面试回电的白人申请者的比例。
处理变量的系数为 3.2%。解释是，作为黑人申请者“导致”接到面试回电的几率减少了 3.2%（6.45% - 9.65%）。
请记住，3.2%是一个很大的幅度，因为它代表了大约 50%的差异。在实际情况下，黑人申请者需要发送 15 份简历才能获得一个面试机会，而白人申请者只需要发送 10 份简历。
处理变量的系数在显著水平（\(\alpha\) = 1%）上也是统计学上显著的。
-4.115 的 t 值是比率：
零假设是：
-4 的 t 值意味着观察值（-3.2%）比均值（\(\beta=0\)）低 4 个标准偏差。得到这个值的概率或概率几乎为 0。因此，我们拒绝零假设，即处理的幅度为 0。
什么定义了一个实验？
处理变量（T）的随机化。
这自动使处理变量（T）独立于其他因素：
在实验中，回归中其他因素的添加不会影响处理变量（\(\beta\)）的估计。如果您看到\(\beta\)有实质性的变化，您可以推断您没有使用实验数据。
请注意，在观察性研究中，您必须始终控制其他因素。否则，您将面临遗漏变量偏差问题。
让我们估计下面的多元回归：
other_factors = ['college', 'yearsexp', 'volunteer', 'military',
'email', 'workinschool', 'honors',
'computerskills', 'specialskills']
multiple_reg = sm.OLS(df['call'],
df[['Intercept', 'Treatment'] + other_factors],
missing='drop').fit()
我们可以看到处理（-3.1%）的系数并没有像预期的那样随着额外的控制变量而发生太大变化。
print(multiple_reg.summary().tables[1])
==================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------
Intercept 0.0547 0.015 3.727 0.000 0.026 0.083
Treatment -0.0311 0.008 -4.032 0.000 -0.046 -0.016
college 0.0068 0.009 0.768 0.443 -0.010 0.024
yearsexp 0.0029 0.001 3.630 0.000 0.001 0.005
volunteer -0.0032 0.011 -0.295 0.768 -0.024 0.018
military -0.0033 0.014 -0.232 0.817 -0.032 0.025
email 0.0143 0.011 1.285 0.199 -0.008 0.036
workinschool 0.0008 0.009 0.093 0.926 -0.016 0.018
honors 0.0642 0.018 3.632 0.000 0.030 0.099
computerskills -0.0202 0.011 -1.877 0.061 -0.041 0.001
specialskills 0.0634 0.009 7.376 0.000 0.047 0.080
==================================================================================
练习
1）在种族歧视的文献中，每个实验研究都有超过 1000 个观察性研究。假设您阅读了 100 个观察性研究，表明种族歧视是真实的。假设您还阅读了 1 个声称没有种族歧视证据的实验研究。您更倾向于接受 100 个观察性研究的结果还是一个实验研究的结果？请解释您的答案。
2）当组间均值差异存在偏差且未捕捉到平均因果效应时？请使用数学方程式来解释您的答案。
3）解释下面列出的列联表的 4 个值。具体地说明含义并比较这些值。
变量‘h’：1 = 更高质量的简历；0 = 更低质量的简历。这个变量是随机化的。
其他变量已经被定义。
contingency_table = pd.crosstab(df['Treatment'], df['h'],
values=df['call'], aggfunc='mean')
contingency_table
| h | 0.0 | 1.0 |
|---|---|---|
| 处理 | ||
| --- | --- | --- |
| 0 | 0.0850 | 0.1079 |
| 1 | 0.0619 | 0.0670 |
4）我创建了一个交互变量'h_Treatment'，它是变量'h'和'treatment'的成对乘法。
如何使用下面回归的系数来获得练习 3 中列联表的值？展示计算过程。
df['h_Treatment'] = df['h']*df['Treatment']
interaction = sm.OLS(df['call'],
df[['Intercept', 'Treatment', 'h', 'h_Treatment'] ],
missing='drop').fit()
print(interaction.summary().tables[1])
===============================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 0.0850 0.008 10.895 0.000 0.070 0.100
Treatment -0.0231 0.011 -2.094 0.036 -0.045 -0.001
h 0.0229 0.011 2.085 0.037 0.001 0.045
h_Treatment -0.0178 0.016 -1.142 0.253 -0.048 0.013
===============================================================================
5）我在没有交互项'h_Treatment'的情况下运行了下面的回归。我能使用下面的系数来获得练习 3 中列联表的值吗？如果可以，展示确切的计算过程。
interaction = sm.OLS(df['call'],
df[['Intercept', 'Treatment', 'h'] ],
missing='drop').fit()
print(interaction.summary().tables[1])
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.0894 0.007 13.250 0.000 0.076 0.103
Treatment -0.0320 0.008 -4.116 0.000 -0.047 -0.017
h 0.0141 0.008 1.806 0.071 -0.001 0.029
==============================================================================
6）编写一个代码以获得下面的列联表：
| 名字 h | 0.0 | 1.0 |
|---|---|---|
| Aisha | 0.010000 | 0.037500 |
| Allison | 0.121739 | 0.068376 |
| … | … | … |
表格内是由简历质量分解的回访率。Kristen 和 Lakisha 的回访率是多少？为什么回访率如此不同？我们能否证明回访率的差异，认为其中一个比另一个更受教育和合格？
7）使用 Bertrand＆Mullainathan（2004）的数据来测试白人和黑人是否具有相同的平均工作经验。陈述零假设。写出测试的数学公式。解释结果。
8）像 Bertrand＆Mullainathan（2004）一样打破常规。提出一个实际的方法来找出蓝眼睛和金发是否会导致更高的薪水。具体说明如何在实践中实施随机化策略。
参考
Bertrand, Marianne, and Sendhil Mullainathan.（2004）。Are Emily and Greg More Employable Than Lakisha and Jamal? A Field Experiment on Labor Market Discrimination。《美国经济评论》，94（4）：991-1013。
三、在伊斯兰或世俗政权下，女性更有可能完成高中吗？
译者：飞龙
电子邮件：econometrics.methods@gmail.com
最后更新：10-3-2020
让我们打开 Meyersson（2014）的数据。每一行代表土耳其的一个市镇。
# Load data from Meyersson (2014)
import numpy as np
import pandas as pd
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "regdata0.dta")
df.head()
| province_pre | ilce_pre | belediye_pre | province_post | ilce_post | belediye_post | vshr_islam2004 | educpop_1520f | hischshr614m | hischshr614f | ... | jhischshr1520m | jhivshr1520f | jhivshr1520m | rpopshr1520 | rpopshr2125 | rpopshr2630 | rpopshr3164 | nonagshr1530f | nonagshr1530m | anyc | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Adana | Aladag | Aladag | Adana | Aladag | Aladag | 0.367583 | 540.0 | 0.0 | 0.0 | ... | 0.478448 | 0.012963 | 0.025862 | 1.116244 | 1.113730 | 0.955681 | 0.954823 | 0.046778 | 0.273176 | 1.0 |
| 1 | Adana | Aladag | Akoren | Adana | Aladag | Akoren | 0.518204 | 342.0 | 0.0 | 0.0 | ... | 0.513089 | 0.023392 | 0.020942 | 1.002742 | 0.993227 | 1.093731 | 1.018202 | 0.020325 | 0.146221 | 0.0 |
| 2 | Adana | Buyuksehir | Buyuksehir | Adana | Buyuksehir | 0.397450 | 76944.0 | 0.0 | 0.0 | ... | 0.348721 | 0.036871 | 0.060343 | 1.006071 | 1.094471 | 1.039968 | 0.990001 | 0.148594 | 0.505949 | 1.0 | |
| 3 | Adana | Ceyhan | Sarimazi | Adana | Ceyhan | Sarimazi | 0.559827 | 318.0 | 0.0 | 0.0 | ... | 0.331343 | 0.022013 | 0.074627 | 1.124591 | 0.891861 | 0.816490 | 0.916518 | 0.040111 | 0.347439 | 0.0 |
| 4 | Adana | Ceyhan | Sagkaya | Adana | Ceyhan | Sagkaya | 0.568675 | 149.0 | 0.0 | 0.0 | ... | 0.503650 | 0.053691 | 0.043796 | 1.079437 | 1.208691 | 1.114033 | 0.979060 | 0.070681 | 0.208333 | 0.0 |
5 行×236 列
变量‘hischshr1520f’是根据 2000 年人口普查数据，15-20 岁女性完成高中的比例。不幸的是，年龄是被聚合的。15 和 16 岁的青少年不太可能在土耳其完成高中。最好是有按年龄分开的数据。由于 15 和 16 岁的人无法从分析中移除，15-20 岁女性完成高中的比例非常低：16.3%。
变量‘i94’是如果一个伊斯兰市长在 1994 年赢得了市镇选举则为 1，如果一个世俗市长赢得了则为 0。伊斯兰党在土耳其执政了 12%的市镇。
# Drop missing values
df = df.dropna(subset=['hischshr1520f', 'i94'])
# Round 2 decimals
pd.set_option('precision', 4)
# Summary Statistics
df.loc[:, ('hischshr1520f', 'i94')].describe()[0:3]
| hischshr1520f | i94 | |
|---|---|---|
| count | 2632.0000 | 2632.0000 |
| mean | 0.1631 | 0.1197 |
| std | 0.0958 | 0.3246 |
在由伊斯兰市长执政的市镇中，15-20 岁女性的平均高中完成率为 14%，而在由世俗市长执政的市镇中为 16.6%。
这是一个天真的比较，因为数据不是来自实验。市长类型并没有被随机化，实际上也不能被随机化。例如，贫困可能导致更高水平的宗教信仰和更低的教育成就。导致高中完成率较低的可能是贫困而不是宗教。
df.loc[:, ('hischshr1520f')].groupby(df['i94']).agg([np.size, np.mean])
| size | mean | |
|---|---|---|
| i94 | ||
| --- | --- | --- |
| 0.0 | 2317.0 | 0.1662 |
| 1.0 | 315.0 | 0.1404 |
图表“天真比较”显示控制组和实验组是基于变量‘iwm94’：伊斯兰赢得的边际。该变量被居中为 0。因此，如果赢得的边际高于 0，伊斯兰市长赢得了选举。另一方面，如果赢得的边际低于 0，伊斯兰市长输掉了选举。
就平均高中学历而言，治疗组（14%）和对照组（16.6%）之间的差异为-2.6%。使用观察数据比较市政结果的问题在于，治疗组与对照组并不相似。因此，混杂因素可能会使结果产生偏差。
import matplotlib.pyplot as plt
# Scatter plot with vertical line
plt.scatter(df['iwm94'], df['hischshr1520f'], alpha=0.2)
plt.vlines(0, 0, 0.8, colors='red', linestyles='dashed')
# Labels
plt.title('Naive Comparison')
plt.xlabel('Islamic win margin')
plt.ylabel('Female aged 15-20 with high school')
# Control vs Treatment
plt.text(-1, 0.7, r'$\bar{y}_{control}=16.6\%$', fontsize=16,
bbox={'facecolor':'yellow', 'alpha':0.2})
plt.text(0.2, 0.7, r'$\bar{y}_{treatment}=14\%$', fontsize=16,
bbox={'facecolor':'yellow', 'alpha':0.2})
plt.show()
由伊斯兰主导的高中学历与世俗主导的高中学历之间的 2.6%差异在 1%的显著水平上是统计学显著的。考虑到高中毕业率的平均值为 16.3%，这个幅度也是相关的。但是，请注意，这是一个天真的比较，可能存在偏见。
# Naive Comparison
df['Intercept'] = 1
import statsmodels.api as sm
naive = sm.OLS(df['hischshr1520f'], df[['Intercept', 'i94']],
missing='drop').fit()
print(naive.summary().tables[1])
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.1662 0.002 83.813 0.000 0.162 0.170
i94 -0.0258 0.006 -4.505 0.000 -0.037 -0.015
==============================================================================
判断天真比较是否可能存在偏见的一种方法是检查由伊斯兰主导的市政是否与由世俗主导的市政不同。
伊斯兰主要赢得的市政当中，1994 年伊斯兰选票比例更高（41% vs 10%），获得选票的政党数量更多（5.9 vs 5.5），对数人口更多（8.3 vs 7.7），19 岁以下人口比例更高（44% vs 40%），家庭规模更大（6.4 vs 5.75），区中心比例更高（39% vs 33%），省中心比例更高（6.6% vs 1.6%）。
df = df.rename(columns={"shhs" : "household",
"merkezi": "district",
"merkezp": "province"})
control = ['vshr_islam1994', 'partycount', 'lpop1994',
'ageshr19', 'household', 'district', 'province']
full = df.loc[:, control].groupby(df['i94']).agg([np.mean]).T
full.index = full.index.get_level_values(0)
full
| i94 | 0.0 | 1.0 |
|---|---|---|
| vshr_islam1994 | 0.1012 | 0.4145 |
| 党派数量 | 5.4907 | 5.8889 |
| 人口 1994 | 7.7745 | 8.3154 |
| ageshr19 | 0.3996 | 0.4453 |
| 家庭 | 5.7515 | 6.4449 |
| 区 | 0.3375 | 0.3937 |
| 省 | 0.0168 | 0.0667 |
使对照组和治疗组相似的一种方法是使用多元回归。对治疗变量'i94'的系数进行解释是ceteris paribus，也就是说，伊斯兰主要对高中学历的影响，考虑其他一切不变。这里的诀窍是“其他一切不变”，这意味着只有在回归中受控的因素。这是一个不完美的解决方案，因为在实际情况下，不可能控制影响结果变量的所有因素。然而，与简单回归相比，多元回归可能更少受到遗漏变量偏差的影响。
多元回归挑战了天真比较的结果。伊斯兰政权对高中毕业率的积极影响比世俗政权高出 1.4%。结果在 5%的显著水平上是统计学显著的。
multiple = sm.OLS(df['hischshr1520f'],
df[['Intercept', 'i94'] + control],
missing='drop').fit()
print(multiple.summary().tables[1])
==================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------
Intercept 0.2626 0.015 17.634 0.000 0.233 0.292
i94 0.0139 0.006 2.355 0.019 0.002 0.026
vshr_islam1994 -0.0894 0.013 -6.886 0.000 -0.115 -0.064
partycount -0.0038 0.001 -3.560 0.000 -0.006 -0.002
lpop1994 0.0159 0.002 7.514 0.000 0.012 0.020
ageshr19 -0.6125 0.021 -29.675 0.000 -0.653 -0.572
household 0.0057 0.001 8.223 0.000 0.004 0.007
district 0.0605 0.004 16.140 0.000 0.053 0.068
province 0.0357 0.010 3.499 0.000 0.016 0.056
==================================================================================
多元回归的结果看起来有些反直觉。治疗变量的符号是如何改变的？
让我们从另一个角度来看数据。图表“天真比较”是所有市政的散点图。每个点代表一个市政。很难看出任何模式或趋势。
让我们绘制相同的图形，但根据高中毕业率的相似性将市政聚合在 29 个箱子中。这些箱子在下图中是蓝色的球。球的大小与用于计算高中毕业率均值的市政数量成比例。
如果你仔细观察截断点附近（垂直红线），当变量伊斯兰胜利边际= 0 时，你会看到高中毕业率水平的不连续或跳跃。
# Library for Regression Discontinuity
!pip install rdd
Collecting rdd
Using cached rdd-0.0.3.tar.gz (4.4 kB)
Requirement already satisfied: pandas in c:\anaconda\envs\textbook\lib\site-packages (from rdd) (1.1.3)
Requirement already satisfied: numpy in c:\anaconda\envs\textbook\lib\site-packages (from rdd) (1.19.2)
Requirement already satisfied: statsmodels in c:\anaconda\envs\textbook\lib\site-packages (from rdd) (0.12.0)
Requirement already satisfied: pytz>=2017.2 in c:\anaconda\envs\textbook\lib\site-packages (from pandas->rdd) (2020.1)
Requirement already satisfied: python-dateutil>=2.7.3 in c:\anaconda\envs\textbook\lib\site-packages (from pandas->rdd) (2.8.1)
Requirement already satisfied: patsy>=0.5 in c:\anaconda\envs\textbook\lib\site-packages (from statsmodels->rdd) (0.5.1)
Requirement already satisfied: scipy>=1.1 in c:\anaconda\envs\textbook\lib\site-packages (from statsmodels->rdd) (1.5.3)
Requirement already satisfied: six>=1.5 in c:\anaconda\envs\textbook\lib\site-packages (from python-dateutil>=2.7.3->pandas->rdd) (1.15.0)
Building wheels for collected packages: rdd
Building wheel for rdd (setup.py): started
Building wheel for rdd (setup.py): finished with status 'done'
Created wheel for rdd: filename=rdd-0.0.3-py3-none-any.whl size=4723 sha256=e114b2f6b2da5c594f5856bd364a95d11da3dc2327c28b419668f2230f7c18c1
Stored in directory: c:\users\vitor kamada\appdata\local\pip\cache\wheels\f0\b8\ed\f7a5bcaa0a1b5d89d33d70db90992fd816ac6cff666020255d
Successfully built rdd
Installing collected packages: rdd
Successfully installed rdd-0.0.3
from rdd import rdd
# Aggregate the data in 29 bins
threshold = 0
data_rdd = rdd.truncated_data(df, 'iwm94', 0.99, cut=threshold)
data_binned = rdd.bin_data(data_rdd, 'hischshr1520f', 'iwm94', 29)
# Labels
plt.title('Comparison using aggregate data (Bins)')
plt.xlabel('Islamic win margin')
plt.ylabel('Female aged 15-20 with high school')
# Scatterplot
plt.scatter(data_binned['iwm94'], data_binned['hischshr1520f'],
s = data_binned['n_obs'], facecolors='none', edgecolors='blue')
# Red Vertical Line
plt.axvline(x=0, color='red')
plt.show()
也许你并不相信存在不连续或跳跃的截断点。让我们用 10 个箱子绘制相同的图形，并限制变量伊斯兰胜利边际的带宽（范围）。与选择任意带宽（h）不同，让我们使用 Imbens＆Kalyanaraman（2012）开发的方法来获得最小化均方误差的最佳带宽。
最佳带宽(\(\hat{h}\))为 0.23，也就是说，让我们在截断上下方创建一个 0.23 的窗口来创建 10 个箱子。
# Optimal Bandwidth based on Imbens & Kalyanaraman (2012)
# This bandwidth minimizes the mean squared error.
bandwidth_opt = rdd.optimal_bandwidth(df['hischshr1520f'],
df['iwm94'], cut=threshold)
bandwidth_opt
0.2398161605552802
以下是 10 个箱子。控制组中有 5 个箱子，伊斯兰赢得的优势<0，处理组中有 5 个箱子，伊斯兰赢得的优势>0。请注意，高中毕业率在索引 4 和 5 之间跳跃，分别为 13.8%和 15.5%，即箱子 5 和 6。13.8%和 15.5%的值是基于分别 141 和 106 个市政府('n_obs')计算的。
# Aggregate the data in 10 bins using Optimal Bandwidth
data_rdd = rdd.truncated_data(df, 'iwm94', bandwidth_opt, cut=threshold)
data_binned = rdd.bin_data(data_rdd, 'hischshr1520f', 'iwm94', 10)
data_binned
| 0 | hischshr1520f | iwm94 | n_obs | |
|---|---|---|---|---|
| 0 | 0.0 | 0.1769 | -0.2159 | 136.0 |
| 1 | 0.0 | 0.1602 | -0.1685 | 142.0 |
| 2 | 0.0 | 0.1696 | -0.1211 | 162.0 |
| 3 | 0.0 | 0.1288 | -0.0737 | 139.0 |
| 4 | 0.0 | 0.1381 | -0.0263 | 141.0 |
| 5 | 0.0 | 0.1554 | 0.0211 | 106.0 |
| 6 | 0.0 | 0.1395 | 0.0685 | 81.0 |
| 7 | 0.0 | 0.1437 | 0.1159 | 58.0 |
| 8 | 0.0 | 0.1408 | 0.1633 | 36.0 |
| 9 | 0.0 | 0.0997 | 0.2107 | 19.0 |
在图表“使用最佳带宽(h = 0.27)进行比较”中，蓝线适用于控制组(5 个箱子)，橙线适用于处理组(5 个箱子)。现在，不连续或跳跃是明显的。我们称这种方法为回归离散度(RD)。红色垂直线(\(\hat{\tau}_{rd}=3.5\)% )是高中毕业率的增加。请注意，这种方法模拟了一个实验。伊斯兰党几乎赢得和几乎输掉的市政府很可能相互类似。直觉是“几乎赢得”和“几乎输掉”是像抛硬币一样的随机过程。选举的相反结果可能是随机发生的。另一方面，很难想象伊斯兰市长会在以 30%的更大优势赢得的市政府中输掉。
# Scatterplot
plt.scatter(data_binned['iwm94'], data_binned['hischshr1520f'],
s = data_binned['n_obs'], facecolors='none', edgecolors='blue')
# Labels
plt.title('Comparison using Optimum Bandwidth (h = 0.27)')
plt.xlabel('Islamic win margin')
plt.ylabel('Female aged 15-20 with high school')
# Regression
x = data_binned['iwm94']
y = data_binned['hischshr1520f']
c_slope , c_intercept = np.polyfit(x[0:5], y[0:5], 1)
plt.plot(x[0:6], c_slope*x[0:6] + c_intercept)
t_slope , t_intercept = np.polyfit(x[5:10], y[5:10], 1)
plt.plot(x[4:10], t_slope*x[4:10] + t_intercept)
# Vertical Line
plt.vlines(0, 0, 0.2, colors='green', alpha =0.5)
plt.vlines(0, c_intercept, t_intercept, colors='red', linestyles='-')
# Plot Black Arrow
plt.axes().arrow(0, (t_intercept + c_intercept)/2,
dx = 0.15, dy =-0.06, head_width=0.02,
head_length=0.01, fc='k', ec='k')
# RD effect
plt.text(0.1, 0.06, r'$\hat{\tau}_{rd}=3.5\%$', fontsize=16,
bbox={'facecolor':'yellow', 'alpha':0.2})
plt.show()
C:\Anaconda\envs\textbook\lib\site-packages\ipykernel_launcher.py:25: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
# RD effect given by the vertical red line
t_intercept - c_intercept
0.03584571077550233
让我们将样本限制在伊斯兰市长以 5%的优势赢得或输掉的市政府。我们可以看到控制组和处理组之间的相似性比本章开头使用全样本进行比较更大。
然而，这种相似性并不接近于“完美实验”。部分原因是控制组和处理组的样本量较小。因此，当我们运行回归离散度时，建议添加控制变量。
# bandwidth (h) = 5%
df5 = df[df['iwm94'] >= -0.05]
df5 = df5[df5['iwm94'] <= 0.05]
sample5 = df5.loc[:, control].groupby(df5['i94']).agg([np.mean]).T
sample5.index = full.index.get_level_values(0)
sample5
| i94 | 0.0 | 1.0 |
|---|---|---|
| vshr_islam1994 | 0.3026 | 0.3558 |
| partycount | 5.9730 | 5.8807 |
| lpop1994 | 8.2408 | 8.2791 |
| ageshr19 | 0.4415 | 0.4422 |
| household | 6.2888 | 6.4254 |
| district | 0.4595 | 0.4037 |
| province | 0.0338 | 0.0826 |
让我们正式阐述回归离散度的理论。
让\(D_r\)成为一个虚拟变量：如果分析单位接受了处理，则为 1，否则为 0。下标\(r\)表示处理(\(D_r\))是运行变量\(r\)的函数。
在 Sharp 回归离散度中，处理(\(D_r\))由运行变量(\(r\))决定。
其中，\(r_0\)是一个任意的截断或阈值。
回归离散度的最基本规范是：
其中\(Y\)是结果变量，\(\beta_0\)是截距，\(\tau\)是处理变量(\(D_r\))的影响，\(\beta_1\)是运行变量(\(r\))的系数，\(\epsilon\)是误差项。
请注意，在实验中，处理是随机的，但在回归离散度中，处理完全由运行变量决定。随机过程的反面是确定性过程。这是反直觉的，但是当决定处理分配的规则(截断)是任意的时，确定性分配具有与随机化相同的效果。
一般来说，观察性研究的可信度非常低，因为遗漏变量偏差（OVB）的基本问题。误差项中的许多未观察因素可能与处理变量相关。因此，在回归框架中的一个大错误是将运行变量留在误差项中。
在所有估计器中，回归离散度可能是最接近黄金标准，随机实验的方法。主要缺点是回归离散度只捕捉局部平均处理效应（LATE）。将结果推广到带宽之外的实体是不合理的。
伊斯兰市长的影响是女性学校完成率高出 4％，使用带宽为 5％的回归离散度。这个结果在 5％的显著水平上是统计上显著的。
# Real RD specification
# Meyersson (2014) doesn't use the interaction term, because
# the results are unstable. In general the coefficient,
# of the interaction term is not statistically significant.
# df['i94_iwm94'] = df['i94']*df['iwm94']
# RD = ['Intercept', 'i94', 'iwm94', 'i94_iwm94']
RD = ['Intercept', 'i94', 'iwm94']
# bandwidth of 5%
df5 = df[df['iwm94'] >= -0.05]
df5 = df5[df5['iwm94'] <= 0.05]
rd5 = sm.OLS(df5['hischshr1520f'],
df5[RD + control],
missing='drop').fit()
print(rd5.summary())
OLS Regression Results
==============================================================================
Dep. Variable: hischshr1520f R-squared: 0.570
Model: OLS Adj. R-squared: 0.554
Method: Least Squares F-statistic: 36.32
Date: Wed, 28 Oct 2020 Prob (F-statistic): 1.67e-40
Time: 17:41:01 Log-Likelihood: 353.09
No. Observations: 257 AIC: -686.2
Df Residuals: 247 BIC: -650.7
Df Model: 9
Covariance Type: nonrobust
==================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------
Intercept 0.3179 0.043 7.314 0.000 0.232 0.403
i94 0.0399 0.016 2.540 0.012 0.009 0.071
iwm94 -0.4059 0.284 -1.427 0.155 -0.966 0.154
vshr_islam1994 -0.0502 0.060 -0.842 0.401 -0.168 0.067
partycount -0.0003 0.003 -0.074 0.941 -0.007 0.007
lpop1994 0.0091 0.005 1.718 0.087 -0.001 0.020
ageshr19 -0.7383 0.065 -11.416 0.000 -0.866 -0.611
household 0.0075 0.002 3.716 0.000 0.004 0.011
district 0.0642 0.010 6.164 0.000 0.044 0.085
province 0.0191 0.019 1.004 0.316 -0.018 0.057
==============================================================================
Omnibus: 17.670 Durbin-Watson: 1.658
Prob(Omnibus): 0.000 Jarque-Bera (JB): 19.403
Skew: 0.615 Prob(JB): 6.12e-05
Kurtosis: 3.546 Cond. No. 898.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
伊斯兰市长的影响是女性学校完成率高出 2.1％，使用回归离散度，最佳带宽为 0.27，根据 Imbens＆Kalyanaraman（2012）计算。这个结果在 1％的显著水平上是统计上显著的。
因此，回归离散度估计器表明，天真的比较在错误的方向上存在偏见。
# bandwidth_opt is 0.2715
df27 = df[df['iwm94'] >= -bandwidth_opt]
df27 = df27[df27['iwm94'] <= bandwidth_opt]
rd27 = sm.OLS(df27['hischshr1520f'],
df27[RD + control],
missing='drop').fit()
print(rd27.summary())
OLS Regression Results
==============================================================================
Dep. Variable: hischshr1520f R-squared: 0.534
Model: OLS Adj. R-squared: 0.530
Method: Least Squares F-statistic: 128.8
Date: Wed, 28 Oct 2020 Prob (F-statistic): 5.95e-161
Time: 17:41:01 Log-Likelihood: 1349.9
No. Observations: 1020 AIC: -2680.
Df Residuals: 1010 BIC: -2630.
Df Model: 9
Covariance Type: nonrobust
==================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------
Intercept 0.2943 0.023 12.910 0.000 0.250 0.339
i94 0.0214 0.008 2.775 0.006 0.006 0.036
iwm94 -0.0343 0.038 -0.899 0.369 -0.109 0.041
vshr_islam1994 -0.0961 0.030 -3.219 0.001 -0.155 -0.038
partycount -0.0026 0.002 -1.543 0.123 -0.006 0.001
lpop1994 0.0135 0.003 4.719 0.000 0.008 0.019
ageshr19 -0.6761 0.032 -20.949 0.000 -0.739 -0.613
household 0.0072 0.001 6.132 0.000 0.005 0.010
district 0.0575 0.006 10.364 0.000 0.047 0.068
province 0.0390 0.010 3.788 0.000 0.019 0.059
==============================================================================
Omnibus: 179.124 Durbin-Watson: 1.610
Prob(Omnibus): 0.000 Jarque-Bera (JB): 373.491
Skew: 1.001 Prob(JB): 7.90e-82
Kurtosis: 5.186 Cond. No. 270.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
练习
1）使用 Meyersson（2014）的数据运行回归离散度：a）使用完整样本，b）另一个带宽为 0.1（两侧均为 10％）。使用本章两个示例的相同规范。解释处理变量的系数。结果“a”或“b”更可信？请解释。
2）下面是变量伊斯兰赢得的边际的直方图。您是否看到截断或异常模式，其中截断= 0？在调查奔跑变量的截断周围是否发生了奇怪的事情的合理性是什么？
import plotly.express as px
fig = px.histogram(df, x="iwm94")
fig.update_layout(shapes=[
dict(
type= 'line',
yref= 'paper', y0 = 0, y1 = 1,
xref= 'x', x0 = 0, x1 = 0)])
fig.show()
3）我修改了“伊斯兰赢得的边际”变量以教育为目的。假设这是 Meyersson（2014）的真实运行变量。请参见下面的直方图。在这种假设情况下，您能从土耳其的选举中推断出什么？在这种情况下运行回归离散度存在问题吗？如果有，您可以采取什么措施来解决问题？
def corrupt(variable):
if variable <= 0 and variable >= -.025:
return 0.025
else:
return variable
df['running'] = df["iwm94"].apply(corrupt)
fig = px.histogram(df, x="running")
fig.update_layout(shapes=[
dict(
type= 'line',
yref= 'paper', y0 = 0, y1 = 1,
xref= 'x', x0 = 0, x1 = 0)])
fig.show()
4）为不熟悉因果推断的机器学习专家解释下面的图形？变量“伊斯兰选票份额”可以用作运行变量吗？推测。
def category(var):
if var <= 0.05 and var >= -.05:
return "5%"
else:
return "rest"
df['margin'] = df["iwm94"].apply(category)
fig = px.scatter(df, x="vshr_islam1994", y="iwm94", color ="margin",
labels={"iwm94": "Islamic win margin",
"vshr_islam1994": "Islamic vote share"})
fig.update_layout(shapes=[
dict(
type= 'line',
yref= 'paper', y0 = 1/2, y1 = 1/2,
xref= 'x', x0 = 0, x1 = 1)])
fig.show()
5）在伊斯兰或世俗政权下，男性更有可能完成高中吗？根据数据和严格的分析来证明您的答案。变量“hischshr1520m”是 15-20 岁男性中高中教育的比例。
参考
Imbens，G.，＆Kalyanaraman，K.（2012）。回归离散度估计器的最佳带宽选择。经济研究评论，79（3），933-959。
Meyersson，Erik。 （2014）。伊斯兰统治与穷人和虔诚者的赋权。经济计量学，82（1），229-269。
四、新教徒是否比天主教徒更喜欢休闲？
原文：
causal-methods.github.io/Book/4%29_Do_Protestants_Prefer_Less_Leisure_than_Catholics.html译者：飞龙
电子邮件：econometrics.methods@gmail.com
最后更新时间：10-4-2020
马克斯·韦伯（1930）认为，新教伦理，尤其是加尔文主义的变体，比天主教更符合资本主义。韦伯（1930）观察到，北欧的新教地区比南欧的天主教地区更发达。他假设新教徒工作更努力，储蓄更多，更依赖自己，对政府的期望更少。所有这些特征都将导致更大的经济繁荣。
也许，宗教并不是经济表现出色的原因。教育是一个混杂因素。从历史上看，新教徒有更高的识字水平，因为他们被激励阅读圣经。
因果效应也可能是相反的。也许，一个工业人更有可能成为新教徒。宗教是一个选择变量。人们自我选择符合他们自己世界观的意识形态。
让我们打开 Basten & Betz（2013）的数据。每一行代表瑞士西部的一个市镇，沃州和弗里堡州。
# Load data from Basten & Betz (2013)
import numpy as np
import pandas as pd
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "finaldata.dta")
df.head(4)
| id | name | district | district_name | foreignpop2000 | prot1980s | cath1980s | noreligion1980s | canton | totalpop2000 | ... | pfl | pfi | pfr | reineink_pc_mean | meanpart | popden2000 | foreignpopshare2000 | sub_hs | super_hs | murten | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2135.0 | Gruyères | 1003.0 | District de la Gruyère | 159 | 0.062251 | 0.907336 | 1.621622 | 10 | 1546.0 | ... | 42.722500 | 52.234196 | 40.143444 | 48.099865 | 40.359428 | 54.455795 | 10.284605 | 85.910339 | 6.770357 | 0.0 |
| 1 | 2128.0 | Châtel-sur-Montsalvens | 1003.0 | District de la Gruyère | 23 | 0.053191 | 0.917526 | 2.061856 | 10 | 205.0 | ... | 49.223751 | 56.793213 | 44.365696 | 42.465569 | 45.434593 | 102.499985 | 11.219512 | 83.832336 | 8.383233 | 0.0 |
| 2 | 2127.0 | Charmey | 1003.0 | District de la Gruyère | 166 | 0.028424 | 0.960818 | 0.255537 | 10 | 1574.0 | ... | 41.087502 | 53.120682 | 39.674942 | 44.451229 | 42.641624 | 20.066286 | 10.546379 | 87.331535 | 6.019766 | 0.0 |
| 3 | 2125.0 | Bulle | 1003.0 | District de la Gruyère | 2863 | 0.053967 | 0.923239 | 1.013825 | 10 | 11149.0 | ... | 47.326248 | 55.033939 | 43.350178 | 50.217991 | 40.885822 | 897.664978 | 25.679434 | 82.203598 | 8.904452 | 0.0 |
4 行×157 列
瑞士在地理和制度方面非常多样化。将生活在阿尔卑斯山的农村天主教徒与苏黎世的城市受过高等教育的新教徒进行比较是不公平的。
从历史上看，城市有不同的激励措施来采纳新教或保持天主教。拥有更强大商人行会的城市更有可能采纳新教；而由贵族统治的城市更有可能保持天主教。
如果我们使用整个国家，混杂因素太多。分析将局限于沃州（历史上是新教徒）和弗里堡（历史上是天主教徒）的州。请参见下面来自 Basten & Betz（2013）的地图。
这片占瑞士 4.5%的 4,883 平方公里地区在制度上和地理上是同质的。1536 年，沃州并不是自愿成为新教徒，而是因为一场战争而被迫。因此，这是一个准实验设置，因为治疗区和对照区在历史事件的影响下相互类似。
来源： Basten & Betz (2013)
在下面的图表中，我们可以看到，市镇中新教徒的比例越高，休闲偏好就越低。蓝色点是历史上的天主教市镇（弗里堡），而红色点是历史上的新教市镇（沃州）。看起来像是不同的子群。我们如何找出是否有因果效应的证据，还是仅仅是相关性？
# Create variable "Region" for the graphic
def category(var):
if var == 1:
return "Protestant (Vaud)"
else:
return "Catholic (Fribourg)"
df['Region'] = df["vaud"].apply(category)
# Rename variables with auto-explanatory names
df = df.rename(columns={"prot1980s": "Share_of_Protestants",
"pfl": "Preference_for_Leisure"})
# Scatter plot
import plotly.express as px
leisure = px.scatter(df,
x="Share_of_Protestants",
y="Preference_for_Leisure",
color="Region")
leisure.show()
让我们深入分析。记住，回归不连续是最接近实验的技术。
在下面的图表中，偏好休闲在边界距离=0 处存在不连续性。边界距离大于 0 的地区包括历史上的新教市镇（沃州）；而边界距离小于 0 的地区包括历史上的天主教市镇（弗里堡）。
运行变量“边界距离”决定了地区，但并不决定新教徒的比例。然而，新教徒的比例随着距离的增加而增加。
df = df.rename(columns={"borderdis": "Border_Distance_in_Km"})
leisure = px.scatter(df,
x="Border_Distance_in_Km",
y="Preference_for_Leisure",
color="Share_of_Protestants",
title="Discontinuity at Distance = 0")
leisure.show()
由于边界是任意的，即由历史事件决定的，靠近边界的市镇很可能彼此更相似，而远离边界的市镇。因此，让我们将分析限制在 5 公里范围内的市镇。
df5 = df[df['Border_Distance_in_Km'] >= -5]
df5 = df5[df5['Border_Distance_in_Km'] <= 5]
简单的均值比较显示，新教市镇的休闲偏好较低（39.5%）与天主教市镇（48.2%）相比。差异为-8.7%。
请注意，新教地区的瑞士法郎平均收入水平较高（47.2K vs 43.7K），而基尼系数捕捉到的不平等程度也较高（0.36 vs 0.30）。
df5 = df5.rename(columns={"reineink_pc_mean": "Mean_Income_(CHF)",
"Ecoplan_gini" : "Gini_1996"})
outcome = ['Preference_for_Leisure', 'Mean_Income_(CHF)', 'Gini_1996']
考虑到只使用了 5 公里范围内的市镇（49 个天主教市镇和 84 个新教市镇），上述比较并不“糟糕”。此外，两个地区在 1980 年的无宗教信仰比例（1.7% vs 2.9%）和海拔高度（642 vs 639 米）方面相似。
然而，一个更可信的方法是使用一个带有运行变量(\(r_i\))的回归不连续框架。
control = ['noreligion1980s', 'altitude']
df5.loc[:, control].groupby(df5['Region']).agg([np.mean, np.std]).T
| 地区 | 天主教（弗里堡） | 新教（沃州） | |
|---|---|---|---|
| 无宗教 1980 年代 | 平均值 | 1.729423 | 2.949534 |
| std | 1.499346 | 2.726086 | |
| 海拔 | 平均值 | 642.591858 | 639.607117 |
| std | 120.230320 | 113.563847 |
在图表“模糊回归不连续”中，我们可以清楚地看到运行变量“边界距离”与处理变量“新教徒比例”之间存在很强的相关性。变量“边界距离”并不决定处理状态，但增加了成为新教徒的概率。因此，这是一个模糊回归不连续的情况，而不是一个尖锐的回归不连续。
让\(D_i\)成为单位\(i\)的处理状态。\(P(D_i=1|r_i)\)是在截断\(r_0\)处处理概率的跳跃：
其中\(f_1(r_i)\)和\(f_0(r_i)\)是可以假定任何值的函数。在尖锐的回归不连续中，\(f_1(r_i)\)为 1，\(f_0(r_i)\)为 0。
fuzzy = px.scatter(df5,
x="Border_Distance_in_Km",
y="Share_of_Protestants",
color="Region",
title="Fuzzy Regression Discontinuity")
fuzzy.show()
在下面的图表中，新教徒比例变量被模拟，以说明尖锐的回归不连续的情况。
def dummy(var):
if var >= 0:
return 1
else:
return 0
df5["Simulated_Share_Protestant"] = df5["Border_Distance_in_Km"].apply(dummy)
sharp = px.scatter(df5,
x="Border_Distance_in_Km",
y="Simulated_Share_Protestant",
color="Region",
title="Sharp Regression Discontinuity")
sharp.show()
假设\(Y\) = 休闲偏好，\(D_r\) = 新教徒比例，\(r\) = 边界距离。现在，我们有一个估计以下方程的问题：
兴趣变量\(D_r\)不再被\(r\)“净化”，也就是说，它不再完全由运行变量\(r\)决定。因此，\(D_r\)很可能与误差项\(\epsilon\)相关：
我们可以使用一个与误差项\(\epsilon\)不相关且在控制其他因素后与\(D_r\)相关的工具变量\(Z\)来解决这个问题：
\(Z\)的自然候选人是变量“vaud”：如果是瓦州的市政当局，则为 1；如果是弗里堡的市政当局，则为 0。没有理由相信这个变量与误差项\(\epsilon\)相关，因为分隔该地区的边界是在 1536 年确定的，当时伯尔尼共和国征服了瓦州。第二个假设也是有效的，因为瓦州地区的新教徒比弗里堡地区更多。
工具变量方法首先是使用\(Z\)“纯化”\(D_r\)：
然后，我们通过进行普通最小二乘（OLS）回归得到\(\hat{D}_r\)的拟合值，并将其代入方程：
逻辑是“纯化”的\(\hat{D}_r\)与误差项\(\epsilon\)不相关。现在，我们可以运行普通最小二乘（OLS）来得到\(\hat{D}_r\)对\(Y\)的孤立影响，即\(\rho\)将是无偏的因果效应。
# Install libray to run Instrumental Variable estimation
!pip install linearmodels
Requirement already satisfied: linearmodels in c:\anaconda\envs\textbook\lib\site-packages (4.17)
Requirement already satisfied: scipy>=1 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.5.3)
Requirement already satisfied: pandas>=0.23 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.1.3)
Requirement already satisfied: property-cached>=1.6.3 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.6.4)
Requirement already satisfied: patsy in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (0.5.1)
Requirement already satisfied: Cython>=0.29.14 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (0.29.21)
Requirement already satisfied: statsmodels>=0.9 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (0.10.2)
Requirement already satisfied: numpy>=1.15 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.19.2)
WARNING: No metadata found in c:\anaconda\envs\textbook\lib\site-packages
ERROR: Could not install packages due to an EnvironmentError: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
计算机自动运行工具变量（IV）程序的两个阶段。我们指出内生变量\(D_r\)是“新教徒比例”，工具变量\(Z\)是“vaud”。我们还添加了变量“t_dist”，即变量“vaud”和“边境距离”的交互。
结果是，“新教徒比例”使休闲偏好减少了 13.4%。
from linearmodels.iv import IV2SLS
iv = 'Preference_for_Leisure ~ 1 + Border_Distance_in_Km + t_dist + [Share_of_Protestants ~ vaud]'
iv_result = IV2SLS.from_formula(iv, df5).fit(cov_type='robust')
print(iv_result)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
IV-2SLS Estimation Summary
==================================================================================
Dep. Variable: Preference_for_Leisure R-squared: 0.4706
Estimator: IV-2SLS Adj. R-squared: 0.4583
No. Observations: 133 F-statistic: 99.425
Date: Fri, Nov 06 2020 P-value (F-stat) 0.0000
Time: 21:12:05 Distribution: chi2(3)
Cov. Estimator: robust
Parameter Estimates
=========================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
-----------------------------------------------------------------------------------------
Intercept 50.528 1.8885 26.755 0.0000 46.826 54.229
Border_Distance_in_Km 0.4380 0.6291 0.6963 0.4862 -0.7949 1.6710
t_dist -0.3636 0.7989 -0.4552 0.6490 -1.9294 1.2022
Share_of_Protestants -13.460 3.1132 -4.3235 0.0000 -19.562 -7.3582
=========================================================================================
Endogenous: Share_of_Protestants
Instruments: vaud
Robust Covariance (Heteroskedastic)
Debiased: False
C:\Anaconda\envs\textbook\lib\site-packages\linearmodels\iv\data.py:25: FutureWarning:
is_categorical is deprecated and will be removed in a future version. Use is_categorical_dtype instead
我们还可以检查第一阶段，看看工具变量“vaud”是否与“新教徒比例”相关联，控制其他因素如“边境距离”和“t_dist”后。瓦州使新教徒比例增加了 67%。“vaud”的 t 值为 20，即在没有任何疑问的情况下具有统计学意义。
因此，我们确信第二阶段的结果比简单的均值比较更可信。13.4%的工具变量影响比 8.7%的简单均值比较更可信。
print(iv_result.first_stage)
First Stage Estimation Results
===============================================
Share_of_Protestants
-----------------------------------------------
R-squared 0.9338
Partial R-squared 0.7095
Shea's R-squared 0.7095
Partial F-statistic 403.04
P-value (Partial F-stat) 0.0000
Partial F-stat Distn chi2(1)
========================== ===========
Intercept 0.1336
(7.7957)
Border_Distance_in_Km 0.0169
(2.8397)
t_dist -0.0057
(-0.4691)
vaud 0.6709
(20.076)
-----------------------------------------------
T-stats reported in parentheses
T-stats use same covariance type as original model
8.7%的简单均值比较结果更接近于下面的天真锐度回归（SRD）的结果为 9%。瓦州地区对休闲的偏好比弗里堡少 9%。我们不能得出新教徒对休闲的偏好比天主教徒少 9%的结论。瓦州地区并不是 100%的新教徒。弗里堡地区也不是 100%的天主教徒。
Fuzz 回归离散度（FRD）使用工具变量（IV）估计，是对天真比较的修正。FRD 隔离了新教徒对休闲偏好的影响。因此，最可信的估计是，新教徒对休闲的偏好比天主教徒少 13.4%。
naive_srd = 'Preference_for_Leisure ~ 1 + vaud + Border_Distance_in_Km + t_dist'
srd = IV2SLS.from_formula(naive_srd, df5).fit(cov_type='robust')
print(srd)
OLS Estimation Summary
==================================================================================
Dep. Variable: Preference_for_Leisure R-squared: 0.3830
Estimator: OLS Adj. R-squared: 0.3686
No. Observations: 133 F-statistic: 91.767
Date: Fri, Nov 06 2020 P-value (F-stat) 0.0000
Time: 21:12:05 Distribution: chi2(3)
Cov. Estimator: robust
Parameter Estimates
=========================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
-----------------------------------------------------------------------------------------
Intercept 48.730 1.5731 30.976 0.0000 45.646 51.813
vaud -9.0303 2.1821 -4.1383 0.0000 -13.307 -4.7534
Border_Distance_in_Km 0.2107 0.6028 0.3496 0.7266 -0.9707 1.3922
t_dist -0.2874 0.8479 -0.3390 0.7346 -1.9493 1.3744
=========================================================================================
练习
1）在估计瑞士（整个国家）的新教徒对天主教徒的经济繁荣的因果影响时，会有哪些混杂因素？详细解释每个混杂因素。
2）有人可能会认为在 1536 年伯尔尼共和国征服瓦州之前，西部瑞士是一个多样化的地区。在这种情况下，弗里堡（天主教徒）对瓦州（新教徒）来说不是一个合理的对照组。在 1536 年之前，可以使用哪些变量来测试该地区的同质性/多样性？指出数据是否存在，或者是否可以收集数据。
3）我复制了 Basten 和 Betz（2013）的主要结果，并注意到他们没有控制教育因素。Becker & Woessmann（2009）认为普鲁士在 19 世纪的新教徒对经济繁荣的影响是由于较高的识字率。Basten 和 Betz（2013）在在线附录中呈现了下表。表中的数字是加强还是削弱了 Basten 和 Betz（2013）的结果？请解释你的推理。
来源： Basten & Betz（2013）
4）休闲偏好是一项调查中的自我报告变量。也许，新教徒比天主教徒更有动机宣称对休闲的偏好较低。在现实世界中，新教徒可能和天主教徒一样喜欢休闲。声明的偏好可能与真实行为不符。相关的研究问题是宗教（新教）是否导致人们实际上努力工作。要有创意地提出一种方法（方法和数据）来解决上述问题。
5）使用巴斯滕和贝茨（2013）的数据，调查新教徒是否导致瑞士西部的平均收入更高。采用您认为最可信的规范来恢复无偏因果效应。解释和证明推理的每一步。解释主要结果。
参考文献
巴斯滕，克里斯托夫和弗兰克·贝茨（2013）。超越职业道德：宗教，个人和政治偏好。《美国经济学杂志：经济政策》5（3）：67-91。在线附录
贝克尔，萨沙 O.和鲁德格·沃斯曼。 （2009）。韦伯错了吗？新教经济史的人力资本理论。《季度经济学杂志》124（2）：531–96。
韦伯，马克思。 （1930）。新教伦理与资本主义精神。纽约：斯克里布纳，（原始出版于 1905 年）。
五、联邦储备是否能够阻止大萧条？
原文：
causal-methods.github.io/Book/5%29_Could_the_Federal_Reserve_Prevent_the_Great_Depression.html译者：飞龙
电子邮件：econometrics.methods@gmail.com
最后更新日期：2020 年 10 月 12 日
关于大萧条，新古典经济学家认为经济活动的下降“导致”了银行倒闭。凯恩斯（1936）则相反地认为，银行破产导致了企业倒闭。
Richardson & Troost（2009）注意到，在大萧条期间，密西西比州被联邦储备（Fed）的不同分支控制，分为圣路易斯和亚特兰大两个地区。
与圣路易斯 Fed 不同，亚特兰大 Fed 采取了凯恩斯主义的贴现贷款和向流动性不足的银行提供紧急流动性的政策，使得借款更加困难。
让我们打开 Ziebarth（2013）的数据。每一行都是 1929 年、1931 年、1933 年和 1935 年的人口普查制造业（CoM）中的一家公司。
# Load data from Ziebarth (2013)
import numpy as np
import pandas as pd
path = "https://github.com/causal-methods/Data/raw/master/"
data = pd.read_stata(path + "MS_data_all_years_regs.dta")
data.head()
| 县 | 平均工资收入人数 | 平均工资 a | 平均工资 b | 平均工资 d | 开始 | 电力容量编号 | 汽油容量编号 | 马车容量编号 | 人口普查年份 | ... | 1933 年未进入 | 1931 年未退出 | 1935 年未进入 | 1933 年未退出 | 1931 年平衡 | 1933 年平衡 | 1935 年平衡 | 产品数量 | rrtsap | delta_indic | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | NaN | NaN | NaN | NaN | 1933 | ... | NaN | NaN | NaN | NaN | 0.0 | 0.0 | 0.0 | 0.0 | NaN | 0.0 | |||||
| 1 | Greene | 12.000000 | 0.706404 | 0.254306 | NaN | 1933 年 1 月 1 日 | 1933 | ... | NaN | NaN | NaN | NaN | 0.0 | 0.0 | 0.0 | 2.0 | 100.797699 | 0.0 | |||
| 2 | Hinds | 4.000000 | 0.242670 | 0.215411 | NaN | 1933 年 4 月 1 日 | 1933 | ... | NaN | NaN | NaN | NaN | 0.0 | 0.0 | 0.0 | 1.0 | 404.526093 | 0.0 | |||
| 3 | Wayne | 36.000000 | 0.064300 | 0.099206 | NaN | 1933 年 1 月 1 日 | 1933 | ... | NaN | NaN | NaN | NaN | 0.0 | 0.0 | 0.0 | 3.0 | 104.497620 | 0.0 | |||
| 4 | Attala | 6.333333 | NaN | 0.437281 | NaN | 1933 年 1 月 1 日 | 1933 | ... | NaN | NaN | NaN | NaN | 0.0 | 0.0 | 0.0 | 1.0 | 148.165237 | 0.0 |
5 行×538 列
首先，我们必须检查 1929 年圣路易斯和亚特兰大两个地区有多相似。所有变量都以对数形式报告。圣路易斯地区公司的平均收入为 10.88；而亚特兰大地区公司的平均收入为 10.78。圣路易斯和亚特兰大的工资收入（4.54 比 4.69）和每个工人的工作时间（4.07 比 4）也相似。
# Round 2 decimals
pd.set_option('precision', 2)
# Restrict the sample to the year: 1929
df1929 = data[data.censusyear.isin([1929])]
vars= ['log_total_output_value', 'log_wage_earners_total',
'log_hours_per_wage_earner']
df1929.loc[:, vars].groupby(df1929["st_louis_fed"]).agg([np.size, np.mean])
| log_total_output_value | log_wage_earners_total | log_hours_per_wage_earner | |
|---|---|---|---|
| 大小 | 平均 | 大小 | |
| --- | --- | --- | --- |
| 圣路易斯 Fed | |||
| --- | --- | --- | --- |
| 0.0 | 424.0 | 10.78 | 424.0 |
| 1.0 | 367.0 | 10.88 | 367.0 |
此外，圣路易斯和亚特兰大的平均价格（1.72 比 1.55）和平均数量（8.63 比 8.83）也相似，如果样本限制在只有 1 种产品的公司。因此，亚特兰大地区是圣路易斯地区的一个合理对照组。
# Restrict sample to firms with 1 product
df1929_1 = df1929[df1929.num_products.isin([1])]
per_unit = ['log_output_price_1', 'log_output_quantity_1']
df1929_1.loc[:, per_unit].groupby(df1929_1["st_louis_fed"]).agg([np.size, np.mean])
| log_output_price_1 | log_output_quantity_1 | |
|---|---|---|
| 大小 | 平均 | |
| --- | --- | --- |
| 圣路易斯 Fed | ||
| --- | --- | --- |
| 0.0 | 221.0 | 1.55 |
| 1.0 | 225.0 | 1.72 |
我们想要看看圣路易斯 Fed 的信贷约束政策是否减少了公司的收入。换句话说，亚特兰大 Fed 是否拯救了濒临破产的公司。
为此，我们必须探索时间维度：比较 1929 年之前和 1931 年之后的公司收入。
让我们将样本限制在 1929 年和 1931 年。然后，让我们删除缺失值。
# Restrict the sample to the years: 1929 and 1931
df = data[data.censusyear.isin([1929, 1931])]
vars = ['firmid', 'censusyear', 'log_total_output_value',
'st_louis_fed', 'industrycode', 'year_1931']
# Drop missing values
df = df.dropna(subset=vars)
现在，我们可以声明一个面板数据结构，即设置分析单位和时间维度。注意表中的变量“firmid”和“censusyear”成为了索引。顺序很重要。第一个变量必须是分析单位，第二个变量必须是时间单位。在表中看到，公司（id = 12）在 1929 年和 1931 年观察到。
注意，在清理数据集后声明了面板数据结构。例如，如果在声明面板数据后删除了缺失值，进行回归的命令可能会返回错误。
df = df.set_index(['firmid', 'censusyear'])
df.head()
| 县 | 平均就业人数 | 平均工资 a | 平均工资 b | 平均工资 d | 开始 | 电力容量编号 | 汽油容量编号 | 马车容量编号 | 城镇或村庄 | ... | 1933 年不进入 | 1931 年不退出 | 1935 年不进入 | 1933 年不退出 | 1931 年平衡 | 1933 年平衡 | 1935 年平衡 | 产品数量 | rrtsap | delta_indic | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| firmid | censusyear | |||||||||||||||||||||
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 12 | 1929 | Oktibbeha | 1.75 | 0.38 | 0.50 | NaN | January 1, 1929 | A & M College | ... | NaN | NaN | NaN | NaN | 1.0 | 0.0 | 1.0 | 1.0 | 328.31 | 0.0 | |||
| 1931 | Oktibbeha | 1.75 | 0.25 | 0.30 | NaN | January 1, 1931 | A and M College | ... | NaN | 1.0 | NaN | NaN | 1.0 | 0.0 | 0.0 | 1.0 | NaN | 0.0 | ||||
| 13 | 1929 | Warren | 7.25 | 0.49 | 0.58 | NaN | January 1, 1929 | Clarksburg | ... | NaN | NaN | NaN | NaN | 1.0 | 0.0 | 1.0 | 2.0 | 670.07 | 1.0 | |||
| 1931 | Warren | 3.50 | NaN | 0.71 | NaN | January 1, 1931 | Vicksburg, | ... | NaN | 1.0 | NaN | NaN | 1.0 | 0.0 | 0.0 | 2.0 | NaN | 1.0 | ||||
| 14 | 1929 | Monroe | 12.92 | 0.17 | 0.23 | NaN | January 1, 1929 | Aberdeen | ... | NaN | NaN | NaN | NaN | 1.0 | 0.0 | 1.0 | 2.0 | 314.90 | 0.0 |
5 行×536 列
让我们解释面板数据相对于横截面数据的优势。后者是某一时间点或时期的快照；而在面板数据中，同一分析单位随时间观察。
让\(Y_{it}\)表示单位\(i\)在时间\(t\)的结果变量。虚拟变量\(d2_{it}\)在第二个时期为 1，在第一个时期为 0。注意解释变量\(X_{it}\)在单位\(i\)和时间\(t\)上变化，但未观察到的因素\(\alpha_i\)在时间上不变。未观察到的因素是一个不可用的变量（数据），可能与感兴趣的变量相关，导致结果偏差。
探索时间变化的优势在于，未观察到的因素\(\alpha_i\)可以通过一阶差分（FD）方法消除。
在第二个时期（\(t=2\)），时间虚拟\(d2=1\)：
在第一个时期（\(t=1\)），时间虚拟\(d2=0\)：
然后：
因此，如果观察到相同的单位随时间变化（面板数据），就不需要担心一个可以被认为在分析的时间内保持不变的因素。我们可以假设公司文化和制度实践在短时间内变化不大。这些因素可能解释了公司之间收入的差异，但如果上述假设是正确的，它们不会对结果产生偏见。
让我们安装可以运行面板数据回归的库。
!pip install linearmodels
Requirement already satisfied: linearmodels in c:\anaconda\envs\textbook\lib\site-packages (4.17)
Requirement already satisfied: numpy>=1.15 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.19.2)
WARNING: No metadata found in c:\anaconda\envs\textbook\lib\site-packages
ERROR: Could not install packages due to an EnvironmentError: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
让我们使用双重差分（DID）方法来估计圣路易斯联邦储备银行政策对公司收入的影响。除了探索时间差异外，还必须使用治疗-对照差异来估计政策的因果影响。
让\(Y\)成为结果变量'log_total_output_value'，\(d2\)成为时间虚拟变量'year_1931'，\(dT\)成为治疗虚拟变量'st_louis_fed'，\(d2 \cdot dT\)成为前两个虚拟变量之间的交互项：
DID 估计量由\(\delta_1\)给出，而不是由\(\beta_1\)或\(\delta_0\)给出。首先，我们计算“治疗组（圣路易斯）”和“对照组（亚特兰大）”之间的差异，然后计算“之后（1931 年）”和“之前（1921 年）”之间的差异。
顺序无关紧要。如果我们首先计算“之后（1931 年）”和“之前（1921 年）”之间的差异，然后计算“治疗组（圣路易斯）”和“对照组（亚特兰大）”之间的差异，结果将是相同的\(\delta_1\)。
让我们正式地展示，我们必须在 DID 估计量\(\delta_0\)中两次取差异：
如果\(d2=0\)且\(dT=0\)，那么\(Y_{0,0}=\beta_0\)。
如果\(d2=1\)且\(dT=0\)，那么\(Y_{1,0}=\beta_0+\delta_0\)。
对于对照组，“之后-之前”的差异是：
让我们对治疗组应用相同的推理：
如果\(d2=0\)且\(dT=1\)，那么\(Y_{0,1}=\beta_0 + \beta_1\)。
如果\(d2=1\)且\(dT=1\)，那么\(Y_{1,1}=\beta_0+\delta_0+ \beta_1+\delta_1\)。
对于治疗组，“之后-之前”的差异是：
然后，如果我们计算“治疗组-对照组”的差异，我们得到：
让我们手动计算“大萧条”期间公司收入图表中的\(\hat{\delta}_1\)。请注意，在双重差分（DID）方法中，根据对照组（亚特兰大）构建了一个反事实。这只是亚特兰大线的平行移位。反事实是治疗组（圣路易斯）的假设结果，如果圣路易斯联邦储备银行遵循了亚特兰大联邦储备银行相同的政策。
圣路易斯联邦储备银行的严格信贷政策使公司的收入减少了约 20%。1931 年底的简单均值比较结果仅为-10%。因此，如果不使用反事实推理，圣路易斯联邦储备银行政策的负面影响将被大大低估。
# Mean Revenue for the Graphic
table = pd.crosstab(df['year_1931'], df['st_louis_fed'],
values=df['log_total_output_value'], aggfunc='mean')
# Build Graphic
import plotly.graph_objects as go
fig = go.Figure()
# x axis
year = [1929, 1931]
# Atlanta Line
fig.add_trace(go.Scatter(x=year, y=table[0],
name='Atlanta (Control)'))
# St. Louis Line
fig.add_trace(go.Scatter(x=year, y=table[1],
name='St. Louis (Treatment)'))
# Counterfactual
end_point = (table[1][0] - table[0][0]) + table[0][1]
counter = [table[1][0], end_point]
fig.add_trace(go.Scatter(x=year, y= counter,
name='Counterfactual',
line=dict(dash='dot') ))
# Difference-in-Differences (DID) estimation
fig.add_trace(go.Scatter(x=[1931, 1931],
y=[table[1][1], end_point],
name='$\delta_1=0.2$',
line=dict(dash='dashdot') ))
# Labels
fig.update_layout(title="Firm's Revenue during the Great Depression",
xaxis_type='category',
xaxis_title='Year',
yaxis_title='Log(Revenue)')
fig.show()
通过回归实施的双重差分（DID）的结果是：
from linearmodels import PanelOLS
Y = df['log_total_output_value']
df['const'] = 1
df['louis_1931'] = df['st_louis_fed']*df['year_1931']
## Difference-in-Differences (DID) specification
dd = ['const', 'st_louis_fed', 'year_1931', 'louis_1931']
dif_in_dif = PanelOLS(Y, df[dd]).fit(cov_type='clustered',
cluster_entity=True)
print(dif_in_dif)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
PanelOLS Estimation Summary
====================================================================================
Dep. Variable: log_total_output_value R-squared: 0.0257
Estimator: PanelOLS R-squared (Between): 0.0145
No. Observations: 1227 R-squared (Within): 0.2381
Date: Fri, Nov 06 2020 R-squared (Overall): 0.0257
Time: 21:12:12 Log-likelihood -2135.5
Cov. Estimator: Clustered
F-statistic: 10.761
Entities: 938 P-value 0.0000
Avg Obs: 1.3081 Distribution: F(3,1223)
Min Obs: 1.0000
Max Obs: 2.0000 F-statistic (robust): 18.780
P-value 0.0000
Time periods: 2 Distribution: F(3,1223)
Avg Obs: 613.50
Min Obs: 575.00
Max Obs: 652.00
Parameter Estimates
================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
--------------------------------------------------------------------------------
const 10.781 0.0723 149.13 0.0000 10.639 10.923
st_louis_fed 0.0945 0.1043 0.9057 0.3653 -0.1102 0.2991
year_1931 -0.3521 0.0853 -4.1285 0.0000 -0.5194 -0.1848
louis_1931 -0.1994 0.1237 -1.6112 0.1074 -0.4422 0.0434
================================================================================
C:\Anaconda\envs\textbook\lib\site-packages\linearmodels\panel\data.py:98: FutureWarning:
is_categorical is deprecated and will be removed in a future version. Use is_categorical_dtype instead
圣路易斯联邦储备银行的政策使公司收入减少了 18%（\(1-e^{-0.1994}\)）。然而，p 值为 0.1074。结果在 10%的显著性水平下并不具有统计学意义。
from math import exp
1 - exp(dif_in_dif.params.louis_1931 )
0.18076309464004925
有人可能会认为公司之间的差异是一个混杂因素。一个或另一个大公司可能会使结果产生偏差。
这个问题可以通过使用固定效应（FE）或 Within Estimator 来解决。该技术类似于 First-Difference（FD），但数据转换不同。时间去均值过程用于消除未观察到的因素\(\alpha_i\)。
然后，我们对每个\(i\)在时间\(t\)上的变量进行平均：
然后，我们取差异，未观察到的因素\(\alpha_i\)消失：
我们可以用更简洁的方式写出上面的方程：
正如我们之前宣布的，公司是这个面板数据集中的分析单位，计算机会自动使用命令“entity_effects=True”实现公司固定效应（FE）。
我们在双重差分（DID）规范中添加了公司固定效应（FE），结果并没有改变太多。直觉是双重差分（DID）技术已经减轻了内生性问题。
圣路易斯联邦政策使公司收入减少了 17%（\(1-e^{-0.1862}\)）。结果在 10%的显著水平上是统计上显著的。
firmFE = PanelOLS(Y, df[dd], entity_effects=True)
print(firmFE.fit(cov_type='clustered', cluster_entity=True))
PanelOLS Estimation Summary
====================================================================================
Dep. Variable: log_total_output_value R-squared: 0.2649
Estimator: PanelOLS R-squared (Between): -0.4554
No. Observations: 1227 R-squared (Within): 0.2649
Date: Fri, Nov 06 2020 R-squared (Overall): -0.4266
Time: 21:12:12 Log-likelihood -202.61
Cov. Estimator: Clustered
F-statistic: 34.361
Entities: 938 P-value 0.0000
Avg Obs: 1.3081 Distribution: F(3,286)
Min Obs: 1.0000
Max Obs: 2.0000 F-statistic (robust): 31.245
P-value 0.0000
Time periods: 2 Distribution: F(3,286)
Avg Obs: 613.50
Min Obs: 575.00
Max Obs: 652.00
Parameter Estimates
================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
--------------------------------------------------------------------------------
const 9.9656 0.4511 22.091 0.0000 9.0777 10.854
st_louis_fed 1.9842 1.0365 1.9144 0.0566 -0.0559 4.0243
year_1931 -0.3666 0.0657 -5.5843 0.0000 -0.4959 -0.2374
louis_1931 -0.1862 0.0982 -1.8965 0.0589 -0.3794 0.0071
================================================================================
F-test for Poolability: 6.8220
P-value: 0.0000
Distribution: F(937,286)
Included effects: Entity
固定效应（FE）可以通过添加虚拟变量手动实现。有不同的固定效应。让我们在双重差分（DID）规范中添加行业固定效应，以排除结果可能受行业特定冲击驱动的可能性。
圣路易斯联邦政策使公司收入减少了 14.2%（\(1-e^{-0.1533}\)）。结果在 10%的显著水平上是统计上显著的。
为什么不同时添加公司和行业固定效应？这是可能的，也是可取的，但由于多重共线性问题，计算机不会返回任何结果。我们每家公司只有两个观察值（2 年）。如果我们为每家公司添加一个虚拟变量，就像运行一个变量比观察值更多的回归。
在他的论文中，Ziebarth（2013）使用了公司和行业固定效应，这是如何可能的？Ziebarth（2013）使用了 Stata 软件。在多重共线性问题的情况下，Stata 会自动删除一些变量并输出结果。尽管这种做法在经济学顶级期刊中很常见，但这并不是“真正”的固定效应。
industryFE = PanelOLS(Y, df[dd + ['industrycode']])
print(industryFE.fit(cov_type='clustered', cluster_entity=True))
PanelOLS Estimation Summary
====================================================================================
Dep. Variable: log_total_output_value R-squared: 0.5498
Estimator: PanelOLS R-squared (Between): 0.5462
No. Observations: 1227 R-squared (Within): 0.3913
Date: Fri, Nov 06 2020 R-squared (Overall): 0.5498
Time: 21:12:12 Log-likelihood -1661.9
Cov. Estimator: Clustered
F-statistic: 29.971
Entities: 938 P-value 0.0000
Avg Obs: 1.3081 Distribution: F(48,1178)
Min Obs: 1.0000
Max Obs: 2.0000 F-statistic (robust): 4.791e+15
P-value 0.0000
Time periods: 2 Distribution: F(48,1178)
Avg Obs: 613.50
Min Obs: 575.00
Max Obs: 652.00
Parameter Estimates
=====================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
-------------------------------------------------------------------------------------
const 10.391 0.2249 46.201 0.0000 9.9497 10.832
st_louis_fed -0.1740 0.0805 -2.1610 0.0309 -0.3319 -0.0160
year_1931 -0.3163 0.0660 -4.7922 0.0000 -0.4458 -0.1868
louis_1931 -0.1533 0.0916 -1.6741 0.0944 -0.3331 0.0264
industrycode.1005 -0.1966 0.3637 -0.5407 0.5888 -0.9101 0.5169
industrycode.101 0.2349 0.2403 0.9775 0.3285 -0.2366 0.7064
industrycode.1014 -0.1667 1.0812 -0.1542 0.8775 -2.2880 1.9546
industrycode.103 1.5898 0.2938 5.4104 0.0000 1.0133 2.1663
industrycode.104 1.1349 0.2818 4.0280 0.0001 0.5821 1.6877
industrycode.105 0.9660 0.4191 2.3049 0.0213 0.1437 1.7882
industrycode.107 0.9374 0.3073 3.0502 0.0023 0.3344 1.5404
industrycode.110 0.2299 0.3447 0.6671 0.5049 -0.4463 0.9062
industrycode.111 2.8271 0.4097 6.9010 0.0000 2.0233 3.6308
industrycode.112 0.2676 0.4185 0.6394 0.5227 -0.5535 1.0888
industrycode.114 2.7559 0.4623 5.9612 0.0000 1.8489 3.6630
industrycode.116 0.4498 0.4351 1.0337 0.3015 -0.4039 1.3035
industrycode.117 -0.1337 0.5436 -0.2461 0.8057 -1.2002 0.9327
industrycode.118 0.1025 0.2483 0.4127 0.6799 -0.3847 0.5896
industrycode.119 -0.3121 0.2336 -1.3360 0.1818 -0.7705 0.1462
industrycode.1204 -0.0102 0.3031 -0.0338 0.9731 -0.6048 0.5844
industrycode.123 0.2379 0.7093 0.3354 0.7374 -1.1537 1.6295
industrycode.126 1.8847 0.4828 3.9039 0.0001 0.9375 2.8319
industrycode.128 -0.0025 0.5277 -0.0047 0.9963 -1.0377 1.0328
industrycode.1303 2.6410 0.2293 11.519 0.0000 2.1912 3.0908
industrycode.136 0.2537 0.2493 1.0178 0.3090 -0.2354 0.7428
industrycode.1410 -1.4091 0.2967 -4.7496 0.0000 -1.9912 -0.8270
industrycode.1502 1.4523 0.3858 3.7644 0.0002 0.6954 2.2093
industrycode.1604 -0.5652 0.4195 -1.3473 0.1781 -1.3882 0.2579
industrycode.1624 -1.0476 1.0108 -1.0365 0.3002 -3.0307 0.9355
industrycode.1640 -0.8320 0.3613 -2.3030 0.0215 -1.5409 -0.1232
industrycode.1676 0.5830 0.2512 2.3205 0.0205 0.0901 1.0759
industrycode.1677 -0.7025 0.2346 -2.9947 0.0028 -1.1627 -0.2422
industrycode.203 -0.6435 0.2241 -2.8712 0.0042 -1.0832 -0.2038
industrycode.210B 1.7535 0.2293 7.6482 0.0000 1.3037 2.2033
industrycode.216 2.7641 0.3549 7.7888 0.0000 2.0679 3.4604
industrycode.234a 1.6871 0.5410 3.1187 0.0019 0.6258 2.7485
industrycode.265 1.2065 0.4193 2.8774 0.0041 0.3838 2.0292
industrycode.266 -0.3244 0.2736 -1.1857 0.2360 -0.8612 0.2124
industrycode.304 1.2156 0.4982 2.4400 0.0148 0.2382 2.1930
industrycode.314 0.8008 0.2411 3.3211 0.0009 0.3277 1.2740
industrycode.317 -0.1470 0.2750 -0.5347 0.5930 -0.6866 0.3925
industrycode.515 -0.4623 0.2862 -1.6153 0.1065 -1.0238 0.0992
industrycode.518 -0.5243 0.2887 -1.8157 0.0697 -1.0908 0.0422
industrycode.520 -0.2234 0.2731 -0.8179 0.4136 -0.7592 0.3124
industrycode.614 1.8029 0.3945 4.5707 0.0000 1.0290 2.5768
industrycode.622 2.9239 0.2430 12.035 0.0000 2.4473 3.4006
industrycode.633 0.0122 1.5690 0.0078 0.9938 -3.0661 3.0906
industrycode.651 1.9170 0.8690 2.2060 0.0276 0.2121 3.6220
industrycode.652 2.5495 0.2211 11.531 0.0000 2.1158 2.9833
=====================================================================================
仅仅出于练习目的，假设未观察到的因素\(\alpha_i\)被忽略。这个假设被称为随机效应（RE）。在这种情况下，\(\alpha_i\)将包含在误差项\(v_{it}\)中，并可能使结果产生偏见。
在一个实验中，处理变量与未观察到的因素\(\alpha_i\)不相关。在这种情况下，随机效应（RE）模型具有产生比固定效应模型更低标准误差的优势。
请注意，如果我们运行简单的随机效应（RE）回归，我们可能错误地得出结论，即圣路易斯联邦政策使公司收入增加了 7%。
from linearmodels import RandomEffects
re = RandomEffects(Y, df[['const', 'st_louis_fed']])
print(re.fit(cov_type='clustered', cluster_entity=True))
RandomEffects Estimation Summary
====================================================================================
Dep. Variable: log_total_output_value R-squared: 0.3689
Estimator: RandomEffects R-squared (Between): -0.0001
No. Observations: 1227 R-squared (Within): 0.0025
Date: Fri, Nov 06 2020 R-squared (Overall): -0.0027
Time: 21:12:13 Log-likelihood -1260.9
Cov. Estimator: Clustered
F-statistic: 716.04
Entities: 938 P-value 0.0000
Avg Obs: 1.3081 Distribution: F(1,1225)
Min Obs: 1.0000
Max Obs: 2.0000 F-statistic (robust): 0.5811
P-value 0.4460
Time periods: 2 Distribution: F(1,1225)
Avg Obs: 613.50
Min Obs: 575.00
Max Obs: 652.00
Parameter Estimates
================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
--------------------------------------------------------------------------------
const 10.518 0.0603 174.36 0.0000 10.399 10.636
st_louis_fed 0.0721 0.0946 0.7623 0.4460 -0.1135 0.2577
================================================================================
练习
假设在非实验设置中，控制组与处理组不同。请解释是否合理使用双重差分（DID）来估计因果效应？在 DID 框架中，你应该修改或添加什么？
假设一项研究声称基于双重差分（DID）方法，联邦储备通过 2008 年的银行纾困避免了大规模的企业倒闭。假设另一项基于回归不连续（RD）的研究声称相反或否认了联邦储备对企业倒闭的影响。你认为哪种是更可信的经验策略，DID 还是 RD 来估计联邦政策的因果影响？请解释你的答案。
在面板数据中，分析单位可以是公司或县，公司级别的结果更可信还是县级别的结果更可信？请解释。
使用 Ziebarth（2013）的数据来估计圣路易斯联邦政策对公司收入的影响。具体来说，运行带有随机效应（RE）的双重差分（DID）。解释结果。关于未观察到的因素\(\alpha_i\)可以推断出什么？
使用 Ziebarth（2013）的数据来估计圣路易斯联邦政策对公司收入的影响。具体来说，运行带有公司固定效应（FE）的双重差分（DID），而不使用命令“entity_effects=True”。提示：你必须为每家公司使用虚拟变量。
参考
凯恩斯，约翰梅纳德。 （1936）。《就业、利息和货币的一般理论》。由美国 Harcourt，Brace 和 Company 出版，并由美国 Polygraphic Company 在纽约印刷。
理查森，加里和威廉特鲁斯特（2009 年）。《大萧条期间货币干预减轻了银行恐慌：来自联邦储备区域边界的准实验证据，1929-1933 年》。《政治经济学杂志》117（6）：1031-73。
齐巴斯，尼古拉斯 L.（2013 年）。《通过大萧条期间密西西比州的自然实验来识别银行倒闭的影响》。《美国经济杂志：宏观经济学》5（1）：81-101。
六、2007-2009 年大衰退期间加拿大就业市场上白人女性名字的溢价
译者：飞龙
电子邮件：econometrics.methods@gmail.com
最近更新：2020 年 9 月 15 日
我重新分析了 Oreopoulos（2011）的实验数据，并发现在 2007-2009 年大萧条期间，在加拿大的就业市场中，拥有白人女性名字是有优势的。白人女性在 2009 年 2 月至 9 月之间的回电率比白人男性高出 8%。考虑到白人男性在不同的回归规范下的回电率约为 10%，这一效应的幅度是相当高的。
Oreopoulos（2011）发现，英文名字的回电率为 15.7%，而印度、巴基斯坦、中国和希腊名字的回电率为 6%。我认为他的主要发现在很大程度上是由白人女性驱动的。我发现，在 2009 年 2 月至 9 月的大萧条最严重时期，拥有白人男性名字与印度、中国和希腊男性名字相比，并没有太多优势。
我使用了 Oreopoulos（2011）的数据集。每一行是发送给多个多伦多和蒙特利尔地区职业的简历。
import numpy as np
import pandas as pd
pd.set_option('precision', 3)
# Data from Oreopoulos (2011)
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "oreopoulos.dta")
df.head(5)
| 公司 ID | 职业类型 | name_ethnicity | 额外证书 | 名字 | 语言技能 | 资格认证 | 参考 | 法律 | 列出的资格认证 | ... | 说话技能 | 社交能力 | 写作能力 | 秋季数据 | 中国人 | 印度人 | 英国人 | 巴基斯坦人 | 中加人 | 相同经验 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | -3 | 行政 | 加拿大 | 0.0 | JillWilson | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 70.0 | 50.0 | 67.0 | 2.0 | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | -3 | 行政 | 印度 | 0.0 | PanavSingh | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 70.0 | 50.0 | 67.0 | 2.0 | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | -3 | 行政 | 印度 | 0.0 | RahulKaur | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | ... | 70.0 | 50.0 | 67.0 | 2.0 | NaN | NaN | NaN | NaN | NaN | NaN |
| 3 | -3 | 行政 | 中国 | 0.0 | 雷丽 | 0.0 | 1.0 | 1.0 | 0.0 | 1.0 | ... | 70.0 | 50.0 | 67.0 | 2.0 | NaN | NaN | NaN | NaN | NaN | NaN |
| 4 | -4 | 行政 | 印度 | 0.0 | MayaKumar | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 80.0 | 70.0 | 65.0 | 2.0 | NaN | NaN | NaN | NaN | NaN | NaN |
5 行×31 列
# Transform the variable of interest in %
df["callback"] = 100*df["callback"]
Oreopoulos（2011）收集的第一批实验数据是在 2008 年 4 月至 8 月之间。这是大萧条的“起始期”。
# Restrict data to April and August 2008
df0 = df[(df.fall_data == 0)]
回电率对工作面试的比例对于加拿大名字（15.84%）要比中国和印度名字（9.23%和 9.53%）高得多。名族是随机的。不可能争辩说加拿大人有更多的教育或经验来证明大约 5%的差异。所有的简历在质量上都是一样的，除了申请人的名字。因此，我们可以得出结论，对移民的歧视是一个真实的现象。
mean = df0.groupby('name_ethnicity').agg([np.mean, np.size])
mean["callback"]
| 平均值 | 大小 | |
|---|---|---|
| name_ethnicity | ||
| --- | --- | --- |
| 加拿大 | 15.845 | 953.0 |
| 中国 | 9.231 | 1430.0 |
| 印度 | 9.534 | 1416.0 |
在 2008 年 4 月至 8 月的样本中，女性名字似乎比男性名字稍微有一些优势，但可以忽略不计，以获得工作面试。这一结果支持了 Oreopoulos（2011）的发现。在他的论文中，女性的系数在大部分回归中都不具有统计学显著性。
prop = pd.crosstab(index= df0['name_ethnicity'], columns=df0['female'],
values=df0['callback'], aggfunc='mean')
prop
| 女性 | 0.0 | 1.0 |
|---|---|---|
| name_ethnicity | ||
| --- | --- | --- |
| 加拿大 | 15.551 | 16.122 |
| 中国 | 9.065 | 9.392 |
| 印度 | 9.490 | 9.577 |
Oreopoulos（2011）收集的第三波实验数据是在 2009 年 2 月至 9 月之间进行的。这是大萧条的最糟糕时期。
# Restrict data to February and September 2009
df2 = df[(df.fall_data == 2)]
加拿大名字的回电率为 14％，而中国名字为 8.96％，英文名和中国姓氏为 7.13％，希腊为 10.11％，印度为 7.9％。
请注意，总体上，这第三波样本中的回电率略低于第一波样本，对于两个样本中的常见种族来说。
mean = df2.groupby('name_ethnicity').agg([np.mean, np.size])
mean["callback"]
| 平均 | 大小 | |
|---|---|---|
| 名字种族 | ||
| --- | --- | --- |
| 加拿大 | 14.080 | 1044.0 |
| 中国 | 8.956 | 1418.0 |
| 中国-加拿大 | 7.128 | 491.0 |
| 希腊 | 10.109 | 366.0 |
| 印度 | 7.899 | 1937.0 |
import plotly.express as px
y = mean["callback"].values[:, 0]
x = mean["callback"].index
fig = px.bar(df2, x, y, color = x,
title="Callback Rate for Interview by Name Ethnicity",
labels={ "y": "Callback Rate (%)",
"x": "Name Ethnicity",
"color": ""} )
fig.update_layout(font_size = 17)
fig.show()
在 2009 年 2 月至 9 月的样本中，白人女性的名字回电率为 18.3％，而白人男性的名字回电率为 10.17％。我们没有看到其他种族有这么大的差异。事实上，对于希腊名字来说，效果是相反的。希腊男性的名字回电率为 10.71％，而希腊女性的名字回电率为 9.6％。白人男性的名字比中国和印度男性的名字有优势，但幅度不像白人男性与白人女性之间的差异那么大。
prop = pd.crosstab(index= df2['name_ethnicity'], columns=df2['female'],
values=df2['callback'], aggfunc='mean')
prop
| 女性 | 0.0 | 1.0 |
|---|---|---|
| 名字种族 | ||
| --- | --- | --- |
| 加拿大 | 10.169 | 18.129 |
| 中国 | 8.715 | 9.177 |
| 中国-加拿大 | 6.410 | 7.782 |
| 希腊 | 10.714 | 9.596 |
| 印度 | 7.252 | 8.559 |
import plotly.graph_objects as go
ethnicity = prop.index
male = prop.values[:,0]
female = prop.values[:,1]
fig = go.Figure(data=[
go.Bar(name='Male', x = ethnicity, y = male),
go.Bar(name='Female', x = ethnicity, y = female) ])
fig.update_layout(barmode='group', font_size = 17,
title = "Callback Rate for Interview by Gender",
yaxis = dict(title='Callback Rate (%)'),
xaxis = dict(title='Name Ethnicity') )
fig.show()
有人可能会争辩说，有一些混杂因素导致了观察到的差异。例如，有人可能会说，在现实世界中，女性比男性更受教育和更合格。请记住，这是实验数据，所有简历都是人为构造的，所有相关维度都是由 Oreopoulos（2011）随机化的。控制变量表明，女性和男性彼此相似。
control = ['additional_credential', 'ba_quality',
'extracurricular_skills', 'language_skills',
'certificate', 'ma', 'same_exp', 'exp_highquality',
'skillspeaking', 'skillsocialper', 'skillwriting']
df2.groupby('female').agg([np.mean])[control]
| 附加证书 | 学士质量 | 课外技能 | 语言技能 | 证书 | 硕士 | 相同经验 | 经验高质量 | 说话技能 | 社交能力 | 写作技能 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 女性 | |||||||||||
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0.0 | 0.059 | 0.640 | 0.598 | 0.306 | 0.007 | 0.170 | NaN | 0.197 | 70.757 | 59.654 | 64.286 |
| 1.0 | 0.054 | 0.655 | 0.595 | 0.322 | 0.008 | 0.185 | NaN | 0.169 | 70.524 | 59.777 | 64.216 |
敏锐的读者可能会争辩说，仅仅表明男性和女性相互之间相似是不够的，以支持我关于白人女性溢价的论点。我必须表明样本中的平均白人女性与平均白人男性相似。对于某些变量，白人女性看起来略微更合格，但对于其他变量，略微不太合格。许多数据维度都是随机的，观察到的差异看起来是抽样变异的产物。总的来说，白人男性和女性看起来很相似。我们可以在回归框架中严格控制所有这些因素。我看到种族之间的变化比性别之间的变化更大。种族之间的变化看起来对实验来说过多。因此，我将按种族分解回归分析，并控制几个因素。
df2.groupby(['female', 'name_ethnicity']).agg([np.mean])[control]
| 附加证书 | 学士质量 | 课外技能 | 语言技能 | 证书 | 硕士 | 相同经验 | 经验高质量 | 说话技能 | 社交能力 | 写作技能 | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | 平均 | |||
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 女性 | 名字种族 | |||||||||||
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0.0 | 加拿大 | 0.056 | 0.746 | 0.623 | 0.343 | 0.004 | 0.209 | NaN | 0.190 | 70.422 | 59.070 | 63.546 |
| 中国人 | 0.062 | 0.607 | 0.592 | 0.326 | 0.007 | 0.165 | NaN | 0.183 | 70.702 | 59.850 | 64.312 | |
| 中-加 | 0.068 | 0.628 | 0.538 | 0.282 | 0.013 | 0.154 | NaN | 0.214 | 71.141 | 59.513 | 63.979 | |
| 希腊人 | 0.065 | 0.774 | 0.631 | 0.321 | 0.012 | 0.214 | NaN | 0.232 | 70.976 | 59.696 | 65.220 | |
| 印度人 | 0.055 | 0.586 | 0.596 | 0.276 | 0.006 | 0.147 | NaN | 0.200 | 70.846 | 59.862 | 64.582 | |
| 1.0 | 加拿大 | 0.055 | 0.789 | 0.577 | 0.327 | 0.004 | 0.228 | NaN | 0.177 | 70.146 | 60.041 | 64.179 |
| 中国人 | 0.055 | 0.590 | 0.596 | 0.300 | 0.011 | 0.189 | NaN | 0.179 | 70.799 | 59.607 | 64.349 | |
| 中-加 | 0.047 | 0.638 | 0.638 | 0.346 | 0.004 | 0.113 | NaN | 0.171 | 70.000 | 59.506 | 63.233 | |
| 希腊人 | 0.045 | 0.808 | 0.566 | 0.308 | 0.010 | 0.217 | NaN | 0.136 | 71.258 | 60.662 | 64.894 | |
| 印度人 | 0.055 | 0.608 | 0.599 | 0.334 | 0.008 | 0.172 | NaN | 0.163 | 70.503 | 59.657 | 64.257 |
让\(y_{rjt}\)是一个虚拟变量，如果简历\(r\)发送到工作\(j\)在时间\(t\)收到回电，则为 1；否则为 0。感兴趣的变量是“女性”虚拟变量和与“简历类型”的交互。
有五种“简历类型”：0) 英文名，加拿大教育和经验；1) 外国名字，加拿大教育和经验；2) 外国名字和教育，加拿大经验；3) 外国名字和教育，混合经验；和 4) 外国名字，教育和经验。
以下的线性概率模型是首选的规范：
其中\(X\)是控制变量的向量，\(\epsilon_{rjt}\)是通常的误差项。所有回归都呈现了对异方差性的稳健标准误差。
对于表 1、2 和 3，我们呈现了 4 个回归，以比较“加拿大人”与特定种族。逻辑是保持一个同质样本，避免可能混淆结果的种族变化。
表 1 呈现了没有交互作用和控制变量的结果。作为女性的优势范围从增加 3.64%到 5.97%的回电率，相对于男性。白人男性的回电率，基准（类型 0），范围从 11.14%到 12.29%。Oreopoulos（2011）提出的类型 0 的估计值从 15.4%到 16%不等，但他的估计捕捉了英文名字的影响，而没有孤立地考虑性别影响。
我们看到一个模式，类型 1、2、3 和 4 的系数都是负数，并且随着“外国”的程度绝对值增加。一个人在名字、教育和经验方面越“外国”，回电率就越低。但仅仅一个外国名字就足以使回电率比英文名字低 3.38%到 5.11%。总体而言，结果在 1%的显著水平上是统计学显著的。一个例外是类型 1 的系数，用于英文名和中国姓氏的回归（3）。这里描述的模式与 Oreopoulos（2011）报告的主要发现相匹配。
import statsmodels.formula.api as smf
# Sample Restriction based on name ethnicity
Canada = df2.name_ethnicity == "Canada"
Indian = df2[(Canada) | (df2.name_ethnicity == "Indian")]
Chinese = df2[(Canada) | (df2.name_ethnicity == "Chinese")]
Chn_Cdn = df2[(Canada) | (df2.name_ethnicity == "Chn-Cdn")]
Greek = df2[(Canada) | (df2.name_ethnicity == "Greek")]
sample = [Indian, Chinese, Chn_Cdn, Greek]
# Run the simple model for each ethnicity
# and save the results
model1 = "callback ~ female + C(type)"
result1 = []
for data in sample:
ols = smf.ols(model1, data).fit(cov_type='HC1')
result1.append(ols)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
# Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer
Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5)
WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer(result1)
stargazer.title('Table 1 - Callback Rates by Resume Type')
names = ['Indian', 'Chinese', 'Chn_Cdn', 'Greek']
stargazer.custom_columns(names, [1, 1, 1, 1])
order = ['female', 'Intercept', 'C(type)[T.1.0]',
'C(type)[T.2.0]', 'C(type)[T.3.0]',
'C(type)[T.4.0]']
stargazer.covariate_order(order)
dict1 = {'C(type)[T.1.0]': '1) Foreign Name, Cdn Educ and Exp',
'C(type)[T.2.0]': '2) Foreign Name and Educ, Cdn exp',
'C(type)[T.3.0]': '3) Foreign Name and Educ, Mixed Exp',
'C(type)[T.4.0]': '4) All Foreign (Name, Educ, and Exp)',
'Intercept': '0) English Name, Cdn Educ and Exp',
'female': 'Female'}
stargazer.rename_covariates(dict1)
stargazer
表 1 - 简历类型的回电率
| 女性 |
| 0) 英文名，加拿大教育和经验 |
| 1) 外国名字，加拿大教育和经验 |
| 2) 外国名字和教育，加拿大经验 |
| 3) 外国名字和教育，混合经验 |
| 4) 所有外国人（姓名，教育和经验） |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
表 2 添加了女性和简历类型的交互项。女性的系数仅捕捉了白人女性的影响，因为外国女性是由女性和类型之间的交互项捕捉的。与白人男性（10.17%的基线）相比，白人女性的回访率增加了 7.96%。
交互项的系数在绝对值上为负，但并非全部统计上显着。这种模式表明，外国女性的回访率与白人女性相比非常低。
有趣的是，类型 1、2、3 和 4 的系数在幅度上较低，并且与表 1 相比在统计上不太显着。这种模式表明，白人男性比印度人和中国姓氏的人有优势，但不包括希腊人或中国人（名字和姓氏）。这两个最后一组的系数在统计上不显着。
model2 = "callback ~ female*C(type)"
result2 = []
for data in sample:
ols = smf.ols(model2, data).fit(cov_type='HC1')
result2.append(ols)
stargazer = Stargazer(result2)
stargazer.title('Table 2 - Callback Rates by Resume Type and Gender')
stargazer.custom_columns(names, [1, 1, 1, 1])
dict2 = {'female:C(type)[T.1.0]':'[Female]x[1]',
'female:C(type)[T.2.0]':'[Female]x[2]',
'female:C(type)[T.3.0]':'[Female]x[3]',
'female:C(type)[T.4.0]':'[Female]x[4]'}
list2 = list(dict2.keys())
dict2.update(dict1)
stargazer.rename_covariates(dict2)
list2 = order + list2
stargazer.covariate_order(list2)
stargazer
表 2 - 简历类型和性别的回访率
| 女性 |
| 0) 英文名，加拿大教育和经验 |
| 1) 外国名字，加拿大教育和经验 |
| 2) 外国名字和教育，加拿大经验 |
| 3) 外国名字和教育，混合经验 |
| 4) 所有外国人（姓名，教育和经验） |
| [女性]x[1] |
| [女性]x[2] |
| [女性]x[3] |
| [女性]x[4] |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
表 3 添加了控制变量作为鲁棒性检查。总体结果与表 2 相似。与表 2 相比，白人女性的影响甚至略有增加。白人女性的巨大溢价仍然超过所有其他类别。白人女性的溢价优于来自世界排名前 200 的大学的学士学位，大型公司的经验，课外活动，流利的法语和其他语言以及加拿大硕士学位的累积影响。
请注意，对于类型 1，只有印度回归的系数在统计上显着。白人男性的名字与中国，中国加拿大和希腊名字没有优势。
control1 = "+ ba_quality + extracurricular_skills + language_skills"
control2 = "+ ma + exp_highquality"
model3 = "callback ~ female*C(type)" + control1 + control2
result3 = []
for data in sample:
ols = smf.ols(model3, data).fit(cov_type='HC1')
result3.append(ols)
stargazer = Stargazer(result3)
stargazer.title('Table 3 - Callback Rates and Robustness Checks')
stargazer.custom_columns(names, [1, 1, 1, 1])
dict3 = {'ba_quality':'Top 200 world ranking university',
'exp_highquality':'High quality work experience',
'extracurricular_skills' :'List extra-curricular activities',
'language_skills':'Fluent in French and other languages',
'ma':'Canadian master’s degree'}
list3 = list(dict3.keys())
dict3.update(dict2)
stargazer.rename_covariates(dict3)
list3 = list2 + list3
stargazer.covariate_order(list3)
stargazer
表 3 - 回访率和鲁棒性检查
| 女性 |
| 0) 英文名字，加拿大教育和经验 |
| 1) 外国名字，加拿大教育和经验 |
| 2) 外国名字和教育，加拿大经验 |
| 3) 外国名字和教育，混合经验 |
| 4) 所有外国人（姓名，教育和经验） |
| [女性]x[1] |
| [女性]x[2] |
| [女性]x[3] |
| [女性]x[4] |
| 世界排名前 200 的大学 |
| 高质量的工作经验 |
| 列出课外活动 |
| 流利的法语和其他语言 |
| 加拿大硕士学位 |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
练习
1）为什么白人女性的溢价出现在 2009 年 2 月至 9 月的大衰退期间，而在 2008 年 4 月和 8 月之前没有出现？推测。
2）招聘人员可能更愿意与白人女性交谈，但不一定会雇佣她们。我如何能确定更高的回电率是否反映在更多的工作提供中。例如，我如何获取数据来检查这种关系？
3）你能从下表推断出什么？你有什么见解要分享吗？
pd.crosstab(index= [df2['name_ethnicity'], df2['female'],
df2['name']], columns=df2['type'],
values=df2['callback'], aggfunc='mean')
| 类型 | 0.0 | 1.0 | 2.0 | 3.0 | 4.0 | ||
|---|---|---|---|---|---|---|---|
| 名字 _ 种族 | 女性 | 名字 | |||||
| --- | --- | --- | --- | --- | --- | --- | --- |
| 加拿大 | 0.0 | 格雷格·约翰逊 | 11.561 | NaN | NaN | NaN | NaN |
| 约翰·马丁 | 8.235 | NaN | NaN | NaN | NaN | ||
| 马修·威尔逊 | 10.638 | NaN | NaN | NaN | NaN | ||
| 1.0 | 艾莉森·约翰逊 | 18.675 | NaN | NaN | NaN | NaN | |
| 凯丽·马丁 | 20.455 | NaN | NaN | NaN | NaN | ||
| 吉尔·威尔逊 | 15.205 | NaN | NaN | NaN | NaN | ||
| 中国人 | 0.0 | 刘东 | NaN | 10.870 | 3.390 | 13.158 | 2.381 |
| 李蕾 | NaN | 14.062 | 9.756 | 8.065 | 11.364 | ||
| 张勇 | NaN | 10.227 | 7.407 | 3.509 | 8.333 | ||
| 1.0 | 刘敏 | NaN | 8.235 | 5.357 | 11.321 | 15.556 | |
| 李娜 | NaN | 10.127 | 12.698 | 5.556 | 2.703 | ||
| 张秀英 | NaN | 11.927 | 6.780 | 9.091 | 7.018 | ||
| 中国-加拿大 | 0.0 | 王艾瑞克 | NaN | 9.677 | 3.390 | 8.889 | 0.000 |
| 1.0 | 王美琪 | NaN | 12.500 | 6.780 | 4.348 | 3.571 | |
| 希腊 | 0.0 | 鲁卡斯·米诺普洛斯 | NaN | 10.714 | NaN | NaN | NaN |
| 1.0 | NicoleMinsopoulos | NaN | 9.596 | NaN | NaN | NaN | |
| 印度人 | 0.0 | 阿尔琼·库马尔 | NaN | 8.642 | 8.333 | 3.125 | 4.762 |
| 帕纳夫·辛格 | NaN | 1.333 | 15.942 | 7.317 | 2.500 | ||
| 拉胡尔·考尔 | NaN | 8.571 | 10.769 | 5.455 | 7.317 | ||
| 萨米尔·夏尔马 | NaN | 5.814 | 6.897 | 10.256 | 6.522 | ||
| 1.0 | MayaKumar | NaN | 14.286 | 6.944 | 3.448 | 5.263 | |
| PriyankaKaur | NaN | 6.481 | 14.815 | 3.636 | 5.128 | ||
| ShreyaSharma | NaN | 13.158 | 8.772 | 4.651 | 5.128 | ||
| TaraSingh | NaN | 14.286 | 8.065 | 9.091 | 4.444 |
4）你能从下表推断出什么？你有什么见解要分享吗？
pd.crosstab(index= df2['occupation_type'],
columns=[df2['name_ethnicity'], df2['female']],
values=df2['callback'], aggfunc='mean')
| 姓名种族 | 加拿大 | 中国 | 中国-加拿大 | 希腊 | 印度 |
|---|---|---|---|---|---|
| 女性 | 0.0 | 1.0 | 0.0 | 1.0 | 0.0 |
| --- | --- | --- | --- | --- | --- |
| 职业类型 | |||||
| --- | --- | --- | --- | --- | --- |
| 会计 | 2.703 | 8.929 | 5.769 | 6.780 | 0.000 |
| 行政 | 7.895 | 23.288 | 10.112 | 8.036 | 6.897 |
| 土木工程师 | 5.556 | 50.000 | 6.250 | 6.250 | 0.000 |
| 文书工作 | 4.444 | 8.140 | 5.172 | 6.000 | 0.000 |
| 电子商务 | 0.000 | 0.000 | 9.091 | 0.000 | 0.000 |
| 电气工程师 | 6.250 | 28.571 | 7.143 | 8.333 | 14.286 |
| 行政助理 | 23.077 | 17.647 | 5.263 | 16.000 | 16.667 |
| 金融 | 16.667 | 26.316 | 5.000 | 12.195 | 0.000 |
| 餐饮服务经理 | 16.667 | 16.667 | 0.000 | 8.333 | 20.000 |
| 人力资源工资 | 20.000 | 18.182 | 0.000 | 0.000 | 0.000 |
| 保险 | 53.846 | 40.000 | 14.286 | 13.636 | 28.571 |
| 市场营销和销售 | 12.791 | 22.222 | 12.409 | 11.966 | 9.091 |
| 生产 | 0.000 | 0.000 | 0.000 | 4.762 | 0.000 |
| 程序员 | 10.256 | 17.391 | 13.462 | 10.526 | 0.000 |
| 零售 | 19.048 | 21.622 | 14.545 | 17.647 | 22.727 |
| 技术 | 0.000 | 16.667 | 3.226 | 4.000 | 0.000 |
5）解释表 4 的结果。重点关注固定效应（职业，姓名和城市）的添加。
FE = "+ C(occupation_type) + C(city) + C(name)"
model4 = "callback ~ female*C(type) " + control1 + control2 + FE
result4 = []
for data in sample:
ols = smf.ols(model4, data).fit(cov_type='HC1')
result4.append(ols)
stargazer = Stargazer(result4)
stargazer.title('Table 4 - Callback Rates and Fixed Effects')
stargazer.custom_columns(names, [1, 1, 1, 1])
stargazer.rename_covariates(dict3)
stargazer.covariate_order(list3)
stargazer.add_line('Fixed Effects', ['', '', '', ''])
stargazer.add_line('Occupation', ['Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('Name', ['Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('City', ['Yes', 'Yes', 'Yes', 'Yes'])
stargazer
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:
covariance of constraints does not have full rank. The number of constraints is 43, but rank is 1
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:
covariance of constraints does not have full rank. The number of constraints is 41, but rank is 39
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:
covariance of constraints does not have full rank. The number of constraints is 37, but rank is 35
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:
covariance of constraints does not have full rank. The number of constraints is 31, but rank is 29
表 4 - 回拨率和固定效应
| 印度 |
| 女性 |
| 0) 英文姓名，加拿大教育和经验 |
| 1) 外国姓名，加拿大教育和经验 |
| 2) 外国姓名和教育，加拿大经验 |
| 3) 外国姓名和教育，混合经验 |
| 4) 所有外国（姓名，教育和经验） |
| [女性]x[1] |
| [女性]x[2] |
| [女性]x[3] |
| [女性]x[4] |
| 世界排名前 200 的大学 |
| 高质量工作经验 |
| 列出课外活动 |
| 流利的法语和其他语言 |
| 加拿大硕士学位 |
| 固定效应 |
| 职业 |
| 名字 |
| 城市 |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
6）Oreopoulos (2011)收集的第二波实验数据是在 2008 年 9 月至 11 月之间。使用这些数据来调查在加拿大就业市场中是否拥有白人女性姓名会有额外的优势。只生成一张专业出版表，并解释主要结果。
7）对于这个问题，要像 Bertrand & Mullainathan (2004)和 Oreopoulos (2011)一样打破常规思维。一些研究认为身高较高的人赚更多钱并不是因为身高的直接影响，而是通过自尊心的间接影响。提出一个可行的研究设计来测试以下因果关系：
a) 身高和薪水。
b) 身高和自尊心。
c) 自尊心和薪水。
参考
Bertrand, Marianne, and Sendhil Mullainathan. (2004). Are Emily and Greg More Employable Than Lakisha and Jamal? A Field Experiment on Labor Market Discrimination. American Economic Review, 94 (4): 991-1013.
Oreopoulos, Philip. (2011). Why Do Skilled Immigrants Struggle in the Labor Market? A Field Experiment with Thirteen Thousand Resumes. American Economic Journal: Economic Policy, 3 (4): 148-71.
七、卖淫合法化对犯罪的影响
原文：
causal-methods.github.io/Book/7%29_The_Impact_of_Legalizing_Prostitution_on_Crime.html译者：飞龙
电子邮件：econometrics.methods@gmail.com
最近更新：2020 年 11 月 2 日
在荷兰，有法定的卖淫区，荷兰称之为 tippelzones。Bisschop 等人（2017）报告称，tippelzone 的开放可以减少大约 30-40%的性虐待和强奸案件。
让我们打开 Bisschop 等人的数据集。每一行是荷兰的一个城市。同一个城市在 1994 年至 2011 年之间被观察到。
import numpy as np
import pandas as pd
pd.set_option('precision', 3)
# Data from Bisschop et al. (2017)
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "CBSregist2015.dta")
df.head(5)
| 城市 | 年份 | 开放 | 关闭 | 城市 1 | logpopdens | openingReg | mayorCDA | mayorCU | mayorD66 | ... | 相似盗窃率 | 聚合盗窃率 | 相似盗窃率的自然对数 | 聚合盗窃率的自然对数 | 盗窃率 | 盗窃率的自然对数 | 职务侵犯率 | 职务暴力率 | 职务侵犯率的自然对数 | 职务暴力率的自然对数 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 阿姆斯特丹 | 1994-01-01 | 0.0 | 0.0 | 1.0 | 8.381 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 59.246 | 57.934 | 8.364 | 8.342 | 117.181 | 9.046 | 7.596 | 5.110 | 6.310 | 5.914 |
| 1 | 阿姆斯特丹 | 1995-01-01 | 0.0 | 0.0 | 1.0 | 8.379 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 50.815 | 43.823 | 8.208 | 8.060 | 94.637 | 8.830 | 7.061 | 4.361 | 6.234 | 5.753 |
| 2 | 阿姆斯特丹 | 1996-01-01 | 1.0 | 0.0 | 1.0 | 8.373 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 42.333 | 37.111 | 8.020 | 7.888 | 79.444 | 8.649 | 7.520 | 5.431 | 6.292 | 5.966 |
| 3 | 阿姆斯特丹 | 1997-01-01 | 1.0 | 0.0 | 1.0 | 8.369 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 46.843 | 32.860 | 8.117 | 7.762 | 79.704 | 8.648 | 6.852 | 4.195 | 6.194 | 5.704 |
| 4 | 阿姆斯特丹 | 1998-01-01 | 1.0 | 0.0 | 1.0 | 8.373 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 45.255 | 33.907 | 8.086 | 7.798 | 79.162 | 8.646 | 6.127 | 4.595 | 6.087 | 5.799 |
5 行×65 列
让我们将城市分成 3 组。大城市和中等城市至少在某一年有 tippelzone，而样本中的其他城市没有 tippelzone。
big_cities = ["Amsterdam", "Rotterdam", "Den Haag"]
medium_cities = ["Utrecht", "Nijmegen", "Groningen",
"Heerlen", "Eindhoven", "Arnhem"]
# Classify cities
def classify(var):
if var in big_cities:
return "Big Cities"
elif var in medium_cities:
return "Medium Cities"
else:
return "No tippelzone"
df['group'] = df["city"].apply(classify)
以下是每 10,000 名居民的年度犯罪报告。总体而言，大城市的犯罪率更高。唯一的例外是与毒品有关的犯罪。
outcome = ['sexassaultpcN', 'rapepcN',
'drugspcN', 'maltreatpcN', 'weaponspcN']
df.groupby('group')[outcome].mean().T
| 组 | 大城市 | 中等城市 | 无 tippelzone |
|---|---|---|---|
| 性侵犯率 | 0.775 | 0.626 | 0.664 |
| 强奸率 | 1.032 | 0.846 | 0.691 |
| 毒品犯罪率 | 14.921 | 15.599 | 12.779 |
| maltreatpcN | 21.259 | 18.665 | 17.864 |
| weaponspcN | 5.635 | 4.385 | 4.207 |
tippelzones 城市的人口和人口密度比没有 tippelzones 的城市更多。平均家庭收入（以 1,000 欧元计）在 3 个组中相似。tippelzones 城市也有受过较高教育的个体。移民比例在大城市中更高（11.4%）。社会保险福利（“insurWWAO_pc”）的份额与 3 个组相似。
demographics = ['popul_100', 'pop_dens_100', 'popmale1565_100',
'inkhh', 'educhpc', 'nondutchpc', 'insurWWAO_pc']
df.groupby('group')[demographics].mean().T
| 组 | 大城市 | 中等城市 | 无 tippelzone |
|---|---|---|---|
| popul_100 | 5974.886 | 1724.191 | 1131.138 |
| pop_dens_100 | 43.258 | 22.977 | 19.560 |
| popmale1565_100 | 2101.446 | 617.019 | 392.255 |
| inkhh | 29.052 | 28.989 | 30.502 |
| educhpc | 0.300 | 0.317 | 0.245 |
| 非荷兰人比例 | 0.114 | 0.059 | 0.052 |
| insurWWAO_pc | 0.074 | 0.081 | 0.078 |
基督教联盟在没有 tippelzone 的城市中拥有更多的市长（31%）。值得一提的是，该党反对开放 tippelzone。
political_party = ['mayorSoc', 'mayorLib', 'mayorChr']
df.groupby('group')[political_party].mean().T
| 组 | 大城市 | 中等城市 | 无 tippelzone |
|---|---|---|---|
| mayorSoc | 0.481 | 0.556 | 0.410 |
| mayorLib | 0.259 | 0.324 | 0.278 |
| mayorChr | 0.259 | 0.120 | 0.312 |
数据集是平衡的面板数据。有必要按顺序声明指数：分析单位和时间单位。
df['year'] = pd.DatetimeIndex(df['year']).year
df['Dyear'] = pd.Categorical(df.year)
# Set Panel Data
# Set city as the unit of analysis
df25 = df.set_index(['city1', 'year'])
设\(Y_{ct}\)为城市\(c\)在年份\(t\)的犯罪率。设\(D_{ct}\) = 1，如果城市\(c\)在年份\(t\)有开放的 tippelzone；否则为 0。让我们估计以下模型：
其中\(\alpha_c\)是城市固定效应，\(X_{ct}\)是控制变量向量，\(\gamma_t\)是年固定效应，\(\epsilon_{ct}\)是通常的误差项。
import statsmodels.formula.api as smf
Ys = ["lnsexassaultN", "lnrapeN", "lndrugsN",
"lnweaponsN", "lnmaltreatN"]
base = "~ 1 + opening"
fe = "+ C(city) + C(Dyear)"
controls = ['logpopmale1565', 'logpopdens', 'inkhh',
'educhpc', 'nondutchpc', 'insurWWAO', 'mayorCDA',
'mayorCU', 'mayorD66', 'mayorVVD']
Xs = ""
for var in controls:
Xs = Xs + '+' + var
columns = []
for Y in Ys:
result = smf.ols(Y + base + fe + Xs, df25).fit(cov_type='cluster',
cov_kwds={'groups': df25['city']})
columns.append(result)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
import pandas.util.testing as tm
# Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer
Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5)
WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
第 1 列表明开放 tippelzone 将性虐待减少 26%（\(e^{-0.302}-1\)）。在其他列中，tippelzone 的系数在统计上不显著。看起来合法化卖淫会减少性虐待，但不会减少其他犯罪，如强奸、攻击、非法武器和毒品相关犯罪。
# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer(columns)
stargazer.title('The Impact of Tippelzone on Crime')
names = ['Sex Abuse', 'Rape', 'Drugs', 'Weapons', 'Assault']
stargazer.custom_columns(names, [1, 1, 1, 1, 1])
stargazer.covariate_order(['opening'])
stargazer.add_line('Covariates', ['Yes', 'Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('City Fixed Effects', ['Yes', 'Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('Year Fixed Effects', ['Yes', 'Yes', 'Yes', 'Yes', 'Yes'])
stargazer
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
'rank is %d' % (J, J_), ValueWarning)
Tippelzone 对犯罪的影响
| 开放 |
| 协变量 |
| 城市固定效应 |
| 年固定效应 |
| 观察 |
| R² |
| 调整 R² |
| 残差标准误差 |
| F 统计量 |
| 注意： |
import math
math.exp(-0.302) - 1
-0.2606619351104681
练习
1）在论文的引言部分，Bisschop 等人（2017: 29）表示：“我们的研究是第一个提供卖淫监管与犯罪之间关联的因果证据之一。”在讨论部分，Bisschop 等人（2017:44）表示：“开放 tippelzone，无论是否有许可制度，都与性虐待和强奸的短期减少 30-40%相关，并且结果在不同规范下都是稳健的。”为什么 Bisschop 等人（2017）在引言部分使用“因果”一词，在讨论部分使用“相关”一词？您是否认为 Bisschop 等人（2017）的主要结果是“因果”还是“相关”？请解释。
2）Bisschop 等人（2017: 29）表示：“我们进行了几项实证测试，以评估 tippelzone 开放周围的内生犯罪趋势。”他们为什么这样做？这其中的逻辑是什么？是否存在内生犯罪趋势？请解释并具体说明您的答案。
3）Bisschop 等人（2017: 36）表示：“...时间趋势\(\mu_t\)是使用年固定效应来建模的。”模拟时间趋势的其他方法是什么？编写不同假设下创建时间趋势的代码片段。提示：记住这是面板数据。在横截面数据中有效的代码将在面板数据结构中创建错误的变量。
4）Bisschop 等人（2017: 36）表示：“我们使用差异中的差异规范来研究 tippelzone 存在对各种犯罪的影响。”部署差异中的差异估计器的关键假设是什么？
5）复制表格“Tippelzone 对犯罪的影响”，不包括阿姆斯特丹、鹿特丹和海牙。此外，用以下四个变量替换变量“opening”：
i) “everopenNoReg”: 如果城市\(c\)在年份\(t\)之前曾开放过没有许可证的 tippelzone，则为 1，否则为 0。
ii) “openingRegP”: 如果城市\(c\)在年份\(t\)之前开放了 tippelzone 并引入了事后许可证，则为 1，否则为 0。
iii) “openingRegA”: 如果城市\(c\)在年份\(t\)之前曾开放过带有许可证的 tippelzone，则为 1，否则为 0。
iv) “closing”: 如果城市\(c\)在年份\(t\)之前关闭 tippelzone，则为 1，否则为 0。
解释结果。
参考
Bisschop, Paul, Stephen Kastoryano, and Bas van der Klaauw. (2017). 街头卖淫区和犯罪. 美国经济学杂志：经济政策，9 (4): 28-63.
八、Airbnb 的主人是否歧视黑人客人？
原文：
causal-methods.github.io/Book/8%29_Do_Hosts_Discriminate_against_Black_Guests_in_Airbnb.html译者：飞龙
电子邮件：econometrics.methods@gmail.com
最后更新时间：11-1-2020
Edelman et al.（2017）发现，黑人名字听起来比白人名字听起来更不可能被 Airbnb 接受为客人，减少了 16%。这个结果不仅仅是相关性。种族变量是随机的。黑人和白人之间唯一的区别是名字。除此之外，黑人和白人客人是一样的。
让我们打开 Edelman 等人（2017）的数据集。每一行是 2015 年 7 月 Airbnb 的一处物业。样本由巴尔的摩、达拉斯、洛杉矶、圣路易斯和华盛顿特区的所有物业组成。
import numpy as np
import pandas as pd
pd.set_option('precision', 3)
# Data from Edelman et al. (2017)
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_csv(path + "Airbnb.csv")
df.head(5)
| 主人回应 | 回应日期 | 消息数量 | 自动编码 | 纬度 | 经度 | 床类型 | 物业类型 | 取消政策 | 客人数量 | ... | 洛杉矶 | 圣路易斯 | 华盛顿特区 | 总客人 | 原始黑人 | 物业黑人 | 任何黑人 | 过去客人合并 | 九月填充 | pr 填充 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 是 | 2015-07-19 08:26:17 | 2.0 | 1.0 | 34.081 | -118.270 | 真正的床 | 房子 | 灵活 | 3.0 | ... | 1 | 0 | 0 | 11.0 | 0.0 | 0.0 | 0.0 | 匹配（3） | 1 | 0.412 |
| 1 | 否或不可用 | 2015-07-14 14:13:39 | NaN | 1.0 | 38.911 | -77.020 | NaN | 房子 | 中等 | 2.0 | ... | 0 | 0 | 1 | 167.0 | 0.0 | 0.0 | 0.0 | 匹配（3） | 1 | 0.686 |
| 2 | 请求更多信息（你能验证吗？有多少... | 2015-07-20 16:24:08 | 2.0 | 0.0 | 34.005 | -118.481 | 拉出沙发 | 公寓 | 严格 | 1.0 | ... | 1 | 0 | 0 | 19.0 | 0.0 | 0.0 | 0.0 | 匹配（3） | 0 | 0.331 |
| 3 | 我会回复你 | 2015-07-20 06:47:38 | NaN | 0.0 | 34.092 | -118.282 | NaN | 房子 | 严格 | 8.0 | ... | 1 | 0 | 0 | 41.0 | 0.0 | 0.0 | 0.0 | 匹配（3） | 0 | 0.536 |
| 4 | 未发送消息 | . | NaN | 1.0 | 38.830 | -76.897 | 真正的床 | 房子 | 严格 | 2.0 | ... | 0 | 0 | 1 | 28.0 | 0.0 | 0.0 | 0.0 | 匹配（3） | 1 | 0.555 |
5 行×104 列
下面的图表显示，黑人客人收到的“是”的回应比白人客人少。有人可能会争辩说 Edelman 等人（2017）的结果是由主人回应的差异驱动的，比如有条件的或非回应。例如，你可以争辩说黑人更有可能有被归类为垃圾邮件的假账户。然而，请注意，歧视结果是由“是”和“否”驱动的，而不是由中间回应驱动的。
# Data for bar chart
count = pd.crosstab(df["graph_bins"], df["guest_black"])
import plotly.graph_objects as go
node = ['Conditional No', 'Conditional Yes', 'No',
'No Response', 'Yes']
fig = go.Figure(data=[
go.Bar(name='Guest is white', x=node, y=count[0]),
go.Bar(name='Guest is African American', x=node, y=count[1]) ])
fig.update_layout(barmode='group',
title_text = 'Host Responses by Race',
font=dict(size=18) )
fig.show()
让我们复制 Edelman 等人（2017）的主要结果。
import statsmodels.api as sm
df['const'] = 1
# Column 1
# The default missing ='drop' of statsmodels doesn't apply
# to the cluster variable. Therefore, it is necessary to drop
# the missing values like below to get the clustered standard
# errors.
df1 = df.dropna(subset=['yes', 'guest_black', 'name_by_city'])
reg1 = sm.OLS(df1['yes'], df1[['const', 'guest_black']])
res1 = reg1.fit(cov_type='cluster',
cov_kwds={'groups': df1['name_by_city']})
# Column 2
vars2 = ['yes', 'guest_black', 'name_by_city',
'host_race_black', 'host_gender_M']
df2 = df.dropna(subset = vars2)
reg2 = sm.OLS(df2['yes'], df2[['const', 'guest_black',
'host_race_black', 'host_gender_M']])
res2 = reg2.fit(cov_type='cluster',
cov_kwds={'groups': df2['name_by_city']})
# Column 3
vars3 = ['yes', 'guest_black', 'name_by_city',
'host_race_black', 'host_gender_M',
'multiple_listings', 'shared_property',
'ten_reviews', 'log_price']
df3 = df.dropna(subset = vars3)
reg3 = sm.OLS(df3['yes'], df3[['const', 'guest_black',
'host_race_black', 'host_gender_M',
'multiple_listings', 'shared_property',
'ten_reviews', 'log_price']])
res3 = reg3.fit(cov_type='cluster',
cov_kwds={'groups': df3['name_by_city']})
columns =[res1, res2, res3]
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
# Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer
Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5)
WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
在第一列中，听起来像白人的名字被接受的概率为 49%；而听起来像黑人的名字被接受的概率大约为 41%。因此，黑人名字带来了 8%的惩罚。这个结果在第 2 列和第 3 列的一组控制变量中非常稳健。
# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer(columns)
stargazer.title('The Impact of Race on Likelihood of Acceptance')
stargazer
种族对接受可能性的影响
| 常数 |
| 黑人客人 |
| 主人性别 M |
| 主人种族黑人 |
| 对数价格 |
| 多个列表 |
| 共享物业 |
| 十次评论 |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
下表显示了有关主人和房产的摘要统计信息。在实验中，控制变量的平均值与分组实验和对照组的平均值相同。
control = ['host_race_white', 'host_race_black', 'host_gender_F',
'host_gender_M', 'price', 'bedrooms', 'bathrooms', 'number_of_reviews',
'multiple_listings', 'any_black', 'tract_listings', 'black_proportion']
df.describe()[control].T
| 计数 | 平均值 | 标准差 | 最小值 | 25% | 50% | 75% | 最大值 | |
|---|---|---|---|---|---|---|---|---|
| host_race_white | 6392.0 | 0.634 | 0.482 | 0.0 | 0.00 | 1.00 | 1.000 | 1.000 |
| host_race_black | 6392.0 | 0.078 | 0.269 | 0.0 | 0.00 | 0.00 | 0.000 | 1.000 |
| host_gender_F | 6392.0 | 0.376 | 0.485 | 0.0 | 0.00 | 0.00 | 1.000 | 1.000 |
| host_gender_M | 6392.0 | 0.298 | 0.457 | 0.0 | 0.00 | 0.00 | 1.000 | 1.000 |
| 价格 | 6302.0 | 181.108 | 1280.228 | 10.0 | 75.00 | 109.00 | 175.000 | 100000.000 |
| 卧室 | 6242.0 | 3.177 | 2.265 | 1.0 | 2.00 | 2.00 | 4.000 | 16.000 |
| 浴室 | 6285.0 | 3.169 | 2.264 | 1.0 | 2.00 | 2.00 | 4.000 | 16.000 |
| 评论数量 | 6390.0 | 30.869 | 72.505 | 0.0 | 2.00 | 9.00 | 29.000 | 1208.000 |
| 多个列表 | 6392.0 | 0.326 | 0.469 | 0.0 | 0.00 | 0.00 | 1.000 | 1.000 |
| 任何黑人 | 6390.0 | 0.282 | 0.450 | 0.0 | 0.00 | 0.00 | 1.000 | 1.000 |
| 地域列表 | 6392.0 | 9.514 | 9.277 | 1.0 | 2.00 | 6.00 | 14.000 | 53.000 |
| 黑人比例 | 6378.0 | 0.140 | 0.203 | 0.0 | 0.03 | 0.05 | 0.142 | 0.984 |
平衡处理测试(t 检验)显示黑人和白人客人是相同的。
result = []
for var in control:
# Do the T-test and save the p-value
pvalue = sm.OLS(df[var], df[['const', 'guest_black']],
missing = 'drop').fit().pvalues[1]
result.append(pvalue)
ttest = df.groupby('guest_black').agg([np.mean])[control].T
ttest['p-value'] = result
ttest
| 黑人客人 | 0.0 | 1.0 | p 值 | |
|---|---|---|---|---|
| host_race_white | 平均值 | 0.643 | 0.626 | 0.154 |
| host_race_black | 平均值 | 0.078 | 0.078 | 0.972 |
| host_gender_F | 平均值 | 0.381 | 0.372 | 0.439 |
| host_gender_M | 平均值 | 0.298 | 0.299 | 0.896 |
| 价格 | 平均值 | 166.429 | 195.815 | 0.362 |
| 卧室 | 平均值 | 3.178 | 3.176 | 0.962 |
| 浴室 | 平均值 | 3.172 | 3.167 | 0.927 |
| 评论数量 | 平均值 | 30.709 | 31.030 | 0.860 |
| 多个列表 | 平均值 | 0.321 | 0.330 | 0.451 |
| 任何黑人 | 平均值 | 0.287 | 0.277 | 0.382 |
| 地域列表 | 平均值 | 9.494 | 9.538 | 0.848 |
| 黑人比例 | 平均值 | 0.141 | 0.140 | 0.919 |
练习
1）据我所知，关于种族歧视的文献中最重要的三篇实证论文是 Bertrand & Mullainathan (2004), Oreopoulos (2011), 和 Edelman et al. (2017)。这三篇论文都使用了实地实验来捕捉因果关系并排除混杂因素。在互联网上搜索并返回一份关于种族歧视的实验论文的参考列表。
2）告诉我一个你热衷的话题。返回一个关于你的话题的实验论文的参考列表。
3）有人认为特定的名字驱动了 Edelman 等人(2017)的结果。在下面的表格中，你可以看到代表黑人和白人的名字并不多。如何反驳这个批评？你可以做什么来证明结果不是由特定的名字驱动的？
female = df['guest_gender']=='female'
df[female].groupby(['guest_race', 'guest_first_name'])['yes'].mean()
guest_race guest_first_name
black Lakisha 0.433
Latonya 0.370
Latoya 0.442
Tamika 0.482
Tanisha 0.413
white Allison 0.500
Anne 0.567
Kristen 0.486
Laurie 0.508
Meredith 0.498
Name: yes, dtype: float64
male = df['guest_gender']=='male'
df[male].groupby(['guest_race', 'guest_first_name'])['yes'].mean()
guest_race guest_first_name
black Darnell 0.412
Jamal 0.354
Jermaine 0.379
Kareem 0.436
Leroy 0.371
Rasheed 0.409
Tyrone 0.377
white Brad 0.419
Brent 0.494
Brett 0.466
Greg 0.467
Jay 0.581
Todd 0.448
Name: yes, dtype: float64
4）根据下表，是否有任何潜在的研究问题可以探讨？请证明。
pd.crosstab(index= [df['host_gender_F'], df['host_race']],
columns=[df['guest_gender'], df['guest_race']],
values=df['yes'], aggfunc='mean')
| 客人性别 | 女性 | 男性 | |
|---|---|---|---|
| 客人种族 | 黑人 | 白人 | |
| --- | --- | --- | --- |
| host_gender_F | host_race | ||
| --- | --- | --- | --- |
| 0 | UU | 0.400 | 0.542 |
| 亚裔 | 0.319 | 0.378 | 0.474 |
| 黑人 | 0.444 | 0.643 | 0.419 |
| 西班牙裔 | 0.464 | 0.571 | 0.375 |
| 多 | 0.568 | 0.727 | 0.408 |
| 不明确 | 0.444 | 0.500 | 0.444 |
| 不明确的三票 | 0.476 | 0.392 | 0.368 |
| 白人 | 0.383 | 0.514 | 0.386 |
| 1 | UU | 0.444 | 0.250 |
| 亚裔 | 0.429 | 0.607 | 0.436 |
| 黑人 | 0.603 | 0.537 | 0.397 |
| 西班牙裔 | 0.391 | 0.667 | 0.292 |
| 不明确 | 0.600 | 0.556 | 0.125 |
| 不明确的三票 | 0.387 | 0.583 | 0.312 |
| 白人 | 0.450 | 0.494 | 0.370 |
5）在 Edelman 等人（2017 年）中，变量“name_by_city”被用来对标准误差进行聚类。变量“name_by_city”是如何基于其他变量创建的？展示代码。
6）使用 Edelman 等人（2017 年）的数据来测试同族偏好假设，即主人可能更喜欢相同种族的客人。使用 Stargazer 库生成一个漂亮的表格。解释结果。
7）总的来说，人们知道社会经济地位与种族有关。Fryer＆Levitt（2004 年）表明，独特的非洲裔美国人名字与较低的社会经济地位相关。Edelman 等人（2017 年：17）明确表示：“我们的发现无法确定歧视是基于种族，社会经济地位，还是这两者的结合。”提出一个实验设计来分离种族和社会经济地位的影响。解释您的假设并详细描述程序。
参考
Bertrand，Marianne 和 Sendhil Mullainathan。 （2004）。艾米丽和格雷格比拉基莎和贾迈尔更受雇用吗？劳动市场歧视的实地实验。《美国经济评论》，94（4）：991-1013。
Edelman，Benjamin，Michael Luca 和 Dan Svirsky。 （2017）。共享经济中的种族歧视：来自实地实验的证据。《美国经济学杂志：应用经济学》，9（2）：1-22。
Fryer，Roland G. Jr.和 Steven D. Levitt。 （2004）。Distinctively Black Names 的原因和后果。《经济学季刊》，119（3）：767–805。
Oreopoulos，Philip。 （2011）。为什么技术移民在劳动市场上挣扎？一项涉及一万三千份简历的实地实验。《美国经济学杂志：经济政策》，3（4）：148-71。
九、股票市场如何缓解以巴冲突？
译者：飞龙
电子邮件：econometrics.methods@gmail.com
最近更新：11-5-2020
“商业是最具破坏性偏见的良药。”（孟德斯鸠，1748 年：第 II 卷，第一章）
Jha & Shayo（2019）随机将 1345 名犹太以色列选民分为金融资产治疗组和对照组。他们报告称，接触股票市场会增加 4-6%的可能性，投票支持主张和平解决冲突的政党。让我们打开 Jha & Shayo（2019）的数据集。每一行都是以色列公民。
import numpy as np
import pandas as pd
pd.set_option('precision', 3)
# Data from Jha & Shayo (2019)
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "replicationdata.dta")
df.head(5)
C:\Anaconda\envs\textbook\lib\site-packages\pandas\io\stata.py:1433: UnicodeWarning:
One or more strings in the dta file could not be decoded using utf-8, and
so the fallback encoding of latin-1 is being used. This can happen when a file
has been incorrectly encoded by Stata or some other software. You should verify
the string values returned are correct.
warnings.warn(msg, UnicodeWarning)
| 统计 | tradestock6 | 愿意承担 1 到 10 的风险 | 用户 ID | 男性 | 关系 | 关系 ID | nafa | 教育 | ses | ... | tradetot_m4 | pricechange_m4 | bought_bm4 | sold_bm4 | active_bm4 | 资产类型 | 上周 | 过去 3 年 | 下周 | facts_0_m4 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 完成 | NaN | 3 | 60814.0 | 0.0 | 犹太教 | 世俗 | 耶路撒冷 | 硕士 | 平均以下 | ... | 3.0 | -0.673 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | 2.0 |
| 1 | 完成 | 1.0 | 2 | 60824.0 | 0.0 | 犹太教 | 世俗 | 特拉维夫 | 文学士 | 平均以上 | ... | 3.0 | 5.323 | 2.0 | 0.0 | 2.0 | 1.0 | 1.0 | 1.0 | 0.0 | 3.0 |
| 2 | 完成 | 0.0 | 3 | 61067.0 | 1.0 | 犹太教 | 世俗 | 中心 | 博士 | 平均以上 | ... | 3.0 | 0.000 | 2.0 | 1.0 | 3.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 3 | 完成 | NaN | 4 | 61095.0 | 1.0 | 犹太教 | 世俗 | 海法 | 文学士学生 | 平均以下 | ... | 3.0 | 5.323 | 3.0 | 2.0 | 3.0 | 0.0 | 0.0 | 1.0 | 0.0 | 1.0 |
| 4 | 完成 | 0.0 | 4 | 61198.0 | 1.0 | 犹太教 | 世俗 | 北部 | 硕士 | 平均以下 | ... | 3.0 | 0.000 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 1.0 |
5 行×526 列
# Drop missing values
df = df.dropna(subset=['left_s3'])
图 1 显示，2013 年选举中，对照组和治疗组的投票情况相似。但在 2015 年，左翼党派在对照组中获得了 24.8%的选票，而在治疗组中获得了 30.9%。相反，右翼党派在对照组中获得了 35.8%的选票，而在治疗组中获得了 31.2%。根据 Jha & Shayo（2019）的说法，右翼和左翼党派在经济政策方面是相似的。主要区别在于左翼党派支持和平进程，而右翼党派认为任何为和平做出的让步都会对以色列国家构成风险。
# Data: Vote Share by year
v2013 = df.groupby('assettreat')['left_2013', 'right_2013'].mean().T
v2015 = df.groupby('assettreat')['left_s3', 'right_s3'].mean().T
prop = v2013.append(v2015)
# Plot Bar Chart
import plotly.graph_objects as go
node = ['Left 2013', 'Right 2013', 'Left 2015', 'Right 2015']
fig = go.Figure(data=[
go.Bar(name='Control', x=node, y = prop[0]),
go.Bar(name='Treatment', x=node, y = prop[1]) ])
fig.update_layout(barmode='group',
title_text = 'Graphic 1 - Elections: 2013 vs 2015 ',
font=dict(size=18) )
fig.update_yaxes(title_text = "Vote Share")
fig.show()
C:\Anaconda\envs\textbook\lib\site-packages\ipykernel_launcher.py:2: FutureWarning: Indexing with multiple keys (implicitly converted to a tuple of keys) will be deprecated, use a list instead.
C:\Anaconda\envs\textbook\lib\site-packages\ipykernel_launcher.py:3: FutureWarning: Indexing with multiple keys (implicitly converted to a tuple of keys) will be deprecated, use a list instead.
This is separate from the ipykernel package so we can avoid doing imports until
下表显示，治疗组与对照组相似。例外情况是年龄和愿意承担风险。治疗组的以色列人比对照组年轻（39.3 岁对 41.5 岁）。治疗组在愿意承担风险方面也更偏好，评估指数从 1 到 10 变化（4.7 对 4.3）。
我们根据 OLS 回归和分层固定效应计算 p 值。
control = ['right_2013', 'left_2013', 'p_index_s1', 'e_index_init',
'tradestock6all', 'male', 'age', 'postsecondary', 'BA_student',
'college_grad', 'married', 'r_sec', 'r_trad', 'r_relig', 'r_ultra',
'g_jerusalem', 'g_north', 'g_haifa', 'g_center', 'g_telaviv', 'g_south',
'g_wb', 'faminc', 'willingrisk1to10', 'patient', 'plitscore']
import statsmodels.formula.api as smf
result = []
for var in control:
# OLS with 104 randomization strata fixed effects
reg = smf.ols(var + "~ 1 + assettreat + C(block13)", df)
# 104 is the last variable: the coefficient of assettreat
pvalue = reg.fit().pvalues[104]
result.append(pvalue)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
table = df.groupby('assettreat')[control].mean().T
table['p-value'] = result
table
| 资产治疗 | 0.0 | 1.0 | p 值 |
|---|---|---|---|
| right_2013 | 0.245 | 0.241 | 0.964 |
| left_2013 | 0.126 | 0.137 | 0.213 |
| p_index_s1 | 0.004 | 0.051 | 0.399 |
| e_index_init | -0.005 | 0.007 | 0.752 |
| tradestock6all | 0.368 | 0.355 | 0.290 |
| 男性 | 0.513 | 0.521 | 0.470 |
| 年龄 | 41.530 | 39.289 | 0.011 |
| 高中后教育 | 0.232 | 0.230 | 0.953 |
| 文学士学生 | 0.152 | 0.148 | 0.834 |
| 大学毕业 | 0.427 | 0.426 | 0.860 |
| 已婚 | 0.629 | 0.598 | 0.295 |
| r_sec | 0.636 | 0.627 | 0.582 |
| r_trad | 0.172 | 0.164 | 0.823 |
| r_relig | 0.119 | 0.124 | 0.780 |
| r_ultra | 0.073 | 0.085 | 0.222 |
| g_jerusalem | 0.096 | 0.091 | 0.800 |
| g_north | 0.089 | 0.097 | 0.595 |
| g_haifa | 0.123 | 0.142 | 0.290 |
| g_center | 0.298 | 0.290 | 0.766 |
| g_telaviv | 0.212 | 0.194 | 0.276 |
| g_south | 0.116 | 0.104 | 0.596 |
| g_wb | 0.066 | 0.081 | 0.341 |
| 家庭收入 | 11162.162 | 10996.970 | 0.511 |
| 愿意承担风险 1 到 10 | 4.344 | 4.716 | 0.009 |
| 患者 | 0.642 | 0.657 | 0.645 |
| plitscore | 69.726 | 70.664 | 0.550 |
# Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer
Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5)
WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
表 1 的第 1 列显示了意向治疗（ITT）的估计值为 6.1%。值得一提的是，这种类型的实验很少有完全的遵从。通常，在需要随时间跟踪个体时，控制组和治疗组都会有一定的流失率。
第 2 列显示 ITT 效应对控制变量和分层固定效应的鲁棒性。
第 3 列使用变量“vote_wgt”作为权重呈现了加权最小二乘法（WLS）。Jha & Shayo（2019）提取了一个随机样本，其中非正统中心选民被过度抽样。逻辑是增加对最有趣的群体：摇摆选民的精度。使用权重可以在不过度抽样的情况下重现结果。
第 4 列显示结果是由接受以色列股票（“isrstock”）和投资券（“cash”）的治疗组个体驱动的。要小心得出巴勒斯坦股票没有影响的结论。在实验期间，以色列股票价格上涨，但巴勒斯坦股票下跌。
ITT1 = smf.ols("left_s3 ~ 1 + assettreat",
df).fit()
Xs = ['right_2013', 'left_2013', 'male', 'age', 'age2',
'postsecondary', 'BA_student', 'college_grad',
'married', 'tradestock6all', 'r_trad', 'r_relig',
'r_ultra', 'g_jerusalem', 'g_north', 'g_haifa',
'g_telaviv', 'g_south', 'g_wb', 'C(newses)',
'willingrisk1to10', 'patient', 'plitscore']
controls = ""
for X in Xs:
controls = controls + '+' + X
ITT2 = smf.ols("left_s3 ~ 1 + assettreat" + controls +
"+C(block13)", df).fit()
WLS = smf.wls("left_s3 ~ 1 + assettreat" + controls +
"+C(block13)", df, weights=df['vote_wgt']).fit()
treatments = "+ isrstock + palstock + cash"
WLS2 = smf.wls("left_s3 ~ 1" + treatments + controls +
"+C(block13)", df, weights=df['vote_wgt']).fit()
# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer([ITT1, ITT2, WLS, WLS2])
stargazer.title('Table 1 - Intent to Treat Estimates of Stock'
+ ' Exposure on Voting for a Left Party')
names = ['ITT', 'ITT', 'WLS', 'WLS']
stargazer.custom_columns(names, [1, 1, 1, 1])
stargazer.covariate_order(['assettreat', 'isrstock',
'palstock', 'cash'])
stargazer.add_line('Strata Fixed Effects', ['No', 'Yes',
'Yes', 'Yes'])
stargazer.add_line('Covariates', ['No', 'Yes', 'Yes', 'Yes'])
stargazer
表 1 - 对左翼党投票的意向治疗估计股票暴露
| 资产处理 |
| isrstock |
| palstock |
| 现金 |
| 分层固定效应 |
| 协变量 |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
表 2 的第 1 列呈现了第一阶段回归。ITT（“资产处理”）是符合分配（“资产 _comp”）的工具变量（IV）。变量“资产 _comp”表示谁实际完成了实验。在控制了几个协变量和分层固定效应之后，ITT 的系数在统计上是显著的。变量“资产处理”是一个完美的 IV。这个变量是随机的。因此，它与误差项不相关。
第 2 列显示了控制函数方法（CF）的结果。对待处理的治疗效应（TOT）估计为 7.3%。接受治疗的个体投票给左翼党的概率增加了 7.3%。在这个框架中，CF 等同于 2SLS。在 CF 中，我们使用第一阶段的残差（\(\hat{u}\)）来控制第二阶段的内生性。注意残差在统计上是显著的。因此，需要进行修正。我们可以得出结论，表 1 低估了金融资产暴露对左翼党投票的影响。
# Fist Stage
FS = smf.ols("asset_comp ~ 1 + assettreat" + controls +
"+C(block13)", df).fit()
# Control Function Approach
df['resid'] = FS.resid
CF = smf.ols("left_s3 ~ 1 + asset_comp + resid" + controls +
"+C(block13)", df).fit()
# Settings for a nice table
stargazer = Stargazer([FS, CF])
stargazer.title('Table 2 - Impact of Stock Exposure'
' on Voting for a Left Party')
names = ['First Stage', 'Control Function']
stargazer.custom_columns(names, [1, 1])
stargazer.covariate_order(['assettreat', 'asset_comp', 'resid'])
stargazer.add_line('Strata Fixed Effects', ['Yes', 'Yes'])
stargazer.add_line('Covariates', ['Yes', 'Yes'])
stargazer
表 2 - 股票暴露对左翼党投票的影响
| 资产处理 |
| 资产 _comp |
| 残差 |
| 分层固定效应 |
| 协变量 |
| 观察 |
| R² |
| 调整后的 R² |
| 残差标准误差 |
| F 统计量 |
| 注： |
练习
1）Jha＆Shayo（2019）的样本由以色列选民组成。推测在巴勒斯坦选民的情况下，结果是否会在质量上相同。证明你的理由。
2）从 Jha＆Shayo（2019）出发，有什么有前途的研究问题？
3）什么是社会期望偏差？描述 Jha＆Shayo（2019）为减轻社会期望偏差所做的工作。
4）使用 2SLS 而不是控制函数方法复制表 2 的结果。不要使用库“linearmodels”。使用库“statsmodels”手动进行 2SLS。
5）暴露于巴勒斯坦股票是否会降低投票右翼政党的概率？运行一些回归来证明你的立场。
参考
Jha，S.和 Shayo，M.（2019）。估价和平：金融市场暴露对选票和政治态度的影响。计量经济学，87：1561-1588。
蒙特斯基耶，C.（1748）。法律精神。伦敦：T. Evans，1777 年，4 卷。第 2 卷。在线自由图书馆。
参考书目
原文：
causal-methods.github.io/Book/Bibliography.html译者：飞龙
作者：Vitor Kamada
电子邮件：econometrics.methods@gmail.com
最后更新日期 11-6-2020
书中使用的论文和数据：Python 因果推断。
巴斯滕，克里斯托夫和弗兰克贝茨（2013）。超越职业道德：宗教，个人和政治偏好。《美国经济学杂志：经济政策》，5（3）：67-91。在线附录。finaldata.dta
贝克尔，萨沙 O.，和 Ludger Woessmann。 （2009）。韦伯错了吗？新教经济史的人力资本理论。《季度经济学杂志》124（2）：531-96。
Bertrand，Marianne 和 Sendhil Mullainathan。 （2004）。艾米丽和格雷格比拉基莎和贾迈尔更受雇用吗？。劳动力市场歧视的实地实验。《美国经济评论》，94（4）：991-1013。lakisha_aer.dta
Bisschop，Paul，Stephen Kastoryano 和 Bas van der Klaauw。 （2017）。街头卖淫区和犯罪。《美国经济学杂志：经济政策》，9（4）：28-63。CBSregist2015.dta
Edelman，Benjamin，Michael Luca 和 Dan Svirsky。 （2017）。共享经济中的种族歧视：来自实地实验的证据。《美国经济学杂志：应用经济学》，9（2）：1-22。Airbnb.csv
Fey，Mark，Richard D. McKelvey 和 Thomas R. Palfrey。 （1996）。恒定总和蜈蚣游戏的实验研究。《国际博弈理论杂志》，25（3）：269-87。
弗莱尔，罗兰 G.，朱尼尔和史蒂文 D.莱维特。 （2004）。独特的黑人名字的原因和后果。《季度经济学杂志》119（3）：767-805。
Imbens，G.，& Kalyanaraman，K.（2012）。回归断点估计器的最佳带宽选择。《经济研究评论》，79（3），933-959。
Jha，S.和 Shayo，M.（2019）。和平的价值：金融市场暴露对选票和政治态度的影响。计量经济学，87：1561-1588。replicationdata.dta
凯恩斯，约翰梅纳德。 （1936）。就业，利息和货币的一般理论。美国哈考特，布雷斯和公司出版，并由美国纽约波利图形公司印刷。
梅耶森，埃里克。 （2014）。伊斯兰统治与穷人和虔诚者的赋权。计量经济学，82（1），229-269。regdata0.dta
Montesquieu，C.（1748）。法的精神。伦敦：T. Evans，1777 年，4 卷。第 2 卷。在线自由图书馆。
Oreopoulos, Philip. (2011). 为什么技术移民在劳动市场上遇到困难？一项涉及一万三千份简历的实地实验. 美国经济学杂志：经济政策, 3 (4): 148-71. oreopoulos.dta
Palacios-Huerta, Ignacio, and Oscar Volij. (2009). 野外蜈蚣. 美国经济评论, 99 (4): 1619-35. Chess.xls
Richardson, Gary, and William Troost. (2009). 货币干预缓解了大萧条期间的银行恐慌：来自联邦储备区域边界的准实验证据，1929-1933. 政治经济学杂志 117 (6): 1031-73。
Weber, Max. (1930). 《新教伦理与资本主义精神》。纽约：斯克里布纳，（原始出版于 1905 年）。
Ziebarth, Nicolas L. (2013). 从密西西比州大萧条期间的自然实验中识别银行倒闭的影响. 美国经济学杂志：宏观经济学, 5 (1): 81-101. MS_data_all_years_regs.dta
浙公网安备 33010602011771号