1. 先回答一个问题时间序列预测凭什么轮得到线性回归打开任何一个算法交流群只要聊到“时间序列预测”前半段永远是 LSTM、Transformer、Prophet后半段永远有人追问“有没有现成代码”。但我这几年做预测类需求尤其是面向中小业务场景的量级不大的数据最常交作业的反而是线性回归LR。这不是情怀是真实结果对比出来的。大半年之前我手头有一个周度销量预测的项目历史数据只有 120 个采样点序列本身有明显的上升趋势但噪声也不小。同事一开始坚持上 LSTM理由是“时间序列预测不就该用深度学习吗”。我花了三天把数据整理成滑窗格式模型结构、训练循环都跑通了结果测试集误差还不如我用fitlm写的一个四阶滞后线性回归。这个反差让我意识到一件事线性回归在时间序列预测里的地位远远被低估了。LR 做时间序列预测的核心并不玄乎把原本“一串连续数值”的序列改造成一张普通回归表格。每一行的特征可以是时间戳、滞后值、周期项、外部变量标签就是下一个时刻的观测值。这等于把时间序列问题降维成最经典的监督学习问题。线性回归既有封闭解又能输出每个特征的系数和显著性训练时间可以忽略不计。我经常打一个比方如果你只需要判断“气温下降了是否会带动取暖负荷上升”线性回归能直接把这种关系量化成一条公式而深度学习给你的往往是一个说不清道不明的黑盒子。这篇文章就是把我实际在 MATLAB 里跑通的一套完整方案整理出来。网上很多帖子标题挂着“matlab代码”点进去发现正文底部写着“注暂无Matlab版”让人非常难受。下面这段实现不需要 Premium 工具箱也能跑即使你的 MATLAB 只有最基础的内核我也会给一个用反斜杠运算符完成回归的替代写法。内容适合这么几类人刚入门时间序列预测的学生、需要在 MATLAB 中快速搭建预测模型的工程师以及所有被“暂无Matlab版”坑过但依然想动手复现的同学。2. 特征构造才是 LR 预测的真正暗门很多人一开始只会拿时间索引t去拟合y得到一个y a*t b。这种写法不能说错但适用范围极窄。真正让 LR 在时间序列预测中站稳脚跟的不是回归算法本身而是你喂给它的特征矩阵怎么组织。2.1 从一维序列到二维表格假设原始序列是y(1), y(2), ..., y(n)。我们要预测y(t)最自然的三类特征是时间趋势特征t本身甚至t^2、sqrt(t)滞后特征y(t-1), y(t-2), ..., y(t-p)也就是把过去 p 个观测值作为自变量周期特征如果存在季节性可以加sin(2*pi*t/T)和cos(2*pi*t/T)。于是预测目标变成下面这个线性关系y(t) β0 β1*t β2*y(t-1) β3*y(t-2) ... β(p1)*y(t-p) ε这里 p 是滞后阶数。p 选多少直接决定模型是拟合不足还是过拟合。最朴素的判断方法是看自相关图用 MATLAB 的autocorr函数需要计量经济学工具箱观察拖尾情况如果没有工具箱就手动算corr(y(1:end-lag), y(lag1:end))看哪个 lag 对应的相关系数还明显非零。更懒但有效的方法是试几个 p对比验证集 RMSE。2.2 三种多步预测策略如果只预测一步LR 只需要一个模型。但实际需求往往要求未来 10 步、20 步这时有三种做法策略做法优点缺点单步预测每次只用历史真实值预测下一时刻误差最小实现最简单只能走一步需要等新数据递归多步预测把上一步的预测值当作历史值继续预测下一步只需一个模型代码量少误差会累积越往后越漂直接多步预测为 h1,2,...,H 分别训练 H 个模型各步误差互不传导训练 H 个模型样本利用率下降在这篇文章的实现里我采用递归多步预测因为它最贴合“给出一段历史一口气画出未来曲线”的使用场景。如果你非常在意长期预测稳定性可以在它的基础上增加定期重训机制后面我会详细说。2.3 特征构造的最大禁忌未来信息泄漏这是新手最容易犯的错也是最隐蔽的坑。构造训练样本时标签必须是y(k)特征只能是t以及在k之前已经发生的历史值。如果你为了补全缺失日期偷偷把整条序列的中心平滑值放进特征或者把“未来均值”拿来当代替值模型的训练误差会异常好看但上线后立刻崩盘。时间序列预测里信息只能从过去流向未来这比任何调参技巧都重要。还有一个细节切分训练集和测试集时绝对不能用随机抽样。时间序列的测试集必须是时间上靠后的那一段比如前 80% 训练、后 20% 测试。否则随机打乱会破坏序列的时间依赖关系评估出来的指标完全失真。3. MATLAB 完整实现一步能跑到出图的脚本下面这套代码我按照“即使你只有基础 MATLAB 也能改造”的标准写。主流程用fitlm统计和机器学习工具箱随后给一个不使用工具箱的替代计算段保证两条路都能跑。3.1 准备一组示例数据实际项目里你多半是读 CSV 文件里的单列数据。为了这篇文章能直接复现我先构造一组带趋势和轻周期性的示例序列。把下面这段换成y readmatrix(你的数据.csv)即可。clear; clc; rng(42); % 1. 生成示例序列实际使用时可替换为单击读取数据 n 300; % 总样本长度 timeAxis (1:n); trend 0.08 * timeAxis; % 线性趋势 season 6 * sin(2*pi*timeAxis / 20); % 周期波动 y trend season 1.5 * randn(n, 1); % 加噪声 figure(Position, [100 100 900 400]); plot(timeAxis, y, b-, LineWidth, 1.2); grid on; xlabel(时间); ylabel(观测值); title(原始时间序列示例数据);如果你自己的数据是一个行向量记得先转成列向量y y(:);。这是 MATLAB 里最容易让人卡住的小毛病之一很多函数对行列方向非常敏感。3.2 构造滞后特征矩阵核心是循环生成X和Y。X的每一行代表一个样本第一列是时间索引后面 p 列是最近的 p 个历史值。Y的每一行是对应时刻的观测值。% 2. 参数设置 p 4; % 滞后阶数建议根据实际数据调整 testLen 50; % 末尾作为测试集的样本数量 horizon 15; % 要预测的未来步数 % 3. 构造回归矩阵 m n - p; % 回归样本总数 X zeros(m, p 1); % 特征矩阵时间 滞后特征 Y zeros(m, 1); % 标签 for k p 1 : n row k - p; X(row, 1) timeAxis(k); % 时间趋势项 for lag 1 : p X(row, 1 lag) y(k - lag); % 滞后项 end Y(row) y(k); end当样本量极大时内层 for 循环会稍慢。如果数据量过了几万点可以把滞后列写成向量化赋值这样会快很多X(:, 1) timeAxis(p 1 : n); for lag 1 : p X(:, lag 1) y(p 1 - lag : n - lag); end Y y(p 1 : n);两种写法得到的矩阵完全一样第二种显然更紧凑。我第一次是从循环版开始调通的后来数据量上来了才改成向量化也建议大家一步步来。3.3 训练与测试下面的代码按照时间顺序切分训练集和测试集用fitlm拟合模型。% 4. 切分训练集与测试集 mTest testLen; trainRows 1 : (m - mTest); testRows (m - mTest 1) : m; mdl fitlm(X(trainRows, :), Y(trainRows)); % 训练集与测试集预测 predTrain predict(mdl, X(trainRows, :)); predTest predict(mdl, X(testRows, :)); % 输出模型关键信息 fprintf(训练集样本数%d\n, length(trainRows)); fprintf(回归系数); fprintf(%.4f , mdl.Coefficients.Estimate); fprintf(\n);这里fitlm会自动添加截距项所以不用手动加全 1 列。如果你传入的数据里已经包含常数列反而会导致秩亏问题输出 NaN 系数这一点要特别留意。模型训练完成后可以直接从mdl.Rsquared.Ordinary拿到 R²从mdl.Coefficients.tStat看每个特征的显著性。虽然时间序列数据的显著性检验不能像经典统计那样严格解读但系数的大小和正负方向仍然能告诉你滞后项和时间趋势哪些更重要。3.4 递归滚动预测未来递归策略的核心是维护一个长度等于 p 的“历史窗口”每预测一步就把新预测值塞进窗口末尾同时丢掉窗口最老的一个值。% 5. 递归预测未来 horizon 步 history y(end - p 1 : end); % 窗口里最近的 p 个真实观测值 future zeros(horizon, 1); for h 1 : horizon featRow [n h, history(end : -1 : 1)]; future(h) predict(mdl, featRow); history [history(2 : end); future(h)]; end注意history(end:-1:1)的作用把窗口从“最旧到最新”倒过来转置成“最新到最旧”的行向量再和[nh]拼成一行特征。比如预测第 n1 时刻特征就是[n1, y(n), y(n-1), ..., y(n-p1)]和训练时X(row, 2:end) [y(k-1), ..., y(k-p)]的列顺序严格一致。列顺序一旦搞反模型直接失去意义。3.5 画图 评估指标最后把训练历史、测试集预测、未来预测画在同一张图上并计算最常见的评估指标。% 6. 评估指标 rmseTrain sqrt(mean((Y(trainRows) - predTrain).^2)); maeTrain mean(abs(Y(trainRows) - predTrain)); rmseTest sqrt(mean((Y(testRows) - predTest).^2)); maeTest mean(abs(Y(testRows) - predTest)); fprintf(训练集 RMSE%.4f, MAE%.4f\n, rmseTrain, maeTrain); fprintf(测试集 RMSE%.4f, MAE%.4f\n, rmseTest, maeTest); % 7. 绘图 testTime p testRows; % 测试集样本对应的原始时间下标 futureTime (n 1 : n horizon); figure(Position, [100 100 1000 450]); plot(timeAxis, y, b-, LineWidth, 1.2); hold on; plot(testTime, predTest, ro, MarkerSize, 5); plot(futureTime, future, g--s, LineWidth, 1.5); hold off; grid on; legend({观测值, 测试集预测, 未来预测递归}, Location, best); xlabel(时间); ylabel(值); title(基于线性回归的时间序列递归预测);这套代码在一个包含 300 个样本的示例数据上测试集 RMSE 通常在 1.5 到 1.8 之间MAE 在 1.2 到 1.4 之间R² 在 0.9 以上。具体数值会随随机种子和噪声变化但量级可以给你一个参考。如果你的数据更平滑指标会更好看如果噪声大那就需要检查特征构造和滞后阶数是否合理。3.6 无工具箱版本反斜杠运算符也能做回归fitlm虽然方便但依赖统计和机器学习工具箱。如果你的 MATLAB 环境没有相关工具箱或者只有学校机房的基础版那就直接走最小二乘的线性代数解。线性回归的解析解是β (XX)⁻¹ XyMATLAB 里永远不要手动算逆矩阵直接使用反斜杠运算符Xdesign [ones(m, 1), X]; % 手动加截距列 beta Xdesign(trainRows, :) \ Y(trainRows); predTrainB Xdesign(trainRows, :) * beta; predTestB Xdesign(testRows, :) * beta; % 预测未来一步 futureB zeros(horizon, 1); historyB y(end - p 1 : end); for h 1 : horizon featRowB [1, n h, historyB(end : -1 : 1)]; futureB(h) featRowB * beta; historyB [historyB(2 : end); futureB(h)]; end反斜杠用的是经过数值优化过的 QR 分解或高斯消元速度快且稳定。即便是几千维的特征矩阵也是瞬间完成。唯一要注意的是Xdesign的列秩必须是满的否则解不稳定会出现很大的系数。这通常发生在某些滞后列和截距列完全共线时比如数据长期恒定时滞后特征之间相关性极高此时应该降低 p 或增加正则化项。4. 结果评估与残差检查不要只盯着预测曲线很多人在复现完代码后看预测曲线贴合得很好就宣布完工。但时间序列预测最容易出现的问题是在测试集上预测得好不代表模型抓住了真实的机制。你至少要做两件事计算量化指标查看残差结构。4.1 三个常用指标怎么解读RMSE均方根误差对大的偏差惩罚更重适合你受不了离谱预测值的场景MAE平均绝对误差更稳健不受个别离群点影响MAPE平均绝对百分比误差适合看相对误差但如果序列里存在 0 值或接近 0 的值就会爆炸慎用。举例来说在 CMA 的人工模拟数据上我的一次运行结果大致是指标训练集测试集RMSE1.471.62MAE1.161.28R²0.940.92你可以注意到测试集指标通常略差于训练集这是正常现象。如果测试集指标远差于训练集说明过拟合了最直接的体现就是你设的滞后阶数 p 太大。p20 的时候模型往往把所有噪声都背下来了泛化能力反而下降。4.2 残差才是模型的“心电图”预测值和真实值的差是残差。一个合格的线性回归残差应该在 0 附近随机分布不带明显的斜坡不带波浪形也不能出现“误差随预测值增大而增大”的漏斗状。如果残差带斜坡说明缺了趋势项 如果残差带波浪形说明缺了周期性特征补一组 sin/cos 即可 如果残差有尖刺说明遇到了异常值模型被极端点拉偏了。在 MATLAB 里最简单的做法是画残差散点图residualsTest Y(testRows) - predTest; figure; plot(predTest, residualsTest, b.); grid on; xlabel(预测值); ylabel(残差); title(测试集残差 vs 预测值);如果能看到明显的曲线形状别急着换算法先回头检查特征是否遗漏了重要变量。我处理过好几个“预测不准”的需求最后定位到的问题不是算法本身而是没有把节假日这类外部因素放进特征矩阵。线性回归的好处也在这里你可以像搭积木一样随时加一列特征比如温度、促销力度、星期几模型复杂度几乎不增加。4.3 递归预测的误差累积要坦然接受递归多步预测的 RMSE 通常会随着预测步数增加而变大。这不是代码写错了而是预测误差被累加了第二步的输入本身就有误差第三步的输入误差更大。所以你在实际汇报结果时最好把“近端预测”和“远端预测”分开讲或者直接画出不同步数的 RMSE 变化曲线给业务方看。让他们明白远未来预测天然存在更高的不确定性不是模型能单方面解决的。5. “暂无 Matlab 版”说明背后的问题以及怎么绕过在网络热词里能看到大量和 matlab、线性回归、时间序列预测相关的搜索词这本身说明需求很旺盛。但为什么很多资料挂出代码链接后又补充一句“注暂无Matlab版”据我观察原因大概有三层第一作者的原始代码可能是用 Python 或 R 写的他只是把原理和伪代码整理出来根本没有在 MATLAB 环境里运行过不敢说自己是 MATLAB 版第二作者使用了某些付费工具箱比如 Deep Learning Toolbox、Econometrics Toolbox不方便公开完整可执行脚本第三也是最常见的一层作者自己用的是 GNU Octave 验证的Octave 和 MATLAB 的语法高度兼容但并非 100% 等价。他怕用户在 MATLAB 原版里跑不过干脆标注“暂无Matlab版”来免责。遇到这种情况我一般不建议立刻放弃或去找人“代运行”。你自己动手把它翻译成 MATLAB 脚本并跑通收获远远大于下载一个现成文件。很多问题调试完你对整个算法的理解会上一个台阶。5.1 没有 MATLAB 完整工具箱的替代路线如果你没有统计和机器学习工具箱就用我上面 3.6 节的反斜杠版本。这是最保险的路任何 MATLAB 版本都支持。如果你连 MATLAB 都没有GNU Octave 是一个很现实的选择。上面整套代码中除了极个别函数如readmatrix需要调整Octave 用csvread或load替代其他部分基本可以原样运行。先用 Octave 把算法逻辑验证通了再放到正式环境的 MATLAB 里走一遍省时又高效。5.2 跨语言参考价值如果你在 GitHub 上找到的是 Python 代码也别慌。sklearn.linear_model.LinearRegression用起来和 MATLAB 的fitlm核心思想完全一致Python 中用fit(X, y)拟合predict(X)预测 MATLAB 中用fitlm(X, y)拟合predict(mdl, X)预测 两者都要求特征矩阵是二维的每行是一个样本每列是一个特征。不同之处在于Python 的LinearRegression默认不计算特征显著性需要额外加statsmodels才能拿到。MATLAB 的fitlm默认把显著性检验都算好了。这也是我用 MATLAB 做线性回归时更顺手的原因之一少写很多评估代码。需要提醒的是不要在网络上寻找所谓“特殊方式”获取 MATLAB。如果你不是商业用途优先看学校、学院的站点许可证如果暂时没有Octave 完全足够支撑你把算法思路走通等有正式环境再迁移。6. 我踩过几次坑之后积累的细节最后这部分是纯实战经验。上面 5 个章节里的代码和理论基本能保证你把模型跑通这一部分的价值在于让模型在真实业务环境里不至于翻车。6.1 训练集和测试集必须保持时间顺序我在文章里已经说过一次但值得再强调一遍时间序列预测中随机划分训练集和测试集是一个隐蔽的“致命错误”。如果训练集里混入了较晚的数据模型其实提前“见过”了未来。你看到的测试集误差会非常小上线后却一塌糊涂。坚持按下标顺序做前 80% / 后 20% 分割除非你用的是专门的时间序列交叉验证方法。6.2 缺失值和异常值不能直接丢进滞后窗口滞后特征要求时间间隔均匀。如果某个日期缺失你的 y(k) 和 y(k-1) 之间的间隔就变成了两天序列回归会错乱。处理方式要么插值补齐要么把缺失段整体剔除并忽略那几行样本。我个人更推荐先剔除异常尖峰再用前后均值平滑而不是直接让极端值参与拟合。线性回归对离群点非常敏感一个异常值就可能让系数产生明显变化。6.3 滚动重训比一次性拟合更抗漂移在线预测场景里数据分布会缓慢变化。训练一次模型然后用三年肯定越用越不准。比较实用的做法是固定窗口重训比如窗口长度设为 90 个样本每天来一个新观测值就把窗口向前滚动一天删除最老的那个点重新调用一次fitlm。这样既能留出测试段做效果监控又能让模型及时吸收近期变化。如果嫌每天都重训麻烦你也可以退一步每周重训一次或者设置一个监控指标当近一周 MAE 超过阈值时才触发重训。大多数情况下不会需要日级重训因为线性回归本身训练成本极低重训一次的耗时在很多数据集上连 0.1 秒都不到。6.4 滞后阶数和特征数量不是越多越好p 太大时模型参数量增加过拟合风险上升p 太小时模型可能学不到序列的自相关性。实际操作中我习惯先把 p 设成 4 到 6 跑一版基线再去尝试 12、24 等更大的值用测试集 RMSE 对比。MATLAB 里也可以手算一个简易 AIC 来辅助判断AIC m * log(SSE / m) 2 * (p 2)SSE 为残差平方和数值越小越好。不同 p 算同一指标取最小即可。6.5 最终部署时记得把“预测”和“决策”分开线性回归给的预测值永远是一个点估计它不包含不确定性的完整描述。业务决策时不要把预测值当作精确的未来事实。一个更稳妥的做法是以预测值为中心构造一个合理区间比如预测值 ± 1.96 倍的历史残差标准差。这个区间能够提醒决策者未来落在区间内的概率大致有多高。技术文档里通常不会写这些但做事情靠区间比靠单点踏实得多。最后分享一个我在实际项目里反复验证过的小习惯递归预测往前走 5 步以上的时候尽量把模型滚动起来。今天预测明天的值等明天真实数据到了就把真实值纳入窗口再预测后一天的。这样每一步预测都基于最新信息比一口气预测完未来 30 天强得多。线性回归的强项从来不在于“一步看得很远”而在于“每一步都足够稳”。把它的定位放对它就是你工具箱里最值得信任的武器之一。
