电工杯 2020 B 范文二：BIPV 板块指数的收益分布诊断、滚动联动演化与蒙特卡洛概率化情景（厚尾压力测试视角）
一、摘要
本篇是本题系列范文的第二篇，与范文一「清洗—静态相关—双法点预测」的主线互补，聚焦三个更深一层的任务：收益率的分布形态诊断、联动的时间维度演化与预测的概率化输出。基于与范文一同一份确定性合成的 5 板块 750 日指数面板（固定种子 SEED=20200418，本文新增蒙特卡洛种子 SEED2=20200419），首先对对数收益率做矩诊断：五板块偏度介于 −0.160+0.030、超额峰度介于 −0.093+0.382，合成面板近高斯、无厚尾——这恰恰提示真实市场的厚尾风险无法从样本自身显现，应在情景引擎中以厚尾创新显式注入。其次以 60 日滚动相关刻画联动的动态结构：BIPV组件均值 0.861（区间 [0.669, 0.924]）、BIPV电站运营均值 0.576（区间 [0.335, 0.762]、极差 0.427），高联动的设备链相关稳定、弱联动的运营端漂移最大，静态相关系数会掩盖「联动阶段性减弱」的风险。再次构建几何布朗运动（GBM）路径引擎对比正态与标准化 t(4) 两种创新分布：30 日终值的 90% 中央区间两者几乎重合（257.3 对 256.2 点、差 −0.4%），但 1 日尺度的 99% 极端区间 t(4) 宽出 27.0%（93.8 对 73.8 点），即便 30 日聚合后 99% 区间仍宽 4.3%——「中央由中心极限定理主导、尾部由幂律律主导」的双面结构得到干净演示。最后把政策与原料冲击参数化为漂移与波动率扰动做三情景万条路径模拟：政策利好情景终值中位 1109.5（较期初 +2.50%）、上涨概率 63.7%，基准 51.4%，硅料涨价冲击情景中位 1055.3（−2.51%）、上涨概率仅 41.2% 且 P95 高达 1276.3——坏情景的本质不是更低的中位数，而是更宽的不确定性带。全文纯标准库、逐位可复现，正文、配图、附录与真源四路数字一致。
二、问题重述与本文切入点
赛题的五项产出中，范文一已覆盖数据清洗、静态相关矩阵、趋势波动率刻画与统计/机器学习两法的滚动一步点预测。尚存三类被浅尝辄止的问题，构成本篇的切入点：
- 分布形态：收益率是厚尾还是高斯？这决定风险度量（分位数、极端损失）的可信度，也决定模拟引擎该用什么创新分布。
- 联动的动态性：静态相关系数是全样本平均，掩盖了「某些阶段联动减弱、分散化失效」的结构变化，需要滚动窗口透视。
- 预测的概率化：单点预测不构成决策依据，「往前若干期的预测值与区间」应来自路径分布的分位数而非残差标准差的简单倍数；外部冲击情景也应给出「涨的概率有多大」而不是一条确定性曲线。
为此本文组织为「诊断（第四、五节）→ 引擎（第六节）→ 情景（第七节）」三段，全部计算与范文一共用同一份数据底座，保证跨篇数字互洽。
三、模型假设与符号
- H1（共享数据底座）：沿用范文一的确定性合成面板（5 板块 × 750 日，共同因子 + 特异噪声结构），本篇新引入的随机性仅来自蒙特卡洛种子 SEED2。
- H2（近高斯样本、厚尾世界）：合成面板的收益率近高斯是生成机制使然；真实板块指数以尖峰厚尾著称，故情景引擎默认采用标准化 t 分布创新做压力测试，并以正态创新作对照。
- H3（GBM 路径法则）：短中期内指数满足 的离散化形式，日增量 ， 为单位方差创新。
- H4（情景=参数扰动）：外部冲击通过漂移 平移与波动率 缩放进入引擎，不改变路径法则本身。
主要符号：日对数收益 ；偏度 、超额峰度 ；60 日滚动相关 ；终值分位 表示路径终值的第 分位数； 为期初点位（本面板 1082.4）。
四、收益率分布诊断（问题①③）
图 2 的收益率序列围绕零轴高频往返、偶有集群式放大，是典型的「带漂移低信噪比」金融序列。表 1 给出五板块的矩诊断。
表 1 对数收益率矩诊断（750 日）
| 板块 | 日均(bp) | 年化σ | 偏度 | 超额峰度 |
|---|---|---|---|---|
| BIPV | 1.1 | 21.06% | −0.160 | +0.382 |
| 硅料 | 2.9 | 17.18% | −0.074 | +0.278 |
| 组件 | 4.1 | 19.10% | −0.030 | −0.093 |
| 逆变器 | 4.5 | 22.48% | −0.112 | +0.058 |
| 电站运营 | 2.1 | 18.57% | +0.030 | −0.088 |
图 3 叠画实测直方与同均值方差的正态曲线，图 4 横向对比超额峰度：五板块峰度均落在 ±0.4 以内、偏度绝对值不超过 0.16，直方形态与正态参考线几乎贴合。结论有两面：其一，在本合成面板上直接套用高斯工具（正态分位区间、参数检验）不会失真；其二，真实市场收益普遍呈 的尖峰厚尾，样本近高斯只说明「厚尾风险藏在样本之外」——这正是第六节坚持用厚尾创新做压力测试的原因。
五、滚动相关与联动结构演化（问题②）
以 60 日窗口滑动计算 BIPV 与各板块的相关系数（图 5），统计见表 2。
表 2 滚动 60 日相关系数的分布
| 板块对 | 均值 | 最小 | 最大 | 极差 |
|---|---|---|---|---|
| BIPV~组件 | 0.861 | 0.669 | 0.924 | 0.255 |
| BIPV~逆变器 | 0.889 | 0.648 | 0.965 | 0.317 |
| BIPV~硅料 | 0.748 | 0.485 | 0.885 | 0.400 |
| BIPV~电站运营 | 0.576 | 0.335 | 0.762 | 0.427 |
两个层次的结构浮现：联动强度分层——设备链（组件、逆变器）全程维持高位（窗口均值 0.860.89），运营端最弱且均值不足 0.6，与范文一静态矩阵的排序完全一致；稳定性反序——联动越弱的板块对极差越大（电站运营 0.427 对组件 0.255），即弱联动并非「稳定的弱」，而是在 0.340.76 之间大幅摇摆。实务含义直接：以静态相关做的分散化或套保比例测算，在与运营端这类弱联动对手方交易时会周期性失效，须按滚动相关设置再平衡触发线。
六、蒙特卡洛预测区间：正态与厚尾创新的对照（问题④）
以训练段估计 （=日均 1.1bp、=日率 1.33%），从 出发各模拟 10000 条 GBM 路径，唯一差别是创新 取正态或标准化 （除以 保持单位方差）。表 3 是核心结果。
表 3 两类创新分布的预测区间宽度对比
| 期限 | 90% 区间宽（正态 / t(4)） | 99% 区间宽（正态 / t(4)） |
|---|---|---|
| 1 日 | 47.2 / 43.7（−7.4%） | 73.8 / 93.8（+27.0%） |
| 30 日 | 257.3 / 256.2（−0.4%） | 410.1 / 427.6（+4.3%） |
这组数字把「尖峰厚尾」教科书式的双面性量化得干干净净：t(4) 标准化后质量同时向中心与极端两端集中、肩部变薄，于是中央 90% 区间反而略窄（1 日 −7.4%）——若只用 90% 区间评估模型，甚至会误以为厚尾假设「更精确」；而 99% 极端区间暴露真相：1 日尺度宽出 27.0%。期限拉长到 30 日后，中心极限定理把增量求和推向正态、90% 区间差异消失（−0.4%），但幂律尾不会被求和平均——99% 区间仍宽 4.3%。结论：风险报告必须同时报告中央与极端两层分位，且尾部资本计提不能依赖中央区间外推。
正态创新的 30 日终值五分位为 P5=960.8、P25=1033.7、P50=1085.2、P75=1138.4、P95=1218.1（图 6），上涨概率 51.4%；t(4) 创新的对应值为 962.2/1032.1/1082.1/1134.5/1218.5、上涨概率 49.8%——中央层面两种引擎几乎可互换，进一步印证差异只在尾部。
七、三情景概率化分析（问题⑤）
以 t(4) 引擎参数化三种外部情景：政策利好（）、基准（历史 ）、硅料涨价冲击（ 且 ），各 10000 条 30 日路径，结果见表 4 与图 7、图 8。
表 4 三情景 30 日终值分布（t(4) 创新，S₀=1082.4）
| 情景 | P5 | P50 | P95 | 中位涨跌幅 | P(终值>期初) |
|---|---|---|---|---|---|
| 政策利好 | 984.9 | 1109.5 | 1250.3 | +2.50% | 63.7% |
| 基准 | 965.8 | 1085.3 | 1219.6 | +0.26% | 51.4% |
| 硅料涨价冲击 | 871.5 | 1055.3 | 1276.3 | −2.51% | 41.2% |
三点解读：① 漂移扰动的量级决定中位位置——每 10bp/日的漂移差经 30 日复利约为 3 个百分点的中位分离（+2.50% 对 −2.51%），线性且可预期；② 波动率放大的效果不对称地体现在两端——冲击情景的 P5 下探至 871.5（较期初 −19.5%）的同时 P95 反而是三情景最高（1276.3），坏情景的正确画像是「更低的中位 + 更宽的区间」，而不是整体下移的一条线；③ 决策口径应落在概率上：即便在冲击情景下仍有 41.2% 的路径收在期初之上，反之政策利好也有 36.3% 的路径收跌——任何「利好必涨」的单点叙事都与路径分布不符。
八、模型检验、优缺点与拓展
检验：① 数据底座与范文一同源同种，跨篇数字互洽（如静态相关排序与滚动均值排序一致）；② 蒙特卡洛种子独立且固定，两次运行输出逐字节一致；③ t(4) 标准化保持单位方差，保证与创新分布无关的可比口径。
优点：① 把「分布诊断—动态联动—概率化情景」串成闭环，恰好补齐赛题产出 2/3/5 的深层要求；② 「90% 对 99% 双指标」的设计直观揭示厚尾双面性，方法论可迁移到任意资产的风险汇报；③ 三情景以概率语言表达，天然衔接决策分析。
局限与拓展：① GBM 忽略波动聚集（GARCH 效应）与跳变强度时变性，可升级为带随机波动的路径法则；② 情景参数（±8bp/日、σ×1.6）为规范设定，正式参赛应以历史事件（政策发布日、硅料急涨段）实证校准；③ 可将滚动相关扩展为 DCC 型动态条件相关，或对终值分布做核密度平滑以给出连续的完整分布而非有限分位。
九、结论
本文以 BIPV 及产业链五板块 750 日确定性面板为底座，完成了分布诊断、动态联动与概率化情景三层深化：① 收益矩诊断显示合成面板近高斯（超额峰度 ≤0.382），据此把厚尾风险以标准化 t(4) 创新显式注入情景引擎而非依赖样本外推；② 滚动 60 日相关揭示「设备链联动强而稳（均值 0.86~0.89）、运营端弱而飘（均值 0.576、极差 0.427）」的双层结构，警示静态相关的分散化陷阱；③ 万条路径对照表明中央区间与创新分布几乎无关（30 日 90% 宽差 −0.4%）而极端区间显著敏感（1 日 99% 宽 +27.0%），尾部度量必须用厚尾引擎；④ 三情景给出可决策的概率口径：政策利好中位 +2.50%、上涨概率 63.7%，硅料冲击中位 −2.51%、上涨概率仅 41.2% 且区间最宽。连同范文一的点预测主线，本题五种产出已形成「清洗—联动—趋势—预测—情景」的完整证据链。
参考文献
[1] Barone-Adesi G, Giannopoulos K, Vosper L. VaR without correlations for portfolios of derivative securities[J]. Journal of Futures Markets, 1999, 19(5): 583–602.
[2] Cont R. Empirical properties of asset returns: stylized facts and statistical issues[J]. Quantitative Finance, 2001, 1(2): 223–236.
[3] Engle R. Dynamic conditional correlation: A simple class of multivariate generalized autoregressive conditional heteroskedasticity models[J]. Journal of Business & Economic Statistics, 2002, 20(3): 339–350.
[4] Glasserman P. Monte Carlo Methods in Financial Engineering[M]. New York: Springer, 2003.
[5] 中国电机工程学会. 全国大学生电工数学建模竞赛历年赛题汇编[M]. 北京: 中国电力出版社, 2020.
附录：核心 Python 实现
以下代码为真源 tools/gen_dgcup2020b_2.py 的等价核心摘录（复用范文一面板、固定种子、纯标准库），可在本文件所在目录直接运行，输出与正文全部数字一致。
# -*- coding: utf-8 -*-
"""收益分布诊断 + 滚动相关 + 蒙特卡洛概率化情景核心（纯标准库）。"""
import math
import random
import os
import sys
_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.abspath(os.path.join(
_HERE, "..", "..", "..", "tools")))
import gen_dgcup2020b as G
SEED2 = 20200419
NSIM = 10000
H = 30
ROLLW = 60
NU = 4
SECTORS = G.SECTORS
def moments(xs):
n = len(xs)
m = sum(xs) / n
m2 = sum((v - m) ** 2 for v in xs) / n
m3 = sum((v - m) ** 3 for v in xs) / n
m4 = sum((v - m) ** 4 for v in xs) / n
sd = math.sqrt(m2)
skew = m3 / sd ** 3 if sd > 0 else 0.0
ekur = m4 / m2 ** 2 - 3.0 if m2 > 0 else 0.0
return m, sd, skew, ekur
def log_rets(levels):
return [math.log(levels[t] / levels[t - 1])
for t in range(1, len(levels))]
def rolling_corr(a, b, w=ROLLW):
return [G.corr(a[i - w + 1:i + 1], b[i - w + 1:i + 1])
for i in range(w - 1, len(a))]
def _bm(rng):
u1 = rng.random(); u2 = rng.random()
if u1 < 1e-12:
u1 = 1e-12
return math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.pi * u2)
def t4_std(rng):
"""标准化 t(4) 创新（单位方差）。"""
z = _bm(rng)
c2 = sum(_bm(rng) ** 2 for _ in range(NU))
return z / math.sqrt(c2 / NU) / math.sqrt(NU / (NU - 2.0))
def mc_terminal(mu, sig, s0, innov, rng, h=H, n=NSIM):
drift = mu - sig * sig / 2.0
terms = []
day_vals = [[] for _ in range(h)]
for _ in range(n):
lv = s0
for k in range(h):
lv *= math.exp(drift + sig * innov(rng))
day_vals[k].append(lv)
terms.append(lv)
med_path = []
for k in range(h):
col = sorted(day_vals[k])
med_path.append(col[n // 2])
return terms, med_path
def quantile(sorted_xs, q):
i = min(int(q * (len(sorted_xs) - 1)), len(sorted_xs) - 1)
return sorted_xs[i]
def run():
D = G.gen_dgcup2020b()
cleaned = D["cleaned"]
Nn = len(cleaned)
print("== 收益率分布诊断（750 日对数收益）==")
stats = {}
for s in SECTORS:
rs = log_rets(D["levels"][s]) if s != "BIPV" \
else log_rets(cleaned)
m, sd, sk, ek = moments(rs)
stats[s] = dict(m=m, sd=sd, skew=sk, ekur=ek)
print("%s: 日均 %.1fbp 年化σ %.2f%% 偏度 %+.3f 超额峰度 %+.3f"
% (s, m * 1e4, sd * math.sqrt(252) * 100, sk, ek))
bi = log_rets(cleaned)
print("== 滚动 60 日相关（BIPV~）==")
roll = {}
for s in ("组件", "逆变器", "硅料", "电站运营"):
rc = rolling_corr(bi, log_rets(D["levels"][s]))
roll[s] = rc
print("~%s: 均值 %.3f 最小 %.3f 最大 %.3f 极差 %.3f"
% (s, sum(rc) / len(rc), min(rc), max(rc),
max(rc) - min(rc)))
print("== 蒙特卡洛预测区间（BIPV，S0=%.1f）==" % cleaned[-1])
mu = sum(bi) / len(bi)
sig = (sum((v - mu) ** 2 for v in bi) / len(bi)) ** 0.5
s0 = cleaned[-1]
for hh in (1, H):
rng_n = random.Random(SEED2)
tg, _ = mc_terminal(mu, sig, s0, _bm, rng_n, h=hh)
tg.sort()
w90g = quantile(tg, 0.95) - quantile(tg, 0.05)
w99g = quantile(tg, 0.995) - quantile(tg, 0.005)
rng_t = random.Random(SEED2 + 1)
tt, _ = mc_terminal(mu, sig, s0, t4_std, rng_t, h=hh)
tt.sort()
w90t = quantile(tt, 0.95) - quantile(tt, 0.05)
w99t = quantile(tt, 0.995) - quantile(tt, 0.005)
print("%d日区间: 90%%宽 正态 %.1f | t(4) %.1f (%+.1f%%) ||"
" 99%%宽 正态 %.1f | t(4) %.1f (%+.1f%%)"
% (hh, w90g, w90t, (w90t / w90g - 1) * 100,
w99g, w99t, (w99t / w99g - 1) * 100))
rng_n = random.Random(SEED2)
tg, pg = mc_terminal(mu, sig, s0, _bm, rng_n)
tg.sort()
qg = [quantile(tg, q) for q in (0.05, 0.25, 0.5, 0.75, 0.95)]
p_up_g = sum(1 for v in tg if v > s0) / len(tg) * 100.0
print("正态创新 30 日五分位: P5 %.1f P25 %.1f P50 %.1f P75 %.1f"
" P95 %.1f" % tuple(qg))
print(" P(终值>期初) %.1f%%" % p_up_g)
rng_t = random.Random(SEED2 + 1)
tt, pt = mc_terminal(mu, sig, s0, t4_std, rng_t)
tt.sort()
qt = [quantile(tt, q) for q in (0.05, 0.25, 0.5, 0.75, 0.95)]
p_up_t = sum(1 for v in tt if v > s0) / len(tt) * 100.0
print("t(4)创新 30 日五分位: P5 %.1f P25 %.1f P50 %.1f P75 %.1f"
" P95 %.1f" % tuple(qt))
print(" P(>期初) %.1f%%" % p_up_t)
print("== 三情景蒙特卡洛（t(4) 创新）==")
specs = [("政策利好", mu + 0.0008, sig, SEED2 + 101),
("基准", mu, sig, SEED2 + 102),
("硅料涨价冲击", mu - 0.0007, sig * 1.6, SEED2 + 103)]
rows_sc = []
for name, m_i, s_i, sd_i in specs:
rr = random.Random(sd_i)
tm, pm = mc_terminal(m_i, s_i, s0, t4_std, rr)
tm.sort()
q5 = quantile(tm, 0.05); q50 = quantile(tm, 0.5)
q95 = quantile(tm, 0.95)
pup = sum(1 for v in tm if v > s0) / len(tm) * 100.0
rows_sc.append(dict(name=name, p5=q5, p50=q50, p95=q95,
pup=pup, path=pm))
print("%s: P5 %.1f P50 %.1f P95 %.1f 涨幅中位 %+.2f%% "
"P(>期初) %.1f%%"
% (name, q5, q50, q95,
(q50 / s0 - 1) * 100, pup))
return dict(stats=stats, roll=roll, bi=bi, mu=mu, sig=sig,
s0=s0, qg=qg, qt=qt, p_up_g=p_up_g,
p_up_t=p_up_t, scen=rows_sc, pg=pg, pt=pt,
cleaned=cleaned)
if __name__ == "__main__":
run()