简介：面向Java开发者与时间序列分析初学者的ARMA、ARIMA模型实践例程，围绕自回归移动平均及差分整合模型的代码实现展开，可用于趋势性、周期性数据的拟合与预测练习。资源包共43个文件，约8.83MB，以12个Java源文件和12个class编译文件为主体，另含5个jar依赖、Eclipse项目配置、示例数据与说明文档，并配套SPSS时间序列分析教程、用户数预测案例及ARMA自相关函数讲义等资料，便于对照理解模型原理。已有1679人学习，说明其在相关学习中具有一定参考价值。代码按ARMA与ARIMA两个工程组织，借助Jama等科学计算库完成参数计算，包含入口程序、数据集和运行说明，读者可据此调试阶数选择、差分处理、模型拟合与预测流程，也能将实现思路迁移到自己的Java项目中，适合作为时间序列建模的入门范例和工程复用素材。
1. 预测曲线并不来自拟合：ARIMA 拟合的是误差的记忆
把一份真实序列丢进 ARMA 模型或 ARIMA 模型的 Java 实现里跑一遍，多数人第一反应是失望：预测线几乎是平的，还带一点点缓慢衰减的弧度。但它对“下一期误差会偏向哪边”的判断，往往比任何曲线外推都准。原因在于，这类模型从一开始就没打算描述序列的形状，它估计的是过去的扰动以什么权重留在当下：自回归项把历史观测值拉回来，滑动平均项把历史误差拉回来，差分只是为了让这两股力有一个稳定的落点。所以同一份数据换个差分阶数，预测能差出几十个点；也所以在那些 Java 实现里，你会看到一堆矩阵递推、数值优化器和参数约束变换——它们都在服务于同一个目标：把不可观测的误差序列反解出来。搞懂这条链路，才算真正看懂这份例程。
2. 平稳性：为什么差分这一步不能省，过差分又为什么翻车
2.1 单位根：系数逼近 1 时方差在悄悄爆炸
先看最简单的一阶自回归。设 $x_t = c + \varphi x_{t-1} + \varepsilon_t$，$\varepsilon_t$ 是方差为 $\sigma^2$ 的白噪声。把式子反复展开，$x_t$ 会写成 $\sum_{j=0}^{t-1}\varphi^j\varepsilon_{t-j}$ 加常数项，于是方差是 $\sigma^2\sum_{j=0}^{t-1}\varphi^{2j}$。$\varphi$ 严格小于 1 时这个和收敛到 $\sigma^2/(1-\varphi^2)$，方差有上界，序列是平稳的；$\varphi = 1$ 时它等于 $t\sigma^2$，方差随时间线性增长，序列成了随机游走，任何基于“历史均值”的估计都会失效。
关键在于“逼近 1”这个过程并不温和。$\varphi = 0.99$ 时理论方差是 $50\sigma^2$，一次冲击的衰减半衰期是 $\ln 0.5 / \ln 0.99 \approx 69$ 期。也就是说，样本量不到 200 的时候，你根本观测不到它回归均值，样本自协方差矩阵接近奇异，最小二乘解出来的系数会把噪声当成信号放大。这就是单位根问题的本质：它属于参数估计的病态，而非假设检验里的概率判据问题。
放到 $p$ 阶的一般情形，平稳性条件写成特征多项式 $\Phi(z) = 1 - \varphi_1 z - \cdots - \varphi_p z^p$ 的所有根都落在单位圆外。Java 侧的实现不需要真的去解复根，工业代码里常用的替代做法是：估计出系数后反解多项式根的模长，只要有模长小于 $1+\epsilon$，就把参数往可行域内收缩，或者直接判定这组参数不可用。这一步在很多手写例程里被省掉，结果就是拟合在训练集上漂亮、外推三期就飞出去。
2.2 d 阶差分的数学含义与过差分的代价
差分的符号写法是 $\nabla x_t = x_t - x_{t-1} = (1-B)x_t$，$B$ 是延迟算子，$B x_t = x_{t-1}$。$d$ 阶差分就是 $(1-B)^d x_t$。它的作用很直接：一阶差分把线性趋势变成常数，二阶差分把二次趋势变成常数，本质上是通过多项式展开把低阶趋势项消掉。ARIMA 的 $d$ 就是在问一个问题——原始序列里有多少阶确定性趋势需要用差分抹平。
但 $d$ 不是越大越好。每多做一次差分，序列就少一个样本，同时把高频噪声放大一次。一份本来只有一阶趋势的数据，你做了二阶差分，等价于在 MA 部分硬塞进一个单位根：$\Theta(z)$ 有一个根落在单位圆上，模型不可逆，误差项无法用历史观测收敛地表达出来。具体表现是预测区间明显变宽、参数估计的方差膨胀、AIC 反而下降但外推效果更差。这就是过差分，它和差分不足一样是常见错误，只是方向相反。
判断 $d$ 的工程做法，通常分两步。先看自相关函数：如果 ACF 衰减到零的速度极慢、几十期还挂在 0.9 以上，基本可以判定存在单位根，先差一次。差分后再看 ACF 是否快速截断或衰减。严格一点的做法是单位根检验，原假设“存在单位根”，统计量拒绝原假设才说明平稳。还要留意一种反向情况：序列过度差分后 ACF 会在滞后 1 期出现显著负值（约 -0.5），这是过差分的典型痕迹，看到它就退回去降一阶。
Java 里的差分实现很朴素，但还原逻辑是个坑。差分丢掉的尾部信息必须在预测阶段补回来，所以要缓存原始序列的最后 $d$ 个值：
public final class Differencer { /** d 阶差分，长度由 n 缩到 n-d */ public static double[] diff(double[] x, int d) { double[] cur = x.clone(); for (int k = 0; k < d; k++) { double[] next = new double[cur.length - 1]; for (int i = 1; i < cur.length; i++) { next[i - 1] = cur[i] - cur[i - 1]; // (1-B) 作用一次 } cur = next; } return cur; } /** * 把差分域上的预测值还原回原尺度 * @param diffForecast 差分域预测，长度 h * @param tail 原始序列最后 d 个值，升序，长度 d * @param d 差分阶数 */ public static double[] integrate(double[] diffForecast, double[] tail, int d) { double[] cur = diffForecast.clone(); for (int k = 0; k < d; k++) { // 每还原一阶，需要上一阶差分序列的最后一个已知值作为起点 double last = tail[tail.length - 1 - k]; double[] next = new double[cur.length]; for (int i = 0; i < cur.length; i++) { last += cur[i]; next[i] = last; } cur = next; } return cur; } }integrate的循环顺序不能反：必须从最低阶差分开始逐级向上累加，每一级都要拿到该级序列在训练集末尾的最后一个真实值作为累加起点。tail就是为此准备的，按原始序列从早到晚取最后 $d$ 个。少存这一个数组，预测结果会整体偏移一个常数，而且偏移量随 $h$ 线性放大，越往后越离谱。
3. ARMA 的算子形式：可逆性决定了预测能不能收敛
3.1 延迟算子、格林函数与冲击的长期影响
ARMA(p,q) 的标准写法是 $x_t = c + \sum_{i=1}^{p}\varphi_i x_{t-i} + \varepsilon_t + \sum_{j=1}^{q}\theta_j \varepsilon_{t-j}$，用算子写就是 $\Phi(B)x_t = c + \Theta(B)\varepsilon_t$，其中 $\Phi(B) = 1 - \varphi_1 B - \cdots - \varphi_p B^p$，$\Theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q$。
平稳的前提下，$\Phi(B)$ 可逆，可以把 $x_t$ 写成无穷阶滑动平均：$x_t = \mu + \Psi(B)\varepsilon_t$，系数 $\psi_j$ 叫格林函数，也就是脉冲响应——第 $t$ 期一个单位的冲击，在第 $t+j$ 期还残留多少。递推关系是把 $\Phi(B)\Psi(B) = \Theta(B)$ 两边按 $B$ 的幂次对齐：
$$\psi_j = \sum_{i=1}^{\min(j,p)} \varphi_i \psi_{j-i} + \theta_j \cdot \mathbb{1}[j \le q]$$
这个数组的用途比看上去多。方差 $\mathrm{Var}(x_t) = \sigma^2\sum_{j\ge0}\psi_j^2$，可以直接校验估计出来的参数是否合理；$h$ 步预测误差的方差是 $\sigma^2\sum_{j=0}^{h-1}\psi_j^2$，预测区间宽度就靠它算；$\psi_j$ 的衰减速度还能反过来验证平稳性——如果 $\psi_j$ 到 50 期还没衰减到 1e-3 以下，说明特征根离单位圆太近，参数不可信。
/** * 由 ARMA 参数递推格林函数，长度 trunc 截断 * psi[0] = 1 */ public static double[] greenFunction(double[] phi, double[] theta, int trunc) { double[] psi = new double[trunc]; psi[0] = 1.0; for (int j = 1; j < trunc; j++) { double v = (j <= theta.length) ? theta[j - 1] : 0.0; // 注意 theta 是 0-based for (int i = 1; i <= Math.min(j, phi.length); i++) { v += phi[i - 1] * psi[j - i]; } psi[j] = v; } return psi; }索引对齐是这段代码最容易写错的地方：phi[i-1]对应 $\varphi_i$，theta[j-1]对应 $\theta_j$，因为数学写法从 1 开始编号而 Java 数组从 0 开始。写错一位，脉冲响应看起来依然“在衰减”，但方差校验会差出一个数量级，排查起来非常费时间。
3.2 可逆性：MA 多项式根跑出单位圆会发生什么
AR 部分的平稳性保证 $\Psi(B)$ 收敛，MA 部分的可逆性保证反过来的方向也收敛——$\varepsilon_t$ 能写成历史观测的无穷级数 $\Pi(B)x_t$。条件同样简单：$\Theta(z) = 1 + \theta_1 z + \cdots + \theta_q z^q$ 的根全在单位圆外。
为什么这一条对实现如此关键？条件最小二乘的残差是用递推算出来的：$\hat\varepsilon_t = x_t - c - \sum\varphi_i x_{t-i} - \sum\theta_j \hat\varepsilon_{t-j}$。这是一个反馈回路，$\hat\varepsilon_t$ 依赖 $\hat\varepsilon_{t-1}$。如果 $|\theta_1| > 1$，这个回路的增益大于 1，初值的一点误差会被逐期放大，几十期后残差就到 $10^{15}$ 量级，目标函数变成 NaN，优化器直接放弃。很多手写例程在这里失败，报出来的现象是“优化不收敛”，实际原因是参数越界后的数值爆炸。
解决办法不是加个 try-catch 兜住，而是在参数化层面做约束。常见做法是把待优化变量 $u$ 映射到可行域，比如 $\theta = 0.98\tanh(u)$ 保证 $|\theta| < 0.98$，$\varphi$ 用同样的方式逐项约束；更严格的做法是约束特征根的模长，用参数化保证根始终在单位圆外。$0.98$ 而不是 $1.0$ 是为了留出数值余量，让 $\theta$ 逼近边界时梯度不至于消失得太快。
还有一个隐蔽的退化情形：AR 部分和 MA 部分发生公因子抵消。比如 $\Phi(B) = (1-0.5B)(1-0.8B)$、$\Theta(B) = (1-0.5B)$，真实模型其实是 AR(1)，但你按 ARMA(2,1) 去估，会得到一组“等价但冗余”的参数，标准误很大且不唯一。检验方法是算两个多项式的根，看有没有接近重合的；工程上更省事的做法是对比 $p$、$q$ 都减一之后的 BIC，如果几乎不变，就选参数更少的那个。
4. Java 里的参数估计：从 Yule-Walker 方程到条件最小二乘
4.1 自协方差、Toeplitz 结构与 Levinson-Durbin 递推
纯 AR 模型的参数有闭式解，这是所有 Java 实现都会先走的路径。把 $x_t$ 两边乘上 $x_{t-h}$ 取期望，利用 $\varepsilon_t$ 与历史不相关，得到 Yule-Walker 方程 $\gamma(h) = \sum_{i=1}^{p}\varphi_i\gamma(h-i)$，写成矩阵形式是 $\Gamma_p\boldsymbol\varphi = \boldsymbol\gamma_p$，$\Gamma_p$ 是 $p \times p$ 的自协方差矩阵，第 $(i,j)$ 项是 $\gamma(|i-j|)$——这就是 Toeplitz 结构，每条对角线上的元素相同。
理论上可以直接做矩阵求逆，但 $O(p^3)$ 且数值不稳。Levinson-Durbin 递推利用 Toeplitz 结构把复杂度降到 $O(p^2)$，而且顺带产出两个极有用的中间量：反射系数 $\kappa_k$ 就是偏自相关系数（PACF），$\sigma_k^2$ 就是 AR($k$) 模型的白噪声方差估计，后面算 AIC 直接拿来用。
样本自协方差的估计式有两种，$\frac{1}{n}\sum_{t}(x_t-\bar x)(x_{t+h}-\bar x)$ 和除以 $n-h$。前者保证 $\Gamma_p$ 半正定，后者是无偏但可能让矩阵不正定。实现里通常用前者，代价是长滞后阶的估计有偏，这个偏差在样本量小的时候会影响定阶，值得在代码注释里留一句。
/** * Levinson-Durbin 递推解 Yule-Walker 方程 * @param gamma 自协方差 gamma[0..p]，gamma[0] 是方差 * @return 长度 p+1 的数组：前 p 个是 phi_1..phi_p，最后一个是 sigma2 */ public static double[] levinsonDurbin(double[] gamma) { int p = gamma.length - 1; double[] phi = new double[p + 1]; // 1-based double[] prev = new double[p + 1]; double sigma2 = gamma[0]; for (int k = 1; k <= p; k++) { double acc = gamma[k]; for (int j = 1; j < k; j++) { acc -= prev[j] * gamma[k - j]; } double kappa = acc / sigma2; // 反射系数 = PACF(k) if (Math.abs(kappa) >= 1.0) { // 数值越界，强行截断 kappa = Math.signum(kappa) * 0.999; } phi[k] = kappa; for (int j = 1; j < k; j++) { phi[j] = prev[j] - kappa * prev[k - j]; } sigma2 *= (1.0 - kappa * kappa); System.arraycopy(phi, 0, prev, 0, p + 1); } double[] out = new double[p + 1]; System.arraycopy(phi, 1, out, 0, p); out[p] = sigma2; return out; }这段递推有个天然优势：循环每走一步就得到一个 AR($k$) 模型的解。所以定阶的时候不用为每个 $p$ 重跑一遍，直接跑一次 $p_{\max}$ 的递推，把每一步的 $\sigma_k^2$ 记下来，AIC 序列就齐了。$\kappa$ 的截断处理是为了防止自协方差矩阵接近奇异时 $\sigma^2$ 被推到负数，进而开根号出 NaN，这是小样本下必然遇到的补丁。
4.2 条件最小二乘：MA 部分没有闭式解，只能上优化器
一旦 $q > 0$，问题变成非线性的，没有闭式解。三种常见路径：
第一种是条件最小二乘（CSS），给定前 $p$ 个观测做初值、前 $q$ 个残差取 0，逐期递推残差，最小化 $\sum\hat\varepsilon_t^2$。目标函数直接、实现简单，缺点是忽略了前 $p$ 个观测的似然贡献，样本短时有偏。
第二种是精确极大似然，用 Kalman 滤波或者协方差矩阵直接算似然，无偏但实现复杂度高一个量级，且每评估一次要算一次 $n \times n$ 协方差，长序列上慢得没法用。
第三种是 Hannan-Rissanen 三步法：先拟合一个长阶 AR（阶数取 $\lceil\sqrt{n}\rceil$ 左右）拿残差，再用残差作为误差项的代理，把 MA 部分转成线性回归，最后用第一步的估计做初值再跑一轮 CSS 精修。工程代码里这个思路性价比极高，因为它把非线性优化变成了“线性回归 + 少量精修”，对初值不敏感。
不管走哪条路，最终都落到一个无约束优化器上。Java 里常用的两种：Nelder-Mead 单纯形法，不需要梯度，代码量小，缺点是维度高时收敛慢；BFGS 需要梯度，可以用数值差分近似，收敛快但每步要评估 $O(k)$ 次目标函数。对于 $p+q \le 10$ 的常规场景，Nelder-Mead 配合参数变换已经够用。
/** * 条件最小二乘目标函数：给定参数返回残差平方和 * x 必须已中心化（均值为 0），否则 c 会和均值混淆 */ public static double css(double[] x, double c, double[] phi, double[] theta) { int n = x.length, p = phi.length, q = theta.length; double[] e = new double[n]; double sse = 0.0; for (int t = p; t < n; t++) { double pred = c; for (int i = 0; i < p; i++) { pred += phi[i] * x[t - 1 - i]; // AR 项用真实观测 } for (int j = 0; j < q; j++) { int idx = t - 1 - j; pred += theta[j] * (idx >= 0 ? e[idx] : 0.0); // 前 q 期残差置 0 } e[t] = x[t] - pred; sse += e[t] * e[t]; } return sse; }三个细节值得说明。循环从t = p开始而不是t = 0，因为前 $p$ 个观测凑不出完整的 AR 项；x必须已中心化，否则常数项 $c$ 和序列均值 $\bar x$ 会互相抵消，估出来的 $c$ 没有物理意义；残差数组e在t < q的位置保持默认的 0，这是条件似然的“条件”所在。
优化器外层还要套一层参数变换。设待优化向量 $u$，映射关系是 $\varphi_i = 0.98\tanh(u_i)$、$\theta_j = 0.98\tanh(u_j)$、$c$ 不做约束。这样无论优化器把 $u$ 推到多大，参数永远落在可行域内，CSS 递推不会炸。代价是目标函数在边界附近变得平坦，梯度消失，所以初值的选择变重要了——用矩估计或者 Hannan-Rissanen 的结果做初值，能少走很多弯路。
4.3 AIC、BIC 与网格搜索：定阶的工程取舍
参数估完了，下一个问题是 $p$ 和 $q$ 取多少。Akaike 信息准则的形式是 $\mathrm{AIC} = n\ln\hat\sigma^2 + 2k$，贝叶斯信息准则是 $\mathrm{BIC} = n\ln\hat\sigma^2 + k\ln n$，$\hat\sigma^2$ 是残差方差，$k$ 是自由参数个数。
$k$ 怎么数是个容易搞错的地方。ARMA($p$,$q$) 带常数项时 $k = p + q + 1$；不带常数项是 $p+q$；如果模型里还有漂移项或者季节项，还要往上涨。数错一个，AIC 就偏向另一侧：少算 $k$ 会系统性地选出高阶模型，多算则选出低阶。常见做法是显式写一个paramCount(p, q, hasConst, hasDrift, seasonalPeriod)方法，把这个逻辑集中在一处，别散在循环里。
AIC 和 BIC 的取舍在于惩罚力度。$\ln n > 2$ 时（$n > 8$）BIC 的惩罚更重，偏好更简洁的模型。样本短、信噪比低的时候 AIC 容易选出高阶模型，在测试集上表现反而更差；样本长、真实模型阶数确实高的时候，BIC 会欠拟合。实用的做法是两个都算，取 AIC 和 BIC 的交叉区间，再留最后 20% 的数据做滚动验证，看哪组参数的外推误差最小。
搜索空间一般限定在 $p \le 5$、$q \le 2$，更准确地说，$p+q$ 的乘积级搜索在 $n$ 不大时没有统计意义，因为高阶模型的参数标准误会大到没法解释。真正需要高 $p$ 的场景（比如长记忆序列），通常应该换成分数阶差分或者谱方法，而不是硬堆 AR 阶数。
5. 预测递推与排查：落地时最容易出问题的几处
5.1 多步预测：递推公式本身不难，容错全在边界条件
$h$ 步预测的公式和 CSS 的残差递推几乎是同一个式子：$\hat x_{T+h} = c + \sum_{i=1}^{p}\varphi_i \hat x_{T+h-i} + \sum_{j=1}^{q}\theta_j \hat\varepsilon_{T+h-j}$。区别只在数据来源：下标 $T+h-i \le T$ 时用真实观测，大于 $T$ 时用前一步的预测值；误差项下标大于 $T$ 的全部取 0，因为未来误差的期望是 0。
/** * 递推多步预测 * @param x 训练序列（已中心化） * @param resid 训练期残差，长度与 x 一致 * @param steps 预测步数 */ public static double[] forecast(double[] x, double c, double[] phi, double[] theta, double[] resid, int steps) { int n = x.length, p = phi.length, q = theta.length; double[] ext = new double[n + steps]; // 拼接历史与预测 System.arraycopy(x, 0, ext, 0, n); double[] e = new double[n + steps]; System.arraycopy(resid, 0, e, 0, n); for (int h = 0; h < steps; h++) { int t = n + h; double pred = c; for (int i = 0; i < p; i++) { pred += phi[i] * ext[t - 1 - i]; // 未来位置已经是预测值 } for (int j = 0; j < q; j++) { int idx = t - 1 - j; pred += theta[j] * (idx < n ? e[idx] : 0.0); // 未来残差为 0 } ext[t] = pred; } double[] out = new double[steps]; System.arraycopy(ext, n, out, 0, steps); return out; }用一个统一数组ext承接历史和预测，是这段代码最省心的地方——AR 项和 MA 项的索引逻辑不用分两套写，边界判断只剩残差那一个idx < n。预测区间宽度按 $\hat\sigma^2\sum_{j=0}^{h-1}\psi_j^2$ 累加，用上一节算好的格林函数直接取前缀和即可，注意 $h$ 步用的是前 $h$ 项而不是全部。
差分还原接在后面：integrate(forecast, tail, d)。这两步之间还差一个反中心化，把估计出的均值加回去。这里经常出现整体偏移，原因下一节细说。
5.2 排查：五个反复出现的失败现场
预测线一路向上发散，第三期开始离谱。现象是差分域预测还挺正常，还原到原尺度后斜率越来越大。原因八成是integrate的起点取错了，或者tail存的是原始序列最后 $d$ 个值但顺序反了。解决方式是把tail的构造和diff的调用放在同一处，写完立刻用d阶差分再d阶还原做一次恒等测试，输出和输入完全一致才算过。
AIC 选出的最优阶数在留出集上明显更差。现象是训练集残差小、测试集误差大，阶数集中在 $p=4$、$q=2$ 这类偏高位置。原因通常是样本量不够而 AIC 惩罚偏轻，或者是差分阶数 $d$ 选小了，模型靠堆 AR 阶数去补趋势。解决方式是换 BIC 重跑一遍，同时检查 ACF 尾部的衰减速度重新确认 $d$；两个准则给出的推荐阶数相差超过 2 时，优先信 BIC。
优化器不收敛，目标函数返回 NaN。现象是迭代几十次后 CSS 值变成NaN或者 $10^{20}$ 量级。原因基本是 MA 系数越界导致残差递推的反馈回路发散，或者初值离真值太远。解决方式是加 $\tanh$ 参数变换把 $\theta$ 约束在 $(-0.98, 0.98)$，同时换成多起点策略：矩估计初值 + 随机扰动各跑一次，取目标函数最小的那组。
预测值和真实值差一个恒定常数。现象是形状对得很准，整体上下平移了一段。原因是序列中心化之后忘了在预测结果上把均值加回来，或者常数项 $c$ 和均值 $\bar x$ 混在一起重复加了一次。解决方式是明确约定：内部计算全部在中心化空间，出结果前统一加回均值，且模型里如果已经有 $c$ 就不再额外加。
长序列上单次拟合要跑几十秒。现象是 $n$ 超过 5000 后，每评估一次目标函数都很慢。原因常常是每次评估都重算了一整遍样本自协方差，而自协方差其实与参数无关，整个优化过程只需要算一次。解决方式是把自协方差、Toeplitz 矩阵的 Cholesky 分解提前算好缓存起来，目标函数只做残差递推，实测能快出一个数量级。
6. 设计取舍：手写 ARIMA 求解器的边界在哪
写完这一套，很少有人会再坚持“所有东西都自己实现”。原因主要落在成本结构上，跟能力无关：ARIMA 的每一部分单看都不难，难点在它们的组合——Yule-Walker 要求自协方差矩阵正定，CSS 要求参数在可行域内，预测还原要求差分和反差分严格对称，任何一处松动都会在另一处以“效果不好”的形式暴露出来，且难以定位。这份例程的价值在于把这条链路拆开摆平，让每个环节的输入输出和边界条件都看得见，这比自己从头造一个更完整的库要划算得多。
一个具体的取舍是优化器的选择。Nelder-Mead 写三十行就能跑，无梯度、对目标函数的光滑性没要求，代价是 $p+q$ 超过 8 之后收敛慢到不可接受。换成 BFGS 要写数值梯度，步长的选取又变成新问题，太小会被浮点噪声淹没，太大则方向不准。多数工程场景里 $p+q \le 5$，这个取舍根本不需要纠结，直接用简单的那个；只有上到季节性模型、参数维度翻几倍的时候，才值得为优化器单独投入。
另一个取舍在似然的精度上。CSS 丢掉前 $p$ 个观测，样本 200 以上时这点损失可以忽略，样本 50 以下时偏差明显。真实业务里更常见的约束是数据量而不是算力，所以与其去实现 Kalman 滤波的精确似然，不如把精力放在数据预处理和阶数选择上——实践里，“$d$ 选对了”带来的收益远超“似然从条件换成精确”。
还有一个常被忽视的扩展成本：从 ARMA 到 ARIMA 只加了一个差分，但到 SARIMA 要引入季节算子 $(1-B^s)$，参数数量乘上季节周期，Yule-Walker 的 Toeplitz 结构不再是简单循环，Levenson-Durbin 得换成块 Toeplitz 版本；再到 ARIMAX，外生变量进入回归项，参数估计要和回归系数联合优化。每加一层，手写代码的维护成本非线性增长。所以合理的做法是把这份例程当作理解工具和对照基准，生产环境里换用成熟库，但保留自己实现的那份当“标准答案”——当库的输出可疑时，拿自己的实现跑同一组数据，差异出现在哪一步一目了然。
我最早写这套东西的时候，踩的坑是直接在原始尺度上拟合，没有中心化，结果常数项估出来是个极小的数，我还以为是优化器没收敛，调了两天步长和多起点，最后发现是均值没扣掉，$c$ 和 $\bar x$ 在最小二乘里互相抵消了。从那以后我每次写时间序列的拟合代码，都强制先跑一遍“中心化 → 拟合 → 加回均值”的闭环校验，再做任何参数解读。希望帮到你。
本文还有配套的精品资源，点击获取
