简介ARMA与ARIMA模型是时间序列分析与预测的核心工具这份Java实现例程面向需要做数据预测、趋势分析的程序员和算法学习者解决从理论到工程落地的衔接问题。压缩包约8.83MB共43个文件12个Java源文件与12个class文件构成两个工程主体5个Jar包含Jama-1.0.2.jar等负责矩阵运算与统计计算README与txt文档给出使用说明docx、pdf及ppt则补充SPSS时间序列分析与自相关函数讲解。资源在CSDN已有1679人浏览学习侧面印证其实用价值。通过运行main函数、阅读源码可以直观掌握ARMA、ARIMA的参数设定、差分处理、拟合与预测流程配套数据文件便于测试验证与实际调参。作者将ARMA-Java--master与ARIMA-master两个工程一并打包目录结构清晰适合作为课程设计、毕业设计及业务预测模块的参考蓝本无论理解原理还是迁移到自有数据集都很有帮助。1. 为什么在 Java 里手写 ARMA/ARIMA这不是复古是现实需求ARMA 与 ARIMA 模型的 Java 实现例程听起来像是教科书作业但真落到生产里它解决的是一个很现实的问题很多公司的核心服务跑在 Java 上实时数据从消息队列进来预测逻辑必须和业务代码同进程不能因为一个时间序列模型就单独开一条 Python 服务。我遇到过几次这种场景最后都是把 ARMA/ARIMA 的核心算法用 Java 重写一遍工程量远没有想象中大。这篇不是复刻某个开源库而是把 ARMA/ARIMA 从公式到 Java 代码的完整路径走一遍先讲清楚两个模型的数学边界再给出可运行的估计例程包括差分还原、自动定阶和预测区间最后是几个我用 Java 实现时踩过的坑。适合两类读者一是要在 Java 服务里嵌入时序预测的开发者二是想搞懂模型内部原理、不想把算法当黑匣子调包的人。2. 先把模型底细说清ARMA/ARIMA 的公式、条件与 Java 选型2.1 四个模型的关系AR、MA 是零件ARMA 是组装ARIMA 是预处理在写代码之前先把数学上的边界划清楚。AR(p) 是自回归当前值由过去 p 个观测值线性组合决定MA(q) 是移动平均当前值由过去 q 个白噪声项的线性组合决定。ARMA(p,q) 把两者拼在一起而 ARIMA(p,d,q) 就是先对原始序列做 d 阶差分让非平稳序列变成平稳序列再套 ARMA。用表把四个模型的关系摆出来更直观模型数学形式适用场景AR(p)x_t c Σφ_i x_{t-i} ε_t序列自相关明显且 ACF 衰减缓慢MA(q)x_t μ Σθ_i ε_{t-i} ε_tACF 在 q 阶后截断短期冲击效应ARMA(p,q)x_t c Σφ_i x_{t-i} Σθ_i ε_{t-i} ε_t平稳序列自相关和短期噪声同时存在ARIMA(p,d,q)对 x_t 做 d 阶差分后套 ARMA非平稳、有趋势或随机游走特征写 Java 实现时我一般会把 AR、MA、ARMA 的估计逻辑收在一个类里差分逻辑单独放一个工具类。因为差分是数据预处理不应该污染模型估计的代码路径。2.2 平稳性与可逆性Java 代码必须校验的边界条件ARMA 模型的估计有一个前提序列必须平稳。所谓平稳直观说就是均值和方差不随时间变化。数学上AR 部分的特征方程 1 - Σφ_i z^i 0 的根必须落在单位圆外MA 部分的可逆性条件则是 1 Σθ_i z^i 0 的根在单位圆外。这个条件不是理论摆设。在用 Java 手写估计时如果解出来的 φ_i 绝对值大于 1预测值会指数爆炸几期之后直接溢出成 NaN。所以我在例程里会加一个校验方法解完 Yule-Walker 方程后检查 AR 系数对应的特征根模长是否大于 1不满足就抛异常或者提示先做差分。判断平稳性还有一种工程上更省事的做法直接看数据ACF 衰减很慢、拖尾拖得很长基本可以判定非平稳先差分一次再说。Java 例程里我习惯写一个简单的isStationaryByACF方法算前 20 阶自相关系数如果有超过 80% 的阶数绝对值大于 0.2就建议差分。2.3 Java 选型为什么我不建议直接引一个重框架Java 生态里不是没有时间序列库但真正落地时选择很少。常见的开源数学库提供的是矩阵运算、概率分布这些底层能力而完整封装好的 ARMA 估计器要么文档不完整要么把估计过程包成黑盒出了问题根本没法查。我的做法是用 Commons Math 提供矩阵分解和正态分布分位数模型估计部分自己写控制在两百行左右。dependency groupIdorg.apache.commons/groupId artifactIdcommons-math3/artifactId version3.6.1/version /dependency这个依赖在整个例程里主要干两件事一是解 Yule-Walker 方程组的线性系统二是预测区间里计算正态分布的 1.96 分位数。如果不想引依赖自己写高斯消元也够用但 Commons Math 的 LU 分解在数值稳定性上更可靠尤其是 p 阶数到 5 以上时手写消元容易因为矩阵病态得到离谱系数。3. 手写 ARMA 估计例程Yule-Walker 求解与 MA 参数近似3.1 第一步用 Yule-Walker 方程解 AR 部分ARMA 参数估计的完整解法是最大似然但对例程来说最稳的切入点是 Yule-Walker 方程。它先把数据去均值然后计算自协方差序列最后解一个线性方程组得到 AR 系数。这个方法的优点是稳定、无迭代、不用调优化参数缺点是估计精度不如最大似然但在样本量大于 200 时差距很小。import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.DecompositionSolver; import org.apache.commons.math3.linear.LUDecomposition; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; private double[] yuleWalker(double[] data, int p) { int n data.length; double mean 0.0; for (double v : data) mean v; mean / n; // 去均值自协方差计算只关心波动不关心水平 double[] x new double[n]; for (int i 0; i n; i) x[i] data[i] - mean; double[] acov new double[p 1]; for (int k 0; k p; k) { double sum 0.0; int cnt 0; for (int i 0; i n - k; i) { sum x[i] * x[i k]; cnt; } acov[k] sum / cnt; } // Yule-Walker 方程R * phi rR 是 Toeplitz 矩阵 RealMatrix A new Array2DRowRealMatrix(p, p); RealVector b new ArrayRealVector(p); for (int i 0; i p; i) { b.setEntry(i, acov[i 1]); for (int j 0; j p; j) { A.setEntry(i, j, acov[Math.abs(i - j)]); } } DecompositionSolver solver new LUDecomposition(A).getSolver(); RealVector phi solver.solve(b); return phi.toArray(); }这段代码里有两个细节值得注意。一是自协方差的分母用n - k而不是n这是有偏估计但小样本下偏差更小Yule-Walker 本来就有偏没必要刻意修正。二是矩阵 A 的构造方式第 i 行第 j 列是acov[|i-j|]这就是 Toeplitz 结构手写时最容易搞错的是索引偏移写错一个下标解出来的系数符号直接反掉。解出来的 phi 数组就是 AR(p) 的系数。如果数据本身是纯 AR 过程到这里参数估计已经完成可以直接拿去做预测。但如果是 ARMA还需要继续估计 MA 部分。3.2 第二步用残差近似估计 MA 参数有了 AR 系数可以先算 AR 过滤后的残差序列再用残差的自相关结构反推 MA 系数。严格做法是最大似然或条件最小二乘但例程里我用不动点迭代把误差控制在可接受范围。private double[] computeResiduals(double[] data, double[] phi) { int n data.length; int p phi.length; double[] resid new double[n]; for (int t p; t n; t) { double pred 0.0; for (int i 0; i p; i) { pred phi[i] * data[t - i - 1]; } resid[t] data[t] - pred; } return resid; } private double[] autocorrelation(double[] resid, int maxLag) { int n resid.length; double mean 0.0; for (double v : resid) mean v; mean / n; double[] ac new double[maxLag 1]; double variance 0.0; for (double v : resid) variance (v - mean) * (v - mean); variance / n; for (int k 0; k maxLag; k) { double sum 0.0; int cnt 0; for (int i 0; i n - k; i) { sum (resid[i] - mean) * (resid[i k] - mean); cnt; } ac[k] sum / cnt / variance; } return ac; } private double[] estimateMA(double[] data, double[] phi, int q) { double[] resid computeResiduals(data, phi); double[] theta new double[q]; double[] ac autocorrelation(resid, q 1); // 初值MA(1) 时 theta 约等于 -ac[1] for (int k 0; k q; k) { theta[k] -ac[k 1]; } // 不动点迭代 5 次让残差自相关趋近于零 for (int iter 0; iter 5; iter) { double[] filtered maFilter(resid, theta); double[] nac autocorrelation(filtered, q 1); for (int k 0; k q; k) { theta[k] - 0.4 * nac[k 1]; } } return theta; } private double[] maFilter(double[] resid, double[] theta) { int n resid.length; double[] filtered new double[n]; for (int t 0; t n; t) { double v resid[t]; for (int k 1; k theta.length t - k 0; k) { v theta[k - 1] * filtered[t - k]; } filtered[t] v; } return filtered; }这段代码的思路是先把数据里 AR 部分剥离干净剩下的残差应该只含 MA 结构然后通过迭代调整 theta让过滤后的序列自相关尽量接近零。迭代步长 0.4 是我试过比较稳的值太大会震荡太小收敛慢。边界要说清楚这个方法对低阶 MAq 不超过 3效果够用q 大了或者数据信噪比低时还是会和最大似然结果有偏差生产环境建议换成梯度优化。3.3 把拟合结果封装起来AIC 与残差方差一起输出估计完参数不能直接结束还要算两个东西残差方差 sigma2 和 AIC 值。前者用于预测区间后者用于后面第 4 章的自动定阶。我习惯用一个简单的内部类把结果攒起来。class FitResult { double[] phi; double[] theta; double sigma2; double aic; int p, q; } private FitResult fitArma(double[] data, int p, int q) { FitResult result new FitResult(); double[] phi yuleWalker(data, p); double[] theta estimateMA(data, phi, q); double[] resid computeResiduals(data, phi); double sigma2 0.0; int n data.length; for (double v : resid) sigma2 v * v; sigma2 / (n - p); result.phi phi; result.theta theta; result.sigma2 sigma2; result.p p; result.q q; result.aic n * Math.log(sigma2) 2 * (p q 1); return result; }AIC 公式里的1是常数项代表均值估计消耗的一个自由度。网格搜索定阶时AIC 越小模型越优但后面第 5 章会专门讲一个 AIC 的坑它偏爱复杂模型小样本下经常选出过大的 p、q。4. ARIMA 落地差分还原、自动定阶与预测区间4.1 差分与还原记住你的数据被改过ARIMA 和 ARMA 的唯一差别就是那层差分。差分本身简单但还原逻辑非常容易写错。下面这段例程处理 d1 的差分和还原d2 时按同样的逻辑递推即可。private double[] difference(double[] data, int d) { double[] current data; for (int k 0; k d; k) { double[] next new double[current.length - 1]; for (int i 1; i current.length; i) { next[i - 1] current[i] - current[i - 1]; } current next; } return current; } private double[] restoreForecast(double[] diffForecast, double[] original) { int n original.length; double last original[n - 1]; double[] restored new double[diffForecast.length]; for (int i 0; i diffForecast.length; i) { if (i 0) { restored[i] diffForecast[i] last; } else { restored[i] diffForecast[i] restored[i - 1]; } } return restored; }还原的要点是预测的是差分序列的下一步要变成原始序列的下一步必须把差分最后一步用到的那个原始值加回来。之后每一步都累加前面已还原的预测值。很多初学者直接在循环里加原始序列最后一位导致从第 2 步预测开始误差叠加越预测越偏。4.2 自动定阶用网格搜索代替肉眼读 ACF/PACF看图定阶是统计学教材教的方法但 Java 例程里更适合写一个网格搜索遍历 p 和 q用 AIC 选最优。这里要注意定阶必须在训练集上做测试集不能参与否则就是数据泄露。private int[] selectOrder(double[] trainData, int maxP, int maxQ, int d) { double[] diffData difference(trainData, d); double bestAic Double.MAX_VALUE; int bestP 0, bestQ 0; for (int p 1; p maxP; p) { for (int q 0; q maxQ; q) { FitResult fit fitArma(diffData, p, q); if (fit.aic bestAic) { bestAic fit.aic; bestP p; bestQ q; } } } return new int[]{bestP, bestQ}; }maxP 和 maxQ 的经验值样本量 500 以下时两者都别超过 5样本量 2000 以上时可以放到 8。搜索范围过大的后果在第 5 章详细说。定阶完成后用选出来的 p、q 对整个训练集重新 fit 一次然后用参数做预测。4.3 预测与置信区间用无穷 MA 表示算误差方差ARMA 模型的 h 步预测误差方差可以通过把模型改写成无穷 MA 形式来计算。这个改写在代码里就是递推计算 psi 权重工程上足够用。private double[] psiWeights(double[] phi, double[] theta, int steps) { int p phi.length; int q theta.length; double[] psi new double[steps]; psi[0] 1.0; for (int j 1; j steps; j) { double s 0.0; for (int i 1; i Math.min(j, p); i) { s phi[i - 1] * psi[j - i]; } if (j q) { s theta[j - 1]; } psi[j] s; } return psi; }这里有个容易写错的地方递推公式右边不是简单地把 phi 和 theta 都卷进去而是利用了 Ψ(B)Φ(B) Θ(B) 的展开关系。theta 项只在前 q 步直接加系数本身之后不再出现。把公式记成“双卷积”的话算出来的 psi 权重会错得很隐晦。有了 psi 权重h 步预测的方差就是 sigma2 乘以 psi_0 到 psi_{h-1} 的平方和。95% 预测区间就是预测值加减 1.96 倍标准差。完整预测流程串起来长这样private void predict(double[] history, FitResult fit, int steps) { double[] psi psiWeights(fit.phi, fit.theta, steps); double[] pred new double[steps]; double[] var new double[steps]; for (int h 1; h steps; h) { pred[h - 1] 0.0; // 例程简化预测值需要按 ARMA 递推计算 double cumVar 0.0; for (int i 0; i h; i) { cumVar psi[i] * psi[i]; } var[h - 1] fit.sigma2 * cumVar; } }严格的多步预测要比上面这段复杂因为每一步都要先算出预测值再作为下一步的输入。例程里我给的是置信区间的计算骨架实际使用时把 ARMA 多步递推补进去即可。核心思想不变误差方差只看 psi 权重预测路径只影响均值不影响区间宽度。5. 避坑六个让 Java 时序例程翻车的细节5.1 估计出来的 AR 系数绝对值大于 1预测发散现象训练时一切正常一跑预测数值快速增大几期之后变成 Infinity 或者负数。原因输入序列没有差分或者差分阶数不够AR 部分特征根落在单位圆内模型本身就是非平稳的。解决在 fit 之后立即做一次平稳性检查判断特征根是否都在单位圆外不满足就提示先增加 d 阶数。不要等到预测结果爆掉才回头查数据。5.2 差分过头预测曲线被强拉回均值现象原始数据有明显趋势d 设成 2预测曲线初始值附近震荡两期后迅速回归均值AIC 还显示模型很“好”。原因过差分把原本的趋势信息也差分掉了模型只剩噪声结构预测自然没有方向性。解决d 的经验值不要超过 2每差一次就看一眼差分序列是否还有明显趋势以及方差是否变大。对大多数业务数据d1 就够。5.3 AIC 网格搜索选出超大 p、q模型过拟合现象maxP 和 maxQ 放到 10网格搜索每次都选 p9、q8 这种组合残差确实白噪声但样本外预测一塌糊涂。原因AIC 对小样本惩罚不足搜索范围越大越容易挑中复杂模型。解决把 maxP 和 maxQ 限制在样本量的十分之一以内另一种办法是搜完之后看系数是否有接近零的有就降阶。前面代码里的1常数项不要省省了 AIC 会更偏好复杂模型。5.4 MA 参数迭代震荡theta 忽正忽负现象estimateMA 方法跑完theta 值每次运行不一样或者迭代 5 次后残差自相关仍然很大。原因MA 估计的初值选得不好0.4 的步长对某些数据太激进或者 q 阶数超过了数据能支撑的复杂度。解决初值从-ac[k1]改成-0.7 * ac[k1]迭代步长降到 0.2同时把 q 上限改小比如 maxQ 不超过 3。如果还是不收敛直接用最大似然优化不要在例程级别死磕。5.5 滚动验证时误差突然爆炸但全样本拟合的误差很小现象用前 80% 数据训练后 20% 数据测试误差很大但用全样本重新拟合后看残差指标很好。原因数据里可能有突变点或结构变化全样本拟合把训练集和测试集的差异都“平均”掉了看似合理实则自欺。解决严格按时间顺序切分训练集永远只包含测试时刻之前的数据遇到误差爆炸就去看测试集在哪个时间点开始偏离结合业务找原因。5.6 预测区间窄得离谱几乎贴着预测线现象95% 置信区间很窄但实际误差经常超出区间。原因sigma2 用的是训练残差方差如果训练集正好避开了一段波动大的时期方差被低估另外 psi 权重递推里 theta 项写错也会导致区间过窄。解决把训练期拉长覆盖至少一个完整波动周期检查 psi 权重前几步是否明显大于 1如果全是 0.x 级别大概率递推公式写错了。6. 验证你的模型残差白噪声与滚动回测的习惯模型拟合完头一件事不是看预测曲线好不好看而是确认残差里没有剩余模式。工程上我习惯做一个简化的白噪声检验算残差前 20 阶自相关超过 1.96 / sqrt(n) 的阶数占比小于 5%就认为残差没有明显自相关。这个阈值来自正态近似虽然不如完整的卡方检验严谨但作为例程的快速体检足够。static boolean residualIsWhiteNoise(double[] resid) { int n resid.length; int maxLag Math.min(20, n / 4); double limit 1.96 / Math.sqrt(n); int exceed 0; double[] ac autocorrelation(resid, maxLag); for (int k 1; k maxLag; k) { if (Math.abs(ac[k]) limit) exceed; } return exceed maxLag * 0.05; }残差通过白噪声检验之后再做滚动回测。滚动回测的做法是固定一个窗口长度比如 500 个点每次用窗口内数据训练预测下一个点然后窗口向后滑动一步重新训练。这样每次预测用的都是“当时”能拿到的数据最能反映真实上线效果。一次性 train/test 切分只能作为快速验证不能作为最终结论。我做过不少次自以为模型没问题、一上线就被真实数据打脸的事最后几乎都是败在验证方式上不是用全样本拟合骗了自己就是预测区间算窄了。现在我的习惯是任何时序模型交付前必须过两道关第一道是残差白噪声第二道是滚动回测下的 MAE 和 RMSE 要稳定不能在某个月份突然翻倍。这两步花了不了一小时但能省掉上线后几天的排查时间。希望这篇例程能帮你在 Java 里把 ARMA/ARIMA 跑通也知道它跑到什么程度算真的能用了。本文还有配套的精品资源点击获取