使用混合方法预测多维时间序列
导言
在时间序列分析和预测中，经典问题仍然是将观测信号分离为确定性和随机分量。 如果确定性部分允许精确的参数描述，那么其余部分通常包含一个相关结构，该结构不会减少到简单的噪声，但也不会用线性方法描述。 将经典分解与余数的神经网络建模相结合的混合方案使得可以利用这两种方法的优势。
在本文中，我们考虑了一种预测多维时间序列的混合方法，该方法包括三个阶段:
- 最小二乘线性趋势去除(OLS);
- 用Nelder-Mead方法选择频率优化的两个主要正弦分量;
- 由循环神经网络lstm预测余数。
使用RMSE和MAE度量评估预测的质量，并且还与没有神经网络块的模型进行比较。
让我们考虑在这项工作中使用的一些定义。:
*趋势是反映其变化的大方向的时间序列的长期，单调的组成部分：增长，下降或停滞。
*季节性是一个系列的周期性成分，它以大致相等的时间间隔以恒定的幅度和相位重复。
*Lstm（长期短期记忆）是一种循环神经网络，能够以较长的时间间隔存储有关序列先前状态的信息。 与经典的循环神经网络不同，LSTM使用一种机制来控制隐藏状态下信息的写入，读取和擦除，这使其能够有效地对时间序列中的长期依赖关系进行建模。
*剩余能量是信号与其模型之间差异的平方和。: . 这是一个标准的OLS目标函数;最小化它意味着找到它最好描述信号的模型参数。 本文将残余能量作为正弦频率优化的标准。
*RMS（均方根，均方根值）–值的均方根: . 对于该系列的其余部分，RMS以与信号相同的单位显示差异的典型幅度。 训练区域中剩余部分的RMS表征该系列的"未解释"部分，然后将其传输到LSTM输入。
*RMSE（均方根误差，均方根误差）-预测质量的度量: . 它以信号单位测量，比MAE更强烈地惩罚大失误。
*MAE（平均绝对误差）-公制: . 它显示了信号单位的平均"典型"预测误差，并且对异常值的敏感度低于RMSE。
*Mnc（最小二乘法）是一种用于估计线性模型的参数的方法，该参数最小化平方残差的总和。 它用于估计正弦曲线的线性趋势和振幅。
*Nelder-Mead方法是一种不需要导数的无条件优化算法;它用于通过残余能量的准则细化一对正弦频率。
加入图书馆
连接必要的库:
Flux -神经网络训练及应用;
Optim -正弦频率的优化;
Statistics-计算余数统计进行归一化;
LinearAlgebra -计算线性趋势和正弦振幅的MNC系统的解决方案;
Random -学习的可重复性。
for pkg in ["Flux", "Optim", "Statistics", "LinearAlgebra", "Random"]
try
eval(Meta.parse("using $pkg"))
catch
Pkg.add(pkg)
eval(Meta.parse("using $pkg"))
end
end
Random.seed!(42)
数据生成
我们将生成三个时间序列 , , 5000点长。
每一行都包含一个总体趋势 和季节性成分 具有相移和恒定 . 为 和 整体趋势的幅度正在缩放，季节性分量的相位正在转移。
让我们绘制所有三行。
t = 0:1:4999
B = 8 .* sin.(0.006 .* t) .+ 0.006 .* t
c = 30
X = B .+ 5 .* sin.(0.016 .* t) .+ c
Y = 0.8 .* B .+ 5 .* sin.(0.016 .* t .+ π/3) .+ c
Z = 1.2 .* B .+ 5 .* sin.(0.016 .* t .+ 2π/3) .+ c
plot(t, [X Y Z], title="时间序列", lw=1.8, label=["X" "Y" "Z"])
分成训练和测试样本
让我们将数据拆分为训练样本（前4000个点）和测试样本（最后1000个点），并形成一个公共数组。 данные_все (5000×3)和时间向量 tt. 让我们分别突出时间轴。 t_обуч 为训练区。
N_обуч = 4000
N_тест = 1000
N = N_обуч + N_тест
данные_все = hcat(X, Y, Z)
данные_обуч = данные_все[1:N_обуч, :]
данные_тест = данные_все[N_обуч+1:end, :]
tt = collect(Float64, 0:N-1)
t_обуч = tt[1:N_обуч]
消除线性趋势
让我们使用每行的最小二乘(OLS)训练样本估计线性趋势，沿着其整个长度绘制趋势并从源数据中减去它。 结果，我们得到了减损的数据。 детренд_все 和 детренд_обуч.
A = hcat(ones(N_обуч), t_обуч)
коэф = A \ данные_обуч
тренд_все = hcat(ones(N), tt) * коэф
детренд_все = данные_все .- тренд_все
детренд_обуч = детренд_все[1:N_обуч, :]
剩余能量函数
定义函数 остаточная_энергия(частоты, сигнал, t)，它基于给定的频率对，建立常数，余弦和正弦的基础，选择OLS的系数并返回残差的平方和。 该函数将作为频率优化的目标函数。
function остаточная_энергия(частоты, сигнал, t)
f1, f2 = частоты
(f1 <= 0 || f2 <= 0) && return 1e30
M = hcat(ones(length(t)),
cos.(2π*f1*t), sin.(2π*f1*t),
cos.(2π*f2*t), sin.(2π*f2*t))
β = M \ сигнал
r = сигнал .- M*β
return sum(abs2, r)
end
两个正弦搜索功能
定义函数 синусоиды(сигнал, t)，它在信号中搜索两个主导频率：首先，它通过网格暴力搜索找到第一个频率，然后，通过减去其贡献，第二个，然后使用NelderMead方法澄清两个频率。 返回找到的对(f1,f2)。
function синусоиды(сигнал, t)
сетка_частот = 1e-5 : 1e-5 : 1.2e-2
лучшая_f1 = 0.0; лучший_r = Inf
for f in сетка_частот
M = hcat(ones(length(t)), cos.(2π*f*t), sin.(2π*f*t))
β = M \ сигнал
r = sum(abs2, сигнал .- M*β)
if r < лучший_r
лучший_r = r; лучшая_f1 = f
end
end
M1 = hcat(ones(length(t)), cos.(2π*лучшая_f1*t), sin.(2π*лучшая_f1*t))
сигнал2 = сигнал .- M1 * (M1 \ сигнал)
лучшая_f2 = 0.0; лучший_r2 = Inf
for f in сетка_частот
abs(f - лучшая_f1) < 5e-4 && continue
M = hcat(ones(length(t)), cos.(2π*f*t), sin.(2π*f*t))
β = M \ сигнал2
r = sum(abs2, сигнал2 .- M*β)
if r < лучший_r2
лучший_r2 = r; лучшая_f2 = f
end
end
целевая = частоты -> остаточная_энергия(частоты, сигнал, t)
res = Optim.optimize(целевая, [лучшая_f1, лучшая_f2],
NelderMead(),
Optim.Options(x_tol=1e-12, f_tol=1e-15,
iterations=N))
f1, f2 = Optim.minimizer(res)
return f1, f2
end
按级数选择正弦曲线
应用函数 синусоиды 对于三个确定的系列中的每一个，将找到的频率与真实值进行比较 和 . 我们将在训练区域中执行振幅的最终MNC选择，并沿着行的整个长度构造正弦分量。 我们将为每行输出RMS。
_, n_сигн = size(данные_все)
синус_все = zeros(N, n_сигн)
частоты_итог = zeros(n_сигн, 2)
for k in 1:n_сигн
сигнал = детренд_обуч[:, k]
f1, f2 = синусоиды(сигнал, t_обуч)
частоты_итог[k, :] = [f1, f2]
f1_истина = 0.006 / (2π)
f2_истина = 0.016 / (2π)
println("Ряд $k: f1 = $(round(f1, digits=10)) (истина $(round(f1_истина, digits=10)), Δ = $(abs(f1-f1_истина)))")
println(" f2 = $(round(f2, digits=10)) (истина $(round(f2_истина, digits=10)), Δ = $(abs(f2-f2_истина)))")
M_обуч = hcat(ones(N_обуч),
cos.(2π*f1*t_обуч), sin.(2π*f1*t_обуч),
cos.(2π*f2*t_обуч), sin.(2π*f2*t_обуч))
β = M_обуч \ сигнал
M_все = hcat(ones(N),
cos.(2π*f1*tt), sin.(2π*f1*tt),
cos.(2π*f2*tt), sin.(2π*f2*tt))
синус_все[:, k] = M_все * β
end
остаток_все = детренд_все .- синус_все
остаток_обуч = остаток_все[1:N_обуч, :]
println("\nRMS:")
for k in 1:n_сигн
println(" Ряд $k: $(round(sqrt(mean(abs2, остаток_обуч[:,k])), sigdigits=3))")
end
平衡正常化
计算平均值 和标准偏差 剩下的训练样本，我们来替换那些太小的 通过1.0，我们得到归一化余数。 остаток_норм 每行的均值和单位方差为零。
μ = mean(остаток_обуч, dims=1)
σ = std(остаток_обуч, dims=1)
σ[σ .< 1e-10] .= 1.0
остаток_норм = (остаток_обуч .- μ) ./ σ
为LSTM创建窗口
让我们为LSTM创建训练序列：窗口长度为l=100，对于三个通道中的每一个。 获取数组 X_посл 表格(3,100,3900)-入口窗户及 Y_посл 表格(3,3900)是目标值。
L = 100
N = N_обуч - L
K = 3
X_посл = zeros(Float32, K, L, N)
Y_посл = zeros(Float32, K, N)
for i in 1:N
X_посл[:, :, i] = Float32.(остаток_норм[i:i+L-1, :]')
Y_посл[:, i] = Float32.(остаток_норм[i+L, :])
end
LSTM培训
让我们定义一个lstm模型（输入3→隐藏状态32→输出3），在输出端有一个完全连接的层。 让我们设置Adam优化器并以每个64个示例的迷你批次中的10个epoch训练模型，最大限度地减少窗口最后一步的预测与真实值之间的MSE。 我们将输出每个epoch的损失函数值。
модель = Chain(LSTM(K => 32), Dense(32 => K))
function потери(модель, x, y)
Flux.reset!(модель)
ŷ = модель(x)
return Flux.mse(ŷ[:, end, :], y)
end
сост_опт = Flux.setup(Flux.Adam(1e-3), модель)
размер_батча = 64
эпох = 10
n_батчей = div(N, размер_батча)
for эпоха in 1:эпох
перест = randperm(N)
потери_эпохи = 0.0
for b in 1:n_батчей
idx = перест[(b-1)*размер_батча+1 : b*размер_батча]
xb = X_посл[:, :, idx]
yb = Y_посл[:, idx]
Flux.reset!(модель)
l, grads = Flux.withgradient(модель) do m
Flux.mse(m(xb)[:, end, :], yb)
end
Flux.update!(сост_опт, модель, grads[1])
потери_эпохи += l
end
println("эпоха $эпоха потери = $(round(потери_эпохи/n_батчей, sigdigits=4))")
end
余额预测
定义函数 прогноз_остатка，它通过递归预测余数 N_тест 前进的步骤：在每一步，长度的最后一个窗口被采取 L，被馈送到LSTM，接收到的预测被添加到缓冲器并用于形成下一个窗口。 然后，我们使用以下方法将预测反规范化回原始比例 和 . 结果，我们得到 остаток_полн 尺寸为5000×3。
function прогноз_остатка(модель, остаток_норм, μ, σ, N_обуч, L, горизонт)
K = 3
всего = N_обуч + горизонт
буфер = zeros(Float32, всего, K)
буфер[1:N_обуч, :] = Float32.(остаток_норм)
Flux.reset!(модель)
for i in 1:горизонт
idx = N_обуч + i
окно = буфер[idx-L:idx-1, :]'
x_вх = reshape(окно, K, L, 1)
ŷ = модель(x_вх)
буфер[idx, :] = vec(ŷ[:, end, :])
end
return Float64.(буфер) .* Float64.(σ) .+ Float64.(μ)
end
остаток_полн = прогноз_остатка(модель, остаток_норм, μ, σ, N_обуч, L, N_тест)
最终预测
让我们把三个部分加起来：线性趋势 тренд_все，正弦分量синус_все 和余额的预测 остаток_полн. 我们会得到完整的预测 y_полн 让我们突出显示测试部分 y_прогноз. 我们将保存真实值以供比较 y_истина.
y_полн = тренд_все .+ синус_все .+ остаток_полн
y_прогноз = y_полн[N_обуч+1:end, :]
y_истина = данные_тест
质量指标
让我们计算测试区域中每三行的RMSE和MAE指标。 此外，我们将计算预测而不考虑lstm进行比较。 它将显示LSTM提高预测质量的程度。
rmse(a, p) = sqrt(mean((a .- p).^2))
mae(a, p) = mean(abs.(a .- p))
имена = ["X", "Y", "Z"]
y_прогноз_безLSTM = (тренд_все .+ синус_все .+ vcat(остаток_обуч, zeros(N_тест,3)))[N_обуч+1:end, :]
println("\nКачество прогноза (t = $N_обуч..$(N_обуч+N_тест-1)):\n")
println(rpad("排", 5), rpad("RMSE（带LSTM）", 16), rpad("RMSE / размах", 16),
rpad("MAE (с LSTM)", 16), rpad("RMSE без LSTM", 16), "Улучшение по RMSE")
for k in 1:3
истина_k = y_истина[:, k]
прогноз_k = y_прогноз[:, k]
безLSTM_k = y_прогноз_безLSTM[:, k]
r_с = rmse(истина_k, прогноз_k)
r_без = rmse(истина_k, безLSTM_k)
m_с = mae(истина_k, прогноз_k)
размах = maximum(истина_k) - minimum(истина_k)
доля = 100 * r_с / размах
улучш = r_без / r_с
println(rpad(имена[k], 5),
rpad(round(r_с, digits=3), 16),
rpad("≈ $(round(доля, digits=1)) %", 16),
rpad(round(m_с, digits=3), 16),
rpad(round(r_без, digits=3), 16),
"в $(round(улучш, digits=1)) раза")
end
预测的可视化
让我们建立三个子图（, , ），在蓝色的每个训练区域上，测试上的真实值以橙色点缀，预测为红色。
p1 = plot(tt[1:N_обуч], X[1:N_обуч], label="", color=:blue, linewidth=2,
title="X", legend=false)
plot!(p1, tt[N_обуч+1:end], X[N_обуч+1:end], label="",
color=:orange, linestyle=:dash, linewidth=3)
plot!(p1, tt[N_обуч+1:end], y_прогноз[:,1], label="",
color=:red, linewidth=3)
vline!(p1, [N_обуч+0.5], color=:black, linestyle=:dash, label="", linewidth=1)
p2 = plot(tt[1:N_обуч], Y[1:N_обуч], label="", color=:blue, linewidth=2,
title="Y", legend=false)
plot!(p2, tt[N_обуч+1:end], Y[N_обуч+1:end], label="",
color=:orange, linestyle=:dash, linewidth=3)
plot!(p2, tt[N_обуч+1:end], y_прогноз[:,2], label="",
color=:red, linewidth=3)
vline!(p2, [N_обуч+0.5], color=:black, linestyle=:dash, label="", linewidth=1)
p3 = plot(tt[1:N_обуч], Z[1:N_обуч], label="", color=:blue, linewidth=2,
title="Z", legend=false)
plot!(p3, tt[N_обуч+1:end], Z[N_обуч+1:end], label="",
color=:orange, linestyle=:dash, linewidth=3)
plot!(p3, tt[N_обуч+1:end], y_прогноз[:,3], label="",
color=:red, linewidth=3)
vline!(p3, [N_обуч+0.5], color=:black, linestyle=:dash, label="", linewidth=1)
display(plot(p1, p2, p3, layout=(3,1), size=(950, 950), legend=false))
结论
本文构建并研究了一种多维时间序列预测的混合方法，将OLS分解与余数的LSTM建模相结合。 所提出的方法表现出高质量：测试区域的预测误差为信号范围的≈2％（RMSE≈0.38–0.52，范围为≈25单位），与没有神经网络单元的模型相比，精度提高了3.5-4.5倍。
混合方案在广泛的应用中有需求：预测电网的能耗和负荷，需求和销售，财务和产品范围，气象和气候指标，技术过程和遥测参数，交通和物流。 经典分解为小样本提供了可解释性和稳定性，而LSTM考虑了其余部分中的非线性和长期依赖性。