手把手教你做事件研究-CAR 计算与显著性检验

计量分析
Author

Tom

Published

February 2, 2023

这篇博客我们会完成事件研究的核心部分，包括正常收益率估计、异常收益率计算、累计异常收益率计算和显著性检验，视频讲解在 b 站。

先回顾一下事件研究法的核心步骤，包括： - 定义事件与事件窗口 - 研究样本选择与数据准备 - 估计正常收益率与计算异常收益率 - 累计异常收益率计算与显著性检验

import datetime
import logging

import pandas as pd
import numpy as np

import statsmodels.api as sm
import tushare as ts
data_2013_2017 = pd.read_csv(
    "./data/日个股回报率文件2013_2017/TRD_Dalyr.csv", dtype={"Stkcd": str}
)
data_2013_2017_1 = pd.read_csv(
    "./data/日个股回报率文件2013_2017/TRD_Dalyr1.csv", dtype={"Stkcd": str}
)
data_2013_2017_2 = pd.read_csv(
    "./data/日个股回报率文件2013_2017/TRD_Dalyr2.csv", dtype={"Stkcd": str}
)
data_2013_2017_3 = pd.read_csv(
    "./data/日个股回报率文件2013_2017/TRD_Dalyr3.csv", dtype={"Stkcd": str}
)

data_2018_2022 = pd.read_csv(
    "./data/日个股回报率文件2018_2022/TRD_Dalyr.csv", dtype={"Stkcd": str}
)
data_2018_2022_1 = pd.read_csv(
    "./data/日个股回报率文件2018_2022/TRD_Dalyr1.csv", dtype={"Stkcd": str}
)
data_2018_2022_2 = pd.read_csv(
    "./data/日个股回报率文件2018_2022/TRD_Dalyr2.csv", dtype={"Stkcd": str}
)
data_2018_2022_3 = pd.read_csv(
    "./data/日个股回报率文件2018_2022/TRD_Dalyr3.csv", dtype={"Stkcd": str}
)
# data:按流通市值加权的A股日收益率数据
data = pd.concat(
    [
        data_2013_2017,
        data_2013_2017_1,
        data_2013_2017_2,
        data_2013_2017_3,
        data_2018_2022,
        data_2018_2022_1,
        data_2018_2022_2,
        data_2018_2022_3,
    ]
)
len(data)
7012621
data.head()
Stkcd Trddt Dretwd
0 000001 2013-01-04 -0.001873
1 000001 2013-01-07 0.019387
2 000001 2013-01-08 -0.018405
3 000001 2013-01-09 -0.008750
4 000001 2013-01-10 0.000631 
data = data.rename(columns={"Trddt": "TradingDate"})  # 更改列名，方便后续做比较合并
# 整理 3因子数据
ff_factor = pd.read_csv(
    "./data/3factor.csv", skiprows=[1, 2], converters={"TradingDate": pd.to_datetime}
)
ff_factor.head(5)
MarkettypeID TradingDate RiskPremium1 SMB1 HML1
0 P9706 1990-12-19 2.473374 NaN NaN
1 P9710 1990-12-19 2.473374 NaN NaN
2 P9709 1990-12-19 2.473374 NaN NaN
3 P9712 1990-12-19 2.473374 NaN NaN
4 P9713 1990-12-19 2.473374 NaN NaN 
# 无风险利率数据
risk_free_rate = pd.read_excel(
    "./data/定期存款利率_3个月.xlsx", dtype={"指标名称": datetime.date}
)
risk_free_rate.head(5)
指标名称 定期存款利率:3个月(月)
0 1989-02-28 7.56
1 1989-03-31 7.56
2 1989-04-30 7.56
3 1989-05-31 7.56
4 1989-06-30 7.56 
r_dict = {}
for t, r in zip(risk_free_rate.iloc[:, 0], risk_free_rate.iloc[:, 1]):
    t_temp = t + datetime.timedelta(
        days=10
    )  # 月底公布下月三月定期存款利率，+10天将日期对应到下月
    r_dict[str(t_temp.year) + str(t_temp.month)] = r / 100

r_list = []
for t in ff_factor.iloc[:, 1]:
    r = r_dict[str(t.year) + str(t.month)]
    r_list.append(r)

ff_factor[["r"]] = np.array(r_list).reshape(
    -1, 1
)  # 将无风险利率添加到 ff 三因子表格内，方便后续回归拟合|
ff_factor.head(5)
MarkettypeID TradingDate RiskPremium1 SMB1 HML1 r
0 P9706 1990-12-19 2.473374 NaN NaN 0.0432
1 P9710 1990-12-19 2.473374 NaN NaN 0.0432
2 P9709 1990-12-19 2.473374 NaN NaN 0.0432
3 P9712 1990-12-19 2.473374 NaN NaN 0.0432
4 P9713 1990-12-19 2.473374 NaN NaN 0.0432 
ff_factor["TradingDate"] = ff_factor["TradingDate"].astype(
    str
)  # 将时间列转化为字符串类型，方便后续做判断合并
sample = pd.read_csv("./data/final.csv", dtype={"secCode": str})
sample.head()
secCode secName orgId announcementTitle announcementTime adjunctUrl announcementContent
0 002152 广电运通 9900003423 广电运通：2016年4月1日投资者关系活动记录表 2016-04-06 http://static.cninfo.com.cn/finalpage/2016-04-... 9、公司提过有储备区块链技术，是否可以介绍？
1 002177 御银股份 9900003781 御银股份：2016年4月21日投资者关系活动记录表 2016-04-21 http://static.cninfo.com.cn/finalpage/2016-04-... 4、公司在区块链技术上未来有什么规划？
2 300386 飞天诚信 9900023058 飞天诚信：2016年4月22日投资者关系活动记录表 2016-04-25 http://static.cninfo.com.cn/finalpage/2016-04-... 公司在区块链部分如何布局?
3 300423 鲁亿通 9900028812 鲁亿通：关于本次交易前12个月内购买、出售资产的说明 2016-06-09 http://static.cninfo.com.cn/finalpage/2016-06-... 同时，上市公司拟采用定价发行的方式向纪法清、孔剑平、孙奇锋、王麒诚和中信建投-数贝泽华人工智...
4 000961 中南建设 gssz0000961 中南建设：关于本公司投资上海承泰信息科技股份有限公司的进展公告 2016-06-23 http://static.cninfo.com.cn/finalpage/2016-06-... ，以区块链技术处理玛娜花园不同版本的个人数据，形成对数据信息的安全保护。 
pro = ts.pro_api()

# 查询当前所有正常上市交易和退市的股票列表
stocks_L = pro.stock_basic(
    exchange="", list_status="L", fields="symbol,name,area,industry,list_date,market"
)
stocks_D = pro.stock_basic(
    exchange="", list_status="D", fields="symbol,name,area,industry,list_date,market"
)
stocks_L.head()
symbol name area industry market list_date
0 000001 平安银行 深圳 银行 主板 19910403
1 000002 万科A 深圳 全国地产 主板 19910129
2 000004 ST国华 深圳 软件服务 主板 19910114
3 000005 ST星源 深圳 环境保护 主板 19901210
4 000006 深振业A 深圳 区域地产 主板 19920427 
stocks_D.head()
symbol name area industry market list_date
0 000003 PT金田A(退) None None 主板 19910703
1 000013 *ST石化A(退) None None 主板 19920506
2 000015 PT中浩A(退) None None 主板 19920625
3 000018 神城A退(退) None None 主板 19920616
4 000024 招商地产(退) None None 主板 19930607 
stocks = pd.concat([stocks_L, stocks_D])
# 构造一个 code:market 字典
cm_dict = {}
for c, m in zip(stocks.loc[:, "symbol"], stocks.loc[:, "market"]):
    cm_dict[c] = m

market_l = []
for code in sample.loc[:, "secCode"]:
    m = cm_dict[code]
    market_l.append(m)

sample[["market"]] = np.array(market_l).reshape(-1, 1)
sample.head()
secCode secName orgId announcementTitle announcementTime adjunctUrl announcementContent market
0 002152 广电运通 9900003423 广电运通：2016年4月1日投资者关系活动记录表 2016-04-06 http://static.cninfo.com.cn/finalpage/2016-04-... 9、公司提过有储备区块链技术，是否可以介绍？ 中小板
1 002177 御银股份 9900003781 御银股份：2016年4月21日投资者关系活动记录表 2016-04-21 http://static.cninfo.com.cn/finalpage/2016-04-... 4、公司在区块链技术上未来有什么规划？ 中小板
2 300386 飞天诚信 9900023058 飞天诚信：2016年4月22日投资者关系活动记录表 2016-04-25 http://static.cninfo.com.cn/finalpage/2016-04-... 公司在区块链部分如何布局? 创业板
3 300423 鲁亿通 9900028812 鲁亿通：关于本次交易前12个月内购买、出售资产的说明 2016-06-09 http://static.cninfo.com.cn/finalpage/2016-06-... 同时，上市公司拟采用定价发行的方式向纪法清、孔剑平、孙奇锋、王麒诚和中信建投-数贝泽华人工智... 创业板
4 000961 中南建设 gssz0000961 中南建设：关于本公司投资上海承泰信息科技股份有限公司的进展公告 2016-06-23 http://static.cninfo.com.cn/finalpage/2016-06-... ，以区块链技术处理玛娜花园不同版本的个人数据，形成对数据信息的安全保护。 主板 
# MarkettypeID [股票市场类型编码] - P9705：创业板;P9706：综合A股市场（不包含科创板、创业板）;
# P9707：综合B股市场;P9709：综合A股和创业板; P9710：综合AB股和创业板；P9711：科创板；
# 构造一个 market:markettypeID 字典
mm_dict = {"中小板": "P9706", "主板": "P9706", "创业板": "P9705", "科创板": "P9711"}

markettypeID = []
for m in sample.loc[:, "market"]:
    m_id = mm_dict[m]
    markettypeID.append(m_id)

sample[["markettypeID"]] = np.array(markettypeID).reshape(-1, 1)
sample
secCode secName orgId announcementTitle announcementTime adjunctUrl announcementContent market markettypeID
0 002152 广电运通 9900003423 广电运通：2016年4月1日投资者关系活动记录表 2016-04-06 http://static.cninfo.com.cn/finalpage/2016-04-... 9、公司提过有储备区块链技术，是否可以介绍？ 中小板 P9706
1 002177 御银股份 9900003781 御银股份：2016年4月21日投资者关系活动记录表 2016-04-21 http://static.cninfo.com.cn/finalpage/2016-04-... 4、公司在区块链技术上未来有什么规划？ 中小板 P9706
2 300386 飞天诚信 9900023058 飞天诚信：2016年4月22日投资者关系活动记录表 2016-04-25 http://static.cninfo.com.cn/finalpage/2016-04-... 公司在区块链部分如何布局? 创业板 P9705
3 300423 鲁亿通 9900028812 鲁亿通：关于本次交易前12个月内购买、出售资产的说明 2016-06-09 http://static.cninfo.com.cn/finalpage/2016-06-... 同时，上市公司拟采用定价发行的方式向纪法清、孔剑平、孙奇锋、王麒诚和中信建投-数贝泽华人工智... 创业板 P9705
4 000961 中南建设 gssz0000961 中南建设：关于本公司投资上海承泰信息科技股份有限公司的进展公告 2016-06-23 http://static.cninfo.com.cn/finalpage/2016-06-... ，以区块链技术处理玛娜花园不同版本的个人数据，形成对数据信息的安全保护。 主板 P9706
... ... ... ... ... ... ... ... ... ...
325 300133 华策影视 9900013428 华策影视：关于全资子公司增资扩股暨关联交易的公告 2021-12-17 http://static.cninfo.com.cn/finalpage/2021-12-... 影业上海基于节目库搭建的正版影视素材库“华策新视界”已进入正式运营，目前正在推进基于区块链技... 创业板 P9705
326 300235 方直科技 9900011849 方直科技：关于对深圳证券交易所关注函的回复公告 2021-12-20 http://static.cninfo.com.cn/finalpage/2021-12-... 元宇宙是整合多种新技术而产生的新型虚实相融的互联网应用和社会形态，它基于扩展现实技术提供沉浸... 创业板 P9705
327 300649 杭州园林 9900023811 杭州园林：关于收购参股公司部分股权的公告 2021-12-24 http://static.cninfo.com.cn/finalpage/2021-12-... 一般项目：网络技术服务；物联网技术服务；技术服务、技术开发、技术咨询、技术交流、技术转让、技... 创业板 P9705
328 300533 冰川网络 9900027730 冰川网络：关于全资子公司设立三家子公司并完成工商注册登记的公告 2021-12-24 http://static.cninfo.com.cn/finalpage/2021-12-... 人工智能通用应用系统；智能控制系统集成；量子计算技术服务；电竞信息科技；云计算装备技术服务；... 创业板 P9705
329 301116 益客食品 9900028946 益客食品：法律意见书及补充法律意见书 2021-12-27 http://static.cninfo.com.cn/finalpage/2021-12-... 不含特种设备制造）；液气密元件及系统制造；工业自动控制系统装置销售；物料搬运装备销售；非居住... 创业板 P9705

330 rows × 9 columns

daily_returns = data.loc[data.iloc[:, 0] == "002152"]
daily_returns
Stkcd TradingDate Dretwd
676981 002152 2013-01-04 -0.019027
676982 002152 2013-01-07 -0.008621
676983 002152 2013-01-08 0.011594
676984 002152 2013-01-09 0.002865
676985 002152 2013-01-10 -0.003571
... ... ... ...
612679 002152 2022-03-10 0.008427
612680 002152 2022-03-11 0.016713
612681 002152 2022-03-14 -0.011872
612682 002152 2022-03-15 -0.046211
612683 002152 2022-03-16 0.023256

2173 rows × 3 columns

def get_window(code, t, m_id):
    """获取事件期和理论收益率估计期

    Args:
        code (str): stock code
        t (str): announcement date
        m_id (str): MarkettypeID

    Returns:
        eventWD (DataFrame): event window
        estimationWD (DataFrame): estimation window
    """
    # 个股日收益率数据
    daily_returns = data.loc[data.iloc[:, 0] == code]
    # 交易日数据
    trading_dates = daily_returns.loc[:, "TradingDate"]

    t_1 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=1)).strftime(
        "%Y-%m-%d"
    )
    t_2 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=2)).strftime(
        "%Y-%m-%d"
    )
    t_3 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=3)).strftime(
        "%Y-%m-%d"
    )
    t_4 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=4)).strftime(
        "%Y-%m-%d"
    )
    t_5 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=5)).strftime(
        "%Y-%m-%d"
    )
    t_6 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=6)).strftime(
        "%Y-%m-%d"
    )
    t_7 = (datetime.date.fromisoformat(t) + datetime.timedelta(days=7)).strftime(
        "%Y-%m-%d"
    )

    # 判断公告日是否为交易日，年假最长 7 天
    if t in list(trading_dates):
        i = list(trading_dates).index(t)
    elif t_1 in list(trading_dates):
        i = list(trading_dates).index(t_1)
    elif t_2 in list(trading_dates):
        i = list(trading_dates).index(t_2)
    elif t_3 in list(trading_dates):
        i = list(trading_dates).index(t_3)
    elif t_4 in list(trading_dates):
        i = list(trading_dates).index(t_4)
    elif t_5 in list(trading_dates):
        i = list(trading_dates).index(t_5)
    elif t_6 in list(trading_dates):
        i = list(trading_dates).index(t_6)
    elif t_7 in list(trading_dates):
        i = list(trading_dates).index(t_7)
    else:
        i = 0

    # 估计期为 200 天
    # 公司上市时间较短时，估计期数据不够
    if i >= 210:
        eventWD = daily_returns.iloc[i - 1 : i + 2, :]
        estimationWD = daily_returns.iloc[i - 210 : i - 10, :]
        eventWD.insert(0, "MarkettypeID", m_id)
        estimationWD.insert(0, "MarkettypeID", m_id)
        return estimationWD, eventWD
    else:
        return
def abnormal_return(estimationWD, eventWD, ff_factor):
    """获取异常收益率

    Args:
        estimationWD (DataFrame): estimation window
        eventWD (DataFrame): event window
        ff_factor (DataFrame): fama-french factors

    Returns:
        ar (array): abnormal returns
    """

    estimation_data = pd.merge(
        estimationWD, ff_factor, how="left", on=["MarkettypeID", "TradingDate"]
    )
    event_data = pd.merge(
        eventWD, ff_factor, how="left", on=["MarkettypeID", "TradingDate"]
    )

    X = estimation_data.loc[:, ["RiskPremium1", "SMB1", "HML1"]].values
    X = sm.add_constant(X)
    y = np.array(
        estimation_data.loc[:, "Dretwd"] - estimation_data.loc[:, "r"]
    ).reshape(-1, 1)
    model = sm.OLS(y, X)
    result = model.fit()

    # rr: 理论收益率
    rr = np.dot(
        sm.add_constant(event_data.loc[:, ["RiskPremium1", "SMB1", "HML1"]].values),
        result.params.reshape(-1, 1),
    ) + event_data.loc[:, "r"].values.reshape(-1, 1)

    # ar: 异常收益率，每只股票在事件期的收益率 array
    ar = event_data.loc[:, "Dretwd"].values.reshape(-1, 1) - rr
    return ar
sample.head()
secCode secName orgId announcementTitle announcementTime adjunctUrl announcementContent market markettypeID
0 002152 广电运通 9900003423 广电运通：2016年4月1日投资者关系活动记录表 2016-04-06 http://static.cninfo.com.cn/finalpage/2016-04-... 9、公司提过有储备区块链技术，是否可以介绍？ 中小板 P9706
1 002177 御银股份 9900003781 御银股份：2016年4月21日投资者关系活动记录表 2016-04-21 http://static.cninfo.com.cn/finalpage/2016-04-... 4、公司在区块链技术上未来有什么规划？ 中小板 P9706
2 300386 飞天诚信 9900023058 飞天诚信：2016年4月22日投资者关系活动记录表 2016-04-25 http://static.cninfo.com.cn/finalpage/2016-04-... 公司在区块链部分如何布局? 创业板 P9705
3 300423 鲁亿通 9900028812 鲁亿通：关于本次交易前12个月内购买、出售资产的说明 2016-06-09 http://static.cninfo.com.cn/finalpage/2016-06-... 同时，上市公司拟采用定价发行的方式向纪法清、孔剑平、孙奇锋、王麒诚和中信建投-数贝泽华人工智... 创业板 P9705
4 000961 中南建设 gssz0000961 中南建设：关于本公司投资上海承泰信息科技股份有限公司的进展公告 2016-06-23 http://static.cninfo.com.cn/finalpage/2016-06-... ，以区块链技术处理玛娜花园不同版本的个人数据，形成对数据信息的安全保护。 主板 P9706 
ar_l = []
delete_code = []
for code, t, m_id in zip(
    sample.loc[:, "secCode"],
    sample.loc[:, "announcementTime"],
    sample.loc[:, "markettypeID"],
):
    try:
        estimationWD, eventWD = get_window(code, t, m_id)
    except:
        delete_code.append(code)
        continue

    ar_array = abnormal_return(estimationWD, eventWD, ff_factor)
    ar = sum(ar_array)[0]
    ar_l.append(ar)
AR = pd.DataFrame(np.array(ar_l).reshape(-1, 1))
AR.describe()
0
count 304.000000
mean 0.014020
std 0.068285
min -0.173427
25% -0.018811
50% 0.003446
75% 0.032451
max 0.337091 
from scipy import stats
stats.ttest_1samp(ar_l, 0)
Ttest_1sampResult(statistic=3.579861379325141, pvalue=0.00040021773286985055)
stats.wilcoxon(ar_l, len(ar_l) * [0], correction=True, alternative="two-sided")
WilcoxonResult(statistic=19144.0, pvalue=0.008515412037149648)
delete_code
['300423',
 '603106',
 '300612',
 '000793',
 '601619',
 '000038',
 '300654',
 '300525',
 '000061',
 '601066',
 '300773',
 '300788',
 '601512',
 '688158',
 '688051',
 '688096',
 '300975',
 '003040',
 '300996',
 '300889',
 '601686',
 '003032',
 '688135',
 '300949',
 '301213',
 '301116']
len(delete_code)
26