生猪养殖生产决策优化（一）：价格周期识别与预测建模

一、问题重述

2014 年国赛 C 题以生猪养殖为背景，要求养殖场依据历史猪价、成本与存栏数据，制定科学的补栏与出栏计划，使一个生产周期内的净回报最大。其本质是一个**"先预测价格周期、再据此做生产决策"**的两阶段问题：若不能准确捕捉猪价的周期性波动，补栏时点就会落在价格低谷，造成亏损；反之，若能在价格高峰前完成育肥并出栏，回报将显著提升。

本题的难点并不在于"养多少头猪"本身，而在于时间序列中隐藏的周期结构。猪价受供需滞后（母猪补栏→育肥→出栏约需 6–10 个月）影响，天然呈现"价高→扩产→过剩→价跌→减产→短缺→价高"的蛛网式循环，周期通常长达 3–4 年。因此，本文第一篇聚焦第一阶段——如何用谐波回归从 60 个月的历史数据中提取猪周期，并对未来 36 个月做出高精度预测。

二、数据理解与预处理

本站提供的练习数据集 data/cumcm2014c.csv 包含 60 个月历史记录（字段：月份、类型、猪价 price 元/kg、仔猪价 piglet_price 元/头、饲料月成本 feed_month 元/头/月），以及后续 36 个月的预测段。

观察历史猪价（图 1）可见明显的振荡形态：价格在 13–19 元/kg 之间起伏，且并非白噪声式的随机波动，而是带有可辨识的波峰与波谷。这说明单纯用移动平均或指数平滑这类"局部平滑"方法会平滑掉周期信息，必须采用能显式刻画周期项的模型。

建模前的常规数据体检亦不可或缺：先检查缺失值与异常值（如录入错误导致的 0 或负值），必要时用邻近月均值插补；再核对量纲（猪价为元/kg、仔猪价为元/头、饲料为元/头/月，三者不可直接相减，必须折算到"每头生猪"口径）；最后观察序列是否含明显趋势与周期。本数据经核对无缺失、量纲一致，且肉眼即可见约 3 年一轮的起伏，与后续谐波回归识别出的 $K=36$ 个月相互印证，说明数据结构健康、建模前提成立。

三、谐波回归模型

本文采用**谐波回归（Harmonic Regression）**对猪价建模。设月份序号为 $t$ ，模型形如：

$p_t = a_0 + a_1 t + b\cos\left(\frac{2\pi t}{K}\right) + c\sin\left(\frac{2\pi t}{K}\right) + \varepsilon_t$

其中 $a_0$ 为基线水平， $a_1 t$ 刻画长期趋势（饲料、通胀等引起的缓慢漂移）， $b\cos(\cdot)+c\sin(\cdot)$ 合成一个周期为 $K$ 的简谐振荡，捕捉猪价的周期性涨跌， $\varepsilon_t$ 为残差。该模型对四个未知参数 $(a_0,a_1,b,c)$ 是线性的，可用最小二乘正规方程直接求解，无需迭代优化，数值稳定且可复现。

3.1 周期 $K$ 的识别

周期 $K$ 是模型的关键超参。本文在 $K\in[24,48]$ 月（即 2–4 年，符合猪周期常识）范围内逐值拟合，以拟合优度 $R^2$ 为准则选取最优周期（图 3）。计算表明，当 $K=36$ 个月时 $R^2$ 达到峰值 0.9944，远高于其他候选周期，说明 3 年一轮的猪周期假设与数据高度吻合。

3.2 拟合与预测

以 $K=36$ 拟合全量历史，得到残差标准差仅 $\sigma=0.257$ 元/kg（图 4），表明模型几乎解释了价格的全部可解释波动。将拟合得到的 $(a_0,a_1,b,c)$ 外推至未来 36 个月，即可得到预测价序列，并以 $\pm1.96\sigma$ 给出 95% 预测区间（图 2）。预测显示未来一年猪价将先探底至约 11.9 元/kg，再回升至 22 元/kg 以上，呈现完整的下行—上行周期。

3.3 参数估计的数值示例

谐波回归的实质是多元线性回归：把原始问题转化为对设计矩阵 $\mathbf{X}=[\mathbf{1},\,\mathbf{t},\,\cos(2\pi\mathbf{t}/K),\,\sin(2\pi\mathbf{t}/K)]$ 的最小二乘求解，即由正规方程 $(\mathbf{X}^\top\mathbf{X})\boldsymbol{\beta}=\mathbf{X}^\top\mathbf{p}$ 得到 $\boldsymbol{\beta}=(a_0,a_1,b,c)$ 。以 $K=36$ 为例，本数据拟合得到趋势斜率 $a_1\approx+0.03$ 元/kg·月（猪价中枢缓慢上行），周期振幅 $\sqrt{b^2+c^2}\approx5$ 元/kg（周期波动幅度），相位由 $b,c$ 的相对大小决定。由于设计矩阵列数仅 4，正规方程规模极小（ $4\times4$ ），用高斯消元即可在毫秒内解出，且解唯一稳定——这比任何迭代式神经网络都更适合"小样本、强结构"的赛题数据。也正因为参数只有 4 个，模型不易过拟合，外推 36 个月的预测依然可靠。

四、预测方法的对比验证

为证明谐波回归的适用性，本文采用留出验证（Hold-out）：用前 48 个月训练，对后 12 个月做预测，以平均绝对百分比误差 MAPE 评价。对比了四种方法：

方法	谐波回归	Holt 线性趋势	移动平均(3)	指数平滑(0.3)
MAPE	1.21%	44.14%	27.75%	25.23%

结果（图 5）显示，谐波回归的误差仅为 1.21%，其余三种方法均超过 25%，Holt 甚至高达 44%。原因在于猪价的核心结构是周期而非单调趋势——Holt 只捕捉线性漂移，完全丢失了振荡信息；移动平均与指数平滑则把波峰波谷平滑成了平缓曲线。这充分说明：对周期性时间序列，必须先识别周期，否则任何"平滑类"预测都会严重失真。

为直观展示谐波预测的精度，图 7 给出留验段前 6 个月的预测价与真实价对照：预测值依次为约 12.48、12.13、11.92、11.85、11.94、12.18 元/kg，而该段真实价恰处于周期下行探底区，二者误差均在 1% 以内。这说明模型不仅抓住了周期相位，也准确刻画了下行斜率。需要强调的是，平滑类方法在拐点处（价格由涨转跌或反之）误差会骤增，而谐波回归因显式建模周期，在拐点附近依然稳健——这正是其 MAPE 远低于对照方法的结构性原因。

五、价格与成本的周期联动

猪价并非孤立变量。仔猪价是未来出栏价的先行指标（当前仔猪贵，往往预示前期母猪补栏多、未来供给足、价格将跌）；饲料成本则受玉米、豆粕等大宗农产品周期影响。图 6 展示三者的联动：仔猪价与猪价高度同向（相关系数高），饲料成本亦呈弱周期。因此，在后续决策模型中，仔猪采购成本与饲料支出都应随预测价同步变动，而不能假设为常数。

六、模型解释力与局限

图 8 将谐波模型的"趋势+周期"拟合曲线叠加到原始序列上，可见二者几乎重合，说明模型成功分离出长期趋势与猪周期两个成分。其局限在于：① 假设周期严格为 36 个月，真实周期可能受疫情、政策等外生冲击而变长或变短；② 外推 36 个月时，远离训练区的预测受趋势项 $a_1 t$ 累积影响，区间逐渐变宽（图 2 右端）。这些不确定性将在第三篇通过情景分析加以讨论。

七、预测精度示例与 ARIMA 对比

在方法选型上，周期性时间序列还有另一常用工具——ARIMA。ARIMA(p,d,q) 通过差分消除趋势、用自回归项捕捉短期依赖，对平稳短相关序列效果良好；但猪价是趋势+强周期的混合体，若直接对原序列建模，ARIMA 的自回归阶数需取得很高才能近似周期项，参数增多易过拟合。谐波回归则把"周期"作为显式结构项，仅用 4 个参数即可表达任意相位、任意幅度的简谐振荡，参数少、可解释、外推稳定。本数据上，谐波回归留验 MAPE 1.21%，而同等信息量下 ARIMA 通常难以低于 5%–8%（其本质仍是对过去值的加权平均，无法像谐波那样"预见"尚未到来的波峰）。因此，本题"先识别周期、再用谐波外推"的路线，在精度与可解释性上均优于黑箱式方法。

进一步看，谐波回归给出的 $a_1$ （趋势项，约 +0.03 元/kg·月）具有明确业务含义：它对应养殖成本长期上行（饲料、人工通胀）推动的猪价中枢缓慢抬升。决策者可据此判断"长期看猪价中枢上移"，但具体买卖时点仍须由周期项 $b\cos+c\sin$ 决定。这种"趋势定性、周期定量"的拆分，比单一预测值更有利于后续的生产计划制定。

八、小结

本篇建立了基于谐波回归的猪价周期预测模型，核心结论为：

猪周期真实存在且周期约为 36 个月，谐波回归对该数据的拟合优度 $R^2=0.9944$ ；
周期建模远优于平滑类方法，留出验证 MAPE 仅 1.21%（对比 Holt 44.14%）；
未来 36 个月预测显示价格将经历一次完整的探底—回升，为下一篇的生产决策提供了可靠输入；
谐波回归仅用 4 个可解释参数即分离出"趋势+周期"，相比 ARIMA 等黑箱方法更透明、外推更稳，特别适合小样本强结构的赛题数据。

综上，准确预测是科学养殖的第一步：只有先"看清周期"，后续才能"踏准节奏"。本篇的方法可平移到其他呈明显周期的价格序列（如禽蛋、饲料原料、部分大宗农产品），具有普遍的建模参考价值。

下一篇将把这一预测结果作为已知条件，建立补栏/出栏的整数规划模型，求解使净回报最大的生产计划。

附录：预测模型复现代码

import sys, os
sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "..", "..", "tools"))
import gen_data as GD

# 1) 谐波回归拟合与周期搜索
price = GD.gen_2014c()["price"]
K, r2 = GD._period_search(price, 24, 48)
beta, fit, sres = GD._harmonic_fit(price, K)
print("识别周期 K =", K, "月,  R² =", round(r2, 4), ", 残差σ =", round(sres, 3))

# 2) 留出验证：前48月训练，后12月测试
def mape(true, pred):
    return 100.0 * sum(abs(true[i] - pred[i]) / true[i] for i in range(len(true))) / len(true)

def harmonic_pred(train, Kk, n_fut):
    bt, _, _ = GD._harmonic_fit(train, Kk)
    nt = len(train)
    return [bt[0] + bt[1] * (nt + k) + bt[2] * GD.math.cos(2 * GD.math.pi * (nt + k) / Kk)
            + bt[3] * GD.math.sin(2 * GD.math.pi * (nt + k) / Kk) for k in range(n_fut)]

trn, tes = price[:48], price[48:]
m_har = mape(tes, harmonic_pred(trn, K, 12))
print("谐波回归 留出 MAPE =", round(m_har, 2), "%   (远优于 Holt/MA/ES)")

运行上述代码可复现本文核心数字：周期 $K=36$ 、拟合优度 $R^2=0.9944$ 、谐波预测 MAPE 约 1.21%，与正文、配图及数据生成工具完全一致。