Прогнозирование многомерных временных рядов гибридным методом
Введение
В задачах анализа и прогнозирования временных рядов классической проблемой остаётся разделение наблюдаемого сигнала на детерминированную и стохастическую компоненты. Если детерминированная часть допускает точное параметрическое описание, то остаток, как правило, содержит коррелированную структуру, которая не сводится к простому шуму, но и не описывается линейными методами. Гибридные схемы, сочетающие классическую декомпозицию с нейросетевым моделированием остатка, позволяют использовать сильные стороны обоих методов.
В настоящей работе рассматривается гибридный метод прогнозирования многомерных временных рядов, включающий три этапа:
- удаление линейного тренда методом наименьших квадратов (МНК);
- подбор двух доминирующих синусоидальных компонент с оптимизацией частот методом Нелдера-Мида;
- прогнозирование остатка рекуррентной нейронной сетью LSTM.
Качество прогноза оценивается метриками RMSE и MAE, а также сопоставляется с моделью без нейросетевого блока.
Рассмотрим некоторые определения, применяемые в данной работе:
- Тренд – это долговременная, монотонная составляющая временного ряда, отражающая общее направление его изменения: рост, спад или стагнацию.
- Сезонность – это периодическая составляющая ряда, повторяющаяся через примерно равные интервалы времени с постоянной амплитудой и фазой.
- LSTM (Long Short-Term Memory, долгая краткосрочная память) – разновидность рекуррентной нейронной сети, способная сохранять информацию о предыдущих состояниях последовательности на длительных интервалах. В отличие от классической рекуррентной нейронная сети, LSTM использует механизм, управляющий записью, чтением и стиранием информации в скрытом состоянии, что позволяет ей эффективно моделировать долгосрочные зависимости во временных рядах.
- Остаточная энергия – сумма квадратов невязок между сигналом и его моделью: . Это стандартная целевая функция МНК; её минимизация означает поиск параметров модели, при которых она наилучшим образом описывает сигнал. В данной работе остаточная энергия используется как критерий при оптимизации частот синусоид.
- RMS (Root Mean Square, среднеквадратичное значение) – корень из среднего квадрата величины: . Для остатка ряда RMS показывает типичную амплитуду невязки в тех же единицах, что и сигнал. RMS остатка на обучающем участке характеризует "необъяснённую" часть ряда, которая далее передаётся на вход LSTM.
- RMSE (Root Mean Square Error, корень из средней квадратичной ошибки) – метрика качества прогноза: . Измеряется в единицах сигнала и сильнее штрафует крупные промахи, чем MAE.
- MAE (Mean Absolute Error, средняя абсолютная ошибка) – метрика: . Показывает среднюю «типичную» ошибку прогноза в единицах сигнала и менее чувствительна к выбросам, чем RMSE.
- МНК (метод наименьших квадратов) – способ оценки параметров линейной модели, минимизирующий сумму квадратов невязок. Применяется для оценки линейного тренда и амплитуд синусоид.
- Метод Нелдера-Мида – алгоритм безусловной оптимизации, не требующий производных; используется для уточнения пары частот синусоид по критерию остаточной энергии.
Присоединение библиотек
Подключим необходимые библиотеки:
Flux – обучение и применение нейронной сети;
Optim – оптимизация частот синусоид;
Statistics– вычисление статистик остатка для нормировки;
LinearAlgebra – решение МНК-систем для вычисления линейного тренда и амплитуд синусоид;
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_обуч]
Удаление линейного тренда
Оценим линейный тренд по обучающей выборке методом наименьших квадратов (МНК) для каждого ряда, построим тренд на всей длине и вычтем его из исходных данных. В результате получим детрендированные данные детренд_все и детренд_обуч.
A = hcat(ones(N_обуч), t_обуч)
коэф = A \ данные_обуч
тренд_все = hcat(ones(N), tt) * коэф
детренд_все = данные_все .- тренд_все
детренд_обуч = детренд_все[1:N_обуч, :]
Функция остаточной энергии
Определим функцию остаточная_энергия(частоты, сигнал, t), которая по заданной паре частот строит базис из константы, косинусов и синусов, подбирает коэффициенты МНК и возвращает сумму квадратов остатков. Эта функция будет использоваться как целевая при оптимизации частот.
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
Подбор синусоид по рядам
Применим функцию синусоиды к каждому из трёх детрендированных рядов, сравним найденные частоты с истинными значениями и , выполним финальный МНК-подбор амплитуд на обучающем участке и построим синусоидальные компоненты на всей длине ряда. Выведем 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 и обучим модель 10 эпох мини-батчами по 64 примера, минимизируя MSE между предсказанием на последнем шаге окна и истинным значением. Выведем значение функции потерь на каждой эпохе.
модель = 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))
Заключение
В данной работе построен и исследован гибридный метод прогнозирования многомерных временных рядов, сочетающий МНК-декомпозицию с LSTM-моделированием остатка. Предложенный метод показал высокое качество: ошибка прогноза на тестовом участке составила ≈ 2 % от размаха сигнала (RMSE ≈ 0.38–0.52 при размахе ≈ 25 ед.), а по сравнению с моделью без нейросетевого блока точность улучшилась в 3.5–4.5 раза.
Гибридные схемы востребованы в широком круге прикладных задач: прогнозирование энергопотребления и нагрузки электросетей, спроса и продаж, финансовых и товарных рядов, метеорологических и климатических показателей, параметров технологических процессов и телеметрии, трафика и логистики. Классическая декомпозиция обеспечивает интерпретируемость и устойчивость к малым выборкам, а LSTM – учёт нелинейных и долгосрочных зависимостей в остатке.