11 生存分析

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.family'] = ['Source Han Serif SC', 'Noto Serif CJK SC', 'SimSun']  # 设置中文字体优先级
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

DATA_DIR = 'C:/qiufei/data' if os.name == 'nt' else '/home/ubuntu/r2_data_mount/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'])  # 去除缺少上市日期的无效记录

observation_end_date = pd.to_datetime('2024-01-01')  # 设定研究观测截止日期
company_survival_records = []  # 初始化列表存储生存记录

for _, company_row in stock_basic_data.iterrows():  # 遍历每家上市公司
    listed_date = company_row['listed_date']  # 获取上市日期
    delisted_date = company_row['de_listed_date']  # 获取退市日期
    if pd.isna(delisted_date) or delisted_date >= observation_end_date:  # 未退市视为删失
        duration_years = (observation_end_date - listed_date).days / 365.25  # 计算存续年数
        event_indicator = 0  # 标记为删失（未退市）
    else:  # 已退市视为事件发生
        duration_years = (delisted_date - listed_date).days / 365.25  # 计算存续年数
        event_indicator = 1  # 标记为事件（已退市）
    if duration_years > 0:  # 排除异常负值记录
        company_survival_records.append({  # 添加生存记录
            'order_book_id': company_row['order_book_id'],  # 股票代码
            'industry': company_row.get('citics_2019_l1_name', 'Unknown'),  # 中信一级行业
            'duration': duration_years,  # 存续时间（年）
            'event': event_indicator  # 事件标记
        })  # 完成记录添加

company_survival_data = pd.DataFrame(company_survival_records)  # 转换为DataFrame
print(f'总样本量: {len(company_survival_data)}')  # 输出总样本量
print(f'退市公司数: {company_survival_data["event"].sum():.0f}')  # 输出退市数量
print(f'退市率: {company_survival_data["event"].mean():.2%}')  # 输出退市比例

financial_statement_path = os.path.join(DATA_DIR, 'stock/financial_statement.h5')  # 财务报表路径
financial_statement_raw = pd.read_hdf(financial_statement_path)  # 读取财务报表数据

financial_statement_raw['report_period_end'] = pd.PeriodIndex(  # 将季度标识转换为报告期末日期
    financial_statement_raw['quarter'], freq='Q'  # 按季度频率
).to_timestamp(how='end')  # 取季度最后一天

financial_groups = financial_statement_raw.groupby('order_book_id')  # 按股票代码分组

initial_leverage_records = []  # 存储每家公司的初始杠杆率

for _, company_row in company_survival_data.iterrows():  # 遍历每家公司
    stock_id = company_row['order_book_id']  # 股票代码
    if stock_id not in financial_groups.groups:  # 无财报数据则跳过
        continue  # 跳过当前循环
    company_financials = financial_groups.get_group(stock_id).sort_values('report_period_end')  # 按日期排序
    basic_info = stock_basic_data[stock_basic_data['order_book_id'] == stock_id]  # 查找基本信息
    if len(basic_info) == 0:  # 无基本信息则跳过
        continue  # 跳过当前循环
    listed_date = pd.to_datetime(basic_info.iloc[0]['listed_date'])  # 上市日期
    post_listing = company_financials[company_financials['report_period_end'] > listed_date]  # 上市后的财报
    if len(post_listing) > 0:  # 存在上市后财报
        first_report = post_listing.iloc[0]  # 取第一份
        total_assets = first_report.get('total_assets', None)  # 总资产
        total_liab = first_report.get('total_liabilities', None)  # 总负债
        if pd.notna(total_assets) and total_assets > 0 and pd.notna(total_liab):  # 数据有效
            leverage = total_liab / total_assets  # 计算资产负债率
            initial_leverage_records.append({'order_book_id': stock_id, 'Initial_Leverage': leverage})  # 记录

initial_leverage_df = pd.DataFrame(initial_leverage_records)  # 转换为DataFrame
company_survival_data = company_survival_data.merge(initial_leverage_df, on='order_book_id', how='left')  # 合并杠杆率
median_leverage = company_survival_data['Initial_Leverage'].median()  # 计算中位数
company_survival_data['Initial_Leverage'] = company_survival_data['Initial_Leverage'].fillna(median_leverage)  # 填充缺失值

province_mapping = stock_basic_data[['order_book_id', 'province']].drop_duplicates('order_book_id')  # 提取省份映射
company_survival_data = company_survival_data.merge(province_mapping, on='order_book_id', how='left')  # 合并省份

COASTAL_PROVINCES = ['广东省', '浙江省', '江苏省', '上海市', '北京市', '福建省', '山东省', '天津市']  # 沿海省份列表
company_survival_data['Region'] = company_survival_data['province'].apply(  # 生成地区分类
    lambda prov: 'Coastal' if prov in COASTAL_PROVINCES else 'Inland'  # 按省份判断沿海/内陆
)  # 完成分类

# 创建行业分组（取前5大行业）
company_survival_data['Industry_Group'] = company_survival_data['industry'].where(  # 条件赋值
    company_survival_data['industry'].isin(top_industry_categories), other='Other'  # 非前5行业标记为Other
)  # 完成行业分组

print(f'有效样本量: {len(company_survival_data)}')  # 输出样本量
print(f'沿海公司占比: {(company_survival_data["Region"]=="Coastal").mean():.1%}')  # 输出沿海占比

from lifelines import CoxPHFitter  # 导入Cox模型拟合器

cox_model_data = company_survival_data[['duration', 'event', 'Industry_Group', 'Region', 'Initial_Leverage']].copy()  # 选取建模变量
cox_model_data['Initial_Leverage'] = cox_model_data['Initial_Leverage'].fillna(  # 填充杠杆率缺失
    cox_model_data['Initial_Leverage'].median() if cox_model_data['Initial_Leverage'].median() == cox_model_data['Initial_Leverage'].median() else 0.5  # 若中位数为NaN则用0.5
)  # 完成填充
cox_model_data['Industry_Group'] = cox_model_data['Industry_Group'].fillna('Other')  # 填充行业缺失
cox_model_data['Region'] = cox_model_data['Region'].fillna('Inland')  # 填充地区缺失
cox_model_data = pd.get_dummies(cox_model_data, drop_first=True, dtype=float)  # 转换为哑变量
cox_model_data = cox_model_data.dropna(subset=['duration', 'event'])  # 删除关键列缺失行

cox_proportional_hazards_model = CoxPHFitter(penalizer=0.1)  # 添加L2正则化防止矩阵奇异
cox_proportional_hazards_model.fit(cox_model_data, duration_col='duration', event_col='event')  # 拟合Cox模型
cox_proportional_hazards_model.print_summary()  # 输出模型完整回归结果

from sklearn.preprocessing import StandardScaler  # 导入标准化工具

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])  # 仅保留数值列

scaler = StandardScaler()  # 实例化标准化器
scaled_predictor_matrix = scaler.fit_transform(predictor_matrix)  # 标准化特征
scaled_predictor_df = pd.DataFrame(scaled_predictor_matrix, columns=predictor_matrix.columns)  # 恢复列名
scaled_predictor_df['duration'] = cox_model_data['duration'].values  # 加回存续时间
scaled_predictor_df['event'] = cox_model_data['event'].values  # 加回事件标记

from sksurv.tree import SurvivalTree  # 导入生存树模型
from sksurv.util import Surv  # 导入生存数据结构化工具
from sklearn.inspection import permutation_importance  # 导入置换重要性

survival_targets = Surv.from_arrays(  # 构造生存目标
    event=cox_model_data['event'].astype(bool),  # 事件标记
    time=cox_model_data['duration']  # 存续时间
)  # 完成构造

tree_predictors = cox_model_data.drop(['duration', 'event'], axis=1)  # 特征矩阵

survival_tree_model = SurvivalTree(max_depth=3, min_samples_split=10, min_samples_leaf=5, random_state=42)  # 限制树深防止过拟合
survival_tree_model.fit(tree_predictors, survival_targets)  # 拟合生存树

perm_result = permutation_importance(  # 计算置换重要性
    survival_tree_model, tree_predictors, survival_targets, n_repeats=10, random_state=42  # 重复10次
)  # 完成计算
feature_importance_df = pd.DataFrame({  # 构建重要性表
    'feature': tree_predictors.columns, 'importance': perm_result.importances_mean  # 特征名和重要性
}).sort_values('importance', ascending=False)  # 按重要性降序排列

stock_price_history = pd.read_hdf(  # 读取后复权股价数据
    os.path.join(DATA_DIR, 'stock/stock_price_post_adjusted.h5')  # 后复权数据路径
).reset_index()  # 重置MultiIndex

haikang_data = stock_price_history[stock_price_history['order_book_id'] == '002415.XSHE'].copy()  # 筛选海康威视
haikang_data = haikang_data.sort_values('date').set_index('date')  # 按日期排序
closing_prices = haikang_data['close']  # 提取收盘价

running_max = closing_prices.cummax()  # 计算历史累计最高价
is_new_high = (closing_prices >= running_max)  # 标记创新高日期
new_high_dates = closing_prices.index[is_new_high]  # 提取新高日期

recovery_durations = []  # 存储回撤持续天数
recovery_events = []  # 存储事件标记

if len(new_high_dates) > 1:  # 至少需要两个新高
    for i in range(len(new_high_dates) - 1):  # 遍历相邻新高
        duration_days = (new_high_dates[i+1] - new_high_dates[i]).days  # 间隔天数
        if duration_days > 5:  # 仅保留显著回撤
            recovery_durations.append(duration_days)  # 记录天数
            recovery_events.append(1)  # 已恢复
    last_gap = (closing_prices.index[-1] - new_high_dates[-1]).days  # 最后一段
    if last_gap > 5:  # 当前仍在回撤中
        recovery_durations.append(last_gap)  # 记录天数
        recovery_events.append(0)  # 删失（未恢复）

drawdown_data = pd.DataFrame({'duration': recovery_durations, 'event': recovery_events})  # 构建数据框
print(f'识别到 {len(drawdown_data)} 次显著回撤周期')  # 输出总数