做时间序列预测很多人一上来就选LSTM、ARIMA这类“看起来高级”的模型。但真到了实际工程和科研论文里线性回归(LR)这个被低估的老模型反而经常被我当成第一根探针。它训练快、可解释、几乎不需要调参还能当复杂模型的baseline。最近看到一个仓库标题是“基于线性回归(LR)算法的时间序列预测”但仓库注释写得很直白暂无Matlab版。这种“Python先行、Matlab缺失”的现象在时间序列代码里特别常见评论区总有人催Matlab实现。我干脆自己动手补了一个完整版本把滑窗特征构造、fitlm训练、滚动预测、指标评估和可视化全流程跑通。这篇文章不会讲太多花哨理论重点放在你能直接复现的Matlab代码和实操经验上。1. 线性回归凭什么能做时间序列预测1.1 把时间序列转换成回归问题核心在“滑窗”时间序列本身是一串有序数据y1, y2, ..., yt。预测yt1LR模型并不会直接处理时间维度而是把过去几期的观测值当作自变量把目标期的观测值当作因变量。这就是经典的“自回归”思路用yt−1, yt−2, ..., yt−p预测yt。数学表达就是yt α β1·yt−1 β2·yt−2 ... βp·yt−p εt翻译成机器学习的话就是构造一个形状为(n−p)行、p列的二维特征矩阵X每一行对应一个历史窗口最后一列是对应的标签Y。这个“滑窗”操作是整个时间序列预测能否用普通回归模型跑通的关键。很多人第一反应是“LR不是做横截面数据吗”这正是没理解时间序列和横截面数据在特征构造上的区别。一旦把时间序列转成滑窗样本后续训练、评估和普通回归没有任何本质区别。1.2 滑动窗口长度、样本数量与信息利用滑动窗口长度p决定了模型能“看到”多少历史信息。p太小模型抓不到周期性p太大不但计算量增加还会引入过多冗余。因为这个原因实际工程里通常会把自相关函数(ACF)和偏自相关函数(PACF)画出来看看阶数大致在什么地方截断。以日频销量数据为例如果ACF在滞后7时仍有较明显峰值就说明有周效应窗口至少应该覆盖到7甚至14。另一个现实问题是滑窗会把原始时间序列压缩n个数据点经过p阶滑窗后只剩n−p个样本数据量越小越要谨慎。我做这个Matlab版本时故意没有把窗口参数写死而是设计成函数参数方便你针对不同数据调整。后面第2节会给出具体的构造代码。这块弄明白了LR时间序列预测就成功了一半。1.3 先想清楚LR适合预测什么、不适合预测什么线性回归擅长捕捉平稳、线性的短期趋势尤其在样本量不大、需要解释模型行为的场景下非常稳。它也有明显短板对剧烈波动、强非平稳、复杂非线性关系比较吃力。用LR预测电负荷、股票收盘价这类数据效果往往不如LSTM但如果你的目标是流量预测、温度趋势、能耗数据这类规律性较强的序列LR完全够用而且训练时间只有深度模型的零头。关键是要摆正心态LR不是万能模型而是基线模型。先用LR跑通完整流程建立评估指标再决定是否升级到更复杂的算法这是我在实战中一直坚持的做法。2. Matlab实现数据准备与特征构造2.1 数据怎么选、怎么切分才不“作弊”时间序列预测最容易犯的错误之一是切分不当。不能像普通分类数据那样随机打乱因为时间序列一旦打乱就破坏了时间顺序训练出的模型等于“偷看”了未来。正确做法是按时间顺序把数据切成训练集和测试集常见比例是8:2或7:3。切分时还有一个容易被忽略的点标准化参数只能从训练集计算。如果你先把整条序列标准化再切分测试集的均值方差已经被模型看到了这就是数据泄露。在Matlab里要用训练集trainset的mean和std去标准化测试集而不是全量zscore。这一点我在第5节还会详细展开。2.2 构造滑窗特征矩阵的Matlab代码我先定义一个通用的滑窗函数。这个函数输入一维序列和滞后阶数lag输出特征矩阵X和标签向量Y。每一行的X(i,:)是data(i:ilag−1)对应的Y(i)是data(ilag)。function [X, Y] createLagFeatures(data, lag) % 将一维时间序列转换为滑窗特征矩阵 % data: 列向量 % lag: 窗口长度使用滞后1~lag n length(data); if n lag error(数据长度必须大于滞后阶数); end samples n - lag; X zeros(samples, lag); Y zeros(samples, 1); for i 1:samples X(i, :) data(i:ilag-1); Y(i) data(ilag); end end这个函数是整篇文章的地基。用循环实现是为了让逻辑更直观数据量上万时性能也足够。如果数据量很大可以考虑用reshape或hankel矩阵优化但不要为了炫技牺牲可读性。注意这里的X每一行内部是按时间顺序排列的千万不要在构造特征时把一行内的顺序打乱那是另一种“时间倒流”。2.3 训练集和测试集的特征构造以我平时用的模拟数据为例用正弦加噪声生成一条时间序列既能覆盖周期成分又带随机扰动非常适合作算法验证。生成并切分数据的代码如下rng(1); time (0:0.05:30); signal sin(0.6*time) 0.3*sin(2.1*time) 0.15*randn(size(time)); % 按时间顺序切分 trainRatio 0.8; splitIdx floor(length(signal) * trainRatio); trainRaw signal(1:splitIdx); testRaw signal(splitIdx1:end); % 构造滑窗特征 lag 6; [trainX, trainY] createLagFeatures(trainRaw, lag); [testX, testY] createLagFeatures(testRaw, lag);这里trainRaw和testRaw都是列向量。用sin和cos组合生成数据只是为了便于复现你完全可以替换成真实业务数据比如CSV里的销量或负载序列。替换时只需要把数据加载成列向量即可。我在实际操作中发现很多人会直接拿原始数据全部样本构造特征再手动切分这样会导致训练集和测试集在窗口上有交叉。假设原始序列有1000条lag6如果先构造全部特征再按行切分测试集第一个样本会用到训练集最后5条真实数据相当于测试时提前看到了部分未来信息。所以正确做法是先切分原始数据再分别构造特征。3. 训练、预测与可视化完整代码3.1 用fitlm训练线性回归模型Matlab里做线性回归有很多选择regress函数偏统计fitlm封装更完整适合快速建模。fitlm会自动处理截距项并给出系数、p值、R²等统计量对后续分析方便。训练和单步预测的代码如下mdl fitlm(trainX, trainY); % 训练集和测试集单步预测 trainPred predict(mdl, trainX); testPred predict(mdl, testX);这里fitlm的默认设置是包含截距项所以不需要手动给X加一列全1。你需要留意的是predict函数要求输入的特征矩阵列数与训练时一致。如果我们用trainX训练那么testX必须同样是lag列否则Matlab会直接报错。训练完成后mdl对象里包含了核心信息。我个人最喜欢看的是mdl.Coefficients中的p值它能快速告诉你哪些滞后项显著。如果某个滞后项的p值很大说明这个历史观测对当前值没有明显解释力可以尝试缩小窗口。这比盲目调参数有意义得多。3.2 评估指标与结果输出的Matlab实现光有预测曲线还不够还得量化误差。时间序列回归最常用的三个指标是RMSE、MAE和R²。RMSE对大误差比较敏感MAE更关注平均偏差R²则看模型解释了多大比例的波动。计算代码如下rmse sqrt(mean((testY - testPred).^2)); mae mean(abs(testY - testPred)); ssRes sum((testY - testPred).^2); ssTot sum((testY - mean(testY)).^2); r2 1 - ssRes / ssTot; fprintf(RMSE %.4f\n, rmse); fprintf(MAE %.4f\n, mae); fprintf(R2 %.4f\n, r2);在我用上述正弦叠噪声数据跑出来的结果里lag6时测试集RMSE大约在0.2附近R²在0.75到0.85之间。具体数值会随随机种子、训练比例和窗口长度变化所以不要太在意单个数字重点看不同配置下的相对变化。3.3 可视化预测曲线与真实曲线对比可视化最重要的不是“画得漂亮”而是能不能看出偏差模式。我会同时画真实数据、训练集预测、测试集预测三条线另外再画一张残差图。残差图能直接暴露系统性的预测偏差这是数字指标很难体现的。figure; plot(trainRaw, k); hold on; plot(trainPred, b); plot(testY, r, LineWidth, 1.2); plot(length(trainRaw) (1:length(testPred)), testPred, m--, LineWidth, 1.2); legend(原始序列, 训练预测, 测试真实, 测试预测, Location, best); title(LR时间序列预测训练与测试对比); xlabel(时间序号); ylabel(数值); grid on;画完图之后我建议先看测试段的预测曲线是否整体滞后。线性回归做单步预测时经常出现“滞后效应”也就是预测值比真实值慢半拍尤其在拐点附近。这种滞后是因为模型把过去信息当作主要输入对突然的转折天然反应迟钝。看到这种图不用慌这是LR的固有属性。3.4 多步滚动预测怎么实现上面代码只做了单步预测也就是用真实的yt, yt−1等预测yt1。真实业务里经常要预测未来5步、30步这时候需要把预测值当作输入继续往后推。实现思路是先用最后一个真实窗口预测第一步然后把得到的预测值加入窗口丢掉窗口最前面的一个值保持窗口长度不变继续预测下一步。这称为递归滚动预测。代码可以先封装成一个函数function predSeq recursivePredict(mdl, recentWindow, steps) % 递归多步预测 % mdl: 线性回归模型对象 % recentWindow: 最近lag个真实值行向量 % steps: 需要预测的步数 predSeq zeros(steps, 1); window recentWindow; for k 1:steps p predict(mdl, window); predSeq(k) p; window [window(2:end), p]; end end注意输入recentWindow必须是一个1×lag的行向量顺序和训练时的特征顺序一致。调用时把测试集最后一个窗口传入即可。滚动预测最大的问题是误差会随着步数累加第一步的小误差会被喂给下一步所以多步预测的评估要和单步分开看。如果你要做多步评估应该把“第k步预测误差”单独统计而不是混成一个RMSE。4. 评估与解读别被R²骗了4.1 评估指标口径的统一时间序列预测的评估口径非常容易踩坑。常见做法有两种一是单步预测统一采用真实历史值作为特征这衡量的是模型的一步拟合能力二是滚动多步预测每一步都用自己的预测值当特征这衡量的才是真实业务中的“预测”能力。两种口径算出的RMSE差异可能会很大写报告时必须明确标注否则自己都会被误导。我一般会给出一张简单表格指标单步预测口径多步滚动口径RMSE0.1830.295MAE0.1410.232R²0.820.56可以看到滚动多步的误差明显更高这是递归式预测的误差累积效应。以后你看到有人只贴一个很高的R²先想想他测的是哪种口径。4.2 残差图是更好的诊断工具数字指标会掩盖问题。预测值只可能偏大或偏小如果RMSE没问题但残差分布有明显规律说明模型漏掉了某些结构。画残差图的方法很简单testResidual testY - testPred; figure; plot(testResidual, o); hold on; plot([1 length(testResidual)], [0 0], k--); title(测试集残差); xlabel(测试样本序号); ylabel(残差值); grid on;正常情况下残差应该围绕0上下波动不存在明显趋势。如果你看到残差在一段时间内几乎全为正接下来全为负那说明模型漏掉了趋势或者季节性。这时候提高窗口长度、增加差分预处理或者加入时间索引特征往往比换模型更有效。4.3 训练集和测试集的表现差距要重视如果训练集R²很高测试集R²断崖式下降模型就是过拟合了。线性回归的过拟合通常由特征过多或共线性引起。lag设得太大、历史数据太少都会导致系数被噪声带偏。处理方式可以是减少窗口长度也可以直接上岭回归后面第6节我会给出方向。反过来如果训练集和测试集误差都很大问题通常不在模型本身而在数据预处理。比如序列里有明显的趋势或突变却直接拿来训练或者测试集分布和训练集不一样。遇到这种情况先回头看看数据而不是急着换LSTM。5. 实战中容易踩的坑与排查技巧5.1 数据泄露比模型选错更致命我遇到过不止一次代码跑出来的R²高达0.98结果发现是数据泄露。时间序列里最典型的数据泄露有三种第一种是前面说的先标准化再切分让测试集的统计信息进入训练第二种是先构造全部样本再进行随机打乱破坏了时间顺序第三种是滑窗重叠导致训练集和测试集之间有重合区间。第三种特别隐蔽因为你把窗口重叠的数据删掉之后训练样本会减少不少但这是必须付出的代价。排查方法很直接把训练样本和测试样本的序号画出来看看有没有重叠。或者干脆在代码里严格遵循“先切分原始序列再分别构造特征”的顺序从源头杜绝。5.2 趋势与季节性线性回归最怕什么线性回归本质是在找输入和输出之间的线性映射它没有内置机制去处理时间结构。如果你的序列有上升趋势比如网站流量每年增长20%LR预测未来时往往会把当前趋势外推导致预测持续偏高或偏低。解决思路有两种一是做一阶差分把原始序列转成增量序列再对增量建模二是直接把时间索引作为一个特征加入模型让模型学到一个全局趋势项。对于有季节性比如周周期、年周期的序列可以把“第几个星期”“第几个月”做成虚拟变量或数值特征。不过要注意一旦引入时间索引你预测未来时就得知道未来的时间序号这在实际预测中是已知信息所以不会有泄露问题。5.3 共线性与特征尺度两个容易被忽略的隐性风险自回归窗口里相邻的滞后项往往是高度相关的这种多重共线性会让系数估计不稳定。今天跑出的β1可能是0.6换一段数据就变成0.1但预测效果看起来还行。这种不稳定会影响模型的可解释性也可能让样本外预测波动很大。用ridge回归或者lasso做正则化是应对共线性最实用的方法。特征尺度方面虽然LR本身对单变量缩放不敏感但如果你的特征里同时有“滞后销量”和“星期几”这种量纲完全不同的变量梯度下降类算法会受影响。Matlab的fitlm基于最小二乘解对尺度相对稳健但为了统一习惯我会在特征进入模型前做标准化。再次强调标准化参数必须在训练集上计算测试集沿用训练集的均值和标准差。5.4 滑窗长度到底怎么选我见过有人把lag设成8有人设成30还有人直接拍脑袋。最靠谱的办法是看自相关图。在Matlab里可以用autocorr函数画出ACF图观察滞后阶数对应的相关性在什么地方显著衰减。以下面这条曲线为例ACF在滞后6到10之间逐渐落到置信区间内那lag取6到10都是合理的。figure; autocorr(signal, 30); title(ACF);另一个实用经验是用不同lag值跑一个快速循环比较测试集RMSE。不要只看训练集指标因为lag增大通常会降低训练误差却可能带来过拟合。我习惯把lag从1扫到20记录测试RMSE画成一条U形曲线选择最低点或相邻几个值。6. 从LR出发还能怎么扩展6.1 从普通LR升级到正则化回归当数据维度变高、特征相关性变强时普通最小二乘的系数会“发散”这时岭回归和LASSO是比普通LR更稳的选择。Matlab里可以直接用lasso函数或fitrlinear配合正则化参数。岭回归适合保留全部特征的场景LASSO则会把不重要的滞后项系数压成0天然做了特征选择。我不建议一上来就上正则化。先跑一个普通LR作为baseline再比较加上L2惩罚后的效果。如果两者差异很小说明模型本身已经够用了不需要增加复杂度。6.2 多变量特征从“单序列自回归”到“外生变量回归”LR的优势之一就是能自由融合不同类型的外生变量。比如预测明天的用电量除了历史用电量还可以加入气温、星期类型、节假日标记。在Matlab里只需把外生变量按时间对齐后拼接到特征矩阵右侧然后用fitlm训练即可。需要注意的是外生变量也必须严格对齐并提前想好“预测时能否拿到这些变量未来的值”。比如气温有天气预报值可以加入当天实时促销活动数据预测明天就未必能拿到这类输入要谨慎。6.3 对比LSTM什么时候该换模型很多人跑完LR之后会问“我是不是该上LSTM”。我的经验是如果LR的测试R²已经到0.8以上换LSTM的收益通常有限还要面临调参、训练时间长、可解释性差的问题。但如果数据量大、非线性特征明显、预测步数长LSTM确实可能做得更好。更合理的做法是把LR当成基准先记录LR的RMSE再看LSTM是否显著超过它。如果没有超过那LSTM的额外复杂度就没有意义。最后再分享一点个人经验。我在实际项目中做时间序列预测从来不会一上来就写深度模型代码。第一步永远是先用线性回归(LR)把数据流跑通有没有数据泄露、指标怎么算、特征怎么构造这些基础问题在简单模型上暴露得最快。把它吃透了再去碰LSTM、Transformer你才知道哪些收益来自模型哪些只是运气。希望这份Matlab补全版本能帮你少走一些弯路。