最近后台一直有人在问时间序列预测该选什么模型做课程设计或者工程项目时既想快速出结果又想把原理讲得明明白白。我给出的答案里反复出现一个名字——线性回归。别看深度学习的声势一年比一年大在时间序列预测这件事上线性回归LR依然是工业界和学术界最常见的baseline没有之一。尤其是用Matlab来做核心求解过程真的可以简洁到夸张一个反斜杠运算符就能把回归系数求出来。这篇内容就把“基于线性回归的时间序列预测”完整讲透覆盖特征构造、滞后期选择、多步预测策略、完整代码以及我实际调试中踩过的坑。网上这类教程铺天盖地都是Python版Matlab版反倒零散这里一次性补齐。1. 线性回归做时间序列预测为什么过了这么多年仍然值得用1.1 深度学习时代为什么还要把LR当回事现在聊时间序列预测大家开口就是LSTM、Transformer、Informer这些似乎不用深度模型就落伍了。但实际做项目的时候你会发现深度学习模型的调参成本和数据需求都远超预期在没有足够数据、算力和时间的情况下复杂模型很容易变成“精美的过拟合机器”。线性回归的价值在于它是一个足够聪明的参照系。第一可解释性强。回归系数直接告诉你过去第几天的数据对当前预测影响有多大是正相关还是负相关这在业务汇报、论文撰写、课程答辩里都非常好讲。第二训练成本几乎为零。Matlab里矩阵求逆加一个反斜杠运算符百万级以下的数据量都是秒出结果。第三它是检验特征工程是否有效的放大镜。如果你加了某个特征之后LR的误差显著下降说明这个特征确实携带了预测信息如果LR怎么调都不行那基本可以断定数据里没有线性可分的模式这时候再上复杂模型才有意义。所以在我的项目流程里拿到一条时间序列的第一件事永远是先跑LR把这个“及格线”画出来。后续无论换什么模型心里都有底。1.2 核心思路把“顺序问题”硬生生变成“表格问题”时间序列预测和普通回归最本质的区别在于普通回归的样本之间是独立的而时间序列的样本不是——它们是按时间顺序排列的彼此之间存在先后依赖。怎么让线性回归这种本来不关心顺序的模型来处理顺序数据答案是滑窗构造特征。假设你有一组按小时记录的电力负荷数据 y1, y2, y3, ..., yn你想预测 yt。线性回归不可能直接理解“t-1时刻在前、t时刻在后”这种顺序关系但它能理解数字之间的函数关系。所以我们要做的事情是把过去 p 个时刻的值作为特征把当前时刻的值作为目标构造出一张“表格”。举个例子假设 p 3那么训练数据长这样特征1yt-3特征2yt-2特征3yt-1目标yty1y2y3y4y2y3y4y5y3y4y5y6这样一构造问题就彻底变了原来是一串按时间排列的序列现在变成一个普通的多元线性回归问题y ≈ beta0 beta1 * yt-3 beta2 * yt-2 beta3 * yt-1。模型完全不关心时间顺序本身它只负责学出这组系数。这就是线性回归做时间序列预测的全部秘密把顺序问题转化为回归问题。后面所有的操作包括代码怎么写、滞后期怎么选、多步预测怎么做都是围绕这张“表格”展开的。1.3 这类模型到底擅长捕捉哪些模式搞清楚LR的适用范围非常重要否则你会被它的失败坑得很惨。线性回归本质上是学一条“直线”或“超平面”去拟合数据所以它对以下三类模式表现不错第一线性趋势。数据整体呈现稳定上升或下降比如月销售额逐步攀升、温室内温度缓慢升高这类趋势LR几乎白捡。第二短期惯性依赖。很多序列在相邻时刻高度相关比如今天的气温、水位、流量和昨天、前天往往强相关。只要这种相关是线性的滞后特征就能很好地捕捉。第三叠加了周期成分的均值回归比如正弦波加上噪声。LR拟合出来的其实是对周期成分的分段线性逼近能大致跟随波峰波谷相位上会有轻微延迟。反过来如果数据存在明显的非线性突变、指数级增长、状态切换比如某个系统突然从正常模式切换到故障模式线性回归就会力不从心。这不是模型写错了而是线性模型本身表达能力的边界。还有一个容易忽略的点LR在“插值”范围之外几乎没有外推能力。测试集的时间点如果超出训练集的时间范围预测本质上是沿着回归超平面外推误差会随着预测距离拉远而快速膨胀。这是多步预测误差累积的根本原因后面细说。2. 建模前的三个决定平稳性、滞后阶数、预测步数2.1 数据拿到手先做平稳性检查很多人拿到数据就直接塞进模型这是时间序列预测里最危险的习惯。线性回归假设特征和目标之间的关系是稳定的也就是说数据的统计性质均值、方差、自相关结构不能随时间发生剧烈变化。如果原始序列带长期趋势或者波动幅度不断增大回归系数会被“带偏”预测后期基本失效。拿Matlab操作来说最简单粗暴的办法是画图观察。plot(1:n, y)之后看均值是否在某个水平附近震荡、方差是否大致恒定。严谨一点用Econometrics Toolbox里的adftest做单位根检验[h, pValue] adftest(y);h 1表示序列平稳h 0表示存在单位根也就是非平稳。如果没有这个工具箱也可以用差分判断对序列做一阶差分 diff(y)如果差分后的序列看起来平稳了那么原始序列大概率是I(1)过程。如果发现数据带趋势最常用的两个对策先做一阶差分对差分序列建模。预测得到差分值后再累加回原始尺度。在特征矩阵里增加一列“时间序号”1, 2, 3, ...让线性回归自己学一个时间趋势项。个人经验带明显上升趋势的温度、流量数据加入时间序号效果立竿见影对于周期性很强的电力负荷数据差分往往更管用。两个都试一下选测试集误差小的方案。2.2 滞后期怎么定别拍脑袋用自相关函数说话滞后阶数 p 是整个方法里最关键的超参数没有之一。p 太小模型学不到足够的“记忆”p 太大特征维度膨胀还容易引入共线性系数变得极不稳定。判断 p 的主要工具是自相关函数ACF。对序列计算滞后 k 阶的自相关系数看它衰减到什么程度。Matlab里可以用autocorr直接绘图figure; autocorr(y, 30);如果0到2阶自相关系数都很高3阶开始明显跌到置信区间以内那 p 取2到3就够了。很多实际场景里p 取3已经能够覆盖绝大多数短周期依赖。这背后的直觉是线性模型对过去信息的利用是“直接叠加”的滞后2阶和滞后3阶之间往往高度相关所以多加几阶的边际收益很小反而容易把噪声放进模型。更系统的做法是网格搜索。在训练集上对 p 从1到10遍历训练模型并在验证集上算RMSE画一条折线选择RMSE开始平台化或开始上升的位置。这套流程在后面的完整代码里我会给出。2.3 预测策略单步、直接多步、滚动多步很多新手做多步预测时的做法是训练好模型之后把测试集的全部特征一次性喂进去得到所有预测值。这在学术评估上勉强说得通但真实场景里完全不成立。真实预测是站在当前时刻没有任何“未来的过去值”可用的。所以多步预测必须明确策略常见有三种单步预测。只预测下一时刻等真实值出来了再滚动。这是最简单、最可靠的方式误差最小但无法提前得到远期预判。直接多步预测。想预测未来h步就分别训练h个模型第1个模型用 yt 的特征预测 yt1第2个模型用 yt 的特征预测 yt2以此类推。每个模型各学各的目标互不干扰不会把预测误差像滚雪球一样传递下去。代价是训练量变成原来的 h 倍。滚动递归多步预测。先用模型预测 yt1然后把这个预测值当作新特征的一部分继续预测 yt2一步一步“自己喂自己”。实现简单但预测误差会叠加尤其是超过5步之后误差常常呈指数放大。实际项目里我用得最多的是直接多步预测和滚动多步预测的结合短期几步用滚动中远期用直接模型。代码部分我会把直接多步预测的实现放进去。2.4 划分训练集和测试集时间序列不能乱切这个坑我见人踩过无数次包括我自己早期也犯过。普通机器学习的训练集和测试集常常是随机划分但时间序列绝对不能这么做。原因很简单你把一条序列顺序打乱之后每个样本的“历史上下文”就被破坏了测试集里可能出现“未来数据帮助预测过去”的荒谬情况最终得到的评估指标虚高得离谱。正确做法是严格按时间顺序切分。比如前80%的时间段作为训练集后20%作为测试集train_len floor(0.8 * size(X, 1)); X_train X(1:train_len, :); y_train y(1:train_len); X_test X(train_len 1:end, :); y_test y(train_len 1:end);另外还要提醒一点如果先对整个序列做归一化再划分训练测试统计量均值、标准差会用到测试集的信息这属于特征泄露。严格讲应该用训练集的均值和标准差去归一化训练集再用同一组参数归一化测试集。小规模实验很多人嫌麻烦直接整段归一化问题不大但你要知道自己做了取巧。3. Matlab完整实现从特征矩阵、回归求解到多步预测3.1 手写滞后特征构造函数Matlab不像Python的pandas那样方便地做shift操作但写一个构造滞后特征矩阵的循环并不复杂。这里给出我自己常用的一版function [X, y] lagged_features(data, p) % 构造滞后特征矩阵 % data: 列向量原始时间序列 % p: 滞后阶数 % X: n-p 行 p 列第 k 列为滞后 k 阶的数据 % y: n-p 行 1 列对应的预测目标 data data(:); n length(data); if n p error(数据长度必须大于滞后阶数); end m n - p; X zeros(m, p); for k 1:p X(:, k) data(p - k 1 : n - k); end y data(p 1 : n); end这个函数的核心逻辑是第 k 列是滞后 k 阶的值即用 yt-k 作为第 k 个特征。循环里 data(p - k 1 : n - k) 的意思是把序列往前平移 k 个位置截取出来和目标 y data(p1:end) 对齐。如果想把时间趋势项也加进去调用之后补一列即可[X, y] lagged_features(y_raw, p); trend (1:size(X, 1)); X [X, trend];3.2 核心一行反斜杠运算符求解线性回归特征矩阵构造好之后线性回归的参数估计就是最小二乘问题beta argmin ||X*beta - y||²。Matlab里面最干净、最地道的写法是beta X_train \ y_train;反斜杠运算符会自动选择合适的算法矩阵规模小时走QR分解大规模稀疏时有专门处理数值稳定性比手写 inv(X*X) * X*y 要好得多。我的原则是永远不要手写正规方程求逆直接用反斜杠。如果你想看到更丰富的统计输出比如系数的置信区间、残差方差、R²等可以用regress函数[Xn_train, ~, mu, sigma] zscore(X_train); Xn_train [ones(size(Xn_train, 1), 1), Xn_train]; [b, bint, r, rint, stats] regress(y_train, Xn_train);不过需要注意regress要求特征矩阵是满秩的如果滞后特征之间严重共线它会直接报矩阵奇异。而反斜杠运算符在这种情况下通常还能给出一个最小范数解虽然结果要打个问号。所以我的建议是先追求简洁用反斜杠出了诡异结果再回去检查特征是否存在共线性。3.3 多步预测的实现以直接多步为例如果只做单步预测预测代码就一行y_hat X_test * beta;如果要做未来 h 步的直接多步预测每个预测步长需要单独训练一个模型。假设滞后期为 p预测步长从1到h那么第 s 步模型的输入特征仍然是 yt-p1 到 yt目标变成 ytsfunction [y_hat_multi, beta_list] direct_multistep(X, y, h, train_len) % 直接多步预测为每一步训练一个独立模型 % X: 滞后特征矩阵不含趋势项 % y: 目标向量 % h: 最大预测步数 % train_len: 训练集长度 [m, p] size(X); beta_list cell(h, 1); y_hat_multi zeros(h, 1); for s 1:h % 第 s 步的目标是 y_{ts} ys y(1 s - 1 : m s - 1); % 这里按函数外部截断方式构造需保证索引不越界 % 通常更安全的做法是在函数外预先扩展目标矩阵 X_s X(1 : m - s 1, :); ys_cut y(s 1 : m 1); % 取前 train_len - s 1 行训练 beta_s X_s(1:train_len - s 1, :) \ ys_cut(1:train_len - s 1); beta_list{s} beta_s; % 用训练集最后一个窗口做递推输入 last_x X(end, :); y_hat_multi(s) last_x * beta_s; end end这段代码我写出来主要是想说明“每个步长单独建模”的结构。实际工程里我会把它整理得更干净预先构造一个目标矩阵 Y_steps第 s 列为 y_{ts}然后对每一列分别训练。滚动多步预测的代码更短核心就是把预测值回填进特征function y_hat rolling_predict(beta, initial_x, h) % initial_x: 1×p 行向量最新的 p 个观测值 % h: 预测步数 % beta: 回归系数向量长度 p 或 p1 p length(initial_x); x_input initial_x; y_hat zeros(h, 1); for s 1:h if length(beta) p 1 y_hat(s) x_input * beta(1:p) beta(p 1); else y_hat(s) x_input * beta; end x_input [x_input(2:end), y_hat(s)]; end end实测下来滚动预测在1到3步内效果尚可5步之后误差增长非常明显。所以不要指望这个极简模型能做超长周期预测它更适合做短期预警。3.4 评估指标四个数看清模型水平预测完必须量化评估。我通常同时看RMSE、MAE、MAPE和R²四者各有侧重rmse sqrt(mean((y_test - y_hat).^2)); mae mean(abs(y_test - y_hat)); mape mean(abs((y_test - y_hat) ./ y_test)) * 100; ss_res sum((y_test - y_hat).^2); ss_tot sum((y_test - mean(y_test)).^2); r2 1 - ss_res / ss_tot;RMSE对大误差敏感如果你想惩罚那些离谱的预测主要盯它MAE反映平均绝对偏差量纲直观MAPE是百分比适合和业务方沟通但y_test里有接近0的值时会爆炸这时候要谨慎使用R²反映模型相比“直接用均值预测”提升了多少接近1说明模型抓住了大部分波动接近0或为负说明预测还不如无脑用历史均值。注意反归一化的问题。如果训练前做了zscore归一化预测结果是归一化尺度必须还原到原始尺度再算指标y_hat_raw y_hat * sigma mu;3.5 完整可运行脚本把上面几块串起来一个完整的单步预测脚本如下%% 1. 准备数据 clear; clc; close all; rng(42); t (1:800); y_raw 0.02 * t 10 * sin(2 * pi * t / 80) 3 * randn(800, 1); %% 2. 数据预处理 mu mean(y_raw(1:600)); sigma std(y_raw(1:600)); y (y_raw - mu) / sigma; %% 3. 构造滞后特征 p 10; [X, y_target] lagged_features(y, p); %% 4. 划分训练集和测试集 train_len floor(0.8 * size(X, 1)); X_train X(1:train_len, :); y_train y_target(1:train_len); X_test X(train_len 1:end, :); y_test y_target(train_len 1:end); %% 5. 训练 beta X_train \ y_train; %% 6. 预测与反归一化 y_hat_norm X_test * beta; y_hat y_hat_norm * sigma mu; y_true y_test * sigma mu; %% 7. 指标 rmse sqrt(mean((y_true - y_hat).^2)); mae mean(abs(y_true - y_hat)); fprintf(RMSE %.4f, MAE %.4f\n, rmse, mae); %% 8. 画图 figure; plot(y_true, b-, LineWidth, 1.2); hold on; plot(y_hat, r--, LineWidth, 1.2); legend(真实值, 预测值); xlabel(测试集样本); ylabel(数值); title(线性回归时间序列预测结果); grid on;这里我故意用了带趋势和周期成分的模拟数据一个p 10的单步LR就能跟上曲线读者可以复制下来直接跑然后换自己的数据试。4. 诊断与可视化预测图、残差图、系数解读4.1 预测图不能只看“贴不贴近”很多人画完预测图看到两条曲线缠在一起就以为大功告成这是典型的自我安慰。正确的做法是先把测试集真实值画出来再把预测值画出来同时标注训练集和测试集的切分位置。这样能一眼看出误差主要集中在哪些区域。我习惯在测试集前段补画一段训练集末尾的预测值用来观察模型是否出现“切换点崩溃”。很多模型在训练集结尾处表现良好一进入测试集立刻偏离说明它过拟合了历史噪声而不是学到了规律。figure; plot(1:train_len, y_target(1:train_len), k-); hold on; plot(train_len1:length(y_target), y_target(train_len1:end), b-); plot(train_len1:length(y_target), y_hat_norm, r--); xline(train_len 0.5, --, 切分点);4.2 残差分析真正的诊断核心残差 真实值 - 预测值。一个合格的线性回归时间序列模型残差应该表现为白噪声均值接近0没有明显的自相关性没有周期性。如果残差序列里还能看出明显的“波浪”说明你的模型漏掉了周期性成分需要把对应周期的特征补进去。比如每24小时一个周期的负荷数据可以在特征里加入“时刻”的哑变量或者加入周期项 sin(2pitime/24)、cos(2pitime/24)。如果残差在某个时间段突然整体偏移说明序列在那个时间点发生了结构变化模型已经不适应新状态这时候建议重新训练或引入状态变量。如果残差方差越来越大像喇叭口一样展开说明数据存在异方差性线性模型无能为力至少要换成加权回归或者对目标做对数变换。4.3 回归系数的业务解读线性回归最香的地方就在这里你能把系数拉出来解释。假设滞后特征是 yt-1、yt-2、yt-3最后学出的系数是 beta [0.62, 0.18, 0.05]含义是当前值主要受上一时刻影响影响力按时间距离快速衰减符合直觉。如果beta[0.7, -0.3, 0.1]说明yt-2对yt有反向作用常见于周期性序列比如波峰过去两个时刻之后开始回落。系数的绝对值大小还能帮助筛选特征。如果一个滞后阶数的系数始终在0附近且换了训练集后符号不稳定这个特征基本可以丢掉了。这比单纯看相关系数更贴近预测任务本身。5. 必踩的坑与排查实录5.1 预测曲线变成一条水平线这是最经典的现象测试集预测值几乎恒等于一个常数真实值还在上下波动。原因通常是滞后特征的系数都很小模型学到的截距项明显大于其他项导致所有预测都向均值回归。解决办法分三步排查检查数据是否被错误归一化。如果训练集和测试集的归一化参数不一致预测值会被压扁。检查滞后期是否太小。p 1时模型只能捕捉“昨天的值直接决定今天”一旦测试集波动形态改变模型就只能输出平均水平。检查训练集是否包含足够多的“极端值”样本。如果训练集里目标值分布非常集中回归模型会把所有预测都推到均值附近。5.2 系数异常大甚至互相抵消滞后特征之间天然强相关yt-1和yt-2的相关系数经常超过0.9这会导致多重共线性。症状是系数数值巨大、符号一正一负、看上去毫无解释力但预测结果又似乎还行。最直接的解决方法是改用岭回归Matlab里一行调用lambda 0.1; beta_ridge ridge(y_train, X_train, lambda); y_hat_ridge [ones(size(X_test, 1), 1), X_test] * beta_ridge;ridge会牺牲一点训练集拟合度来换取系数的稳定性实际效果常常比普通最小二乘更稳。lambda可以用crossval网格搜索我个人经验是从0.01开始每10倍往上加看测试集RMSE的“浴缸曲线”找到最优值。5.3 序列有趋势却忘了处理如果你的序列整体在上升而你没加趋势项也没做差分线性回归学出来的其实是一条“平均趋势直线”对局部波动的预测基本无效。我的建议是先差分再建模把趋势彻底去掉。差分之后的预测值是“变化量”最后累加回原始值。这套组合在金融、气象、流量数据上都比直接对原始值建模稳定。还有一种混合方案对差分序列和趋势项同时建模相当于让模型自己决定是用“惯性”还是用“趋势”。5.4 训练集完美测试集崩盘这种情况十有八九是过拟合。滞后阶数取得太大模型把训练集里的噪声当作规律一遇到测试集的新噪声就乱套。一个重要的实操原则滞后期宁少勿多。我在调试中反复验证很多数据p取2到6就足够超过10通常只是把问题复杂化。先把p设成3跑一遍基线再逐步增加观察测试集误差是否真的在下降。若误差不降反升果断回退。5.5 多步预测误差像滚雪球滚动多步预测中的误差累积是数学上注定的结果不是你的代码写错了。每预测一步误差会作为“输入噪声”进入下一步特征形成自我放大的恶性循环。应对策略预测步数控制在3步以内把任务设计成短期预警而不是长期预报。改用直接多步预测每个步长单独建模切断误差链。每步都做一次特征更新真实观测一旦到手立即用真实值替换预测值重新构造特征窗口再预测下一步。这种方法叫做“校准的滚动预测”在工程上非常实用。以我个人做负荷预测的经验来说LR模型在1到3步内能给出非常可靠的参考值超过6步就开始明显偏差。如果你需要的是更长周期的预测我的建议是放弃LR转用带时间注意力机制的模型或者在LR基础上引入外生变量天气、节假日、事件标记把预测问题从“纯时间序列外推”变成“带辅助信息的回归问题”效果往往比换个复杂模型更立竿见影。最后分享一个我自己的习惯任何时间序列项目无论最终决定用多高级的模型第一版永远是线性回归。它像一个诚实的老伙计不会给你惊喜但也不会给你幻觉。把这条基准线跑明白你对自己数据的理解会提升一个档次。往后加特性、换模型也都有一个清晰的比较对象。这套Matlab代码从头到尾都是可以逐行跑的你可以把它当模板换成自己的数据跑通之后再按上面的思路做诊断和优化。