1. 项目概述：从业务痛点出发的量化风险洞察

在金融、供应链、企业集团等复杂商业生态中，风险从来不是孤立存在的。一个核心客户的经营波动，可能通过担保链、应收账款、股权关联或市场情绪，像多米诺骨牌一样传导至整个集团乃至产业链上下游。过去，我们依赖经验判断和定性分析来评估这种“风险传递”，但面对海量、动态的关联数据，人脑的局限性日益凸显。定性分析往往只能给出“风险较高”或“需重点关注”这类模糊结论，难以量化具体的影响范围和冲击强度，更无法进行压力测试和前瞻性预警。

这正是“集团客户风险传递的数学建模”要解决的核心问题。这个项目的本质，是构建一个数字化的“风险传导模拟器”。它不再满足于“感觉有风险”，而是要精确回答：如果客户A出现偿付困难，其风险有多大比例会传递到关联方B、C？整个集团体系的脆弱点在哪里？在最坏的情况下，可能的损失边界是多少？为了实现这种从定性到定量的跨越，数学模型成为了不可或缺的工具。通过构建网络模型、设置传导规则、输入实际数据，我们可以将隐性的风险关联显性化、动态化。

我选择使用北太天元来完成这个模型的构建与求解，并非偶然。在金融工程和风险管理领域，我们常常面临一个困境：MATLAB功能强大但授权费用高昂，Python生态丰富但在某些数值计算和矩阵运算上需要额外的库和调试。北太天元作为一款新兴的科学计算软件，提供了与MATLAB高度兼容的语法和强大的内置数学库，同时兼顾了国产化和易用性。对于这类涉及矩阵运算（如邻接矩阵）、线性方程组求解、随机模拟的数学模型，北太天元的向量化操作和丰富的金融工具箱函数能极大提升开发效率。接下来，我将完整拆解这个模型的构建思路、数学原理，并附上可直接运行的北太天元代码，让你不仅能理解概念，更能亲手复现一个可用的风险传递分析工具。

2. 风险传递建模的核心思路与框架选择

构建风险传递模型，第一步是抽象现实。我们需要将复杂的集团客户关系，抽象为一种可计算的结构。最直观且有效的框架是网络图模型。在这个模型中，每个客户（或子公司）是图中的一个“节点”，客户之间的风险传导渠道（如担保、借贷、持股、交易依赖）则是连接节点的“边”。边的方向代表风险传导方向，边的权重代表传导的强度或概率。

2.1 模型框架选型：为什么是传染病模型与熵权法结合？

在众多网络动力学模型中，我选择了结合SIR传染病模型的思想和熵权法来确定边权重。这是经过实际项目验证的有效组合。

• SIR模型的适应性：经典的SIR模型将人群分为易感者(S)、感染者(I)、移除者(R)，描述了疾病通过接触传播的过程。这与风险传递高度相似：一个“健康”的客户（节点）可能因关联方“违约”而“感染”风险；风险在其关联网络中传播；部分客户在承受损失后可能“移除”（如破产清算），不再传导风险但也无法恢复。这个类比为我们提供了成熟的微分方程框架来描述风险状态的动态变化。
• 熵权法的客观性：风险传导的强度（边的权重）不能主观设定。我们需要依据客观数据来计算。熵权法是一种根据各项指标数据提供的信息量大小来确定权重的方法。在这里，我们可以选取多个衡量客户间关联紧密度的指标（如担保金额占比、交易额占比、持股比例等），为每对关联客户计算一个综合的“关联强度”权重。信息熵越小，说明该指标值的变异程度越大，提供的信息量越多，其权重也应越大。这保证了权重确定的客观性和数据驱动性。

2.2 模型核心输入与输出定义

在动手写代码前，必须明确模型的输入和输出是什么。

核心输入数据：

1. 客户清单：集团内所有需要监控的客户或子公司列表，编号为1, 2, ..., N。
2. 关联关系矩阵：一个N×N的矩阵，初始时，矩阵元素R(i,j)表示客户i对客户j的风险暴露度。例如，这可以是j对i的担保金额、i对j的应收账款、i持有j的股份比例等。R(i,i)通常为0（自身不构成外部风险）。
3. 客户个体风险指标：每个客户自身的财务状况指标，如资产负债率、流动比率、净利润增长率等，用于计算其初始的“感染概率”或“脆弱性”。
4. 风险传导参数：包括风险传导率（类比传染病感染率）、风险恢复率（类比疾病治愈率）等，这些参数可以通过历史数据校准或专家经验设定。

模型核心输出：

1. 风险传递路径与范围：模拟结果显示，从某个风险源出发，风险依次影响了哪些客户。
2. 风险冲击强度：每个客户最终承受的风险值或损失估计。
3. 系统脆弱性评估：识别出那些一旦出事会对整个网络产生最大影响的“关键节点”。
4. 压力测试结果：模拟极端情景（如某个大客户突然违约）下，整个集团体系的损失分布。

注意：初始关联矩阵的构建是关键，也是难点。数据可能分散在财务系统、合同管理系统、股权系统中。在实际项目中，往往需要数据团队的支持，进行大量的数据清洗、关联和归一化处理，才能形成一份干净、可用的关联矩阵。这是模型能否成功的“地基”。

3. 模型构建的详细步骤与北太天元实现

下面，我们进入实操环节，一步步用北太天元构建模型。我将假设一个包含5个客户的简化集团为例。

3.1 数据准备与关联矩阵构建

首先，我们定义客户和模拟他们之间的关联关系。这里，我们用“应收账款占比”作为关联强度的代理变量，即客户B的营收中有多大比例来自客户A，那么A出问题就会影响B。

% 假设有5个客户，编号1-5 num_clients = 5; % 初始化一个5x5的关联矩阵R，R(i,j)表示客户j对客户i的依赖程度（风险传导方向：j -> i） % 例如，R(2,1)=0.3表示客户1的营收占客户2总营收的30%。 R = zeros(num_clients, num_clients); % 手动设定关联关系（在实际中从数据库读取） R(2,1) = 0.3; % 客户1是客户2的大客户 R(3,1) = 0.15; R(4,2) = 0.4; % 客户2是客户4的大客户 R(5,3) = 0.25; R(5,4) = 0.2; % R(1,5) = 0.1; % 可以设定循环依赖，增加复杂性 % 可视化关联矩阵（非必须，但有助于理解） disp('客户关联矩阵 R (行i受列j影响的程度):'); disp(R); % 计算每个客户的“入度”风险暴露（即受他人影响的总和） risk_exposure_in = sum(R, 2); % 对每一行求和 % 计算每个客户的“出度”风险影响（即影响他人的总和） risk_impact_out = sum(R, 1)'; % 对每一列求和，转置成列向量 disp('每个客户的总风险暴露（受他人影响）:'); disp([(1:num_clients)', risk_exposure_in]); disp('每个客户的总风险影响（影响他人）:'); disp([(1:num_clients)', risk_impact_out]);

3.2 基于熵权法计算综合关联权重

假设我们有三个指标来衡量客户间关联：应收账款占比(ind1)、担保金额占比(ind2)、交易频率标准化值(ind3)。我们将为每对有关联的客户计算一个综合权重。

% 假设我们为有关联的边收集了三个指标数据 % 这里用随机数模拟，实际中替换为真实数据 % 假设有m条边 edge_list = find(R > 0); % 找出R中所有非零元素的位置（线性索引） [m, ~] = size(edge_list); m = length(edge_list); % 为这m条边生成3个指标的模拟数据矩阵 X (m行，3列) X = rand(m, 3); % 随机生成0-1之间的数 % 熵权法计算权重 % 1. 数据标准化（正向指标） X_normalized = zeros(size(X)); for col = 1:size(X, 2) X_normalized(:, col) = (X(:, col) - min(X(:, col))) / (max(X(:, col)) - min(X(:, col)) + eps); end % 2. 计算第j项指标下，第i条边的特征比重 p_ij p = X_normalized ./ sum(X_normalized, 1); % 3. 计算第j项指标的熵值 e_j k = 1 / log(m); e = -k * sum(p .* log(p + eps), 1); % 加eps防止log(0) % 4. 计算信息熵冗余度 d_j d = 1 - e; % 5. 计算各指标权重 w_j w = d / sum(d); disp('三个指标（应收、担保、交易频）的熵权法权重:'); disp(w); % 6. 计算每条边的综合关联强度 % 首先将权重应用到原始标准化数据上，然后加权平均（或加权和） edge_weight = X_normalized * w'; % 7. 将计算出的综合权重填回关联矩阵W（加权后的矩阵） W = zeros(size(R)); W(R > 0) = edge_weight; % 将计算出的权重赋给有关联的边 % 8. 对W进行归一化，使得每个客户受到的所有影响权重之和为1（行归一化） % 这是为了满足概率转移矩阵的性质，用于后续的模拟。 W_normalized = zeros(size(W)); for i = 1:num_clients row_sum = sum(W(i, :)); if row_sum > 0 W_normalized(i, :) = W(i, :) / row_sum; end end disp('加权并归一化后的风险传导权重矩阵 W_normalized:'); disp(W_normalized);

3.3 风险传递的SIR动力学模拟

现在我们有了传导权重矩阵W_normalized，可以开始模拟风险的动态传递过程了。我们将每个客户的状态定义为一个三维向量[S, I, R]，分别代表“健康”、“感染风险”、“风险出清”的概率或比例。这里简化为每个客户处于单一状态，用状态值（0健康，1感染）来模拟。

% 模拟参数设置 beta = 0.6; % 风险传导率，表示连接带来的感染概率强度 gamma = 0.1; % 风险恢复率，表示单位时间内从“感染”状态恢复的概率 T = 50; % 模拟总时间步长 dt = 1; % 时间步长 % 初始化状态：假设客户1初始爆发风险（“感染”） state = zeros(num_clients, 1); % 0表示健康，1表示感染风险 state(1) = 1; % 记录每个时间点感染客户的数量 infected_count = zeros(T+1, 1); infected_count(1) = sum(state); % 记录状态历史，用于绘图 state_history = zeros(num_clients, T+1); state_history(:, 1) = state; % 开始时间迭代模拟 for t = 1:T new_state = state; % 复制当前状态 for i = 1:num_clients if state(i) == 1 % 如果客户i已感染 % 它有 gamma 的概率恢复 if rand() < gamma new_state(i) = 0; % 恢复健康（风险出清） end else % 如果客户i健康 % 计算它被所有已感染邻居传染的风险 infection_force = 0; for j = 1:num_clients if state(j) == 1 % 传导强度 = 传导率beta * 权重W_normalized(i,j) infection_force = infection_force + beta * W_normalized(i, j); end end % 根据总感染力计算被感染的概率 prob_infect = 1 - exp(-infection_force * dt); % 常用的一种概率映射方式 if rand() < prob_infect new_state(i) = 1; end end end state = new_state; state_history(:, t+1) = state; infected_count(t+1) = sum(state); end % 可视化模拟结果 figure; subplot(2,1,1); plot(0:T, infected_count, 'b-o', 'LineWidth', 1.5); xlabel('时间步'); ylabel('感染风险客户数量'); title('风险传递模拟：感染客户数随时间变化'); grid on; subplot(2,1,2); imagesc(state_history); colorbar; xlabel('时间步'); ylabel('客户编号'); title('客户风险状态演化（黄色：感染，蓝色：健康）'); yticks(1:num_clients);

3.4 关键节点识别：基于网络中心性指标

模拟结束后，我们想知道哪些客户在风险网络中最重要。我们可以计算几个经典的中心性指标。

% 将权重矩阵W视为邻接矩阵，构建有向加权图 % 计算度中心性（加权） weighted_out_degree = sum(W, 1)'; % 加权出度 weighted_in_degree = sum(W, 2); % 加权入度 % 计算特征向量中心性（衡量一个节点与重要节点相连的程度） % 使用MATLAB兼容的eig函数 [V, D] = eig(W_normalized'); [~, idx] = max(abs(diag(D))); % 找到主特征值 eigenvector_centrality = abs(V(:, idx)); % 取主特征向量 eigenvector_centrality = eigenvector_centrality / sum(eigenvector_centrality); % 归一化 % 将结果汇总成表 node_metrics = table((1:num_clients)', weighted_in_degree, weighted_out_degree, eigenvector_centrality, ... 'VariableNames', {'客户编号', '加权入度中心性', '加权出度中心性', '特征向量中心性'}); disp('客户网络中心性指标:'); disp(node_metrics); % 找出关键节点（例如，特征向量中心性最高的前2个） [~, sorted_idx] = sort(eigenvector_centrality, 'descend'); key_nodes = sorted_idx(1:2); fprintf('\n基于特征向量中心性识别的关键节点（风险枢纽）是: 客户 %d 和 客户 %d\n', key_nodes(1), key_nodes(2));

4. 模型应用、调优与常见问题

构建出模型只是第一步，更重要的是应用和迭代。这个模型可以用于多种场景：

1. 定期风险扫描：每月或每季度运行一次模型，输入最新的关联数据和财务指标，生成集团风险热力图，识别出风险暴露上升最快的客户和关联路径。
2. 新客户/新业务准入评估：当集团拟新增一个重要客户或投资一家新公司时，可以将其虚拟加入现有网络，模拟其可能带来的风险传导效应，作为决策参考。
3. 压力测试与应急预案：设定极端情景（如某个核心客户违约率上升50%），运行模型，评估对集团整体资产质量的影响，并提前制定针对关键传导路径的应急预案。

4.1 模型参数校准与验证

模型中的关键参数，如风险传导率beta和恢复率gamma，不能永远靠猜。校准是让模型贴合现实的关键一步。

• 历史数据回溯：如果集团历史上发生过风险事件，且有详细的传导记录，可以用历史数据来反推最优的beta和gamma。例如，寻找一组参数，使得模型模拟的传导范围和速度与历史情况最吻合。这可以转化为一个优化问题，使用北太天元的fminsearch或fmincon函数求解。
• 专家经验结合：在没有足够历史数据时，可以邀请风控专家对几组预设参数下的模拟结果进行评价，选择他们认为最合理的一组。
• 敏感性分析：对beta和gamma在一定范围内进行扰动（例如±20%），观察模型输出（如最终感染客户数、峰值时间）的变化程度。如果输出变化剧烈，说明模型对该参数敏感，需要更谨慎地确定其值。

% 一个简单的参数敏感性分析示例 beta_range = 0.3:0.1:0.9; gamma_range = 0.05:0.05:0.2; peak_infected = zeros(length(beta_range), length(gamma_range)); for i = 1:length(beta_range) for j = 1:length(gamma_range) % 调用上面定义的模拟函数，传入不同的beta和gamma % 这里假设有一个封装好的模拟函数 simulate_risk_spread(beta, gamma) % [~, infected_ts] = simulate_risk_spread(beta_range(i), gamma_range(j), W_normalized); % peak_infected(i, j) = max(infected_ts); % 为示例，我们简单计算一个替代指标 peak_infected(i, j) = beta_range(i) / (gamma_range(j) + eps); % 简化的R0近似值 end end figure; imagesc(gamma_range, beta_range, peak_infected); colorbar; xlabel('风险恢复率 (gamma)'); ylabel('风险传导率 (beta)'); title('参数敏感性分析（颜色代表理论传播强度）'); set(gca, 'YDir', 'normal');

4.2 常见问题与排查技巧

在实际运行模型时，你可能会遇到以下问题：

1. 模拟结果不稳定，每次运行差异很大

   • 原因：这通常是因为模型中包含了概率性过程（rand()函数），而风险传导力又设置得比较临界。
   • 解决：进行蒙特卡洛模拟。不要只运行一次，而是运行成百上千次，然后对结果取统计量（如平均值、中位数、95%分位数）。这能给出更稳健的预期范围和风险分布。

   num_simulations = 1000; final_infected_counts = zeros(num_simulations, 1); for sim = 1:num_simulations % 运行一次模拟，获取最终感染数 % state_final = run_one_simulation(...); % final_infected_counts(sim) = sum(state_final); final_infected_counts(sim) = poissrnd(3); % 用泊松分布随机数示例 end fprintf('模拟%d次，平均最终感染客户数: %.2f\n', num_simulations, mean(final_infected_counts)); fprintf('风险价值(VaR) 95%%: %.2f\n', prctile(final_infected_counts, 95));

2. 关联矩阵非常稀疏，风险传不下去

   • 原因：真实网络中，很多客户可能只有少数几个强连接，导致风险传导路径很容易中断。
   • 解决：考虑间接关联。例如，引入二阶、三阶关联的影响。在数学上，这可以通过计算关联矩阵的幂（W^2,W^3）来评估间接暴露。或者，在模拟时，允许风险以较低的概率“跳跃”到非直接关联但同处一个脆弱子群的客户。

3. 模型输出难以向业务部门解释

   • 原因：直接展示矩阵和时序图，业务人员可能看不懂。
   • 解决：进行可视化翻译。将关键输出转化为业务语言和图表：
     • 风险传导图：使用北太天元的绘图功能或导出数据到专业工具，绘制动态风险传导网络图，用节点大小表示风险值，用箭头粗细表示传导强度。
     • 风险仪表盘：将“关键节点列表”、“最大潜在损失”、“风险传导Top 3路径”等核心结论，以卡片、排行榜、桑基图的形式展示。
     • 故事线叙述：不要只说“客户A是关键节点”，要说“如果客户A发生违约，根据模型，其风险将通过担保链在3个月内波及B、C、D三家子公司，预计最大损失约为X万元，其中对子公司B的影响占比最高，建议对A的授信加强监控，并为B准备流动性预案。”

4. 计算速度慢，客户数量多时效率低

   • 原因：使用多层for循环进行模拟，时间复杂度高。
   • 解决：向量化运算。北太天元擅长矩阵运算。尽量将循环操作改写为矩阵和向量操作。例如，计算感染力的循环可以用矩阵乘法替代。

   % 向量化计算感染力的示例（假设state是列向量） infection_force_matrix = beta * W_normalized; % 传导率乘以权重矩阵 % 当前感染节点的状态向量 infected_vec = (state == 1); % 每个健康节点受到的总感染力 = 矩阵对应行 * 感染节点向量 force_on_susceptible = infection_force_matrix * infected_vec; % 然后向量化计算感染概率 prob_infect = 1 - exp(-force_on_susceptible * dt); % 向量化判断是否感染 new_infections = (rand(num_clients,1) < prob_infect) & (state == 0);

   对于超大规模网络（成千上万个节点），可能需要考虑稀疏矩阵存储（sparse）和更高效的算法，甚至需要借助高性能计算或分布式计算框架，北太天元可能就需要与其他平台结合使用了。

这个模型不是一个一劳永逸的“黑箱”，而是一个需要与业务持续对话、用数据不断喂养和校准的“活工具”。从简单的5客户demo扩展到真实成百上千客户的系统，挑战主要在于数据的质量、计算的效率和结果的可解释性。但无论如何，迈出数学建模这一步，就已经将风险管控从模糊的经验主义，带向了清晰的数字驱动决策的新阶段。