简介这份资源是2024年美国大学生数学建模竞赛ICM Problem E的完整参赛作品围绕财产保险可持续性这一现实议题用Python构建了从数据预处理到建模预测的全流程方案。内容适合数学建模初学者、进阶学习者以及需要完成课程设计、大作业或毕业设计的学生参考也可作为保险精算与风险评估方向的入门实践素材。压缩包共43个文件约18.43MB包含7个Python脚本、3个Excel数据表、16张PNG与7张SVG图表、2份PDF文档及LaTeX源码等覆盖灰色预测、层次分析、模糊综合评价、灰色关联与支持向量机等多类模型实现。资源中附有ROC曲线、PCA降维对比、混淆矩阵、灵敏度分析等可视化结果以及SVM模型文件与数据集便于读者复现建模思路、理解参数调优与结果验证过程。目前已有111人学习下载适合希望系统掌握数学建模完整流程与Python实现技巧的读者。1. 从一道美赛题说起财产保险的可持续性到底在算什么2024年美国大学生数学建模竞赛的这道题表面上是保险精算骨子里是一道典型的多目标决策与时间序列预测题。财产保险的可持续性说白了就是保险公司收上来的保费能不能覆盖未来因极端天气事件导致的赔付同时还要让公司活下去、让投保人愿意继续买。这个平衡点一旦被打破要么保费高到没人买要么赔付多到公司破产。适合谁看正在准备数学建模竞赛、想用 Python 把保险定价和气候风险量化跑通的人以及做金融风控、想了解巨灾模型落地路径的工程师。热搜里“python数据分析与可视化”“数学建模优秀论文”这些词恰好对应了这道题的两个核心动作用 Python 处理历史赔付数据用可视化把风险敞口讲清楚。接下来我不谈虚的直接拆这道题从数据到模型再到论文图表的完整链路。2. 拆解题目财产保险可持续性的三个量化维度2.1 保费充足率与赔付率的动态平衡保险可持续性的第一层是保费收入与赔付支出之间的时间错配。财产险尤其是巨灾险赔付不是均匀发生的一场飓风可能吃掉十年的利润。常见做法是构建一个赔付率时间序列用历史数据估计其分布再结合保费增长率判断未来若干年是否会出现累计赤字。这里的关键参数是目标偿付能力充足率通常监管要求不低于100%但实际运营中保险公司会留到150%甚至200%的缓冲。用 Python 实现时我一般会先算滚动12个月的赔付率再做蒙特卡洛模拟看第5年、第10年的破产概率。这一步不需要复杂模型pandas 的 rolling 加 numpy 的随机抽样就能跑出可信结果。2.2 极端天气事件的频率与强度建模财产险可持续性被打破往往不是因为日常小赔案而是极端天气。2024年美赛题里明确提到了天气相关的巨灾风险。对频率常用泊松分布或负二项分布拟合每年发生次数对强度用对数正态或广义帕累托分布拟合单次事件的损失金额。这里有个容易翻车的地方历史数据里极端值太少直接拟合会低估尾部风险。我一般会引入极值理论设定一个阈值对超过阈值的损失单独用广义帕累托分布建模。Python 里 scipy.stats 的 genpareto 和 poisson 可以直接调用但阈值选多少需要看平均超额图不能拍脑袋。2.3 再保险与资本约束下的可持续性判据保险公司不是独自扛下所有风险再保险是维持可持续性的关键工具。题目里通常会给再保险的层结构比如自留额、分出比例、赔付上限。可持续性判据可以定义为在给定再保险安排下未来T年内累计盈余小于零的概率低于某个阈值比如1%。这个判据把保费定价、天气风险和资本管理串在一起。用 Python 做的时候我会把再保险的赔付函数写成一个分段函数对每次模拟的年度总损失计算分出部分和自留部分再累加盈余。参数上自留额和分出比例是决策变量可以通过网格搜索找最优组合。3. 用 Python 跑通数据清洗与特征工程3.1 读取与合并多源赔付数据美赛题通常会提供多个 CSV 或 Excel 文件包括历史赔付记录、保单信息、天气事件表。第一步是把它们按年份和地区对齐。我一般用 pandas 的 read_csv 读入然后用 merge 按公共键合并。注意日期格式不统一是常态比如有的文件是“2020/1/1”有的是“01-01-2020”统一用 pd.to_datetime 加 format 参数处理不要依赖自动推断否则会出玄学错误。import pandas as pd import numpy as np # 读取三个数据源假设文件在当前目录 claims pd.read_csv(claims.csv, parse_dates[claim_date]) policies pd.read_csv(policies.csv, parse_dates[start_date, end_date]) weather pd.read_csv(weather_events.csv, parse_dates[event_date]) # 统一日期格式如果自动解析失败手动指定 claims[claim_date] pd.to_datetime(claims[claim_date], format%Y-%m-%d, errorscoerce) # 合并保单信息到赔付记录按保单号 df claims.merge(policies, onpolicy_id, howleft) # 按年月聚合赔付总额 df[year_month] df[claim_date].dt.to_period(M) monthly_claims df.groupby(year_month)[claim_amount].sum().reset_index() print(monthly_claims.head())这段代码的逻辑是先读入三个表把日期列强制转成统一格式再按保单号左连接最后按月聚合赔付金额。参数说明parse_dates 指定需要解析的列errorscoerce 让无法解析的日期变成 NaT避免报错中断。合并时用 howleft 保留所有赔付记录即使保单信息缺失也不丢数据。聚合后的 monthly_claims 是后续时间序列建模的基础。3.2 构造天气相关的特征变量天气事件表里通常有事件类型、风速、降雨量、影响区域。要把这些变成模型可用的特征需要做几件事按地区统计每年极端天气次数按事件强度分箱再和赔付数据按地区和年份合并。我一般会构造“年度巨灾次数”“最大单次风速”“累计降雨量”三个特征。注意地区编码要统一有的表用 FIPS 码有的用州名缩写得先映射。# 假设 weather 表有 event_type, wind_speed, rainfall, region, event_date weather[year] weather[event_date].dt.year # 只保留飓风、洪水、野火三类极端事件 extreme weather[weather[event_type].isin([hurricane, flood, wildfire])] # 按地区和年份聚合 weather_features extreme.groupby([region, year]).agg( event_count(event_type, count), max_wind(wind_speed, max), total_rain(rainfall, sum) ).reset_index() # 将地区映射到赔付数据中的地区列 region_map {AL: Alabama, CA: California, FL: Florida, TX: Texas} weather_features[region_full] weather_features[region].map(region_map) # 合并到赔付数据 df[year] df[claim_date].dt.year df df.merge(weather_features, left_on[region, year], right_on[region_full, year], howleft)逻辑说明先提取极端天气事件按地区和年份做聚合生成三个特征。region_map 是示例映射实际要根据数据里的地区编码调整。合并时用左连接保证赔付记录不丢。参数上agg 里的 count、max、sum 分别对应次数、最大风速和累计降雨。这一步做完每条赔付记录就带上了当年的天气背景可以进入建模阶段。3.3 缺失值与异常值的处理边界保险数据里缺失值很常见尤其是小额赔付的天气关联字段。我的原则是赔付金额缺失且无法从其他字段推算的直接删天气特征缺失的用同地区同年份的中位数填充并加一个缺失指示列。异常值方面赔付金额超过99.9分位数的不要直接删那可能是真实巨灾应该单独标记为“极端事件”并保留。用 pandas 的 quantile 和 clip 要小心clip 会改变分布我一般只做标记不做截断。# 标记极端赔付 threshold df[claim_amount].quantile(0.999) df[is_extreme] (df[claim_amount] threshold).astype(int) # 天气特征缺失用同地区同年份中位数填充 df[max_wind] df.groupby([region, year])[max_wind].transform( lambda x: x.fillna(x.median()) ) # 如果整组都是缺失用全局中位数兜底 df[max_wind] df[max_wind].fillna(df[max_wind].median()) # 删除赔付金额缺失的行 df df.dropna(subset[claim_amount])这段代码先算99.9分位数作为极端阈值生成标记列。然后用 groupby 加 transform 做分组填充transform 保证返回的索引和原表一致。最后兜底填充和删除缺失赔付行。参数说明quantile(0.999) 可根据数据量调整数据少时用0.99。is_extreme 列后续可以放进模型作为特征也可以用来分层建模。4. 建模与模拟从泊松-伽马到蒙特卡洛破产概率4.1 频率-强度模型的参数估计财产险巨灾模型的标准框架是频率-强度分离频率用泊松分布拟合每年事件次数强度用伽马或对数正态分布拟合单次损失。用 Python 的 scipy.stats 做最大似然估计。注意泊松的 lambda 估计就是样本均值但伽马分布的形状和尺度参数需要用 fit 方法。我一般会先画 QQ 图看拟合效果再决定是否换分布。from scipy import stats # 假设 annual_events 是每年极端事件次数序列 lambda_est annual_events.mean() # 单次损失序列只取极端事件对应的赔付 single_losses df[df[is_extreme] 1][claim_amount].values # 拟合伽马分布 shape, loc, scale stats.gamma.fit(single_losses, floc0) print(f泊松lambda: {lambda_est:.2f}, 伽马shape: {shape:.2f}, scale: {scale:.2f}) # 拟合对数正态作为对比 mu, sigma stats.lognorm.fit(single_losses, floc0)[1:3] print(f对数正态mu: {mu:.2f}, sigma: {sigma:.2f})逻辑说明lambda_est 是年事件频率gamma.fit 返回形状、位置、尺度三个参数floc0 强制位置为0因为损失不能为负。对数正态拟合返回 mu 和 sigma。参数说明shape 越大分布越集中scale 是尺度参数。实际选哪个分布看 AIC 或 KS 检验 p 值我一般两个都跑选拟合优度高的。4.2 蒙特卡洛模拟年度总损失有了频率和强度分布就可以模拟未来每年的总损失。步骤是先抽泊松决定当年事件次数再对每次事件抽伽马得到单次损失求和得到年度总损失。重复一万次得到年度损失的分布。这一步用 numpy 的随机数生成器设置种子保证可复现。np.random.seed(42) n_sim 10000 n_years 10 annual_losses np.zeros((n_sim, n_years)) for i in range(n_sim): for y in range(n_years): n_events np.random.poisson(lambda_est) if n_events 0: losses np.random.gamma(shape, scale, n_events) annual_losses[i, y] losses.sum() else: annual_losses[i, y] 0 # 计算每年损失的均值和95%分位数 mean_loss annual_losses.mean(axis0) var_95 np.percentile(annual_losses, 95, axis0) print(年度损失均值:, mean_loss) print(95%分位数:, var_95)逻辑说明双重循环外层模拟次数内层年份。每年先抽事件次数再抽损失金额并求和。参数说明n_sim 一万次足够稳定n_years 根据题目要求设通常5到10年。var_95 是95%分位数对应巨灾情景。注意这里假设年份之间独立实际可能有趋势但美赛题通常接受独立假设。4.3 再保险结构下的盈余模拟与破产概率加入再保险后保险公司的自留损失是年度总损失的一个分段函数。假设自留额为A分出比例为r赔付上限为L则自留损失 min(年度总损失, A) r * max(0, min(年度总损失, L) - A)。盈余 初始资本 保费收入 - 自留损失。累计盈余小于零即破产。用模拟结果算破产概率。initial_capital 1e7 # 初始资本 premium 2e6 # 年保费 retention 5e5 # 自留额 ceded_ratio 0.8 # 分出比例 limit 5e6 # 再保险赔付上限 def net_loss(gross_loss): layer1 np.minimum(gross_loss, retention) layer2 np.maximum(0, np.minimum(gross_loss, limit) - retention) return layer1 (1 - ceded_ratio) * layer2 surplus np.zeros((n_sim, n_years)) for i in range(n_sim): cum initial_capital for y in range(n_years): net net_loss(annual_losses[i, y]) cum cum premium - net surplus[i, y] cum ruin_prob (surplus 0).any(axis1).mean() print(f破产概率: {ruin_prob:.4f})逻辑说明net_loss 函数实现再保险分层layer1 是自留额内全赔layer2 是超出部分按分出比例赔。盈余逐年累加。参数说明initial_capital、premium、retention、ceded_ratio、limit 都是决策变量可以通过改变它们看破产概率变化。ruin_prob 是至少有一年盈余为负的模拟比例。5. 避坑与排查财产险建模里最容易翻车的五件事5.1 现象模拟破产概率为零但实际赔付率很高原因再保险参数设得太保守比如自留额极低、分出比例极高导致自留损失被压到很小。或者初始资本设得过大掩盖了风险。解决检查参数是否合理自留额通常与公司资本规模挂钩不能随意设。用敏感性分析画出破产概率随自留额变化的曲线找到拐点。5.2 现象伽马分布拟合报错提示形状参数无效原因单次损失数据里有零或负值。伽马分布定义在正实数上零和负数会导致拟合失败。解决先检查 single_losses 是否全为正如果有零说明极端事件标记有问题或者赔付金额字段有误。用single_losses single_losses[single_losses 0]过滤但要在论文里说明过滤了多少条。5.3 现象蒙特卡洛结果每次跑都不一样论文数据无法复现原因没有设随机种子。numpy 的随机数生成器默认从系统时间取种子。解决在模拟开始前加np.random.seed(42)或者用np.random.default_rng(42)创建独立生成器。论文里要写明种子值方便评委复现。5.4 现象合并数据后行数暴增原因合并键不唯一。比如保单表里一个保单号对应多条记录赔付表里也有多条merge 后产生笛卡尔积。解决合并前先检查键的唯一性用df.duplicated(subset[policy_id]).sum()看重复情况。如果确实需要多对多先聚合到同一粒度再合并。5.5 现象极端值标记后模型完全被极端事件主导原因is_extreme 标记的比例过高比如用了0.99分位数数据量又小导致10%的记录被标为极端。解决调整分位数阈值或者用绝对金额阈值而不是分位数。我一般会同时看分位数和业务含义比如赔付超过100万的才算巨灾而不是机械地用0.999。6. 把模型变成论文图表三个让评委一眼看懂的可视化技巧6.1 用累积分布图展示尾部风险财产险可持续性的核心是尾部风险但直方图看不出尾部。我一般画累积分布函数图横轴是年度总损失纵轴是累积概率在95%和99%分位处画竖线标注。这样评委一眼能看到“有5%的概率年度损失超过X”。用 matplotlib 的plt.plot(sorted_losses, np.linspace(0,1,len(sorted_losses)))即可比 seaborn 的 ecdfplot 更可控。import matplotlib.pyplot as plt sorted_losses np.sort(annual_losses[:, -1]) # 取最后一年 cdf np.arange(1, len(sorted_losses)1) / len(sorted_losses) plt.figure(figsize(8,5)) plt.plot(sorted_losses, cdf, labelCDF of Annual Loss) plt.axvline(np.percentile(sorted_losses, 95), colororange, linestyle--, label95% VaR) plt.axvline(np.percentile(sorted_losses, 99), colorred, linestyle--, label99% VaR) plt.xlabel(Annual Loss (USD)) plt.ylabel(Cumulative Probability) plt.legend() plt.title(Tail Risk of Property Insurance Losses) plt.show()逻辑说明先排序损失再算累积概率画线。两条竖线标出VaR。参数说明percentile 的95和99对应置信水平。这张图放在论文里比任何文字都直观。6.2 用热力图展示破产概率对再保险参数的敏感性再保险参数有两个关键变量自留额和分出比例。我一般会做一个网格每个格点跑一次模拟算破产概率然后用 seaborn 的 heatmap 画出来。这样能直接看到哪个区域破产概率低哪个区域高。注意模拟次数可以降到1000次以加快速度但论文里要说明。import seaborn as sns retentions np.linspace(1e5, 1e6, 10) ceded_ratios np.linspace(0.5, 0.95, 10) ruin_matrix np.zeros((len(retentions), len(ceded_ratios))) for i, r in enumerate(retentions): for j, c in enumerate(ceded_ratios): # 简化模拟只跑1000次 # 这里省略模拟细节假设有函数 calc_ruin(r, c) ruin_matrix[i, j] calc_ruin(r, c) sns.heatmap(ruin_matrix, xticklabelsnp.round(ceded_ratios,2), yticklabelsnp.round(retentions,0), cmapYlOrRd) plt.xlabel(Ceded Ratio) plt.ylabel(Retention) plt.title(Ruin Probability under Different Reinsurance Structures) plt.show()逻辑说明双重循环遍历参数组合每个组合算破产概率存进矩阵。heatmap 用颜色深浅表示概率高低。参数说明retentions 和 ceded_ratios 的范围根据实际资本调整。这张图能帮你在论文里论证最优再保险方案。6.3 用时间序列分解图讲清趋势与季节性如果数据按月或按季度可以画分解图把趋势、季节性和残差分开。statsmodels 的 seasonal_decompose 一行搞定。但保险数据季节性往往不明显趋势更重要。我一般会画滚动12个月赔付率曲线叠加原始月度数据用半透明线表示原始粗线表示滚动平均。这样既能看出波动又能看出趋势。from statsmodels.tsa.seasonal import seasonal_decompose # 假设 monthly_claims 是月度赔付序列索引为时间 monthly_claims.set_index(year_month, inplaceTrue) result seasonal_decompose(monthly_claims[claim_amount], modeladditive, period12) result.plot() plt.show()逻辑说明seasonal_decompose 把序列拆成趋势、季节、残差。参数说明modeladditive 适用于波动幅度不随趋势增大的情况否则用multiplicative。period12 表示年度周期。这张图放在论文里能展示你对数据结构的理解。最后说个我自己的习惯每次跑完模拟先把随机种子、参数组合、破产概率记在一个单独的 CSV 里不要只存在内存。论文写到一半发现某个参数要改回头找不到原始结果那种血泪经验一次就够了。这个方案值不值得做如果你在准备数学建模竞赛或者想入门保险精算的 Python 实现它是一条从数据到决策的完整链路跑通一次后面换数据换场景都能复用。希望帮到你。本文还有配套的精品资源点击获取