本章我们将探讨生存分析和删失生存数据。这些方法源于一种独特的结局变量类型:事件发生时间 :
学习闭环
先修与可观察目标
完成本章后,学生应能:
从事件定义、时间原点、观察终点构造 \((Y,\delta)\) ,并逐行审计右删失编码,错误率为零。
对 10 个对象的小表写出每个风险集并手算 Kaplan–Meier 曲线,概率误差不超过 \(10^{-3}\) 。
写出 Cox 偏似然、score 与信息量,解释基线风险为何消去,并区分风险比与绝对风险。
诊断比例风险与独立删失假设;若不成立,提出时间交互、分层或替代模型。
在时间锁定的生存任务中报告删失率、事件数、C-index 与失败条件;将未惩罚与 Lasso-Cox 的稳定性比较列为拓展目标。
入口检查(5 分钟)
某公司在研究开始后第 3 年退市,另一公司到第 5 年末仍上市。写出两者的 \((Y,\delta)\) 。
展开入口评分、补救与异形复测
答案分别为 \((3,1)\) 与 \((5,0)\) 。若把第二家公司写成事件,先回到 式 11.1 与 式 11.2 ,完成“提前还款/研究结束/失访”三例分类后重测。
无提示检索、渐隐链与迁移
闭书写出“风险集 \(r_j\) —事件数 \(d_j\) —条件生存率—累乘”链。第一轮对照正文小表;第二轮只给事件/删失时点;第三轮把第 5 月事件改到第 6 月并自行更新风险集。陌生迁移任务是为订阅客户定义“从首次付费到流失”的原点、事件、删失和竞争风险,说明为什么不能把静态状态当成事件日期。
目标—评价映射
1—2
练习 1、2、6、7、9
风险集和时点逐行可复算
3
练习 10—12
核心题 10 给出偏似然与 score/info 骨架;11—12 为拓展推导
4
练习 3—5
能指出假设破坏及修复
5
练习 10 与时间队列案例
核心出口含删失率、事件数、时间外 C-index 和失败条件;练习 8 的惩罚路径为拓展
本章练习均标为 项目:无;生存分析不属于贯穿项目里程碑。但 小节 11.10 是两条标准路线都必交的独立审计 A11:60 分钟、10 分章内证据、占课程成绩 5%;它不能替代任何 M02—M09、M12 或 M13 里程碑,也不能用 A10 替换。
核心概念:生存分析的应用场景
生存分析虽然名称源于医学研究,但其应用远不止于医学领域:
医学研究 :患者生存时间、疾病复发时间
商业分析 :客户流失(churn)时间、产品使用寿命
工程技术 :设备故障时间、系统可靠性
金融领域 :信用卡违约时间、贷款提前还款时间
生存时间与删失时间
对于每个个体,存在两个真实时间:
生存时间(Survival Time) \(T\) :事件发生的真实时间(也称为失败时间或事件时间)
删失时间(Censoring Time) \(C\) :删失发生的时刻
我们实际观测到的是:
\[ Y = \min(T, C) \tag{11.1}\]
同时观测到一个状态指示变量:
\[ \delta = \begin{cases}
1 & \text{如果 } T \leq C \text{(观测到事件)} \\
0 & \text{如果 } T > C \text{(被删失)}
\end{cases} \tag{11.2}\]
直观理解:删失的含义
假设你进行一项为期 5 年的客户流失研究:
如果客户在第3年取消服务,你观测到完整的流失时间
如果客户在第 5 年研究结束时仍然活跃,你只知道其流失时间至少为 5 年,但不知道确切值
这第二种情况就是右删失(right censoring)
删失机制的重要假设
为了分析生存数据,我们需要对删失机制 做出关键假设:
独立删失假设 :在给定特征条件下,事件时间 \(T\) 与删失时间\(C\) 独立。
违反独立删失假设的例子
病情导致的失访 :如果病情严重的患者更倾向于退出研究,会导致对平均生存时间的高估
选择性失访 :如果男性重症患者比女性重症患者更容易失访,可能导致错误的性别生存时间比较
这些情况都违反了独立删失假设,分析结果会有偏差。
删失的类型
右删失(Right Censoring) :在本章 \(\delta=I(T\le C)\) 的约定下,删失个体满足 \(Y=C\) 且真实事件时间 \(T>Y\) ;\(T=Y\) 会被编码为事件而不是删失
左删失(Left Censoring) :真实事件时间\(T \leq Y\)
区间删失(Interval Censoring) :只知道事件发生在某个时间区间内
本章重点讨论右删失
Kaplan-Meier 生存曲线
生存函数(Survival Function)定义为:
\[ S(t) = \Pr(T > t) \tag{11.3}\]
这是时间 \(t\) 的递减函数,表示生存超过时间\(t\) 的概率。
Kaplan-Meier 估计量
令 \(t_1 < t_2 < \cdots < t_K\) 为不同的事件时刻,\(d_k\) 为 \(t_k\) 的事件数,\(q_k\) 为该时刻事件发生后退出观察的删失数,\(r_k\) 为 \(t_k\) 之前仍在风险集中的个体数。
Kaplan–Meier 乘积极限估计以各事件时点的条件生存概率连乘来处理右删失 (Kaplan 和 Meier 1958年 ) :
\[ \widehat{S}(t_k) = \prod_{j=1}^{k} \frac{r_j - d_j}{r_j} \tag{11.4}\]
对于介于 \(t_k\) 和 \(t_{k+1}\) 之间的时间 \(t\) ,设 \(\widehat{S}(t)=\widehat{S}(t_k)\) 。删失不会在当时降低生存概率,但会减少后续风险集。
Kaplan-Meier 的直观解释
考虑 10 个产品。第 2 个月有 2 个失效,随后 1 个删失;第 5 个月又有 1 个失效。第 2 月事件发生前 \(r_1=10,d_1=2\) ,所以 \(\widehat S(2)=8/10=0.8\) 。删失后剩 7 个进入后续风险集,因此第 5 月 \(r_2=7,d_2=1\) , \[
\widehat S(5)=\frac{8}{10}\times\frac{6}{7}\approx0.686.
\] 这张小时间线同时展示了事件降低生存概率、删失只改变后续分母的区别。
中国实际案例:A股上市公司退市生存分析
让我们使得A 股上市公司的基本数据,来进行一个关于公司退市(Delisting) 的生存分析。随着中国推行全面注册制并常态化退市机制,上市公司的“存活状态”成了市场和监管层关注的重点焦点。
时间起点 :公司上市日的(listed_date)。
事件 :de_listed_date 有记录且不晚于研究截止日期的“任何原因退市”;截止日当日事件按 \(T=C\) 编码为事件。本数据没有可靠区分主动退市、并购吸收、监管退市等竞争事件,因此结果不能直接解释成财务困境或破产风险。
删失 :截至数据截止,公司仍在正常交易(右删失)。
下面用 stock_basic_data 演示任何原因退市的 Kaplan–Meier 描述。索引所列核心动作的本地 setup 是 表 11.1 :运行前设置 BOOK_DATA_DIR,随后按源码执行至 列表 11.1 ;Cox consumer 另需 表 11.3 的两个 BOOK_M11_* 环境变量,再执行至 表 11.6 。核心报告标签分别为 ch11_survival_sample_audit 与 ch11_cox_results;缺任一输入时应停止,不能借用交互环境中的旧对象。曲线反映本样本、当前截止日和混合退市原因下的存续分布;它不能单独解释行业因果机制、财务困境或投资价值。
表 11.2 用三个确定性教学行锁定 式 11.2 的边界约定:\(T=C\) 必须归入事件,只有 \(T>C\) 才是右删失。该小表不代表经验样本。
import pandas as pd # 允许边界fixture在fresh kernel独立运行
censoring_boundary_fixture = pd.DataFrame({'case' : ['T<C' , 'T=C' , 'T>C' ], 'event_time' : [2.0 , 3.0 , 5.0 ], 'censor_time' : [4.0 , 3.0 , 4.0 ]}) # 构造三个确定性边界情形
censoring_boundary_fixture['observed_time' ] = censoring_boundary_fixture[['event_time' , 'censor_time' ]].min (axis= 1 ) # 按Y=min(T,C)构造观测时间
censoring_boundary_fixture['event' ] = censoring_boundary_fixture['event_time' ].le(censoring_boundary_fixture['censor_time' ]).astype(int ) # 按T<=C编码事件
assert censoring_boundary_fixture['event' ].tolist() == [1 , 1 , 0 ] # T=C必须与T<C同属事件
assert censoring_boundary_fixture['observed_time' ].tolist() == [2.0 , 3.0 , 4.0 ] # 逐行验证观测时间
censoring_boundary_fixture # 输出可人工审计的三行证据
数据加载完成后,我们需要遍历每一家上市公司,根据其上市日期与退市日期计算存续时间,构建生存分析所需的事件-时间数据集。
observation_end_date = pd.to_datetime('2024-01-01' ) # 设定研究观测截止日期
survival_records = [] # 初始化列表,用于收集每家公司的生存记录
for _, row in stock_basic_data.iterrows(): # 遍历每一家上市公司
listed_date = row['listed_date' ] # 获取上市日期
delisted_date = row['de_listed_date' ] # 获取退市日期
if pd.isna(delisted_date) or delisted_date > observation_end_date: # 缺失或截止日后退市才右删失
duration = (observation_end_date - listed_date).days / 365.25 # 从上市到截止日的存续年数
event = 0 # 标记为删失(未观测到退市事件)
else : # 截止日当日或之前退市均为事件
duration = (delisted_date - listed_date).days / 365.25 # 从上市到退市的存续年数
event = 1 # 标记为事件(退市)
if duration > 0 : # 排除存续时间为负的异常记录
survival_records.append({ # 将计算结果追加到列表
'order_book_id' : row['order_book_id' ], # 股票代码
'duration' : duration, # 存续年数
'event' : event, # 事件标记(0=删失,1=退市)
'industry' : row.get('citics_2019_l1_name' , '未知' ) # 中信一级行业分类
}) # 完成构建
company_survival_data = pd.DataFrame(survival_records) # 将列表转换为DataFrame
top_industry_categories = company_survival_data['industry' ].value_counts().nlargest(5 ).index # 提取样本量前五的行业
company_survival_data['Industry_Group' ] = company_survival_data['industry' ].apply ( # 对每行/列应用函数
lambda x: x if x in top_industry_categories else 'Other' # 非前五行业统一归类为Other
) # 完成构建
生存数据构建完毕后,我们首先绘制 A 股市场整体的 Kaplan-Meier 生存曲线,展示所有上市公司随时间推移的”存活概率”变化趋势。
kaplan_meier_fitter = KaplanMeierFitter() # 实例化KM生存函数估计器
kaplan_meier_fitter.fit( # 拟合整体生存曲线
company_survival_data['duration' ], # 存续时间
company_survival_data['event' ], # 事件标记
label= 'A股整体生存率 (未退市)' # 图例标签
) # 完成构建
fig, ax = plt.subplots(figsize= (10 , 6 )) # 创建10×6英寸画布
kaplan_meier_fitter.plot(ax= ax) # 绘制KM生存曲线
ax.set_title('A股上市公司生存曲线 (Time to Delisting)' , fontsize= 14 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置标题
ax.set_xlabel('上市年数' , fontsize= 12 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置 x 轴标签
ax.set_ylabel('生存概率 (未退市)' , fontsize= 12 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置 y 轴标签
ax.grid(True , alpha= 0.3 ) # 添加半透明网格线
plt.tight_layout() # 自动调整布局
plt.show() # 显示图形
请从 图 11.1 的当次运行输出读取曲线水平与置信区间。只能将它解释为本数据覆盖、观察截止日与“任何原因退市”定义下的描述性存续分布。曲线不能自动识别制度执行、“壳价值”或退市门槛的因果作用;尾部置信区间变宽则提醒风险集中样本减少。
整体生存曲线展示了 A 股市场的”宏观韧性”。接下来,我们进一步按行业分解,比较前五大行业的退市风险差异。
fig, ax = plt.subplots(figsize= (10 , 6 )) # 创建10×6英寸画布
for ind in top_industry_categories: # 遍历前五大行业
if ind != '未知' : # 跳过行业信息缺失的记录
mask = company_survival_data['industry' ] == ind # 构建当前行业的布尔掩码
kaplan_meier_fitter.fit( # 拟合当前行业的KM生存曲线
company_survival_data.loc[mask, 'duration' ], # 当前行业的存续时间
company_survival_data.loc[mask, 'event' ], # 当前行业的事件标记
label= ind # 使用行业名称作为图例
) # 完成构建
kaplan_meier_fitter.plot(ax= ax) # 将曲线绘制到同一画布
ax.set_title('不同行业A股的退市生存曲线比较' , fontsize= 14 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置标题
ax.set_xlabel('上市年数' , fontsize= 12 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置 x 轴标签
ax.grid(True , alpha= 0.3 ) # 添加网格线
plt.legend(prop= {'family' : 'Source Han Serif SC' }) # 使用已安装中文字体设置图例
plt.tight_layout() # 自动调整布局
plt.show() # 显示图形
图 11.2 的行业、曲线顺序和区间宽度都由当次样本决定。报告时应点名实际出现的行业、给出风险集人数和事件数,并用下一节的预先规定检验补充目测。即使曲线分离,也只是当前混合退市结局的组间关联,不是行业因果效应。
对数秩检验:行业间比较
我们使用对数秩检验来比较不同行业的“仍维持上市”生存曲线是否存在显著差异。
当我们用肉眼观察上一节中不同行业的Kaplan-Meier 生存曲线时,可能会觉得它们之间似乎存在明显的高低差异。但作为严谨的量化分析师,我们绝不能仅仅依靠视觉来做判断——由于各个行业在不同历史时期的上市节奏(即进入风险集的队列)以及最终被删失的比例大相径庭,那些看似分离的曲线,在统计学上是否真的存在不能被随机波动所解释的本质区别? 对数秩检验在各事件时点比较两组的观察事件数与零假设下期望事件数。它给出当前分组下的统计证据,不预设 p 值大小,也不能把组间差异解释为行业“基因”或因果机制。
from lifelines.statistics import logrank_test # 导入对数秩检验用于比较生存曲线差异
# 假设 company_survival_data 已经存在 (来自上一节)
# 我们选取两个主要行业进行比较,例如"电子" (Tech) 和"房地产" (Real Estate)
# 注意: 具体的行业名称取决于数据中的 industry_name
# 这里我们先打印一下行业名称看看(假设前两个Top行业)
industry_first = top_industry_categories[0 ] # 定义industry_first变量
industry_second = top_industry_categories[1 ] # 定义industry_second变量
print (f"比较行业: { industry_first} vs { industry_second} " ) # 输出结果到控制台
mask_industry_first = company_survival_data['industry' ] == industry_first # 定义mask_industry_first变量
mask_industry_second = company_survival_data['industry' ] == industry_second # 定义mask_industry_second变量
rank_test_results = logrank_test( # 提取测试集数据
company_survival_data.loc[mask_industry_first, 'duration' ], # 执行数据处理操作
company_survival_data.loc[mask_industry_second, 'duration' ], # 执行数据处理操作
company_survival_data.loc[mask_industry_first, 'event' ], # 执行数据处理操作
company_survival_data.loc[mask_industry_second, 'event' ] # 执行数据处理操作
) # 完成构建
print (f'对数秩检验统计量: { rank_test_results. test_statistic:.4f} ' ) # 输出结果到控制台
print (f'p值: { rank_test_results. p_value:.4f} ' ) # 输出结果到控制台
fig, ax = plt.subplots(figsize= (10 , 6 )) # 创建子图布局
# 训练/拟合模型
kaplan_meier_fitter.fit(company_survival_data.loc[mask_industry_first, 'duration' ], company_survival_data.loc[mask_industry_first, 'event' ], label= industry_first)
kaplan_meier_fitter.plot(ax= ax) # 绘制图形
# 训练/拟合模型
kaplan_meier_fitter.fit(company_survival_data.loc[mask_industry_second, 'duration' ], company_survival_data.loc[mask_industry_second, 'event' ], label= industry_second)
kaplan_meier_fitter.plot(ax= ax) # 绘制图形
# 设置子图标题
ax.set_title(f'行业生存曲线比较: { industry_first} vs { industry_second} ' , fontsize= 14 , fontproperties= 'Source Han Serif SC' )
ax.set_xlabel('上市年数' , fontsize= 12 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置子图 X 轴标签
ax.grid(True , alpha= 0.3 ) # 添加子图网格线
plt.legend(prop= {'family' : 'Source Han Serif SC' }) # 使用已安装中文字体添加图例
plt.show() # 显示图形
比较行业: 机械 vs 基础化工
对数秩检验统计量: 0.9421
p值: 0.3317
图 11.3 上方输出给出本次动态选择行业后的检验统计量和 p 值。应按预先规定的显著性水平报告“拒绝”或“未拒绝”零假设,并同时说明组别是按样本量选择、事件包含所有退市原因。目测曲线分离不能替代检验,未拒绝也不等于两组完全相同。
风险函数的Cox 比例风险模型
除了比较曲线,我们更希望量化各个因素(如初始财务状况、所在地域)对遭遇退市风险的影响。
这里先定义 Cox 模型使用的风险函数:
\[ h(t\mid x)=\lim_{\Delta t\to0}\frac{\Pr(t<T\le t+\Delta t\mid T>t,x)}{\Delta t}. \]
它是已存续到 \(t\) 时事件在紧接着的瞬间发生的条件率,不是普通概率。本节使用这一定义建模,后文再详述它与生存函数、累计风险的关系。
数学推导:Cox 比例风险模型与偏似然函数 (Partial Likelihood)
Cox 比例风险模型 (Proportional Hazards Model) 的核心思想是将风险函数 \(h(t|x_i)\) 分解为两部分。 \[ h(t|x_i) = h_0(t) \exp(x_i^T \beta) \] 其中 \(h_0(t)\) 是基线风险函数 (Baseline Hazard Function),表示当所有协变量为0 时的潜在基准风险。\(\exp(x_i^T \beta)\) 是相对基准的风险乘数 :
Cox 比例风险模型通过偏似然估计回归参数 \(\beta\) ,不必为基准风险函数 \(h_0(t)\) 预先指定参数形式 (Cox 1972年 ) ;这项性质不意味着绝对风险已经由偏似然单独识别。
假设所有样本中的\(K\) 个非删失的独特事件(例如退市)发生时间 \(t_1 < t_2 < \dots < t_K\) 。令 \(\mathcal{R}(t_i)\) 为在时间 \(t_i\) 刚好处于风险集 (Risk Set)中的所有个体的集合(含义是:那些存活到\(t_i\) 并且尚未发生事件、也未被删失的全部个体)。
给定在时的\(t_i\) 确切地一个人 发生了事件(暂不考虑 Tied Times 打平),那么这个人恰巧是受试的\(j \in \mathcal{R}(t_i)\) 的条件概率依据风险律直接写出为: \[ P(\text{个体 } j \text{ 在} t_i \text{ 发生事件} \mid \text{在} t_i \text{ 有一个事件发生}) = \frac{h(t_i|x_j)}{\sum_{k \in \mathcal{R}(t_i)} h(t_i|x_k)} \] 清爽地代入Cox 模型的风险函数乘性公式: \[ \frac{h_0(t_i) \exp(x_j^T \beta)}{\sum_{k \in \mathcal{R}(t_i)} h_0(t_i) \exp(x_k^T \beta)} = \frac{\exp(x_j^T \beta)}{\sum_{k \in \mathcal{R}(t_i)} \exp(x_k^T \beta)} \]
奇妙的数学分离在这里发生了:未知的非参数基准常数的\(h_0(t_i)\) 在多项式分子和分母中被平顺地消去了!
然后将所有发生独立事件的时刻 \(t_1, \ldots, t_K\) 的这些离散的条件概率在时间轴上累积相乘,即可得到全样本的偏似然函数: \[ L(\beta) = \prod_{i=1}^K \frac{\exp(x_{(i)}^T \beta)}{\sum_{k \in \mathcal{R}(t_i)} \exp(x_k^T \beta)} \] 其中 \(x_{(i)}\) 就是命运指针抽中在时的\(t_i\) 发生失败(如退市)的那个“倒霉”候选人的特征向量。
对该函数取对数,得到对数偏似然函数 \(\ell(\beta)=\log L(\beta)\) ,可以用牛顿—拉弗森等数值方法估计 \(\hat{\beta}\) 。半参数设计避免预先指定基线风险形状,但其适用性仍依赖比例风险、删失和风险集定义;是否优于其他生存模型必须由相同数据与评价设计决定。
准备数据:point-in-time landmark 设计
Cox 模型允许加入连续变量,但报告期末不是信息公开日。本节与第 12—13 章统一采用 point-in-time 版本合同:课程运行器必须提供含可靠 info_date 的冻结财务切片、合同版本、内容哈希和冻结日;若缺 info_date、哈希不匹配或任一记录的可得日晚于其 landmark,流程必须结构化停止。不得再用“报告期末 + 固定天数”替代披露日,也不得把后续修订回填到较早 landmark。
下面把逐公司记录转换为稳定表,并单独执行 point-in-time 断言;lst-cox-landmark-extraction 与 lst-cox-landmark-validation 必须在同一 fresh kernel 中依次运行。
上方现场输出给出实际样本量、协变量分布和缺失处理结果。这里 Landmark_Date 就是所选记录的可靠 info_date,代码已经断言 available_date <= landmark(取等号);解释 Cox 系数前还应检查极端杠杆值与修订版本分布。正文不固定记录可能随切片版本变化的样本数。
拟合 Cox 比例风险模型
在完成数据预处理后,我们使用 CoxPHFitter 拟合带 L2 惩罚的 Cox 模型。事件是任何原因退市,不是财务困境;惩罚模型的常规 p 值只作描述,不能按未惩罚模型的精确推断方式解释。行业、地区和杠杆的关系是条件关联,不是因果效应。
Cox 表报告 landmark 样本中的条件关联。带惩罚估计的系数和常规 p 值不应机械套用 5% 阈值;风险比只描述比例风险假设成立时、其他已纳入协变量相同的相对退市率。若诊断拒绝比例风险假设,应考虑分层或时变效应,并停止给出全时段不变的风险比。
结果解释
Cox 模型的回归系数可以通过森林图展示。图中的系数描述协变量与“任何原因退市”风险率的条件关联;惩罚拟合下的区间与 p 值只作描述,不应称为财务困境因果效应。
图 11.4 用于描述样本中的条件关联,不能揭示因果“核心因子”。行业、地区和杠杆系数还可能受进入 landmark 风险集的选择、退市原因混合及遗漏变量影响;具体数值以当前执行结果为准。
风险函数与累计风险(Hazard Function)
风险函数(也称为风险率)定义为:
\[ h(t) = \lim_{\Delta t \to 0} \frac{\Pr(t < T \leq t + \Delta t | T > t)}{\Delta t} \tag{11.5}\]
这是在给定生存到时间 \(t\) 的条件下,在紧接着的瞬间发生事件的瞬时率 。
风险函数 vs 生存函数
生存函数 \(S(t)\) :回答生存到时间\(t\) 的概率是多少。
风险函数 \(h(t)\) :回答如果已经生存到时间\(t\) ,在下一瞬间发生事件的概率强度是多少。
关键关系。 \[ S(t) = \exp\left(-\int_0^t h(u) du\right) \tag{11.6}\]
中国案例:杠杆分组的累计风险
如果说生存函数是在俯瞰全局的“存活概率”,那么风险函数(Hazard Function) 就是在凝视每一个当下的“死亡威胁”。 下面使用 NelsonAalenFitter 估计累计风险 \(H(t)\) ,并按 landmark 杠杆率分组。累计风险的跳跃不是瞬时风险率 \(h(t)\) ;若要估计 \(h(t)\) ,需要额外平滑与带宽选择。 Nelson–Aalen 曲线是累计风险的阶梯估计,跳跃表示累计事件强度增加,并不是某一时点的平滑瞬时风险。曲线本身不能识别制度改革、估值或公司基本面机制;这些解释需要预先规定的协变量、时间模型和额外证据。
图 11.5 展示各组从 landmark 起累积的退市事件强度。曲线越高表示截至该时点累计风险越大,但它既不是瞬时风险率,也不能单凭组间差异作因果解释。若曲线差异随时间明显改变,应进一步用残差诊断或时间交互项检验比例风险假设。
Cox 模型的正则化
当特征数量较多时,我们可以使用正则化 Cox 模型进行变量选择。
当我们赋予Cox 比例风险模型处理数百个财务指标的能力时,多重共线性和过拟合的阴影便会随之降临。幸运的是,我们在第 6 章学过的 Lasso 正则化同样可以完美嫁接到生存分析中。 下面的这段“抗噪压力测试”代码非常有趣:我们故意在原本真实的上市财务数据中,人为地注入了 5 列完全由随机数构成的纯噪音特征(Noise_0 了Noise_4)。接着,我们利了lifelines 库开启了 L1 惩罚(l1_ratio=1.0),让惩罚力度(Lambda)从极小逐渐呈指数级放大,并记录下每一个特征对应系数的挣扎轨迹。 路径图仅展示惩罚增大时各系数如何收缩。在有限样本、相关特征与一次随机噪声注入下,噪声变量不保证比业务变量更早归零,未归零变量也不等于真实风险因子。该图是惩罚机制演示;若要做筛选,还需在时间外验证中选惩罚强度,并用重抽样入选频率评估稳定性。
from sklearn.preprocessing import StandardScaler # 导入标准化工具用于特征缩放
from lifelines import CoxPHFitter # 导入Cox比例风险模型拟合器
# 增加一些噪声特征来测试 Lasso 的筛选能力
np.random.seed(1 ) # 设定随机种子确保可复现性
for i in range (5 ): # 循环生成5个噪声特征列
cox_model_data[f'Noise_ { i} ' ] = np.random.normal(0 , 1 , len (cox_model_data)) # 生成标准正态分布随机噪声
predictor_matrix = cox_model_data.drop(['duration' , 'event' ], axis= 1 ) # 提取协变量矩阵(排除响应变量)
predictor_matrix = predictor_matrix.select_dtypes(include= [np.number]) # 仅保留数值型列
# 标准化特征(使各变量均值为0、标准差为1)
scaler = StandardScaler() # 实例化标准化器
scaled_predictor_matrix = scaler.fit_transform(predictor_matrix) # 拟合并转换特征矩阵
scaled_predictor_df = pd.DataFrame(scaled_predictor_matrix, columns= predictor_matrix.columns) # 转回DataFrame保留列名
scaled_predictor_df['duration' ] = cox_model_data['duration' ].values # 将存续时间列加回
scaled_predictor_df['event' ] = cox_model_data['event' ].values # 将事件标记列加回
以下代码使用带L1惩罚的Cox回归(Lasso正则化),通过逐步增大惩罚参数来观察各协变量系数的衰减路径。
图 11.6 展示一次样本、一次噪声注入下的惩罚路径。随着惩罚增大,系数趋向零;但某个噪声变量较早归零并不能证明选择一致性,也不能把其余变量认证为“真实信号”。稳定筛选需要重复生成或重抽样并报告入选频率。
生存树模型
生存树可以自动发现非线性关系和交互作用。
生存树可以用分割表示非线性与交互,但仍需要规定深度、叶节点样本量并做样本外验证。下面用 scikit-survival 拟合一棵限深树,再用置换重要性描述当前样本中特征与预测表现的关系。重要性没有方向,相关特征会分享或替代重要性,且本练习的训练样本重用结果不是因果“话语权”或外部预测证据。
# 导入生存树模型和数据结构化工具
from sksurv.tree import SurvivalTree # 导入生存树模型
from sksurv.util import Surv # 导入生存数据结构化工具
# 准备生存树所需的结构化目标变量
survival_targets = Surv.from_arrays( # 构造生存目标数组
event= cox_model_data['event' ].astype(bool ), # 事件标记转为布尔型
time= cox_model_data['duration' ] # 存续时间
) # 完成构建
tree_predictor_matrix = cox_model_data.drop(['duration' , 'event' ], axis= 1 ) # 提取特征矩阵
# 拟合生存树模型(限制树深为3以防止过拟合)
survival_tree_model = SurvivalTree(max_depth= 3 , min_samples_split= 10 , min_samples_leaf= 5 , random_state= 42 ) # 配置生存树参数
survival_tree_model.fit(tree_predictor_matrix, survival_targets) # 拟合模型
# 计算并排序特征重要性(使用置换重要性替代不可用的impurity-based方法)
from sklearn.inspection import permutation_importance # 导入置换重要性工具
perm_importance_result = permutation_importance( # 计算置换重要性
survival_tree_model, tree_predictor_matrix, survival_targets, # 执行数据处理操作
n_repeats= 10 , random_state= 42 # 重复10次取平均
) # 完成构建
calculated_feature_importance = pd.DataFrame({ # 构建DataFrame存储特征重要性
'feature' : tree_predictor_matrix.columns, # 特征名称
'importance' : perm_importance_result.importances_mean # 置换重要性均值
}).sort_values('importance' , ascending= False ) # 按重要性降序排列
在完成生存树模型的拟合和特征重要性计算后,下面输出特征重要性排名并以水平条形图可视化展示,直观比较各财务指标对企业生存时间的预测贡献。
图 11.7 的重要性是当前训练样本和模型设定下的描述量。排序应从现场输出读取,并通过未来测试期置换重要性检查;它不能单独证明某变量具有经济主导性、因果作用或稳定预测能力。
模型评估:C-index
使用 C-index 评估模型对“谁先退市”的风险排序能力。
构建完这些花哨的生存模型后,我们该如何评价它们的好坏呢?传统的准确率由于“删失”特征的存在而彻底失效,于是统计学家祭出了专属的裁判指标:C-index(一致性指数,Concordance Index) 。 C-index 比较可比样本对的风险排序;0.5 附近表示排序能力接近随机。为避免训练内评价的乐观偏差,下面按 landmark 日期划分较早训练队列与较晚测试队列,并只在后者报告结果。该指标评价的是当前事件定义下的排序,不代表概率校准、财务困境预测或因果效应。
请以代码实际输出为准,并同时报告切分日期与测试队列规模。单次时间队列切分仍可能不稳定;严谨研究还应预先固定特征、事件口径和评价窗口,并用多个后续队列或外部样本复核。
案例:股价回撤恢复时间分析
在投资风险管理中的最大回撤 (Max Drawdown)是一个关键指标。我们不仅关心回撤的深度,还关心回撤恢复时间 (Time to Recovery),即股价从跌破高点到重创新高所需的时间。这可以看作是一个生存分析问题。
事件 :股价创新高(恢复)。
时间 :从上一次创新高到下一次创新高的持续时间。
删失 :截至目前仍未创新高。
我们将分析海康威视(002415.XSHE)的历史回撤恢复情况。
在这个章节的最后,我们要打破思维的桎梏:谁说生存分析只能用来研究“死亡”与“退市”?在充满奇迹的交易桌上,一次灾难性的暴跌也可以被视为一次“死亡”,而股价重新跌回并突破前高,则是一场漫长的“复活”。 下面的代码,将生存分析的矛头直指向了一个令无数股民夜不能寐的金融难题:“买在最高点被套牢,到底需要熬多少天才解套?”(即最大回撤恢复时间分析) : 我们以智能安防龙头——海康威视的全历史日线数据为手术刀。代码首先计算了股价的历史阻力位与实时回撤幅度,然后用近乎苛刻的条件找出了每一次超了5 天的显著“套牢解套”周期,这构成了我们的生存(持续套牢)时间,而股价重新突破前高就是那激动人心的最终“事件”。同时,对于当前依然被套牢在山顶无法解脱的区间,代码极其专业地将其标记为了“右删失”。 最后输出的这张生存曲线以及底部的硬核文字报告,不再是一份冰冷的医学死亡率,而是一份实打实的交易心理学防线评估:它不仅告诉你海康威视中位数的解套天数,更无情地向你宣告了如果你在错误的时间冲进去,在接下来的半年(180天)内还要继续承受无期徒刑般煎熬的概率。
# 在真实 HDF 存储层直接选择海康威视,避免为单股恢复期示例载入全市场行情
drawdown_price_path = Path(os.environ['BOOK_DATA_DIR' ]).expanduser().resolve() / 'stock' / 'stock_price_post_adjusted.h5' # 独立定位后复权股价文件
haikang_data = pd.read_hdf( # 选择性读取海康威视的全部真实历史行情
drawdown_price_path, # 指定后复权股价 HDF 文件
where= "order_book_id='002415.XSHE'" , # 在存储层限定唯一教学标的
columns= ['date' , 'order_book_id' , 'close' ], # 只读取恢复时间所需字段
).reset_index() # 恢复日期与证券索引列并重建唯一行索引
haikang_data = haikang_data.sort_values('date' ).set_index('date' ) # 按日期排序并设为索引
closing_prices = haikang_data['close' ] # 提取收盘价序列
# 计算滚动最高价与回撤幅度
running_maximum_prices = closing_prices.cummax() # 计算历史累计最高价
drawdown_series = (closing_prices - running_maximum_prices) / running_maximum_prices # 计算回撤百分比
# 找出所有创新高的日期
is_new_high = (closing_prices >= running_maximum_prices) # 标记当日是否创历史新高
new_high_dates = closing_prices.index[is_new_high] # 提取所有创新高的日期
基于上述回撤计算结果,我们接下来识别每一次显著的回撤-恢复周期(持续时间超过5个交易日的回撤),并计算其持续天数:
# 遍历所有创新高的相邻日期对,计算回撤恢复时间
recovery_durations = [] # 存储每次回撤的持续天数
recovery_events = [] # 存储事件标记(1=已恢复, 0=删失/未恢复)
if len (new_high_dates) > 1 : # 至少需要两个新高日期才能计算间隔
for i in range (len (new_high_dates) - 1 ): # 遍历每对相邻新高日期
start_date = new_high_dates[i] # 本次创新高日期
end_date = new_high_dates[i+ 1 ] # 下次创新高日期
duration_days = (end_date - start_date).days # 计算间隔天数
if duration_days > 5 : # 仅保留超过5天的显著回撤
recovery_durations.append(duration_days) # 记录回撤持续天数
recovery_events.append(1 ) # 标记为已恢复事件
# 处理最后一段(当前是否仍在回撤中)
last_high_date = new_high_dates[- 1 ] # 最后一次创新高日期
last_recorded_date = closing_prices.index[- 1 ] # 数据最后一天
if last_recorded_date > last_high_date: # 若最后一天不是新高,说明当前仍在回撤中
duration_days = (last_recorded_date - last_high_date).days # 计算当前回撤持续天数
if duration_days > 5 : # 仅保留显著回撤
recovery_durations.append(duration_days) # 记录当前回撤天数
recovery_events.append(0 ) # 标记为删失(尚未恢复)
drawdown_recovery_data = pd.DataFrame({'duration' : recovery_durations, 'event' : recovery_events}) # 构建回撤恢复数据框
print (f'识别到 { len (drawdown_recovery_data)} 次显著回撤周期' ) # 输出回撤周期总数
recovered_mean_days = drawdown_recovery_data.loc[drawdown_recovery_data['event' ]== 1 , 'duration' ].mean() # 计算已恢复周期的平均天数
print (f'平均恢复时间: { recovered_mean_days:.1f} 天' ) # 输出平均恢复天数
识别到 38 次显著回撤周期
平均恢复时间: 100.7 天
回撤周期数量与恢复时间以当前数据快照的代码输出为准。这里的周期由“创新高—再次创新高”规则构造,彼此并非天然独立,且价格路径、复权方式和样本终点都会改变结果;因此它是描述性练习,不是对未来解套时间的承诺。
我们已识别出所有显著回撤周期的生存数据。下面使用 Kaplan-Meier 方法估计回撤恢复时间的生存曲线,并标注中位恢复时间:
# 拟合 Kaplan-Meier 生存曲线
kaplan_meier_fitter = KaplanMeierFitter() # 实例化KM估计器
kaplan_meier_fitter.fit(drawdown_recovery_data['duration' ], drawdown_recovery_data['event' ], label= '回撤恢复' ) # 拟合模型
fig, ax = plt.subplots(figsize= (10 , 6 )) # 创建10×6英寸画布
kaplan_meier_fitter.plot(ax= ax) # 绘制KM生存曲线
# 添加中位恢复时间参考线
median_recovery_time = kaplan_meier_fitter.median_survival_time_ # 获取中位恢复时间
plt.axvline(median_recovery_time, color= 'r' , linestyle= '--' , label= f'中位恢复时间: { median_recovery_time:.0f} 天' ) # 绘制竖直参考线
ax.set_title('股价回撤恢复时间生存曲线 (Time to Recovery)' , fontsize= 14 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置图标题
ax.set_xlabel('天数' , fontsize= 12 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置 x 轴标签
ax.set_ylabel('未恢复概率' , fontsize= 12 , fontproperties= 'Source Han Serif SC' ) # 使用已安装中文字体设置 y 轴标签
ax.grid(True , alpha= 0.3 ) # 添加网格线
plt.legend(prop= {'family' : 'Source Han Serif SC' }) # 使用已安装中文字体设置图例
plt.tight_layout() # 自动调整布局
plt.show() # 显示图形
# 预测特定时间点的未恢复概率
prob_not_recovered_30d = kaplan_meier_fitter.predict(30 ) # 30天内未恢复的概率
prob_not_recovered_180d = kaplan_meier_fitter.predict(180 ) # 180天内未恢复的概率
print (f'30天内无法恢复的概率: { prob_not_recovered_30d:.1%} ' ) # 输出30天未恢复概率
print (f'半年(180天)内无法恢复的概率: { prob_not_recovered_180d:.1%} ' ) # 输出180天未恢复概率
30天内无法恢复的概率: 44.7%
半年(180天)内无法恢复的概率: 13.2%
图 11.8 描述当前样本与所选回撤定义下“尚未恢复”的经验分布。30 天、180 天及中位恢复时间应直接从本次拟合对象读取;删失、重叠回撤和单一标的限制使其不能直接转化为止损规则、系统性风险归因或基本面判断。
回撤恢复支线审计任务
本任务为 60 分钟、10 分的可独立描述性课堂支线。它使用与正文海康威视案例不同的 TEACHING_SECURITY_B、8 个营业日纳入阈值和 20/60 个营业日评价窗口;fixture 是无随机数的人工教学轨迹,不代表任何真实证券、交易所行情或投资结果。
consumer 与独立证据
独立运行时先执行 列表 11.4 。consumer 只能接收 ch11_audit_fixture,先断言版本、列集合、唯一 (security_id,date) 键、日期严格递增、价格为正和 source_status == 'synthetic_teaching_fixture';随后用 adjusted_close.cummax() 找创新高,以相邻创新高间隔至少 8 个营业日定义已恢复周期,末端未恢复区间编码为右删失。不得读取 closing_prices、drawdown_recovery_data 或正文的 kaplan_meier_fitter。
0—15 分钟提交 ch11_event_contract,登记时间原点、恢复事件、行政删失、8 日阈值和 20/60 日窗口;15—30 分钟提交 ch11_period_audit,逐周期报告稳定 period_id、起止日、duration_business_days、event 与每个事件时点风险集;30—45 分钟提交 ch11_km_evidence,包括周期数、事件数、删失数、KM 表及可定义性;45—60 分钟提交 ch11_window_summary 与 ch11_limitations,报告中位数、20/60 日未恢复概率,以及相邻周期依赖、单一人工轨迹和非交易结论。
评分为:事件/删失合同 3 分,逐周期与风险集 3 分,KM 独立输出 2 分,窗口摘要、限制与非真实数据边界 2 分。缺版本/哈希/schema 守卫或事件合同,结构化失败且最高 4 分;风险集不可复算,第二项为 0 分;把曲线写成真实证券发现、未来承诺或交易建议,结论项为 0 分。此支线的证据必须由 ch11-recovery-audit-v1 生成,正文 图 11.8 的图形和数值不能充当提交物。
本章小结
本章我们系统学习了生存分析的核心方法:
Kaplan-Meier 估计量 :非参数估计生存曲线
对数秩检验 :比较两组或多组生存曲线
Cox 比例风险模型 :半参数回归模型,评估协变量对生存的影响
模型评估 :使用C-index 评估预测性能
正则化 :对 Cox 模型应用 Lasso/Ridge 正则化
生存分析的应用建议
医学研究 :患者生存时间、无病生存时间
商业分析 :客户流失时间、产品使用寿命
工程质量 :设备故障时间、系统可靠性
人力资源 :员工离职时间、晋升时间
关键假设 :
独立删失(数据收集时要仔细设计)
比例风险(Cox 模型需要检验)
软件工具 :
Python: lifelines, scikit-survival
R: survival, survminer
理论来源与前沿
生存分析最早在可靠性工程 和医学随访研究 中系统化发展:前者关注设备故障时间,后者关注患者生存时间。Kaplan-Meier 的乘积极限估计提供了在删失存在时对生存函数\(S(t)\) 的非参数估计;Cox 模型则用偏似然把协变量效应与基线风险函数 \(h_0(t)\) 分离,使得回归分析在不指定\(h_0(t)\) 的情况下仍然可行。
近十年的研究前沿主要集中在三类问题上:
高维协变量与正则化 :当特征数远大于样本量时,需要用 Lasso/Elastic Net 等正则化构建稀疏Cox 模型,并配套稳定的变量选择与不确定性量化。
时间变化效应与非比例风险 :在真实商业场景(例如客户流失)中,某些协变量的影响往往随时间变化,此时需要扩展Cox 模型(时间交互项、分段比例风险、AFT 模型等),并用Schoenfeld 残差等方法诊断。
机器学习与因果推断的结合 :将生存模型与树模型、Boosting、深度学习以及因果推断框架结合,用于处理复杂非线性、异质性处理效应与动态干预策略。
练习
概念题
[核心|难度:1|时间:8分钟|分值:4|项目:无] 解释什么是右删失、左删失和区间删失,并各举一个实际例子。
[核心|难度:2|时间:10分钟|分值:5|项目:无] Kaplan-Meier 估计量的核心思想是什么?为什么它比简单地计算生存比例更合理?
[核心|难度:2|时间:10分钟|分值:5|项目:无] Cox 比例风险模型中的“比例风险”假设是什么意思?如何检验这个假设?
[核心|难度:2|时间:10分钟|分值:5|项目:无] 解释风险函数 \(h(t)\) 和生存函数\(S(t)\) 之间的关系。
[拓展|难度:2|时间:10分钟|分值:5|项目:无] 在什么情况下独立删失假设会被违反?这会对分析结果造成什么影响?
应用题
[核心|难度:2|时间:15分钟|分值:8|项目:无] 内嵌小表的 Kaplan–Meier 手算(必做) :10 个产品在第 2 月发生 2 个事件,随后 1 个删失;第 5 月发生 1 个事件;第 7 月有 2 个删失;其余产品在第 10 月研究结束时删失。列出每个时点的风险集、事件数和删失数,计算 \(\widehat S(2)\) 、\(\widehat S(5)\) ,并说明第 7 月为何不产生生存曲线下降。
lung 或 colon 数据可作为选做扩展,但不是完成本题的外部依赖。 ### 应用于(金融场景)
[拓展|难度:3|时间:40分钟|分值:15|项目:无] 使用 financial_statement.h5 数据,选取两个行业,构建“从首份可得财报到首次亏损”的生存数据集。
要求:
绘制两个行业的Kaplan-Meier 生存曲线。
使用对数秩检验判断两条曲线是否有显著差异。
解释结果:哪个行业的“财务安全期”更长?
[拓展|难度:3|时间:40分钟|分值:15|项目:无] 模拟信贷违约数据 : 假设你是一家网贷平台的风控分析师,请模拟一份贷款数据(包含借款人年龄、收入、信用分、借款金额等特征)。
要求:
设定真实的违约机制(例如低信用分、高杠杆导致风险指数增加)。
添加删失(例如贷款尚未到期)。
使用带Lasso 正则化的 Cox 模型筛选关键风险因子。
[拓展|难度:3|时间:40分钟|分值:15|项目:无] 综合案例:A股新股生存分析 利用 stock_basic_data.h5 中的 listed_date 和 de_listed_date(退市日期)。
要求:
定义事件:公司退市。
定义删失:公司截至目前仍在交易。
比较不同板块(主板vs 创业板vs 科创板)的退市风险(生存曲线)。
若退市事件过少,应报告估计不稳定并合并教学分组;不得改用缺少生效日期的静态 ST 状态替换事件。
理论题
[核心|难度:3|时间:20分钟|分值:10|项目:无] 给定三个无并列事件时点及其风险集,写出 Cox 偏似然的一般乘积、score 与观察信息的作用;再说明为什么基线风险相消、风险比不等于绝对风险,以及时间外 C-index 在零事件或全删失测试集为何不可定义。
[拓展|难度:3|时间:25分钟|分值:12|项目:无] 推导 Cox 模型的偏似然函数,并解释为什么不需要估计基线风险\(h_0(t)\) 。
[拓展|难度:3|时间:25分钟|分值:12|项目:无] 证明在单变量二元协变量的情况下,Cox 模型的score test 等价于对数秩检验。
练习参考解答
展开完整解答、评分点与常见失败模式
评分以题面分值为准;概率与风险集手算允许 \(10^{-3}\) 绝对误差。应用题若事件、时间原点、删失或观察终点任一未登记,最高得 50%;若把风险比写成绝对概率、静态状态写成事件日期或忽略零事件组,解释项不得分。
概念题解答
右删失、左删失与区间删失
右删失 (right censoring) :只知道事件时间 \(T\) 大于观测结束时刻 \(C\) 。例:研究结束时公司仍未退市。
左删失 (left censoring) :只知道事件时间 \(T\) 小于某个观测起点。例:只知道公司在数据收集前已经违约,但不知确切日期。
区间删失 (interval censoring) :只知道 \(T\) 落在区间 \((L, R]\) 。例:只知道违约发生在两次财报披露之间。
Kaplan-Meier 的核心思想
将总体生存函数分解为一连串条件生存概率的乘积: \[ S(t) = \prod_{t_j \le t} P(T > t_j \mid T \ge t_j) \] 在每个事件时刻 \(t_j\) ,用风险集大小 \(n_j\) 与事件数 \(d_j\) 估计条件概率为 \(1 - d_j/n_j\) 。
比例风险假设与检验
Cox 模型假设两组个体的风险比与时间无关。常用检验:Schoenfeld 残差与时间的相关性检验。若存在显著相关性,说明违反比例风险假设,应考虑含时间依赖系数的 Cox 模型。
风险函数与生存函数
\[ S(t) = \exp\Big(-\int_0^t h(u)\,du\Big) \]
独立删失
若删失机制包含有关 \(T\) 的信息(例如财务状况恶化导致的数据缺失),则不再独立,标准方法失效。
应用题解答
内嵌小表手算
第 2 月事件前 \(r_1=10,d_1=2\) ,所以 \(\widehat S(2)=8/10=0.8\) ;同月随后 1 个删失,故第 5 月事件前 \(r_2=7,d_2=1\) ,\(\widehat S(5)=0.8(6/7)\approx0.686\) 。第 7 月只有删失,条件生存概率没有事件因子,曲线保持 0.686,但第 7 月后风险集减少 2。第 10 月其余对象行政删失,曲线仍不下降。
行业首次亏损比较(实现契约)
在公司内按可靠 info_date 排序,将首次 net_profit < 0 的可得日定义为事件日;从首份可得财报开始计时,研究截止仍未亏损者右删失。下面复用 表 11.3 的冻结切片与哈希,不允许固定日数代理;selected_industries 是事前指定的比较族,不按结果选择。
import os
import hashlib
import pandas as pd
from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test
from pathlib import Path # 使用跨平台路径对象解析显式切片
financial_path = Path(os.environ['BOOK_M11_FINANCIAL_SLICE' ]).expanduser().resolve() # 读取同一point-in-time切片接口
expected_hash = os.environ['BOOK_M11_FINANCIAL_SLICE_SHA256' ].strip().lower() # 读取冻结内容哈希
if hashlib.sha256(financial_path.read_bytes()).hexdigest() != expected_hash: # 读取HDF前校验版本内容
raise RuntimeError ({'status' : 'stopped' , 'reason' : 'input_sha256_mismatch' }) # 哈希不符即停止
financial = pd.read_hdf(financial_path).copy() # 读取已验证的小切片
DATA_DIR = Path(os.environ['BOOK_DATA_DIR' ]).expanduser().resolve() # 基本信息仍由显式本地根提供
stock_basic_data = pd.read_hdf(DATA_DIR / 'stock' / 'stock_basic_data.h5' ).copy()
if 'info_date' not in financial.columns: # 禁止缺失披露日时继续
raise RuntimeError ({'status' : 'stopped' , 'reason' : 'missing_info_date' }) # 缺时点即停止
financial['available_date' ] = pd.to_datetime(financial['info_date' ], errors= 'coerce' ) # 以可靠披露或修订日为可得日
if financial['available_date' ].isna().any (): # 禁止不可解析时点进入风险集
raise RuntimeError ({'status' : 'stopped' , 'reason' : 'invalid_info_date' }) # 无法解析即停止
observation_end_date = financial['available_date' ].max ().normalize()
financial = financial.sort_values(['order_book_id' , 'available_date' ])
first_available = financial.groupby('order_book_id' )['available_date' ].min ().rename('origin' )
first_loss = (financial[financial['net_profit' ].lt(0 )]
.groupby('order_book_id' )['available_date' ].min ().rename('event_date' ))
图 11.9 的行业曲线必须从现场输出读取;它依赖当前冻结切片与 info_date 版本,不能证明行业因果差异。检验输出单独见 表 11.9 。
补充审计:旧的任何原因退市行业模板(不构成第 7 或第 9 题答案)
#| output: false
import pandas as pd # 导入pandas用于数据处理
from lifelines import KaplanMeierFitter # 导入Kaplan-Meier生存函数估计器
from lifelines.statistics import logrank_test # 导入对数秩检验用于比较生存曲线差异
import matplotlib.pyplot as plt # 导入matplotlib用于数据可视化
import os # 导入操作系统模块用于跨平台路径处理
# ---------- 1. 数据准备(复用正文 Cell 1 已生成的 company_survival_data) ----------
# company_survival_data 包含 'industry'、'duration'、'event' 列
# 'industry' 列的值来自 citics_2019_l1_name(中信一级行业分类)
# 动态获取样本量前两大的行业,而非硬编码行业名称
industry_value_counts = company_survival_data['industry' ].value_counts() # 统计各行业样本数
industry_first_name = industry_value_counts.index[0 ] # 样本量最多的行业
industry_second_name = industry_value_counts.index[1 ] # 样本量第二多的行业
print (f'选取的两个行业: { industry_first_name} , { industry_second_name} ' ) # 打印实际行业名称
# ---------- 2. 构建行业筛选掩码 ----------
mask_first_industry = company_survival_data['industry' ] == industry_first_name # 第一个行业的布尔掩码
mask_second_industry = company_survival_data['industry' ] == industry_second_name # 第二个行业的布尔掩码
# ---------- 3. 绘制两个行业的 KM 生存曲线并进行 Log-Rank 检验 ----------
fig_exercise, ax_exercise = plt.subplots(figsize= (10 , 6 )) # 创建画布
kaplan_meier_fitter_exercise = KaplanMeierFitter() # 实例化 KM 估计器
# 拟合并绘制第一个行业
kaplan_meier_fitter_exercise.fit( # 训练/拟合模型
company_survival_data.loc[mask_first_industry, 'duration' ], # 第一行业的存续时间
company_survival_data.loc[mask_first_industry, 'event' ], # 第一行业的事件标识
label= industry_first_name # 图例标签
) # 完成构建
kaplan_meier_fitter_exercise.plot(ax= ax_exercise) # 绘制到画布上
完成第一个行业的 Kaplan-Meier 曲线绘制后,接下来对第二个行业执行相同的拟合与绘制流程,并通过 Log-Rank 检验比较两个行业生存曲线之间的差异是否具有统计学显著性。
# 拟合并绘制第二个行业
kaplan_meier_fitter_exercise.fit( # 训练/拟合模型
company_survival_data.loc[mask_second_industry, 'duration' ], # 第二行业的存续时间
company_survival_data.loc[mask_second_industry, 'event' ], # 第二行业的事件标识
label= industry_second_name # 图例标签
) # 完成构建
kaplan_meier_fitter_exercise.plot(ax= ax_exercise) # 绘制到画布上
# 进行 Log-Rank 检验,比较两个行业生存曲线的差异是否显著
rank_test_result = logrank_test( # 提取测试集数据
company_survival_data.loc[mask_first_industry, 'duration' ], # 第一行业持续时间
company_survival_data.loc[mask_second_industry, 'duration' ], # 第二行业持续时间
company_survival_data.loc[mask_first_industry, 'event' ], # 第一行业事件
company_survival_data.loc[mask_second_industry, 'event' ] # 第二行业事件
) # 完成构建
ax_exercise.set_title(f'行业生存曲线比较: { industry_first_name} vs { industry_second_name} ' ) # 设置标题
ax_exercise.set_xlabel('存续年数' ) # 设置 x 轴标签
ax_exercise.set_ylabel('生存概率' ) # 设置 y 轴标签
plt.tight_layout() # 自动调整布局
plt.show() # 展示图形
print (f'Log-Rank 检验 p 值: { rank_test_result. p_value:.4f} ' ) # 输出检验 p 值
上述代码动态选取两个行业并输出 Log-Rank 检验。是否存在统计证据必须按现场 p 值和预先阈值判断;它是任何原因退市的辅助练习,不是“首次亏损”答案。
信贷违约模拟与 Lasso-Cox(完整实现)
本题明确允许模拟,因为目标是检查已知真机制下的变量筛选。设信用分与收入降低风险,杠杆提高风险,年龄与三个噪声变量的真系数为零;独立行政删失并不依赖事件时间。
按样本编号锁定前 70% 为开发集、后 30% 为验证集。这里只在验证集选择惩罚强度;验证分数相同到 \(10^{-4}\) 时选择更大的惩罚以偏向稀疏。
评分证据包括真系数表、删失率、惩罚网格、验证 C-index 和全部系数。允许估计误差,但应检查三项:信用分/收入方向为负、杠杆方向为正、无关项整体收缩更强。若验证集事件少于 20、只含一个事件类别或模型不收敛,应扩大预注册样本而不是报告筛选成功。
A股板块任何原因退市风险
# ---------- A股不同上市板块的退市风险比较:数据准备 ----------
import os # 导入操作系统模块用于跨平台路径处理
import pandas as pd # 导入pandas用于数据处理
from lifelines import KaplanMeierFitter # 导入KM估计器
from lifelines.statistics import logrank_test # 导入对数秩检验
import matplotlib.pyplot as plt # 导入matplotlib用于可视化
from pathlib import Path # 使用跨平台路径对象解析显式数据根
DATA_DIR = Path(os.environ['BOOK_DATA_DIR' ]).expanduser().resolve() # 从必需环境变量取得数据根
if not DATA_DIR.is_dir(): # 在读取前验证数据根
raise FileNotFoundError (f'BOOK_DATA_DIR 不存在或不是目录: { DATA_DIR} ' ) # 失败即停止,不回退固定路径
stock_basic_data_ex = pd.read_hdf(DATA_DIR / 'stock' / 'stock_basic_data.h5' ) # 加载上市公司基本信息
# 将上市和退市日期转换为 datetime 格式
stock_basic_data_ex['listed_date' ] = pd.to_datetime(stock_basic_data_ex['listed_date' ], errors= 'coerce' ) # 上市日期
stock_basic_data_ex['de_listed_date' ] = pd.to_datetime(stock_basic_data_ex['de_listed_date' ], errors= 'coerce' ) # 退市日期
stock_basic_data_ex = stock_basic_data_ex.dropna(subset= ['listed_date' ]) # 去除缺少上市日期的记录
observation_end = pd.to_datetime('2024-01-01' ) # 研究观测截止日期
在加载基本信息数据后,下一步根据股票代码后缀规则将各股票归入对应的上市板块(主板、创业板、科创板),以便后续按板块维度进行生存分析比较。
# 根据股票代码后缀判断上市板块
def classify_board(order_book_id): # 定义函数classify_board
"""根据股票代码判断所属板块""" # 执行数据处理操作
code = str (order_book_id).split('.' )[0 ] # 提取纯数字代码
if code.startswith('688' ): # 科创板
return '科创板' # 返回结果
elif code.startswith('300' ): # 创业板
return '创业板' # 返回结果
elif code.startswith('60' ): # 上交所主板
return '主板(沪)' # 返回结果
elif code.startswith('00' ): # 深交所主板
return '主板(深)' # 返回结果
else : # 默认分支
return '其他' # 返回结果
stock_basic_data_ex['board_type' ] = stock_basic_data_ex['order_book_id' ].apply (classify_board) # 添加板块分类列
stock_basic_data_ex = stock_basic_data_ex[stock_basic_data_ex['board_type' ] != '其他' ] # 排除无法分类的股票
基于上述板块分类结果,接下来构建各股票的生存数据并绘制 Kaplan-Meier 生存曲线:
# ---------- 构建生存数据并绘制各板块KM曲线 ----------
board_survival_records = [] # 存储各股票的生存记录
for _, row in stock_basic_data_ex.iterrows(): # 遍历迭代
listed = row['listed_date' ] # 上市日期
delisted = row['de_listed_date' ] # 退市日期
if pd.isna(delisted) or delisted > observation_end: # 缺失或截止日后退市才删失
dur = (observation_end - listed).days / 365.25 # 存续年数
evt = 0 # 事件标记为删失
else : # 截止日当日或之前退市均为事件
dur = (delisted - listed).days / 365.25 # 存续年数
evt = 1 # 事件标记为退市
if dur > 0 : # 排除无效记录
# 将计算结果追加到列表
board_survival_records.append({'board' : row['board_type' ], 'duration' : dur, 'event' : evt})
board_survival_df = pd.DataFrame(board_survival_records) # 转换为DataFrame
完成生存数据构建和 KM 估计器初始化后,下面遍历各板块分别拟合 Kaplan-Meier 模型并在同一画布上绘制生存曲线,从而直观比较不同板块上市公司的退市风险差异。
fig_board, ax_board = plt.subplots(figsize= (10 , 6 )) # 为板块曲线创建独立画布
km_board = KaplanMeierFitter() # 实例化板块曲线共用的KM估计器
for board_name in ['主板(沪)' , '主板(深)' , '创业板' , '科创板' ]: # 遍历各板块
mask_board = board_survival_df['board' ] == board_name # 筛选该板块数据
if mask_board.sum () < 2 : # 跳过样本量不足的板块
continue # 跳过本次循环
km_board.fit( # 拟合 KM 模型
board_survival_df.loc[mask_board, 'duration' ], # 存续时间
board_survival_df.loc[mask_board, 'event' ], # 事件标识
label= f' { board_name} (n= { mask_board. sum ()} )' # 图例含样本量
) # 完成构建
km_board.plot(ax= ax_board) # 绘制生存曲线
ax_board.set_title('A股不同板块上市公司生存曲线比较' ) # 设置标题
ax_board.set_xlabel('存续年数' ) # x 轴标签
ax_board.set_ylabel('生存概率' ) # y 轴标签
plt.tight_layout() # 自动调整布局
plt.show() # 展示图形
图 11.10 只展示当前快照下任何原因退市的描述性曲线;板块样本量与事件率另由 表 11.10 报告,避免把表格输出混入图形对象。
上述代码根据股票代码将 A 股公司分为主板(沪、深)、创业板和科创板四个板块。解释 表 11.10 时必须同时报告事件数;静态退市比例不等于生存概率,也不能替代 图 11.10 的风险集校正。
理论题解答
Cox 偏似然与时间外评价骨架
无并列事件时,若 \(i(j)\) 在 \(t_j\) 发生事件、风险集为 \(R_j\) ,则
\[L_p(\beta)=\prod_j\frac{\exp(x_{i(j)}^\top\beta)}{\sum_{k\in R_j}\exp(x_k^\top\beta)}.\]
score 是 \(U(\beta)=\partial\log L_p/\partial\beta\) ,观察信息是 \(I(\beta)=-\partial^2\log L_p/\partial\beta\partial\beta^\top\) ;前者给出估计方程,后者刻画局部曲率并用于标准误。条件在“风险集中谁发生事件”后,公共的 \(h_0(t_j)\) 从分子分母相消,因此偏似然识别相对风险参数,但不直接给出绝对生存概率。时间外 C-index 需要可比较对以及结果差异;测试期零事件、全删失或没有可比较对时必须报告不可定义,而不能填成 0.5。
Cox 偏似然的连续推导
设 \(h_i(t)=h_0(t)\exp(x_i^\top\beta)\) ,在无并列事件的 \(t_j\) ,已知风险集 \(R_j\) 中恰有一人发生事件,则事件人为 \(i(j)\) 的条件概率为
\[
\frac{h_{i(j)}(t_j)}{\sum_{k\in R_j}h_k(t_j)}
=\frac{h_0(t_j)e^{x_{i(j)}^\top\beta}}
{\sum_{k\in R_j}h_0(t_j)e^{x_k^\top\beta}}
=\frac{e^{x_{i(j)}^\top\beta}}{\sum_{k\in R_j}e^{x_k^\top\beta}}.
\]
基线风险在条件概率分子分母中相消。跨事件时点相乘得偏似然
\[
L_p(\beta)=\prod_{j=1}^{D}\frac{e^{x_{i(j)}^\top\beta}}{\sum_{k\in R_j}e^{x_k^\top\beta}},
\]
其对数为 \(\ell_p(\beta)=\sum_j[x_{i(j)}^\top\beta-\log\sum_{k\in R_j}e^{x_k^\top\beta}]\) 。因此可以先估计相对风险参数 \(\beta\) ;若要绝对生存概率,仍需另估累计基线风险。
二元协变量下 score test 与 log-rank 的等价
令 \(x_i\in\{0,1\}\) 。对上式求导:
\[
U(\beta)=\sum_j\left[x_{i(j)}-
\frac{\sum_{k\in R_j}x_ke^{\beta x_k}}{\sum_{k\in R_j}e^{\beta x_k}}\right].
\]
在 \(H_0:\beta=0\) 下,若 \(Y_{1j}=\sum_{k\in R_j}x_k\) 、\(Y_j=|R_j|\) ,且无并列事件,则
\[
U(0)=\sum_j\left[d_{1j}-\frac{Y_{1j}}{Y_j}\right]
=\sum_j(O_{1j}-E_{1j}),
\]
正是 log-rank 分子。负二阶导给出信息量;允许每时点 \(d_j\) 个事件时,其超几何方差为
\[
I(0)=\sum_j\frac{Y_{1j}Y_{0j}d_j(Y_j-d_j)}{Y_j^2(Y_j-1)}.
\]
所以 Cox score 统计量 \(U(0)^2/I(0)\) 与 log-rank 的卡方统计量相同(相同风险集与并列事件处理下)。若使用不同的 tie 近似、权重或时间变化效应,等价关系不再原样成立。
章末闭环
逐项目标自检 :定义生存时间、事件与删失;从风险集写出 Kaplan–Meier 乘积极限;解释 Cox 偏似然为何消去基线风险;检查比例风险和 landmark 时点;用 C-index 或预注册时间点损失比较模型。事件/删失合同、landmark 先于结果窗口及时间外评价可定义性是 must-pass;三项全过且其余两项至少一项有证据,才继续学习。
禁用情境 :事件起点、删失机制或风险集无法审计时,不应报告生存概率;比例风险明显不成立且未建模时间变化效应时,不应把单一 Cox 风险比作为全时期摘要。常见误区是把删失当普通缺失删除,或把风险比直接解释为概率差、因果效应或行动收益。
无提示检索 :1)删失时点是否在 Kaplan–Meier 乘积中产生事件因子?2)Cox 偏似然消去了什么、没有估计什么?3)训练和验证的 landmark 为什么必须先于各自结果窗口?
展开检索反馈与学习决策
1)不产生事件因子,但会改变后续风险集人数。2)条件化后基线风险相消,偏似然估计相对风险参数;绝对生存概率仍需累计基线风险。3)否则特征或入组会使用结果发生后的信息。第 1 题错:补修 小节 11.3 ,把第 7 月删失改到第 4 月并重算两次风险集;第 2 题错:补修 小节 11.4 ,用两个风险集写出分子分母并指出被约去项;第 3 题错:补修 小节 11.8 ,为“首次逾期”画 landmark、特征截止和测试窗口,要求严格不交叉。三项异形复测均通过才从 小节 11.11 重入。
陌生迁移 :为一个贷款提前还款、供应链恢复或客户流失任务,提交时间起点、事件、删失、landmark、风险集和失败条件。下一章转向没有响应变量时的结构探索。
Cox, D. R. 1972年.
《Regression Models and Life-Tables》 .
Journal of the Royal Statistical Society: Series B (Methodological) 34 (2): 187~220.
https://doi.org/10.1111/j.2517-6161.1972.tb00899.x .
Kaplan, E. L., 和 Paul Meier. 1958年.
《Nonparametric Estimation from Incomplete Observations》 .
Journal of the American Statistical Association 53 (282): 457~81.
https://doi.org/10.1080/01621459.1958.10501452 .