17 事件研究法
17.1 它是如何工作的？
17.1.1 古老如时间本身的研究方法
事件研究可能是最古老、最简单的因果推断研究设计。它早于实验方法，早于统计学，甚至可能早于人类语言的出现，或许比人类本身的历史还要悠久。其核心思想是：在某个特定时间点发生了一个事件，导致某种处理（干预）在该时间点被实施。通过比较事件发生前后的变化，即可得出该处理的效果。与\(第16章\)介绍的固定效应方法类似，事件研究仅利用”内部变异\(（within \space variation）\)“进行分析。但不同之处在于，事件研究专注于处理从”关闭”到”开启”的转变，且（通常）不使用面板数据，而是采用时间序列数据——即只追踪单个”个体”在多个时期的变化。
所以，当某人在深夜喝了一杯啤酒后立刻睡着，并在第二天得出结论认为啤酒会让人犯困时——这就是一个事件研究。当你在门上挂了一个”禁止推销”的牌子，之后发现推销的人变少了，并得出结论认为牌子有效时——这也是一个事件研究。当你的狗饿了，找到你并呜呜叫，随后被喂饱，于是认为呜呜叫能换来食物时——这同样是一个事件研究。当一只公鸡认为太阳升起是它的功劳，因为每天早晨它一打鸣太阳就升起时——这还是一个事件研究。
事件研究的设计适用于因果关系图如 \(图17.1\) 所示的情况。在某个特定时间段，我们从”事件前”进入”事件后”，此时某项处理开始生效1。此处的”后”并非指”处理来了又走”，而是指”事件发生后的时间段”，因此处理可能仍在持续生效。通过比较处理实施后的结果与未实施时的结果，我们希望能估计出”\(处理\rightarrow结果\)“的效应。
当然，这里还存在一个需要解决的后门路径：
\[ \text{Treatment} \leftarrow \text{AfterEvnet} \leftarrow \text{Time} \rightarrow \text{Outcome} \]
此处”时间”变量实质上是”所有随时间变化因素”的集合变量。我们之所以能判定公鸡打鸣与日出无关，是因为太阳无论如何都会升起2。要使事件研究有效，必须通过某种方式阻断后门路径。实际上，我们无法直接控制”事件后”变量——这样做等同于同时控制了”处理”变量，因为二者在此情境下具有完全共线性：要么处于”事件前+未处理”状态，要么处于”事件后+已处理”状态。图中将二者分离仅为表述清晰。但由于”时间”仅通过”事件后”影响”处理”，理论上可以通过控制”时间”的某些成分来阻断后门路径，同时保留”事件后”的有效变异。
能否严谨分析 \(图17.1\) 的因果结构，正是区分优质与劣质事件研究的关键。此刻您脑海中可能已浮现出劣质案例——事实上，本书中诸多”明显错误的因果推断”示例，本质上都是失败的事件研究。许多迷信观念同样源于此类缺陷：
友人穿着绿色运动裤后，其支持球队获胜，便认定裤子带来好运
电视剧因收视率连年下滑引入新角色，不久后停播，观众归咎于该角色
这些本质上都是事件研究（尽管是错误的设计）。
但这些都属于设计不当的事件研究！严谨分析因果图意味着需要做到两点：
深入思考该因果图是否适用于当前场景——从”事件前”到”事件后”，处理变量是否是唯一发生变化的因素？是否存在其他同时发生的混淆变量可能导致观测结果？事件研究本身并非缺陷方法，但其适用情境具有特定性：您的研究设定是否符合这些条件？
周密考虑如何阻断通过时间变量的后门路径
电视剧被取消的粉丝们在 \(第2步\)上存在疏漏——他们本应控制”收视率持续下滑”这一趋势变量，才能准确评估新角色的影响。而那位穿运动裤的体育迷的问题更多源于 估计方法缺陷 而非因果图设计——在仅有单一重复观测值的情况下，如何正确处理 抽样变异 本身就是一个方法论难题。
17.1.2 预测与偏离
事件研究真正棘手之处在于 后门路径 的干扰。该研究设计的核心逻辑依赖于时间维度上的单一变化——事件发生导致处理生效，进而通过对比事件前后差异来识别处理效应。但这一方法的 关键前提 是：处理变量必须是唯一发生变化的因素。若结果变量本身会因其他未控制的干扰因素产生趋势性变化，则该事件研究将产生严重偏误。
这正是方法论的核心难点——时间维度上永远存在多重变化。问题的关键在于：如何确保在控制所有时间相关干扰因素的同时，精准保留”事件前后临界点”这一关键变异来源？
简言之，我们必须构建反事实预测机制。利用事件前的丰富数据，在满足趋势延续性假设（即若无处理干预，事件前趋势将保持恒定）的前提下，可通过以下步骤识别处理效应：
基于事件前数据预测 无处理状态下 的预期结果
计算实际观测值与预测值的 系统性偏离
将该偏离量确认为 处理效应估计量
例如，参见\(图17.2\)。此处展示了两组时序数据图。每幅图中，黑色实线 代表实际观测到的时间序列数据。在事件发生日期之前，数据序列均呈现平稳走势；而当事件发生时，时序路径出现显著偏离——\(图17.2（a）\)表现为陡然上跃，\(图17.2（b）\)则呈现断崖式下跌。
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(extrafont)## Registering fonts with Rlibrary(ggpubr)
library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:ggpubr':
##
## get_legend
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(lubridate)
library(tidyquant)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## ── Attaching core tidyquant packages ─────────────────────── tidyquant 1.0.11 ──
## ✔ PerformanceAnalytics 2.0.8 ✔ TTR 0.24.4
## ✔ quantmod 0.4.28 ✔ xts 0.14.1── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
## ✖ zoo::as.Date() masks base::as.Date()
## ✖ zoo::as.Date.numeric() masks base::as.Date.numeric()
## ✖ dplyr::filter() masks stats::filter()
## ✖ xts::first() masks dplyr::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ xts::last() masks dplyr::last()
## ✖ PerformanceAnalytics::legend() masks graphics::legend()
## ✖ cowplot::stamp() masks lubridate::stamp()
## ✖ quantmod::summary() masks base::summary()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsd <- list()
set.seed(1000)
d[[1]] <- tibble(Date = ymd('2018-01-01') + days(0:199), y = c(runif(150,0,.5), runif(50,0,.5) + 1/sqrt(1:50)),
yalt = c(rep(NA_real_,150),runif(50,0,.5))) %>%
pivot_longer(cols = c('y','yalt')) %>%
mutate(name = factor(case_when(name == 'y' ~ 'What We See', name == 'yalt' ~ 'What We Expect'),
levels = c('What We See','What We Expect')))
p1 <- ggplot(na.omit(d[[1]]), aes(x = Date, y = value, color = name)) +
geom_line() +
geom_vline(aes(xintercept = ymd('2018-05-30')), linetype = 'dashed') +
scale_color_manual(values = c('black','gray')) +
theme_pubr() +
labs(y = 'Outcome', color = NULL, title = '(a) Flat Beforehand') +
theme(text = element_text(family = 'sans', size = 10),
legend.position = c(.3, .8),
legend.background = element_rect(color = 'black'),
legend.text = element_text(size = 10))## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.d[[2]] <- tibble(Date = ymd('2018-01-01') + days(0:199),y = c(rnorm(150,0,.5) + (1:150)/75, rnorm(50,0,.5) + (151:200)/75 - 1),
yalt = c(rep(NA_real_,150),rnorm(50,0,.5) +(151:200)/75)) %>%
pivot_longer(cols = c('y','yalt')) %>%
mutate(name = factor(case_when(name == 'y' ~ 'What We See',name == 'yalt' ~ 'What We Expect'),
levels = c('What We See','What We Expect')))
p2 <- ggplot(na.omit(d[[2]]), aes(x = Date, y = value, color = name)) +
geom_line() +
geom_vline(aes(xintercept = ymd('2018-05-30')), linetype = 'dashed') +
scale_color_manual(values = c('black','gray')) +
theme_pubr() +
guides(color = 'none') +
labs(y = 'Outcome',title = '(b) Prior Trend') +
theme(text = element_text(family = 'sans', size = 10),
legend.position = c(.2, .8),
legend.background = element_rect(color = 'black'),
legend.text = element_text(size = 10))
plot_grid(p1,p2)
ggsave('EventStudies/event_study_examples.pdf', width = 8, height = 3.5,device=cairo_pdf)要使事件研究有效，我们必须合理推断若无事件发生时的潜在结果。通常我们假设事件前的趋势模式将会延续：以 \(图17.2（a）\)为例，虽然时序数据每日波动，但既无整体上升/下降趋势，也无结构性变化——若事件未发生，该模式很可能持续。据此，灰色虚线 表示趋势延续的预测值，而事件研究估计量即通过对比事件后的黑色实线（实际数据）与灰色虚线（反事实预测）得出。结果显示：处理实施后黑色线立即显著高于灰色线，表明处理存在短期正向效应；但随着时间推移，两条线逐渐收敛，说明效应具有衰减特征。
该原理同样适用于\(图17.2（b）\)。此处事件前存在 显著上升趋势，故假设趋势将持续具有合理性。因此，灰色虚线 沿原上升趋势延伸，代表预期发展路径；而 黑色实线 的实际数据却陡然下降。此时我们可估计出处理对结果存在持续负向效应——与 \(图17.2（b）\)不同，该效应呈现持久性特征：直至数据期末端，黑色线始终稳定低于灰色线。
那么，如何构建这些预测？虽然确定图中灰色反事实预测线涉及大量技术细节与估计方法，但从 方法论框架 层面而言，主要存在三种核心路径：
方法1：完全忽略。没错！您可以直接忽略时间趋势变化，仅简单对比事件前后数据。这竟是实践中惊人常见的做法——但显然，这是个糟糕的方法，不是吗？
该方法的有效性本质上取决于 研究情境。诚然，“完全忽略法”的某些应用纯属一厢情愿；但若确有理论依据表明时间效应可忽略不计，此时采用该方法反而构成合理的策略性简化。
该方法在两种情境下具有合理性：
情境一：当数据特征近似 \(图17.2（a）\)3——事件前时序数据呈现绝对平稳性（无趋势性变化、无结构变动、无异常波动），仅围绕固定值小幅随机波动（且信噪比理想）。若同时能确认事件发生时 除处理变量外无其他干扰因素，则可合理推定原数据模式将持续延续。
情境二：分析极微观时间尺度的数据。该思路常见于拥有高频数据的金融研究（如分钟级/秒级序列）。虽然公司股价长期可能呈现趋势性或异常波动，但当分析窗口压缩至60秒级时，此类宏观趋势的影响可被有效忽略。
方法2：基于事件前数据的预测建模。该方法通过分析事件前的结果变量数据，提取其内在模式以预测事件后的潜在结果。如 \(图17.2（a）\)与\(（b）\)所示：
当数据存在 明显趋势特征 时，可通过 趋势外推法 获得反事实预测值
最终通过 预测值与实际值的系统性偏离 识别处理效应
原理简明——至少概念层面如此。实际操作中，此类预测存在多种建模方法，其效果参差不齐。本章后续将展开具体讨论。
方法3：基于事件后数据的预测建模。该方法与方法2高度相似：如同方法2，您首先利用事件前数据确立数据的内在规律模式，随后将该模式外推至事件后阶段，从而获得预测结果。
核心差异在于：本方法不仅捕捉结果变量的趋势特征，还需识别事件前期 变量与其他协变量的关联模式。在事件后期，通过观测这些协变量的变化轨迹，可进一步增强预测的稳健性。
该方法在 股市分析 中应用典型：研究者对\(X\)股票进行事件研究时，首先基于事件前数据建立\(X\)股价与大盘指数的关联模型。例如，可能发现市场指数每上涨\(1\%\)，\(X\)股价相应上涨\(0.5\%\)的量化关系。
在事件后期，研究者除延续方法2对\(X\)股价自身趋势的外推外，还需监测大盘指数变动：若事件后大盘指数上涨 \(5\%\)，则根据前期模型应预期\(X\)股价自然上涨\(2.5\%\)（与处理效应无关）。最终，需从实际股价涨幅中 扣除这\(2.5\%\)的市场关联成分，方可分离出纯净的处理效应。
这三种方法共同揭示了一个关键原则：事件研究通常仅适用于 短期后事件窗口（事件研究术语称为”短视界”）。随着时间推移，基于事件前的预测效力会持续衰减，时间后门路径 的干扰风险将显著增加——即便已实施控制变量。因此，除非满足以下严格条件之一：
能确证长期中除处理变量外无其他干扰因素
具备长期趋势变化的有效建模能力
Sagar P Kothari and Jerold B Warner. Econometrics of event studies. In Handbook of Empirical Corporate Finance, pages 3–36. Elsevier, 2007.
否则应保持短视界分析。虽然长视界分析在技术上可行，但此时时间序列设计提供的因果识别力大幅削弱，研究者必须：
强化研究设计效度
严密监控无关趋势干扰
持续排除后门路径与噪声干扰
核心要旨大抵如此。尽管具体操作远比上述复杂，但事件研究能成为古老而持久的方法，正源于其概念层面的简洁性。真正的复杂性出现在估计环节——我将在”操作实施”章节详细展开。
library(tidyverse)
library(lubridate)sp500 <- read_csv('EventStudies/SP500_Historical.csv',show_col_types = FALSE) %>%
mutate(Date = mdy(Date)) %>%
arrange(Date) %>%
mutate(SP500_Return = `Close/Last`/lag(`Close/Last`) - 1)
# 读取数据并处理
both <- read_csv('EventStudies/GOOG.csv', show_col_types = FALSE) %>%
dplyr::filter(Date >= '2015-05-01' & Date <= '2015-08-31') %>%
# 转换时间戳为日期
mutate(Date = as.Date(Date)) %>%
# 计算收益率
mutate(Google_Return = as.double(Close) / lag(as.double(Close)) - 1) %>%
# 选择需要的列
select(Date, Google_Return) %>%
# 合并SP500数据
inner_join(sp500 %>% select(Date, SP500_Return),by = join_by(Date)) %>%
# 移除第一行(因为有NA收益率)
slice(-1)
# 保存结果
write_csv(both, 'EventStudies/google_stock_data.csv')
head(both)## # A tibble: 6 × 3
## Date Google_Return SP500_Return
## <date> <dbl> <dbl>
## 1 2015-05-04 0.00535 0.00294
## 2 2015-05-05 -0.0185 -0.0118
## 3 2015-05-06 -0.0124 -0.00446
## 4 2015-05-07 0.0124 0.00377
## 5 2015-05-08 0.0142 0.0135
## 6 2015-05-11 -0.00468 -0.0050917.1.3 术语说明
本章将探讨我称之为“事件研究”的方法论。然而，正如对这种直观方法所预期的那样，不同领域独立发展出了各自的研究范式，导致术语体系存在差异——有时甚至伴随细微的内涵变异。
例如在公共卫生领域，您可能遇到”统计过程控制”方法——该方法仅适用于 无显著事件前趋势 的情境（即前文所述 \(方法1\)）。本书将此类方法及其衍生变体统称为事件研究的范式分支。
另一易混淆术语是“断点时间序列”——理论上可泛指任何考虑事件前趋势的事件研究方法，但实践中通常特指以下两类：
基于时间序列计量经济学的事件研究
分段回归法（分别为事件前/后拟合独立回归线）
需警惕：当该术语被混用时，两类方法的应用者往往产生概念恐慌4。
事件研究还存在控制组变体——通过引入未受事件影响的对照组来消除时间因素干扰。本章虽不展开讨论该设计5，但其核心思想（即“事件前后比较+控制组”的双重识别策略）将在后续章节系统阐述：
双重差分法（第\(18\)章）
合成控制法（第\(21\)章）
需注意的是，实践中标榜”含对照组的事件研究”或”含控制组的断点时间序列”的研究，往往比标准双重差分/合成控制方法包含更多高阶时间序列技术（参见”专业实践”章节）。
更令人困惑的是，部分经济学家将“事件研究”术语严格限定于多处理组设计（各组接受处理的时间点不同，无论是否存在对照组）。此类方法我们同样将在第\(18\)章（双重差分法）系统探讨。
当前章节仅聚焦经典前后对照式事件研究设计。
17.2 操作实施方法
探讨”事件研究的操作实施”实则存在一定概念张力——因其作为方法论框架的高度普适性，加之各学科衍生出的变体体系，导致具体操作流程呈现显著异质性。
为避免冗赘地遍历所有事件研究变体，本节将聚焦三大主流范式：
金融计量应用：通过显式计算事件后各时点的预测值，并与实际数据对比
回归框架构建6：基于参数化模型的分析路径
专业实践方法（见专门章节）：更严格处理数据的时间序列特性
17.2.1 股市事件研究法
事件研究在 股市分析 中应用广泛，尤其适用于 股票收益率 作为结果变量的场景。其理论根基在于：任何有效且流动性充足的市场中，股价理应已反映所有公开信息。尽管收益率存在日常波动，或随大盘整体上涨，但股价出现 剧烈变动的唯一诱因 应是公司相关 未预期信息 的释放——这正是”事件”的本质所在。
股价的 信息即时反映 特性带来双重优势：
短期窗口可行性：支持“无同期干扰”假设；精准定位事件时点7
预期效应捕捉：
以企业并购公告为例，虽实际合并可能滞后数月/年，但投资者基于 预期交易 的买卖行为会立即使股价反应，从而避免长期观测导致的 时间后门路径 干扰。
标准事件研究流程包含以下步骤：
步骤1：时段划分，估计窗口：选取事件前足够长的清洁期；观测窗口：覆盖事件临界点及后续效应持续期
步骤2：预期收益率建模
使用估计窗口数据构建收益率预测模型\(\hat R\)，常用方法包括：
均值调整模型：直接采用估计窗口内 股票收益率 的算术平均值：\(\hat R=\bar R\)（适用于无系统性风险暴露的标的）
市场调整模型：以 市场收益率 作为预测基准：\(\hat R=R_m\)（隐含完全市场同步假设）
风险调整模型：估计阶段，通过 市场模型回归 量化风险暴露：\(R=\alpha+\beta R_m+\epsilon\)；预测阶段，利用实际\(R_m\)计算预期收益率：\(\hat R=\hat \alpha+\hat\beta R_m\)（需满足市场模型稳定性假设8）
步骤3：异常收益率计算，在观测窗口内，通过实际收益率与预测值的差值计算异常收益率\(（Abnormal \space Return, AR）\)：\(AR=R-\hat R\)
步骤4：在观测期内分析异常收益率\(AR\)：事件发生前出现非零\(AR\)值表明市场已预期到该事件；事件后出现的非零\(AR\)值则反映了该事件对股票收益的实际影响。正是基于同样的市场有效性原理（这也是股票收益率适合用于事件研究的根本原因），事件效应通常会在短期内快速达到峰值，随后逐渐消退。
让我们通过一个实例来说明这个方法。金融领域事件研究的一个优势在于其广泛适用性——几乎任何场景都能应用。因此，我不打算引用已发表论文中的案例，而是自己构建一个研究示例。
\(2015年8月10日\)，谷歌宣布将重组公司架构：原有”谷歌”作为母公司控股\(FitBit\)、\(Nest\)等子公司的模式将被取代，转而成立名为”\(Alphabet\)“的新控股公司，统一管理谷歌及其他业务板块。
市场对此作何反应？为探明这一点，我收集了 \(2015年5月至8月底\) 期间的数据：
标的股票：\(GOOG\)股价数据（该代码在重组前代表\(Google\)，重组后代表\(Alphabet\)）
市场基准：标普\(500\)指数价格（作为市场整体表现的参照）
随后，我分别计算了\(GOOG\)股票和标普\(500\)指数的日收益率（计算公式：\(当日价格/前1日价格-1\)）。为何使用收益率而非原始价格？主要原因包括：
信息反应隔离：收益率能直接反映新信息对价格的边际影响
跨股票可比性：消除绝对价格差异带来的比较偏差
library(ggpubr)
library(tidyverse); library(lubridate)
goog <- causaldata::google_stock
event <- ymd("2015-08-10")
# Create estimation data set
est_data <- goog %>% filter(Date >= ymd('2015-05-01') & Date <= ymd('2015-07-31'))
# And observation data
obs_data <- goog %>% filter(Date >= event - days(4) & Date <= event + days(14))
# Estimate a model predicting stock price with market return
m <- lm(Google_Return ~ SP500_Return, data = est_data)
# Get AR
obs_data <- obs_data %>%
# Using mean of estimation return
mutate(AR_mean = Google_Return - mean(est_data$Google_Return),
# Then comparing to market return
AR_market = Google_Return - SP500_Return,
# Then using model fit with estimation data
risk_predict = predict(m, newdata = obs_data),
AR_risk = Google_Return - risk_predict)# Graph the results
ggplot(obs_data, aes(x = Date, y = AR_risk)) +
geom_line() +
geom_vline(aes(xintercept = event), linetype = 'dashed') +
geom_hline(aes(yintercept = 0))
rdata <- goog %>%
rename(Google = Google_Return, `S&P 500` = SP500_Return) %>%
pivot_longer(cols = c('Google','S&P 500'))
ggplot(rdata, aes(x = Date, y = value, linetype = name)) +
geom_line() +
geom_text(data = rdata %>% filter(Date == max(Date)),
mapping = aes(x = Date, y = value, label = name),
hjust = -.2,
family = 'sans', size = 14/.pt) +
geom_vline(aes(xintercept = event), linetype = 'dashed') +
geom_hline(aes(yintercept = 0)) +
guides(linetype = FALSE) +
scale_x_date(limits = c(min(rdata$Date),max(rdata$Date) + days(18))) +
scale_y_continuous(labels = function(x) scales::percent(x, accuracy = 1)) +
labs(x = 'Date', y = 'Abnormal Return') +
theme_pubr() +
theme(text = element_text(family = 'sans', size = 13))## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave('FixedEffects/google_sp500_data.pdf', width = 6, height = 5,device=cairo_pdf)从 \(图17.3\) 的日收益率曲线可以观察到两个关键特征：
市场联动性：谷歌收益率与标普\(500\)指数收益率呈现同步波动，这正是三种预测方法中两种需引入市场指数作为参照的原因——需区分个股异动与市场整体波动。
常态平稳性：收益率多数时间保持平稳（如\(7\)月中旬的显著峰值属例外），这种”事件驱动型突变+常态平稳”的模式正是事件研究的基础前提。
library(tidyverse)
library(lubridate)
goog <- read_csv('EventStudies/google_stock_data.csv', show_col_types = FALSE)
event <- ymd("2015-08-10")
# Create estimation data set
est_data <- goog %>% filter(Date >= ymd('2015-05-01') & Date <= ymd('2015-07-31'))
# And observation data
obs_data <- goog %>% filter(Date >= event - days(4) & Date <= event + days(14))
# Estimate a model predicting stock price with market return
m <- lm(Google_Return ~ SP500_Return, data = est_data)
# Get AR
obs_data <- obs_data %>%
# Using mean of estimation return
mutate(AR_mean = Google_Return - mean(est_data$Google_Return),
# Then comparing to market return
AR_market = Google_Return - SP500_Return,
# Then using model fit with estimation data
risk_predict = predict(m, newdata = obs_data),
AR_risk = Google_Return - risk_predict)
goog <- goog %>% mutate(AR_risk = Google_Return - predict(m, newdata = goog))
sd(goog$AR_risk)## [1] 0.02050506让我们开始进行事件研究。首先需要确定估计窗口和观测窗口：估计窗口选择 \(5月至7月\)9，观测窗口则从公告日前几天\(（8月6日）\)开始，持续至\(（8月24日）\)。
接下来，我们将利用观测窗口内的数据构建预测模型并进行比较。具体实现代码如下：
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import seaborn.objects as so
import matplotlib.pyplot as plt
from causaldata import google_stock
goog = google_stock.load_pandas().data
# Create estimation data set
goog['Date'] = pd.to_datetime(goog['Date'])
est_data = goog.loc[(goog['Date'] >= '2015-05-01') & (goog['Date'] < '2015-07-31')]
# And observation data
obs_data = goog.loc[(goog['Date'] >= '2015-08-06') & (goog['Date'] < '2015-08-24')]
# Estimate a model predicting stock price with market return
m = smf.ols('Google_Return ~ SP500_Return', data = est_data).fit()
# Get AR
goog_return = np.mean(est_data['Google_Return'])
obs_data = (obs_data
# Using mean of estimation return
.assign(AR_mean = lambda x: x['Google_Return'] - goog_return,
# Then comparing to market return
AR_market = lambda x: x['Google_Return'] - x['SP500_Return'],
# Then using model fit with estimation data
risk_pred = lambda x: m.predict(x),
AR_risk = lambda x: x['Google_Return'] - x['risk_pred']))
# Graph the results. The fig parts are necessary to add ref lines
fig = plt.figure()
so.Plot(obs_data, x = 'Date', y = 'AR_risk').add(so.Line()).on(fig).plot()## <seaborn._core.plot.Plotter object at 0x0000021C82D15BE0>fig.axes[0].axvline(pd.Timestamp(2015,8,10), linestyle = '--')
fig.axes[0].axhline(0)
causaldata google_stock.dta, use clear download
* Get average return in estimation window and calculate AR
summ google_return if date <= date("2015-07-31","YMD")
g AR_mean = google_return - r(mean)
* Get difference between Google and the market
g AR_market = google_return - sp500_return
* Estimate a model predicting stock price with market return
reg google_return sp500_return if date <= date("2015-07-31","YMD")
* Get the market prediction for AR_risk
predict risk_predict
g AR_risk = google_return - risk_predict
* Graph the results in the observation window
g obs = date >= date("2015-08-06", "YMD") & date <= date("2015-08-24", "YMD")
* Format the date nicely for graphing
format date %td
* Use the tsline function for time series plots
tsset date
tsline AR_risk if obs, tline(10aug2015)将\(3\)条预测线绘制在同一图表中，得到 \(图17.4\)。我们可以观察到：不论使用哪种方法，\(8月10日\) 后都出现了明显的短期峰值，但这个峰值仅持续了\(1\)天。正如有效市场理论所预期的，谷歌（现为\(Alphabet\)）的股价迅速调整至包含重组信息的新价位，随后日收益率回落至\(0\)附近（尽管股价本身仍高于事件前水平）。此外，事件前并未出现显著的市场异动，这表明此次重组并未被市场充分预期。
library(ggplot2)
gdata <- obs_data %>%
rename(Mean = AR_mean, Market = AR_market, Risk = AR_risk) %>%
pivot_longer(cols = c('Mean','Market','Risk'))
ggplot(gdata, aes(x = Date, y = value, linetype = name)) +
geom_line() +
geom_text(data = gdata %>% filter(Date == max(Date)),
mapping = aes(x = Date, y = value, label = name),
hjust = -.2,
family = 'sans',
size = 10/.pt) +
geom_vline(aes(xintercept = event), linetype = 'dashed') +
geom_hline(aes(yintercept = 0)) +
guides(linetype = FALSE) +
scale_x_date(limits = c(min(gdata$Date),max(gdata$Date) + days(3))) +
scale_y_continuous(labels = function(x) scales::percent(x, accuracy = 1)) +
labs(x = 'Date', y = 'Abnormal Return') +
theme_pubr() +
theme(text = element_text(family = 'sans', size = 10))
ggsave('EventStudies/google_stock_event.pdf', width = 6, height = 5,device=cairo_pdf)观察 \(图17.4\) 右侧数据，我们可以发现该方法论的几个特点，这些提醒我们需要保持警惕：
首先，注意 \(8月20日\) 左右，“均值调整线”出现明显下跌，而”市场调整线”和”风险调整线”仍维持在\(0\)值附近？这表明谷歌股价确实下跌，但市场整体也在下跌。这很可能与\(Alphabet\)重组无关，而是整个市场的波动。这正是我们需要采用考虑 市场整体走势 的方法的重要原因。
\(图17.4\) 中 \(8月20日左右\) 的波动数据提醒我们：必须严格控制事件窗口期。\(Alphabet\)重组仅产生单日效应便立即消退，因此一个半星期后的股价变动很可能与重组无关——但这些波动仍被误判为事件效应！在瞬息万变的股票市场中，一个半星期足以发生其他重大变化，若采用过长的观测窗口，这些无关波动就会被错误归因为重组效应。
获得日度效应数据后，我们可进行以下分析：我们首先可以 通过对特定日期的异常收益率\(AR\)进行显著性检验，来判断该日是否存在显著效应。最常用的检验方法非常简单：只需计算整个时间窗口（而不仅是观测窗口）的\(AR\)值，然后求得这些\(AR\)的标准差——这就是标准误。接下来，用目标日期的\(AR\)值除以该标准误，得到\(t\)统计量，最后检查该\(t\)值是否超过临界值（通常 \(95\%\) 显著性水平下为 \(\pm1.96\)）即可完成检验。根据计算，\(GOOG\)股票在整个时间窗口内的风险调整异常收益率\(AR_{risk}\)标准差为 \(0.021\)。而公告次日的\(AR_{risk}\)达到 \(0.054\)，由此得出的\(t\)统计量为 \(0.054/0.021=2.571\)。该值显著高于\(1.96\)的临界值（\(95\%\)置信水平），因此我们可以得出结论：\(Alphabet\)重组公告确实在当天显著提升了\(GOOG\)股票的收益率。
我们可能还关注事件对日均异常收益率或累计异常收益率的影响。计算方法非常直接：
日均异常收益率：观测窗口（事件后）内所有\(AR\)的算术平均值
累计异常收益率：观测窗口内所有\(AR\)的代数和
不过，要为这些综合指标计算标准误或进行显著性检验，则需要更复杂的方法。针对这些聚合指标（如日均/累计异常收益率）的统计检验存在多种不同方法，其核心思路通常是：通过计算多只股票的平均或累计异常收益率，并利用股票间的横截面标准差来构建检验统计量。
最基础的是”横截面检验”方法：先计算每家公司平均或累计异常收益率\(AR\)，再分别用这些均值或累计值的标准差作为分母，乘以公司数量的平方根，即可构建\(t\)统计量——操作十分简便。
当然，这种简单的方法忽略了各种各样的因素 — 时间序列相关性、波动性增加、企业之间的相关性等等。正如我提到的，有无数种替代方法。其中包括许多源自\(Patell（1976）\)的检验方法。该检验考虑了收益率在观察窗口中出现得早或晚的可变性。调整后\(Patell\)还考虑了企业之间收益率的相关性。还有标准化横截面检验，它也考虑了时间序列的序列相关性以及事件增加波动性的可能性。而调整后的标准化横截面检验则增加了考虑企业之间相关性的能力。而且这还没涉及基于异常回报排名的非参数检验呢。哎呀，选项太多了。
Patell, James M. 1976. “Corporate Forecasts of Earnings Per Share and Stock Price Behavior: Empirical Test.” Journal of Accounting Research 14 (2): 246–76.这些额外的检验可以在 \(R\) 语言的 \(estudy2\) 包中实现，也可以在 \(Stata\) 的 \(eventstudy2\) 包中实现。在 \(Python\) 中，你可能得自己想办法（自己动手实现）。
17.2.2 基于回归分析的事件研究法
金融领域的方法专为 短期脉冲效应（快速冲高后消退）设计，这在该领域十分常见。但其他学科的研究者若关注 持久性结构变化 时该如何处理？此时，基于回归框架的事件研究设计将成为更合适的工具。
基于回归的事件研究核心思路
预处理模型：使用事件前数据建立基准回归模型
后处理模型：使用事件后数据建立对比回归模型
效应评估：通过模型参数差异识别事件影响
该方法可通过简单的交互项（事件研究术语中称为”分段回归”）实现：
\[ \text{Outcome} = \beta_0 + \beta_1 t + \beta_2 \text{After} + \beta_3 \times t \times \text{After} + \varepsilon \quad (17.1) \]
其中，\(t\)代表时间变量，\(After\)为二元虚拟变量（事件后取值为\(1\)）。这当然只是一个基础的线性模型设定，强制约束了时间趋势的线性特征。如需更灵活的形状拟合，可参照第\(13\)章讨论的方法，通过添加多项式项来扩展模型。
这种方法具有明显的优点和局限性。其优势在于，相较于逐日（或逐秒等高频数据）分析，它能提供更精确的时间趋势估计。但缺点在于难以捕捉效应的具体形态——比如前文提到的谷歌股价单日跳涨这类精细变化。此外，若对整体时间趋势的形态设定错误，将导致严重偏误。不过，若能通过添加适当的多项式项准确刻画事件两侧的时间序列特征，就能有效解决这一问题。
这种方法还存在一个重大缺陷：在进行统计显著性检验时需要格外谨慎。若数据存在自相关性（即每期数值依赖于前期数值），即使事件实际上毫无影响，也极可能得出统计显著的结果10。这是因为时间序列数据具有”粘性”特征——相邻时点的数值会持续聚类，导致回归模型对趋势判断的置信度虚高。在自相关性较强的数据中，随机选取一个虚假”事件日”进行分段回归检验时，约半数情况会显示统计显著性（远超过真实零效应下\(95\%\)置信水平应有的\(5\%\)误判率）。
不过，自相关问题仍有解决之道。使用 \(第13章\) 讨论的异方差-自相关稳健标准误 \(HAC\) 能有效改善推断结果。更优的方法是采用能够直接刻画自相关结构的模型，具体技术细节可参考本章”专业操作”部分。
分段回归方法无需特殊代码，只需将结果变量对时间变量的多项式（含\(After\)虚拟变量的交互项）进行回归即可。由于处理的是时间序列数据，务必使用异方差-自相关稳健标准误 \(HAC\)。关于交互项和 \(HAC\) 标准误的具体说明，请参阅 \(第13章\) 内容。
让我们通过实际案例来演示该方法的应用。\(Taljaard \space et.al（2014）\)采用基于回归的事件研究法，评估了政策干预对健康结果的影响。具体而言，他们研究了英国\(2010年\)中期实施的一项政策——该政策旨在通过提升救护车转运期间的医疗质量，来降低患者心脏病发作和中风的风险。政策内容包括组建质量改进团队、促进多方协作、共享创新方案以及加强医护人员培训。那么，这项政策是否真的改善了心脏病救治效果？
Monica Taljaard,Joanne E McKenzie,Craig R Ramsay, and Jeremy M Grimshaw. The use of segmented regression in analysing interrupted time series studies: An example in pre-hospital ambulance care. Implementation Science, 9(1):1–4,2014.\(Taljaard \space et \space al.（2014）\)的研究采用以下回归模型设定：以心脏病救治绩效（\(AMI\)，急性心肌梗死救治指标）为因变量，自变量包括：
时间基准变量：\(Week-27\)（以政策实施\(第27周\)为原点进行中心化处理，使截距项直接反映政策突变效应）
政策虚拟变量：\(After\)（政策实施后时期取值为\(1\)）
交互项：\((Week-27) \times After\)11
\[ \text{AMI} = \beta_0 + \beta_1(\text{Week} - 27) + \beta_2 \text{After} + \beta_3 (\text{Week} - 27) \times \text{After} + \varepsilon \quad (17.2) \]
\(Taljaard \space et \space al.\) 针对心脏病救治的研究结果可参见\(图17.5\)。如图所示，事件起始期\(（第27周）\)左右两侧分别拟合了两条回归线，这正是交互项作用的直观体现：
左侧线段（政策前）：\(\beta_0 + \beta_1(Week-27)\)
右侧线段（政策后）： \((\beta_0 + \beta_2) + (\beta_1 + \beta_3)(Week-27)\)
从分析结果中可以得出什么结论？首先，在政策临界点\(（第27周）\)并未观察到任何跳跃式变化，表明政策 未产生即时效果——考虑到此类需要多方协作的行为干预政策本就需时见效，这一结果并不意外。政策效应量可通过\(第27周\)时的预测值变化来衡量：当\(Week=27\)时，左侧线段预测值为\(\beta_0\)，右侧线段预测值为 \(\beta_0+\beta_2\)（此时\(Week-27=0\)，两条线中的\(\beta_1\)和\(\beta_3\)项均被消去）。
两者的差值为：\(\beta_0 + \beta_2 - \beta_0 = \beta_2\)。\(\beta_2\) 即代表政策生效时点的预测值跃升幅度，也就是我们所要估计的事件研究效应量。
# Ambulance study
hosp <- read.csv(text = 'X, Y
0.9411764705882355, 0.48863636363636365
1.8235294117647047, 0.4545454545454546
2.882352941176471, 0.4386363636363636
3.9411764705882355, 0.38181818181818183
5, 0.42272727272727273
5.882352941176469, 0.4113636363636364
6.9411764705882355, 0.15227272727272734
8, 0.38636363636363646
8.882352941176471, 0.34090909090909094
9.941176470588236, 0.42727272727272736
10.823529411764707, 0.4136363636363637
11.882352941176471, 0.3772727272727273
12.941176470588236, 0.35
14, 0.49090909090909096
14.882352941176471, 0.44090909090909103
15.941176470588236, 0.3568181818181818
17.000000000000004, 0.5318181818181817
18.058823529411768, 0.4931818181818183
19.117647058823533, 0.44999999999999996
19.823529411764707, 0.4568181818181818
21.058823529411764, 0.44772727272727275
21.941176470588236, 0.5181818181818182
23.000000000000004, 0.4181818181818182
24.058823529411764, 0.4
24.941176470588243, 0.415909090909091
25.823529411764707, 0.36363636363636365
27.05882352941177, 0.44090909090909103
27.941176470588236, 0.38636363636363646
29.000000000000004, 0.4113636363636364
30.058823529411764, 0.49090909090909096
31.117647058823533, 0.3318181818181819
32, 0.3772727272727273
32.88235294117648, 0.44772727272727275
33.94117647058823, 0.49090909090909096
35, 0.4363636363636364
36.05882352941177, 0.3545454545454545
37.117647058823536, 0.42499999999999993
38, 0.4795454545454546
38.882352941176464, 0.46818181818181825
39.94117647058823, 0.42954545454545456
41, 0.5227272727272727
41.88235294117648, 0.4977272727272727
43.117647058823536, 0.4772727272727273
44, 0.5545454545454545
44.88235294117648, 0.5636363636363637
45.941176470588246, 0.55
47, 0.5340909090909091
48.05882352941177, 0.5590909090909091
48.941176470588246, 0.5977272727272728
50.176470588235304, 0.5568181818181819
50.882352941176464, 0.4136363636363637
51.94117647058823, 0.4636363636363636
53, 0.6136363636363636
54.05882352941177, 0.625
55.11764705882352, 0.5545454545454545
56, 0.6022727272727273
57.05882352941177, 0.515909090909091
58.11764705882352, 0.5840909090909091
59, 0.5136363636363637
60.05882352941177, 0.6204545454545455
61.117647058823536, 0.5636363636363637
62, 0.6590909090909092
63.05882352941177, 0.6136363636363636
64.11764705882354, 0.5636363636363637
65, 0.5977272727272728
66.05882352941177, 0.5863636363636364
67.11764705882354, 0.5818181818181818
68.1764705882353, 0.5681818181818181
69.05882352941177, 0.5363636363636364
69.94117647058823, 0.5545454545454545
71, 0.6227272727272728
72.05882352941177, 0.5636363636363637
73.11764705882354, 0.5613636363636363
74, 0.6204545454545455
75.05882352941177, 0.6386363636363637
75.94117647058825, 0.6000000000000001
77, 0.625
78.05882352941177, 0.6590909090909092
79.11764705882354, 0.6272727272727272
80.1764705882353, 0.665909090909091
80.88235294117646, 0.6090909090909091
82.11764705882354, 0.5863636363636364
83, 0.6045454545454545
84.23529411764707, 0.6090909090909091
85.11764705882354, 0.5954545454545455
86.1764705882353, 0.5954545454545455
87.05882352941177, 0.675
88.11764705882354, 0.6454545454545455
89, 0.6590909090909092
89.88235294117646, 0.6795454545454545
90.94117647058825, 0.65
92.17647058823529, 0.6227272727272728
93.23529411764706, 0.6022727272727273
94.11764705882352, 0.6363636363636364
95.00000000000001, 0.6318181818181818
96.05882352941177, 0.6840909090909091
96.94117647058825, 0.6363636363636364
98, 0.6954545454545455
99.05882352941178, 0.6863636363636363
100.11764705882354, 0.6409090909090909
101.00000000000001, 0.6909090909090909
102.23529411764706, 0.7159090909090909
103.11764705882352, 0.6681818181818182
104, 0.7
105.23529411764706, 0.7
106.11764705882354, 0.7568181818181818
107.17647058823529, 0.6363636363636364
108.05882352941177, 0.6977272727272728
109.11764705882352, 0.7318181818181818
110.1764705882353, 0.6704545454545454
111.05882352941177, 0.7318181818181818
111.94117647058825, 0.6795454545454545')
hosp <- hosp %>%
mutate(X = 1:112) %>%
mutate(After = X >= 27) %>%
rename(Weeks = X, AMI = Y) %>%
mutate(Weeks_Centered = Weeks - 27)
ggplot(hosp, aes(x = Weeks, y = AMI)) +
geom_line(alpha = .5, color = 'gray') + geom_point() +
geom_smooth(formula='y~x', mapping = aes(group = After), method = 'lm', color = 'black', se = FALSE) +
geom_vline(aes(xintercept = 27), linetype = 'dashed') +
scale_y_continuous(limits = c(0,1)) +
labs(y = 'AMI Performance') +
theme_pubr() +
theme(text = element_text(family = 'sans', size = 10))
ggsave('EventStudies/ambulances.pdf', width = 6, height = 5,device=cairo_pdf)正如图表中缺乏明显跃升所预示的，\(表17.1\) 中对应的 \(Logit\) 回归结果显示：\(After\)变量的系数\(\beta_1\)与\(1\)无显著差异。与多数医学研究一致，该系数以比值比\(OR\)形式呈现——这意味着效应是乘积性而非加性的。此时零效应基准值为\(1\)（而非\(0\)），因为乘以\(1\)不会改变原有数值。
不过，\(图17.5\)确实显示出斜率的变化趋势——这正是回归模型擅长捕捉、而金融方法难以识别的效应类型。虽然政策未产生立竿见影的效果，但或许随着时间推移逐步改善了救治水平？表面数据似乎支持这一推断，然而回归结果表明：斜率变化量\(\beta_3\)（即交互项\((Week-27) \times After\)的系数）同样不显著。具体而言，斜率从\(\beta_1\)变为\(\beta_1+\beta_3\)，其差异\((\beta_1 + \beta_3) - \beta_1 = \beta_3\)在 \(表17.1\) 中未通过显著性检验。看来这项政策的效果确实乏善可陈！
17.3 专业实践方法
17.3.1 多处理组事件研究法
本章截至目前所述方法均基于 单一处理组 的研究设定——即事件仅影响特定群体（或我们仅关注该群体的效应）。例如，我们研究\(Alphabet\)重组公告对谷歌股价的影响时，既未预期也不关注其对其他股票的潜在溢出效应。
然而，许多事件会同时影响多个群体。以\(Alphabet\)重组为例，其影响仅局限于谷歌股票；但若涉及像《金融工具市场指令》(\(MiFID\))这类欧盟证券市场监管政策时，所有在该地区交易的股票都会受到影响。这种情况下该如何进行事件研究？
针对这种情况，我们有以下几种分析路径可供选择：
我们可以采用一种简化处理方式：虽然实际存在多个受影响群体，但将各组时间序列数据横向聚合——若有\(10\)个可能受影响的组别，只需计算各时间点上所有组别的均值，生成一条综合时间序列，再对其应用标准事件研究方法即可。
这个方法其实没有太多需要额外解释的——本质上就是把多组时间序列数据简单平均，然后套用单组分析流程。但这种简易处理存在明显缺陷：数据聚合过程会导致大量信息丢失，使得最终估计结果的精确度显著低于直接利用原始数据全部变异信息的方法。
我们可以采用分组独立分析法：既然已掌握单组事件研究方法，且拥有多组数据，那么只需为每个组别单独进行事件研究，即可获得各组差异化的处理效应估计12。
通过这种方法，您将获得每个组别的独立效应估计——可以是事件后单个时点的即时效应、回归模型中 \(After\) 变量的系数、事件后整个时期的平均效应或累计效应。在此基础上，您可以：
分析效应分布：考察所有组别的效应异质性
计算综合效应：聚合得出整体即时/平均/累计效应
我们可以采用回归聚合方法。前文”操作实施”章节介绍的回归式事件研究框架天然支持每个时间点包含多个观测值，因此只需将各组事件研究数据纵向堆叠即可。事实上，\(Taljaard \space et \space al.\)的救护车研究正是这样操作的——他们整合了多家医院的运营数据进行联合分析。
该回归模型设定与单时间序列分析几乎完全相同：
\[ \text{Outcome} = \beta_i + \beta_1 t + \beta_2 \text{After} + \beta_3 t \times \text{After} + \varepsilon \quad (17.3) \]
与之前类似，模型在事件前拟合的直线为\(\beta_i + \beta_1 t\)，事件后则为\((\beta_i + \beta_2) + (\beta_1 + \beta_3)t\)。唯一区别在于将\(\beta_0\)替换为组别固定效应\(\beta_i\)（如\(第16章\)所述）。虽然并非必须，但当各组别的结果变量基线水平存在差异时，固定效应能有效控制这一特征。此外，该方程同样允许通过添加多项式项来拟合非线性趋势13。
这种回归方法得到的结果与先将所有数据取平均值再进行回归分析的结果非常相似。但由于它使用了所有个体数据点，因此能够准确识别每个时间段的变异程度，从而更精确地估计标准误。
当然，这种回归方法将我们限制在特定的函数形式中。如果我们研究的事件效应仅持续一天，或需要完整刻画效应动态轨迹时该怎么办？此时可以采用另一种回归框架14：
\[ \text{Outcome} = \beta_0 + \beta_t + \varepsilon \quad (17.4) \]
这个简化方程有何含义？它实质上是将结果变量对 时间段固定效应 进行回归。这里我们采用了一种特殊设定：以事件产生效应前的最后\(1\)期作为参照组（注意：当存在一组二元指示变量时，必须留出一个参照组）。
特定时间段的固定效应系数，实际上表示该时间段结果变量的均值相对于事件发生前最后一期的估计差异。如果将这些时间段固定效应绘制出来，会形成一条类似于我们自行对各期数据取均值后合成的单一时间序列。唯一的区别在于：\(（1）\)现在每个时期都配有相应的标准误；\(（2）\)所有效应值均已相对于事件前基准期进行了标准化处理。
通过将事件发生前最后一期（如\(第3期\)）设为参照组，我们可以更清晰地识别事件效应：
若事件后 \(第4期\) 的估计系数 \(\hat {\beta_4} = 0.3\)，则表明事件导致结果变量从\(第3期\)到\(第4期\)增加了\(0.3\)，即单日效应为\(0.3\)
若 \(第51期\) 系数 \(\hat{\beta_5}=0.2\)，则表示\(第3期\)至\(第5期\)的累计效应为\(0.2\)
此外，通过对特定时段的系数取平均或求和，还可计算平均效应或长期累计效应。
以下是使用模拟数据进行该分析的代码实现步骤：
library(tidyverse)
library(fixest)
set.seed(10)
# Create data with 10 groups and 10 time periods
df <- crossing(id = 1:10, t = 1:10) %>%
# Add an event in period 6 with a one-period positive effect
mutate(Y = rnorm(n()) + 1*(t == 6))
# Use i() in feols to include time dummies,specifying that we want to drop t = 5 as the reference
m <- feols(Y ~ i(t, ref = 5), data = df, cluster = 'id')
# Plot the results, except for the intercept,and add a line joining them and a space and line for the reference group
coefplot(m, drop = '(Intercept)',pt.join = TRUE, ref = c('t:5' = 6), ref.line = TRUE)
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import seaborn.objects as so
import matplotlib.pyplot as plt
# Use a seed to make the results consistent
rng = np.random.default_rng(10)
# Ten groups with ten periods each
id = pd.DataFrame({'id': range(0,10), 'key': 1})
t = pd.DataFrame({'t': range(1,11), 'key': 1})
d = id.merge(t, on = 'key')
# Add an event in period 6 with a one-period effect
d['Y'] = rng.normal(0,1,100) + 1*(d['t'] == 6)
# Estimate our model using time 5 as reference
m = smf.ols('Y~C(t, Treatment(reference = 5))', data = d)
# Fit with SEs clustered at the group level
m = m.fit(cov_type = 'cluster',cov_kwds={'groups': d['id']})
# Get coefficients and CIs
# The original table will have an intercept up top
# But we'll overwrite it with our 5 reference
p = pd.DataFrame({'t': [5,1,2,3,4,6,7,8,9,10],
'b': m.params, 'se': m.bse,
'ci_top': list(m.conf_int().iloc[:,1]),
'ci_bottom': list(m.conf_int().iloc[:,0])})
# Add our period-5 zero
p.iloc[0] = [5, 0, 0, 0, 0]
p.head(10)## t b se ci_top ci_bottom
## Intercept 5 0.000000 0.000000 0.000000 0.000000
## C(t, Treatment(reference=5))[T.1] 1 -0.185602 0.304436 0.411082 -0.782286
## C(t, Treatment(reference=5))[T.2] 2 -0.161321 0.297543 0.421853 -0.744495
## C(t, Treatment(reference=5))[T.3] 3 0.312131 0.271123 0.843523 -0.219260
## C(t, Treatment(reference=5))[T.4] 4 0.572810 0.274704 1.111220 0.034400
## C(t, Treatment(reference=5))[T.6] 6 0.957727 0.549068 2.033881 -0.118426
## C(t, Treatment(reference=5))[T.7] 7 0.465978 0.360241 1.172037 -0.240081
## C(t, Treatment(reference=5))[T.8] 8 0.031695 0.291827 0.603666 -0.540276
## C(t, Treatment(reference=5))[T.9] 9 0.702904 0.203525 1.101806 0.304002
## C(t, Treatment(reference=5))[T.10] 10 0.548795 0.342223 1.219538 -0.121949# Plot the estimates as connected lines with error bars
fig = plt.figure()
so.Plot(p, x = 't', y = 'b', ymax = 'ci_top', ymin = 'ci_bottom').on(fig).add(so.Line()).add(so.Range()).plot()## <seaborn._core.plot.Plotter object at 0x0000021C82E28CD0>fig.axes[0].axvline(5, linestyle = '--')
fig.axes[0].axhline(0)
library(RStata)
options("RStata.StataPath"='D:/ProgramData/stata15/StataSE-64.exe')
options("RStata.StataVersion"=15.1)
dotext <- '
* Ten groups with ten time periods each
clear
set seed 10
set obs 100
g group = floor((_n-1)/10)
g t = mod((_n-1),10)+1
* Add an event in period 6 with a one-period effect
g Y = rnormal() + 1*(t == 6)
* ib#. lets us specify which # should be the reference
* We want 5 to be the reference, right before the event
* And exclude the regular constant (nocons)
reg Y ib5.t, vce(cluster group) nocons
* Plot the coefficients
coefplot, drop(_cons) base vertical recast(line)
'
stata(dotext, data.in = NULL, data.out = TRUE)## .
## . * Ten groups with ten time periods each
## . clear
## . set seed 10
## . set obs 100
## number of observations (_N) was 0, now 100
## . g group = floor((_n-1)/10)
## . g t = mod((_n-1),10)+1
## .
## . * Add an event in period 6 with a one-period effect
## . g Y = rnormal() + 1*(t == 6)
## .
## . * ib#. lets us specify which # should be the reference
## . * We want 5 to be the reference, right before the event
## . * And exclude the regular constant (nocons)
## . reg Y ib5.t, vce(cluster group) nocons
##
## Linear regression Number of obs = 100
## F(9, 9) = 379.53
## Prob > F = 0.0000
## R-squared = 0.1357
## Root MSE = 1.0559
##
## (Std. Err. adjusted for 10 clusters in group)
## ------------------------------------------------------------------------------
## | Robust
## Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## t |
## 1 | .0089498 .337251 0.03 0.979 -.7539649 .7718645
## 2 | -.1907903 .2160079 -0.88 0.400 -.679434 .2978535
## 3 | -.5842148 .2252918 -2.59 0.029 -1.09386 -.0745693
## 4 | .3265403 .3086343 1.06 0.318 -.3716391 1.02472
## 6 | 1.037317 .4517982 2.30 0.047 .0152786 2.059356
## 7 | -.1208886 .3356411 -0.36 0.727 -.8801614 .6383842
## 8 | .0754189 .3308649 0.23 0.825 -.6730495 .8238873
## 9 | .0702192 .3787089 0.19 0.857 -.7864797 .9269182
## 10 | -.0808192 .3330625 -0.24 0.814 -.834259 .6726206
## ------------------------------------------------------------------------------
## .
## . * Plot the coefficients
## . coefplot, drop(_cons) base vertical recast(line)由于随机数生成的影响，每次分析结果会略有差异，但\(R\)语言输出的结果如 \(图17.6\) 所示15。可以看到，我们已将事件前最后一期（第\(5\)期）的效应固定为\(0\)，使所有效应均以此为基准。结果确实在第\(6\)期（真实效应所在期）出现了预期跃升——置信区间完全偏离\(0\)值，表明该日效应具有统计显著性。不过，在第\(2\)期和第\(4\)期（实际无真实效应）也出现了显著结果，这正是小样本分析难以避免的假阳性问题。
library(tidyverse)
library(fixest)
set.seed(10)
# Create data with 10 groups and 10 time periods
df <- crossing(id = 1:10, t = 1:10) %>%
# And an event in period 6 with a one-period positive effect
mutate(Y = rnorm(n()) + 1*(t == 6))
# Use i() in feols to include time dummies,specifying that we want to drop t = 5 as the referencen and no intercept (-1)
m <- feols(Y ~ -1 + i(t, ref = 5), data = df, cluster = 'id')
d <- broom::tidy(m) %>% mutate(ci.bot = estimate - std.error, ci.top = estimate + std.error) %>% mutate(t = c(5,1:4,6:9))
d[1,] <- list('t::5', 0, 0, 0, 0,0,0,5)
d %>%
arrange(t) %>%
ggplot(aes(x = t, y = estimate)) +
geom_line() + geom_point() +
geom_errorbar(aes(ymin = ci.bot, ymax = ci.top), width = .2) +
geom_hline(aes(yintercept = 0)) +
geom_vline(aes(xintercept = 5), linetype = 'dashed') +
scale_x_continuous(breaks = c(1,3,5,7,9)) +
scale_y_continuous(breaks = c(-.5, 0, .5, 1, 1.5)) +
labs(x = 't', y = 'Estimate') +
theme_pubr() +
theme(text= element_text(family = 'sans', size = 10))
ggsave('EventStudies/period_regression.pdf', width = 6, height = 5,device=cairo_pdf)采用这种方法进行事件研究时，必须明确其核心逻辑：该模型未纳入任何时间趋势——它将每个时期视为独立单元，各期的反事实预测完全基于\(第5期\)（所有效应均以此为基准）。这种方法仅适用于不存在显著时间趋势的情况16。即便如此，仍需严格检验事件前各期的伪效应。
该方法默认事件前各期的效应值应为零——若非如此，则暗示模型存在潜在问题。这实际上带来一个额外优势：通过检验事件前效应的显著性，可有效识别研究设计缺陷。虽然可能仅是抽样变异所致（如本例中已知的数据生成过程），但同样可能意味着存在未被妥善控制的时变混杂因素。理想情况下，若发现事件前存在非零效应，它们应呈现某种规律性趋势（可通过建模加以控制）；但若呈现无序波动，则可能意味着研究设计本身存在根本性缺陷。
在多组别事件研究中，无论采用何种方法，我们在构建反事实预测时都存在一定局限。那些依赖其他组别信息来构建反事实的方法（例如利用整体市场趋势预测处理组若无干预时的表现）在多数组别均受处理时将失效。不过，不依赖跨组别信息的方法仍可适用——包括：
回归方法：通过时间趋势建模
均值调整法：以事件前收益率的均值作为事件后预测基准（如股票收益率分析）
当然，我们还可以寻找 未受处理的额外组别 作为对照组——但这属于其他章节的讨论范畴，包括：
第\(18\)章：双重差分法\(（DID）\)
第\(21\)章：合成控制法
17.3.2 基于时间序列方法的规范预测
事件研究的核心逻辑在于：通过对比事件后的实际观测值与反事实预测值来识别处理效应。该预测始终基于一个关键假设——若无事件发生，事件前的趋势模式将会延续。换言之，我们正利用事件前的时序数据（及可能的其他预测变量）进行 外推预测 ——这正是时间序列预测方法的本质功能。
严格来说，本章目前介绍的所有方法都属于时间序列预测的范畴——但仅限于形式上的归类。真正专业的时间序列预测方法会严谨处理数据的时序特性：时序数据具有跨期关联性，某一时期的变化不仅影响当期，还会持续影响后续所有观测路径。这种特性会通过无数难以量化的方式干扰统计估计和预测结果。
严格来说，本章目前介绍的所有方法都属于时间序列预测的范畴——但仅限于形式上的归类。真正专业的时间序列预测方法会严谨处理数据的时序特性：时序数据具有跨期关联性，某一时期的变化不仅影响当期，还会持续影响后续所有观测路径。这种特性会通过无数难以量化的方式干扰统计估计和预测结果。
因此，若需严谨构建反事实预测，且研究对象非股票收益率这类理想场景时，就必须采用真正考虑时序特性的预测模型。
需要说明的是，时间序列估计与预测不仅是统计学中的一个完整领域，更是统计学、计量经济学、机器学习及金融学中的庞大分支——甚至形成了独立产业体系，还包括完整的贝叶斯方法体系。鉴于其复杂性，本文不试图全面涵盖，仅以\(ARMA\)模型作为标准时间序列预测方法的示例。如需系统学习，推荐参考 \(Montgomery、Jennings和Kulahci（2015）\)的经典著作。
Douglas C Montgomery, Cheryl L Jennings, and Murat Kulahci. Introduction to Time Series Analysis and Forecasting. John Wiley & Sons, 2015\(ARMA\)模型名称源自时间序列数据的两个关键特征：
\(AR\)（自回归）：当期数值直接影响下一期。例如今日储蓄账户存入\(1\)美元，会使今日余额增加\(1\)美元，明日余额也增加\(1\)美元（加上利息），后日继续增加（利息更多），这种持续影响就是自回归。储蓄账户属于 \(AR(1)\) 过程——“\(1\)”表示只有最近一期的\(Y\)值对当期预测有用。
预测特性：
预测今日余额时，昨日余额最关键
若已知昨日余额，前日数据无增量信息（否则需用\(AR(2)\)模型）
\(AR(1)\)模型的数学表达式为：
\[ Y_t = \beta_0 + \beta_1 Y_{t-1} + \varepsilon_t \quad (17.5) \]
其中下标\(t\)和\(t-1\)表示用上一期的\(Y\)值（即\(Y\)的滞后项）预测当期\(Y\)值。
\(MA\) 是”移动平均\(（moving-average）\)“的缩写。移动平均指的是，某事物受到的暂时性影响会持续一小段时间，之后逐渐消失。比如，我从自己的储蓄账户里借了 \(10\) 美元给朋友，而她在接下来的几周里分批次把钱还给我 —— 那么今天我的账户会先减少 \(10\) 美元；到了下周，她还回一部分钱，账户的减少额就会变少一些；再下一周，她还得更多，账户的减少额会进一步降低；等她把钱全部还清后，我的账户余额就会恢复正常，仿佛这件事从未发生过。这就是一个移动平均过程。
\(MA(1)\)过程——即仅最近一期暂时性冲击产生影响，其模型形式为：
\[ Y_t = \beta_0 + (\varepsilon_t + \theta \varepsilon_{t-1}) \quad (17.6) \]
（注：由于无法直接观测扰动项\(\varepsilon_t\)和\(\varepsilon_{t-1}\)，无法通过常规回归估计\(\theta\)参数，需采用最大似然估计等替代方法）
根据上述特征描述，储蓄账户余额的时间序列明显同时具备自回归\(AR\)和移动平均\(MA\)特性，因此可以判定其遵循\(ARMA\)过程，需采用\(ARMA\)模型进行未来余额预测。以\(ARMA(2,1)\)模型为例，该模型包含两个自回归项和一个移动平均项，其当前值取决于前两期的历史值，以及前一期的暂时性冲击影响。
\[ Y_t = \beta_0 + \beta_1 Y_{t-1} + \beta_2 Y_{t-2} + (\varepsilon_t + \theta \varepsilon_{t-1}) \quad (17.7) \]
通过 \(图17.7\) 可以直观理解这些模型的差异：\(（a）\)图展示完全随机的白噪声；\(（b）\)图显示\(MA(1)\)过程——相邻时点的数值相关性开始增强；\(（c）\)图呈现\(AR(1)\)过程，序列明显更平滑，因为当期值直接受前期值影响；最后的\(ARMA(1,1)\)过程则表现出更强的时序关联性。需特别注意：这四组数据源自相同的初始随机序列，只是时序结构被逐步强化。
library(cowplot)
# ARMA data
d <- data.frame(t = 1:100, Y = rnorm(100)) %>%
mutate(Y_AR = Y,Y_MA = Y,Y_ARMA = Y)
for (i in 2:100) {
d[i, 'Y_AR'] <- d[i,'Y_AR'] + .8*d[i-1,'Y_AR']
d[i, 'Y_MA'] <- d[i, 'Y_MA'] + .9*d[i-1, 'Y']
d[i, 'Y_ARMA'] <- d[i, 'Y_ARMA']+ .8*d[i-1,'Y_ARMA'] + .9*d[i-1, 'Y']
}
names <- c('(a) Totally Random', '(b) MA(1) Process', '(c) AR(1) Process', '(d) ARMA(1,1) Process')
vars <- c('Y','Y_MA','Y_AR','Y_ARMA')
p <- list()
for (i in 1:4) {
p[[i]] <- ggplot(d, aes_string(x = 't', y = vars[i])) +
geom_line() +
labs(x = 'Time', y = 'Value', title= names[i]) +
theme_pubr() +
theme(text= element_text(family = 'sans', size = 8))
}## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plot_grid(p[[1]],p[[2]],p[[3]],p[[4]], ncol = 2)
ggsave('EventStudies/different_processes.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)\(AR\)和\(MA\)绝非时间序列模型仅有的要素——分析师常通过叠加各类扩展项构建复杂模型17：
\(ARIMA\)：加入差分项\(（I）\)，将绝对量转为日变动量分析
\(ARIMAX\)：引入时变外生变量\(X\)（如月工资收入）
\(SMARIMAX\)：纳入季节性 \(S\) 成分（如圣诞季储蓄波动）18
时间序列预测的核心步骤是：\(（1）\)选择合适模型 \(\rightarrow\)\(（2）\)拟合参数 \(\rightarrow\) \(（3）\)进行预测。这一思路虽与常规回归模型类似，但关注点截然不同——此时不再聚焦变量间的函数形式或控制变量选择，而是专注于单一时序本身的特性：是否存在季节性？需差分处理吗？自相关阶数多少？移动平均项又该取几阶？这些问题才是建模关键。
正如前文所述，解答这些问题本身就是一个完整的学科领域。虽然本文介绍了基础概念，但远不足以支撑完整的实操能力——为避免提供”残缺的代码示例”造成误导，建议系统学习时间序列专业教材。
最后提供几个关键提示：首先注意多数软件包不单独提供 \(ARMA\) 函数，而是通过 \(ARIMA\) 实现（当差分阶数\(I\)设为\(0\)时即为\(ARMA\)）。在\(R\)语言中，可通过\(fable\)包轻松实现\(ARIMA\)等核心模型的估计与预测；\(Stata\)内置了\(arima\)等标准模型命令，使用\(predict\)配合\(dynamic()\)选项即可进行预测；\(Python\)的\(statsmodels.tsa\)模块则提供多种时间序列模型，均配备\(.predict()\)预测方法。
17.3.3 联合检验问题
事件研究设计高度依赖于反事实预测的准确性19——该预测是事件研究有效性的唯一保障。例如，若你假设股票收益率在无干预时会保持不变，但实际应上涨\(0.1\%\)，则所有效应估计将系统性偏差\(0.1\%\)。对于日频投资组合收益而言，这种偏差足以导致重大误判。
这意味着事件研究的结果本质上是双重混合体——既包含真实的事件效应，也反映了所用反事实预测模型的准确性。当你对效应进行显著性检验时，检验的不仅是真实效应，还隐含了对预测模型正确性的联合检验。二者无法分离！这种情况令人担忧——因为我们永远无法确知模型是否正确。即便所选模型能完美拟合事件前数据，但反事实的本质决定了其不可观测性：很可能模型在事件发生的同时突然失效。
面对联合检验问题，我们能采取哪些措施？该问题本质上无法彻底消除，但可通过优化反事实预测模型来降低风险。不过，确实存在一种方法能部分剥离这个联合检验中的成分。
事件研究中的安慰剂检验应用由来已久，其方法体系在\(Brown和Warner（1985）\)的研究中已基本确立规范化20。
Stephen J Brown and Jerold B Warner. Using daily stock returns: The case of event studies. Journal of Financial Economics, 14(1):3–31, 1985.其核心逻辑是：当理论上应存在效应时，我们无法区分检测结果中哪些部分来自真实效应，哪些源于反事实模型的设定错误。但若在已知\(0\)效应的情境下进行检验，任何估计结果都只能反映模型缺陷——此时若仍检出”效应”，便可确证模型存在问题并需调整方法。
因此，需要选取大量不同的时间序列进行测试，特别是那些 不受事件影响 的样本。以股票收益为例，不仅应包含特定日期发布重大公告的公司，还需纳入同期无公告 的其他公司股票作为对照。
接着，需要重复进行事件研究——每次随机选择一个时间序列，并随机指定一个”伪事件日”作为分析基准21。
理论上，这些随机选取的”伪事件研究”不应存在真实效 应——尽管偶尔会碰巧选中确实发生异常波动的时点，但同样可能选中下跌或毫无变动的时段，最终效应均值应趋近于\(0\)。若大量检验后发现均值显著非零，则证明反事实模型存在系统性偏差，必须更换建模策略。
某些事件研究会区分三个时间段——事前、事后以及事件发生期间（或称”事件日/期”）。当预期处理效应持续时间极短时，研究者会采用这种方法。↩︎
这个顿悟正是1991年电影《摇滚大笨鸡》（Rock-a-Doodle）中那只公鸡的认知突破——现在您的教授终于有借口在课堂上播放这部电影，代替一天的正式教学了。↩︎
当研究者完全掌握数据生成过程，并能确证结果变量不存在任何潜在变化路径时（即不存在后门路径），该方法同样成立。↩︎
若出现此类概念混淆，不妨建议他们深呼吸，并默念研究假设必然成立——此乃缓解方法论焦虑的应急良方。↩︎
方法3的应用场景除外——因其本质上已部分实现了这类控制组功能。↩︎
经济学家通常将此称为“断点时间序列”（但需注意前文所述术语差异的说明）。↩︎
事件时点界定常面临方法论挑战。例如研究”某人当选公职”的影响时，事件时点可能存在多重定义：民调显示其将胜选时？实际胜选时刻？正式就职日？首个法案通过时？↩︎
两项技术说明：
参数符号惯例：此处沿用金融学传统，使用\(\alpha\)和\(\beta\)而非\(\beta_0\)和\(\beta_1\)表示回归系数。
模型扩展性：可轻松引入其他投资组合（如构建Fama-French三因子模型等高阶定价模型）。
有人可能会提出合理的质疑：考虑到谷歌在7月出现的巨大波动峰值，或许应该仅使用5月和6月数据作为估计窗口，而剔除7月数据。这一观点很可能成立。请仔细思考你是否认同这一批评意见。↩︎
关于自相关性的深入讨论，请参阅第13章及本章”专业实践”部分。↩︎
研究者实际采用了Logit回归（而非普通最小二乘法），但最终以比值比（Odds Ratio）形式呈现结果。这种处理虽保持趋势线的线性特征，但较难直观解释——至少对像我这样的大脑而言颇具挑战。↩︎
该方法的显著优势在于能获取组别特异性效应。正如第10章关于处理效应的讨论所指出的，不同组别间很可能存在真实的效应差异，而这种方法能清晰呈现这种异质性。↩︎
甚至可以采用第4章讨论的非参数方法——尽情发挥吧！↩︎
在估计时应当采用能够反映数据时序特性的标准误计算方法，例如在组别层面进行聚类（允许同一组别内不同时间点的观测值存在相关性）。↩︎
若不设置随机种子而多次运行分析，结果会出现明显波动——事件研究本质上存在较大随机误差，除非满足以下三个条件之一：足够多的组别、足够长的时间跨度，或者极其规整的数据质量。↩︎
至少该方法本身无法自动处理时间趋势。但我们可以先通过回归估计时间趋势（或与其他变量的关联），然后在应用事件研究方法前，从原始数据中减去这部分趋势成分。↩︎
差分处理的目的是获得平稳过程。简言之，由于复利效应，储蓄账户今日存入1美元对明日余额的影响会大于对当日的影响——若放任不管，这种自我强化机制将导致余额长期趋向无穷大！显然，统计模型难以处理这种无限循环的反馈机制。因此通过差分转换，期望新序列能消除这种发散性趋势。↩︎
此外，别忘了同样存在处理条件异方差（CH）的进阶模型链——当误差项的方差随时间波动时，可用ARCH、GARCH、EGARCH等模型捕捉这种特性。↩︎
从本质上说，所有因果推断设计都极度依赖准确的反事实预测。正因如此，类似此处的安慰剂检验方法也常见于其他章节的因果分析中。↩︎
在因果推断框架下，安慰剂检验是指在 理论效应为零 的情境下进行的检验——若检测到显著效应，则表明研究设计存在缺陷。↩︎
实现该分析所需的完整代码详见第15章《模拟分析》。↩︎