阿片成瘾的时间序列建模与短期预测

摘要：本文针对美赛 2019-C「阿片成瘾」的第一层问题——时间序列刻画与外推预测。以某州 7 年（84 个月，月度）的阿片过量死亡率 $y(t)$ （单位：每 10 万人，全样本均值 34.22、区间 27.80–40.29）为对象，建立"趋势 + 季节 + 干预阶跃"加性模型 $y(t)=a+bt+c\cdot\mathrm{ramp}(t)+\gamma_{t\bmod 12}+\varepsilon$ 。第 5 年初（ $t_0=54$ ）启动一项综合干预，其效果以 12 个月爬坡（ $\mathrm{ramp}$ ）的阶跃项刻画。分解显示：死亡率存在显著冬季偏高、夏季偏低的季节项（12 月指数 $+2.375$ 、7 月 $-2.253$ ，振幅约 4.6/10 万）；拟合得 $a=29.847,\;b=0.165,\;c=-8.790$ ，即干预前以 +0.165/月 上升、干预使水平下移约 8.79/10 万。以末尾 12 个月做留出校验，本文模型 MAPE=1.335%，远优于"忽略干预的趋势外推"（26.13%）与"季节朴素"（4.40%）；据此外推未来 12 个月（第 85–96 月）死亡率落在 33.5–38.5/10 万、95% 区间 [32.5, 39.5]。结论：该时间序列兼具趋势性、季节性与结构性突变（干预），朴素外推会严重高估未来负担，必须显式建模干预效应（图1—图8）。

关键词：时间序列；加性分解；干预哑变量；留出校验；MAPE

一、问题重述

赛题给出成瘾治疗与复发的纵向数据，要求建立可解释、可外推的模型描述阿片相关危害的演化，并支撑后续的风险分层与干预评估。第一层问题天然是一个单变量时间序列问题：把"月度阿片过量死亡率"当作观测序列，识别其趋势、季节与结构性变化，对未来做出点预测与区间预测，并量化预测不确定性。本文聚焦这一层——把序列形式化、可计算化，为第二层（个体风险分类）与第三层（干预因果评估）提供时间维度的基准锚点。

从公共卫生视角，阿片危机在美国呈明显的"三波"演化：处方阿片、 heroin、合成阿片（芬太尼）先后主导，宏观上表现为死亡率的阶段性抬升；同时，寒冷季节的社交孤立与用药风险使月度数据呈现冬季高峰。因此序列并非平稳，任何"拿历史均值当未来"或"纯线性外推"的做法都会失真——这正是本模型引入季节项与干预阶跃的根本动机。

从运行层面，州卫生系统真正需要的是"未来三到六个月大约会发生多少过量事件"这类粗粒度但可靠的预判，用以排程急诊容量、预置纳洛酮与戒毒床位。纯描述性结论（如"今年高于去年"）无法支撑这类排程，必须给出带不确定区间的量化外推——这正是本文模型的落脚点。

二、模型假设

观测为月度频率，共 $T=84$ 个月，时间索引 $t=0,\dots,83$ ；
序列由四部分叠加：长期趋势、周期为 12 的加性季节项、第 $t_0=54$ 月起的干预阶跃、零均值噪声；
干预效果在启动后 12 个月内线性爬坡至满额，之后保持，即 $\mathrm{ramp}(t)=\min((t-t_0)/12,1)$ （ $t<t_0$ 时为 0）；
季节项以月为周期、每年同形（无季节漂移）；
噪声独立同分布、方差恒定；所有参数由种子 2019 确定性给定，可一键复现。

三、符号说明

符号	含义	单位
$t$	月份索引	月
$y(t)$	月度阿片过量死亡率	每 10 万人
$a,b,c$	截距 / 趋势斜率 / 干预阶跃系数	— / (月) $^{-1}$ / (10万) $^{-1}$
$\gamma_m$	第 $m$ 月季节指数（ $m=0..11$ ）	每 10 万人
$\mathrm{ramp}(t)$	干预爬坡函数（0→1）	—
MAPE	平均绝对百分比误差	%

四、模型的建立

4.1 加性结构

将数据写为

$y(t)=a+bt+c\cdot\mathrm{ramp}(t)+\gamma_{t\bmod 12}+\varepsilon(t)$

其中趋势 $a+bt$ 捕捉长期演化，季节 $\gamma_{t\bmod 12}$ 捕捉年内周期，阶跃 $c\cdot\mathrm{ramp}(t)$ 捕捉干预带来的水平位移。该结构比"纯乘法"更利于分离干预效应，也便于后续做反事实（令 $c=0$ ）。相比乘法结构，加性假设"季节波动幅度不随水平放大"，对本文死亡率这类量级稳定、无指数增长的序列更贴切，也避免了对数为负的边界问题。

4.2 分解算法

先以窗宽 12 的中心移动平均提取趋势 $\hat\tau(t)$ （内部点 $t=6..77$ ）；去趋势后按月份取均值并归一化得季节指数 $\hat\gamma_m=\frac1{|S_m|}\sum_{t\in S_m}(y(t)-\hat\tau(t))-\bar\gamma$ ，使 12 个月指数和为 0。残差 $\hat\varepsilon(t)=y(t)-\hat\tau(t)-\hat\gamma_{t\bmod12}$ 用于估计噪声标准差 $\hat\sigma$ （图2）。中心移动平均窗宽取 12（一个完整周期）可无偏提取趋势；窗宽过窄会残留季节、过宽会过度平滑突变。序列两端（前/后 6 个月）缺失趋势估计，故首尾预测依赖线性外推，若近端恰逢结构突变会略偏，但留出校验已证明整体外推可靠。

4.3 拟合与预测

季节固定后，对残差 $r(t)=y(t)-\hat\gamma_{t\bmod12}$ 用最小二乘拟合 $r(t)=a+bt+c\cdot\mathrm{ramp}(t)$ ，解析解三元线性回归（正规方程）。预测

$\hat y(t)=\hat a+\hat b\,t+\hat c\cdot\mathrm{ramp}(t)+\hat\gamma_{t\bmod12}$

未来 $h$ 步的 95% 区间取 $\hat y(t)\pm1.96\hat\sigma$ 。

图1 阿片过量死亡率月度序列（红虚线=干预启动时点 t=54）

图2 加性分解：趋势 / 季节 / 残差

4.4 为何必须建模干预

若不显式建模干预（ $c=0$ ，仅趋势+季节外推），则模型假定后干预期仍沿 +0.165/月继续攀升，对末尾 12 月的 MAPE 高达 26.13%；而"季节朴素"（用 $y(t-12)$ 作预测）为 4.40%。本文模型因嵌入阶跃项，MAPE 降至 1.335%——差距本身即量化了"忽略干预"的代价（图7）。

五、模型求解与结果

趋势与干预： $\hat a=29.847,\;\hat b=0.165,\;\hat c=-8.790$ 。干预前序列以约 +0.165/月上行；干预使稳态水平下移约 8.79/10 万。
季节性（图3）：冬季偏高、夏季偏低。12 月指数 $+2.375$ 、1 月 $+2.145$ 、2 月 $+2.145$ 量级，7 月最低 $-2.253$ 、6 月 $-2.223$ ；峰谷振幅约 4.6/10 万，说明季节波动不可忽略。
留出校验（图5）：末尾 12 个月实际值（如 35.08、35.62、34.64、33.85）与预测（35.29、35.22、34.84、33.35）高度贴合，MAPE=1.335%，证明结构设定正确。
未来预测（图6）：第 85–96 月预测中值 37.26→38.51/10 万 间波动（季节驱动），95% 区间约 [32.5, 39.5]。

图3 月度季节指数（±，单位：每 10 万人）

图4 加性模型拟合 vs 实测（训练段，含干预哑变量）

六、结果分析与灵敏度

季节真实存在：峰（12 月/1 月）与谷（7 月）相差约 4.6/10 万，占均值 13%——任何忽略季节的模型都会系统性误判冬春高峰，对资源调度有害。
干预是结构性突变：阶跃系数 $|c|=8.79$ 远大于季节振幅，说明政策干预对水平的冲击比年内波动更剧烈；这正是"趋势外推"失效的根源（MAPE 26% vs 1.3%）。
拟合优度：残差标准差约 0.55/10 万，相对均值 34.2 仅 1.6%，说明加性结构已解释绝大部分变异。
预测区间合理：未来 12 月 95% 区间半宽约 ±3.5/10 万，主要来自噪声不确定；区间未覆盖"二次突变"，若再出新药（如更强合成阿片）则需重新估计 $c$ 。
模型比较结论（图7）：本模型 1.34% < 季节朴素 4.40% < 趋势外推 26.13%，层次清晰，佐证干预项不可或缺。
与后两问衔接：本层给出的 $t_0,c$ 是第三问反事实评估的支点；本层的"月度率"亦是第二问个体风险在人群层面的聚合表现。
冬季高峰的可操作性：季节振幅约 4.6/10 万，意味着每年 12—2 月会规律出现一波可预测的升高；公共卫生部门可据此在秋末前置纳洛酮库存与急诊值守，把"事后救火"转为"事前布防"。
作为早期预警监测器：本模型把最新月度数据代入即可更新 $\hat c$ 与残差；若某月残差持续偏离 $1.96\hat\sigma$ 带，提示出现了模型未涵盖的新冲击（如新型合成阿片流入），可作为预警触发人工复核。

图5 末尾 12 月留出校验：实测 vs 预测（MAPE=1.34%）

图6 未来 12 个月预测（蓝带=95% 区间）

七、模型评价

优点：(1) 加性结构可解释、各分量物理意义明确（趋势/季节/干预各司其职）；(2) 嵌入 $\mathrm{ramp}$ 阶跃显式分离政策效应，可直接做反事实；(3) 全程纯标准库实现、种子固定，可一键复现，满足竞赛可复现要求；(4) 留出校验 + 区间预测给出误差与不确定性双重度量。同时，加性形式天然支持"成分分解展示"，政策沟通时能把"本月升高是季节、趋势还是干预"一目了然地拆给决策者看，降低误读。

局限：(1) 季节项假设每年同形，未建模季节漂移；(2) 干预以固定爬坡刻画，未区分干预组成（戒毒床位、纳洛酮可及性等）；(3) 单变量，未引入失业率、处方量等外生协变量；(4) 未来区间仅含噪声不确定，不含"结构突变"尾部风险。

八、结论

该死亡率序列是一个"趋势上升 + 强季节 + 干预突变"的三元叠加过程。朴素外推因忽略干预会高估未来负担达 26% 之巨；本文加性模型以 1.335% 的留出误差准确刻画了结构，并给出未来 12 月 [32.5, 39.5]/10 万的区间。核心管理启示：政策干预对水平的冲击（≈8.8/10万）远大于季节波动（≈4.6/10万），评估与预测必须把干预"放进模型"，否则将系统性误判危机走向。此外，模型可滚动更新、充当早期预警监测器：当新观测持续逸出 95% 区间即提示结构性突变，应启动人工复核。

图7 留出 MAPE 模型对比（越低越好）

图8 本研究技术路线（Q1 时间序列）

九、管理建议（落地指引）

基于上述模型，给出五条可操作的监测与决策建议。第一，将加性模型嵌入公共卫生部门的月度监测看板，把"趋势、季节、干预"三项分解结果同步上墙，使决策者同时看见长期走向、周期波动和政策效果，避免被单一环比数字误导。第二，任何干预政策出台后，必须在模型里保留其效应项，并持续观察残差是否收敛到零附近；若残差长期偏离，说明干预强度已发生变化，需重新估计。第三，将 95% 预测区间作为早期预警触发器：当连续两个月实测值逸出区间上界，即自动告警，提示可能出现了新的合成阿片品种或供应冲击等结构性突变。第四，针对冬季高峰，在每年十月至次年二月前预先部署额外的纳洛酮储备与急诊人力，平抑季节性超额死亡。第五，模型应每月滚动更新，用新到数据重估参数，形成"监测—预警—复盘"的闭环。这五条建议的共同内核是：把预测从一次性的"算个数"升级为持续的"看过程"。

附录：核心 Python 实现

# 附录：核心 Python 实现（独立可运行，复现本文权威数字）
import os, sys
_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.abspath(os.path.join(_HERE, "..", "..", "..", "tools")))
import gen_mcm2019c as G

D = G.gen_mcm2019c()
f = D["fit"]
print("序列长度 T =", D["T"], " 干预启动 t0 =", D["t0"], " 爬坡月 =", D["ramp_len"])
print("y 区间 [%.2f, %.2f], 均值 %.2f" % (
    min(v for _, v in D["series"]), max(v for _, v in D["series"]),
    sum(v for _, v in D["series"]) / D["T"]))
print("季节指数(月1-12) =", [round(s, 3) for s in D["decomp"]["seasonal"]])
print("拟合 a=%.3f b=%.3f c=%.3f" % (f["a"], f["b"], f["c"]))
print("留出 MAPE = %.3f%%" % f["mape"])
print("未来12月预测 =", [round(v, 2) for v in f["fut_p"]])
print("未来12月 95% 区间下界 =", [round(v, 2) for v in f["fut_lo"]])
print("未来12月 95% 区间上界 =", [round(v, 2) for v in f["fut_hi"]])