import scipy.integrate # 导入scipy积分模块以添加向后兼容性补丁
scipy.integrate.trapz = scipy.integrate.trapezoid # 为旧版lifelines添加trapz别名（scipy 1.14+已移除trapz）
import pandas as pd # 导入pandas用于数据处理
import numpy as np # 导入numpy用于数值计算
import matplotlib.pyplot as plt # 导入matplotlib用于可视化
import os # 导入os模块用于跨平台路径
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC'] # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
DATA_DIR = 'C:/qiufei/data' if os.name == 'nt' else '/home/ubuntu/r2_data_mount/qiufei/data' # 跨平台数据路径
stock_basic_data = pd.read_hdf(os.path.join(DATA_DIR, 'stock/stock_basic_data.h5')) # 读取上市公司基本信息
stock_basic_data['listed_date'] = pd.to_datetime(stock_basic_data['listed_date'], errors='coerce') # 转换上市日期格式
stock_basic_data['de_listed_date'] = pd.to_datetime(stock_basic_data['de_listed_date'], errors='coerce') # 转换退市日期格式
stock_basic_data = stock_basic_data.dropna(subset=['listed_date']) # 去除缺少上市日期的无效记录11 生存分析
A股整体退市生存曲线
Code
from lifelines import KaplanMeierFitter # 导入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)) # 创建画布
kaplan_meier_fitter.plot(ax=ax) # 绘制KM生存曲线
ax.set_title('A股上市公司退市生存曲线', fontsize=14) # 设置图标题
ax.set_xlabel('上市年数', fontsize=12) # 设置x轴标签
ax.set_ylabel('生存概率', fontsize=12) # 设置y轴标签
ax.grid(True, alpha=0.3) # 添加网格线
plt.tight_layout() # 自动调整布局
plt.show() # 显示图形
按行业分组的生存曲线
Code
industry_counts = company_survival_data['industry'].value_counts() # 统计各行业样本量
top_industry_categories = industry_counts.head(5).index.tolist() # 选取前5大行业
fig, ax = plt.subplots(figsize=(10, 6)) # 创建画布
for ind in top_industry_categories: # 遍历前5大行业
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) # 设置图标题
ax.set_xlabel('上市年数', fontsize=12) # x轴标签
ax.grid(True, alpha=0.3) # 添加网格线
plt.legend(prop={'family': 'Source Han Serif SC'}) # 设置图例中文字体
plt.tight_layout() # 自动调整布局
plt.show() # 显示图形
Cox模型结果：森林图
不同杠杆率分组的风险函数
Code
from lifelines import NelsonAalenFitter # 导入Nelson-Aalen累积风险估计器
company_survival_data['leverage_group'] = pd.cut( # 按杠杆率分组
company_survival_data['Initial_Leverage'], # 杠杆率变量
bins=[0, 0.3, 0.6, float('inf')], # 分组边界
labels=['低杠杆(0-30%)', '中杠杆(30-60%)', '高杠杆(>60%)'] # 分组标签
) # 完成分组
fig, ax = plt.subplots(figsize=(10, 6)) # 创建画布
nelson_aalen_fitter = NelsonAalenFitter() # 实例化NA估计器
for leverage_group in company_survival_data['leverage_group'].dropna().unique(): # 遍历各分组
mask = company_survival_data['leverage_group'] == leverage_group # 筛选当前分组
if mask.sum() < 2: # 样本不足则跳过
continue # 跳过
nelson_aalen_fitter.fit( # 拟合累积风险
company_survival_data.loc[mask, 'duration'], # 存续时间
company_survival_data.loc[mask, 'event'] # 事件标记
) # 完成拟合
cumulative_hazard = nelson_aalen_fitter.cumulative_hazard_ # 获取累积风险
hazard_estimate = cumulative_hazard.diff().fillna(0) # 差分得到风险函数估计
hazard_estimate.plot(ax=ax, label=leverage_group) # 绘制到画布
ax.set_title('不同初始杠杆率分组的退市风险函数', fontsize=14) # 设置标题
ax.set_xlabel('生存时间（年）', fontsize=12) # x轴标签
ax.set_ylabel('风险函数 h(t)', fontsize=12) # y轴标签
ax.grid(True, alpha=0.3) # 添加网格线
plt.legend(title='杠杆率分组', prop={'family': 'Source Han Serif SC'}) # 图例
plt.tight_layout() # 自动调整布局
plt.show() # 显示图形
Lasso正则化路径
Code
penalty_parameters = np.logspace(-4, 0, 20) # 生成20个惩罚参数
coefficient_path = [] # 记录系数变化路径
successful_penalties = [] # 记录成功的惩罚参数
for lambda_val in penalty_parameters: # 遍历每个惩罚参数
try: # 捕获可能的收敛错误
penalized_model = CoxPHFitter(penalizer=lambda_val, l1_ratio=1.0) # Lasso正则化
penalized_model.fit(scaled_predictor_df, duration_col='duration', event_col='event') # 拟合模型
coefficient_path.append(penalized_model.params_.values) # 记录系数
successful_penalties.append(lambda_val) # 记录成功参数
except Exception: # 收敛失败则跳过
continue # 跳过
penalty_parameters = np.array(successful_penalties) # 转为数组
coefficient_path = np.array(coefficient_path) # 转为二维数组
fig, ax = plt.subplots(figsize=(12, 7)) # 创建画布
for i, col in enumerate(predictor_matrix.columns): # 遍历每个特征
ax.plot(penalty_parameters, coefficient_path[:, i], label=col) # 绘制系数路径
ax.set_xscale('log') # x轴对数刻度
ax.set_xlabel('Lambda (正则化参数)', fontsize=12) # x轴标签
ax.set_ylabel('系数值', fontsize=12) # y轴标签
ax.set_title('Cox模型的Lasso正则化路径', fontsize=14) # 标题
ax.grid(True, alpha=0.3) # 网格线
plt.legend(prop={'family': 'Source Han Serif SC'}, fontsize=8) # 图例
plt.tight_layout() # 调整布局
plt.show() # 显示
生存树特征重要性
Code
fig, ax = plt.subplots(figsize=(10, 6)) # 创建画布
top_features = feature_importance_df.head(10) # 取前10重要特征
ax.barh(range(len(top_features)), top_features['importance']) # 水平条形图
ax.set_yticks(range(len(top_features))) # y轴刻度
ax.set_yticklabels(top_features['feature']) # y轴标签
ax.set_xlabel('特征重要性', fontsize=12) # x轴标签
ax.set_title('生存树：特征重要性排名', fontsize=14) # 标题
ax.grid(True, alpha=0.3, axis='x') # x方向网格线
plt.tight_layout() # 调整布局
plt.show() # 显示
回撤恢复时间的KM生存曲线
Code
km_recovery = KaplanMeierFitter() # 实例化KM估计器
km_recovery.fit(drawdown_data['duration'], drawdown_data['event'], label='回撤恢复') # 拟合生存曲线
fig, ax = plt.subplots(figsize=(10, 6)) # 创建画布
km_recovery.plot(ax=ax) # 绘制生存曲线
median_recovery = km_recovery.median_survival_time_ # 获取中位恢复时间
plt.axvline(median_recovery, color='r', linestyle='--', label=f'中位恢复时间: {median_recovery:.0f}天') # 参考线
ax.set_title('海康威视：股价回撤恢复时间生存曲线', fontsize=14) # 标题
ax.set_xlabel('天数', fontsize=12) # x轴标签
ax.set_ylabel('未恢复概率', fontsize=12) # y轴标签
ax.grid(True, alpha=0.3) # 网格线
plt.legend(prop={'family': 'Source Han Serif SC'}) # 图例
plt.tight_layout() # 调整布局
plt.show() # 显示
prob_30d = km_recovery.predict(30) # 30天未恢复概率
prob_180d = km_recovery.predict(180) # 180天未恢复概率
print(f'30天内无法恢复的概率: {prob_30d:.1%}') # 输出30天概率
print(f'半年(180天)内无法恢复的概率: {prob_180d:.1%}') # 输出180天概率
30天内无法恢复的概率: 44.7%
半年(180天)内无法恢复的概率: 13.2%