线性回归实现时间序列预测：Matlab完整代码与实战

做时间序列预测，很多人一上来就选LSTM、ARIMA这类“看起来高级”的模型。但真到了实际工程和科研论文里，线性回归(LR)这个被低估的老模型反而经常被我当成第一根探针。它训练快、可解释、几乎不需要调参，还能当复杂模型的baseline。最近看到一个仓库，标题是“基于线性回归(LR)算法的时间序列预测”，但仓库注释写得很直白：暂无Matlab版。这种“Python先行、Matlab缺失”的现象在时间序列代码里特别常见，评论区总有人催Matlab实现。我干脆自己动手补了一个完整版本，把滑窗特征构造、fitlm训练、滚动预测、指标评估和可视化全流程跑通。这篇文章不会讲太多花哨理论，重点放在你能直接复现的Matlab代码和实操经验上。

1. 线性回归凭什么能做时间序列预测

1.1 把时间序列转换成回归问题，核心在“滑窗”

时间序列本身是一串有序数据，y1, y2, ..., yt。预测yt+1，LR模型并不会直接处理时间维度，而是把过去几期的观测值当作自变量，把目标期的观测值当作因变量。这就是经典的“自回归”思路：用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:i+lag−1)，对应的Y(i)是data(i+lag)。

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:i+lag-1)'; Y(i) = data(i+lag); 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(splitIdx+1:end); % 构造滑窗特征 lag = 6; [trainX, trainY] = createLagFeatures(trainRaw, lag); [testX, testY] = createLagFeatures(testRaw, lag);

这里trainRaw和testRaw都是列向量。用sin和cos组合生成数据只是为了便于复现，你完全可以替换成真实业务数据，比如CSV里的销量或负载序列。替换时只需要把数据加载成列向量即可。

我在实际操作中发现，很多人会直接拿原始数据全部样本构造特征，再手动切分，这样会导致训练集和测试集在窗口上有交叉。假设原始序列有1000条，lag=6，如果先构造全部特征再按行切分，测试集第一个样本会用到训练集最后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);

在我用上述正弦叠噪声数据跑出来的结果里，lag=6时测试集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等预测yt+1。真实业务里经常要预测未来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差异可能会很大，写报告时必须明确标注，否则自己都会被误导。

我一般会给出一张简单表格：

可以看到滚动多步的误差明显更高，这是递归式预测的误差累积效应。以后你看到有人只贴一个很高的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补全版本能帮你少走一些弯路。