Неувядающая классика или «чёрный ящик»: кто кого в битве за прогноз. Глава вторая. Завершение

В прошлых частях второй главы мы с вами определили оптимальные значения параметров (p, d, q) статистических моделей семейства АРПСС по одноимённой методологии, и выполнили две подходящие модели, включая сезонную модель. В этой, завершающей, части мы будем использовать временной ряд со значениями температуры в качестве сигнала и применим к нему дискретное преобразование Фурье, чтобы выявить сезонные компоненты и разложить их на составляющие гармоники, сумму которых будем использовать в качестве экзогенной переменной. В конце выясним, сможет ли это улучшить точность предсказаний моделей.

1. Описание набора климатических данных и его предварительная подготовка.

2. Статистический анализ данных, включая оценку сезонной компоненты, для разработки моделей из семейства ARIMA с целью прогнозирования температуры: начало, продолжение и завершение.

3. Разработка моделей глубокого обучения Keras с использованием слоя LSTM для аналогичной задачи.

4. Сравнение и оценка эффективности результатов прогнозирования.

5. Рассмотрение подхода прогнозирования смещённой температуры с помощью модели глубокого обучения Keras с использованием слоёв LSTM.

Перед тем как перейти к основному материалу, давайте определим, что такое сигнал, и рассмотрим несколько примеров различных сигналов, включая их спектральные характеристики и гармоническую структуру — совокупность синусоид, из которых они формируются. Затем перейдём к анализу временного ряда температурных данных.

Сигнал, как правило, содержит информацию, которая может быть представлена в виде функции времени или частотного спектра. Обычно сигналы классифицируют по двум признакам: по времени – являются ли они непрерывными, то есть определёнными для любого момента времени, либо дискретными, когда они определены только для некоторых моментов; а также по повторяемости – являются ли они периодическими, то есть повторяющими свои значения через определённый интервал времени, или непериодическими (апериодическими). В качестве сигнала мы будем рассматривать температурный временной ряд, который сочетает в себе как повторяющиеся значения, свидетельствующие о наличие сезонных колебаний, так и непериодические компоненты, отражающие нерегулярные изменения. Временной ряд состоит из последовательных наблюдений, сделанных в определённые моменты времени, что делает его дискретным по своей природе.

Основными компонентами, на которые может быть разложен сигнал, являются: амплитуда, частота и фаза.

Амплитуда – это максимальное значение смещения или изменения переменной величины от среднего значения сигнала при его колебательном или волновом движении. Как правило безразмерная величина, которая показывает, насколько сильно проявляется сигнал.

Частота – это количество повторяющихся колебаний или циклов сигнала, совершаемых за единицу времени. Измеряется в герцах.

Фаза – это сдвиг сигнала относительно начала отсчёта времени или другого сигнала, указывающий в какой момент времени начинается одно полное колебание или цикл. Измеряется в градусах или радианах.

Примером периодического сигнала является синусоидальный сигнал. В общем виде синусоидальный сигнал можно представить с помощью следующей формулы:

xt = A * np.sin(2 * np.pi * f * time_vec + phi)

где А – амплитуда сигнала, ƒ – частота повторов (циклов) в герцах, time_vec – массив с моментами времени, а phi – фаза сигнала в радианах (именно в радианах, поскольку используется компонента “2π“).

С помощью следующего кода сгенерируем простой синусоидальный сигнал по заданным параметрам:

import numpy as np

from matplotlib import pyplot
from scipy.fft import fft, fftfreq, ifft
from scipy.signal import periodogram
from pandas import DataFrame, read_parquet, Series
from sklearn.metrics import mean_absolute_percentage_error
from statsmodels.nonparametric.smoothers_lowess import lowess
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import Holt

pyplot.rcParams["figure.autolayout"]=True
pyplot.rcParams["figure.figsize"]=(14, 4)

np.random.seed(1)

n = 1000
time_vec = np.arange(0, 1, 1/n)
A = 1
f = 10
phi = 0

xt = A * np.sin(2 * np.pi * f * time_vec + phi)

pyplot.figure()
pyplot.plot(time_vec, xt, marker='.')
pyplot.axvline(x=time_vec[100], color='orange', linestyle='dashed')
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
text_1 = "Простой периодический синусоидальный сигнал с периодом 0.1 сек., "
text_2 = "амплитудой 1, частотой 10 Гц и фазой 0 радиан"
pyplot.title(text_1 + text_2)
pyplot.grid(True)
pyplot.show()

Обратите внимание, как в данном исследовании будет определяется временная ось. Сначала с помощью параметра n задаётся общее количество точек сигнала (=длина сигнала), а затем с помощью функции np.arange(0, 1, 1/n) создаётся массив моментов времени time_vec с равномерно распределёнными между 0 и 1 секундой шагами 1/n. В первом рассматриваемом примере сигнала массив time_vec содержит моменты времени, равномерно распределённые между 0 и 1 секундой с шагом 0.001. Все рассматриваемые далее сигналы будут находиться в интервале от 0 до 1 секунды (их длительность будет считаться равной одной секунде). Добавление «дискретные» к моментам времени, как по мне, было бы избыточным, поскольку при цифровой обработке сигналов уже подразумевается дискретизированная структура сигнала.

Итак, выполнив код, мы получили простой периодический синусоидальный сигнал, состоящий из 1000 точек. Видно, что это периодический сигнал: имеется 10 полных повтора (или цикла); длительность каждого цикла, называемая периодом, составляет 0.1 секунду. Частота 10 Гц означает, что в одну секунду в сигнале имеется 10 таких периодов. Амплитуда сигнала равна 1. А его фаза равна 0, что означает что сигнал начинается точно в начале координат.

Теперь поговорим о частотном спектре. Благодаря трудам французского математика Жана-Батиста Жозефа Фурье стало возможным представить любой периодический сигнал в виде комбинации гармонических колебаний. Фурье в своей работе «Аналитическая теория тепла» (1822) утверждал, что любой периодический сигнал может быть разложен в ряд синусоидальных функций с различными частотами, амплитудами и фазами. Совокупность амплитуд и частот гармоник в разложении Фурье называют спектром сигнала.

Рассмотрим разложение периодических функций (сигналов) в ряд Фурье с помощью дискретного преобразования Фурье (далее – ДПФ). ДПФ – это математический метод, который используется для анализа частотного состава сигналов. Он преобразует сигнал из временной области в частотную, представляя его набором комплексных чисел. Каждое комплексное число соответствует одной частоте и содержит информацию о её амплитуде и фазе. Обычно количество комплексных чисел равно количеству точек сигнала, что указывает на совпадение числа частот с числом точек сигнала. Исследователи часто фокусируются на частотах, которые связаны с синусоидальными компонентами исходного сигнала. Для этого надо правильно интерпретировать результаты ДПФ. Далее мы подробно рассмотрим, с помощью каких действий это можно выполнить.

При обработке временных рядов ДПФ полезен тем, что позволяет выявлять периодические компоненты и шумы, которые затем можно отфильтровать. Вычисления ДПФ основаны на алгоритме под названием «быстрое преобразование Фурье» (БПФ).

Python-библиотека SciPy предоставляет ряд функций для работы с сигналами, включая функции для выполнения ДПФ (scipy.fft). В данном исследовании будут рассмотрены три из них — fft(), fftfreq() и ifft(). Рассмотрим каждую из них по отдельности.

• Функция scipy.fft.fft() [ссылка] применяется для вычисления ДПФ сигнала, используя алгоритм БПФ (в англоязычной терминологии – Fast Fourier Transform, FFT).

  Основными аргументами функции, на которые хотелось обратить внимание, являются:

  – x: входной массив данных, который требуется преобразовать.

  – n: количество точек в выходном преобразовании.

  – norm: параметр, определяющий нормализацию результата. По умолчанию None (без нормализации).

• Функция scipy.fft.fftfreq() [ссылка] используется для вычисления частот, соответствующих результатам ДПФ.

  Одними из основных аргументов, которые принимает данная функция, являются:

  – n: количество точек в преобразовании Фурье (или размер окна, как в описании по документации).

  – d: временной шаг. Или расстояние между образцами (обратное частоте дискретизации), как в описании по документации. По умолчанию d считается равным 1.

  Понятие «частота дискретизации» будет затронуто при анализе последнего из приводимых в данной главе примеров сигнала: его рассмотрение (а равно и применение) в простых примерах опущено, чтобы не усложнять материал.

• Функция scipy.fft.ifft() [ссылка] используется для вычисления обратного ДПФ (в англоязычной терминологии – Inverse Fast Fourier Transform, IFFT). Эта функция позволяет преобразовывать данные из частотной области обратно во временную область.

  Обязательным аргументом функции является:

  – x: массив комплексных чисел, представляющих частотные компоненты, которые необходимо преобразовать обратно во временную область. Это может быть массив любого размерности. О комплексных числах будет рассказано далее.

  Необязательный аргумент:

  – n: количество точек для обратного преобразования.

Отдельно отмечу, что функция scipy.fft.fft() вычисляет ДПФ, но сама по себе не учитывает частоту выборки данных. Поэтому для интерпретации результатов ДПФ необходимо использовать функцию scipy.fft.fftfreq(), чтобы получить соответствующие частоты для каждого элемента выходного спектра.

Как видно, у всех трёх функций есть важный аргумент – параметр n, который указывает, сколько точек будет использоваться в выполняемом действии, что позволяет фиксировать размер выходного массива и разрешение частотного спектра. Применение этого параметра позволяет упростить обработку данных и эффективнее использовать память. Но поскольку рассматриваемые примеры сигналов и временной ряд с температурными данными статичны и не прожорливы в плане использования памяти, а также в целях облегчения объяснения выполняемых действий, данный параметр у вышеописанных функций будет применяться со значением по умолчанию. Однако для отображения фактических значений амплитуд и частот на графиках, а также при фильтрации частот этот параметр сыграет ключевую роль, как будет показано далее и в своё время.

Про комплексные числа. Результатом ДПФ является набор комплексных чисел, которые могут быть использованы для расчёта амплитуд и фаз синусоидальных компонентов, составляющих исходный сигнал.

(В книге «Цифровая обработка сигналов на языке Python» (М.:ДМК Пресс, 2017. – 160 с.: ил.) хорошо объясняются комплексные числа применительно к обработке сигналов, а также такое понятие, как биение сигнала.

При этом раз речь зашла о литературе, то рекомендую книгу Питера Блумфилда Fourier Analysis of Time Series, 2nd edition. Изложенный в ней материал последовательный – от простых концепций до сложных. И хотя автор представляет книгу в качестве «вводного текста по методам Фурье, не обременённого обилием математических, вероятностных или статистических подробностей», однако в ней всё же присутствует достаточное количество математического материала, чтобы обеспечить глубокое понимание анализа Фурье применительно к временным рядам. Именно в этой книге я находил ответы на вопросы, которые появлялись при изучении спектрального анализа.)

В python комплексное число имеет тип complex и состоит из двух чисел с плавающей точкой, представляющих его действительную и мнимую части. Для представления мнимой части используется j (число, оканчивающееся на j, считается мнимым). При этом у комплексных чисел есть атрибуты real и imag для получения действительной или мнимой части соответственно. Действительная и мнимая части комплексного числа в результате ДПФ не являются амплитудой и фазой напрямую, но они могут быть использованы для расчёта этих величин: для отдельного получения модуля (амплитуды) используется функция np.abs(), для отдельного получения угла – np.angle(). Демонстрационный код представлен ниже.

>>> x = complex(1, 2)
>>> print(x)
# (1+2j)
>>> y = complex(3, 4)
>>> print(y)
# (3+4j)
>>> z = x + y
>>> print(z)
# (4+6j)
>>> print(z.real)
# 4.0
>>> print(z.imag)
# 6.0
>>> print(np.abs(z))
# 7.211102550927979
>>> print(np.angle(z))
# 0.982793723247329

Итак, приступим к выполнению дискретного преобразования Фурье. Ниже представлена шпаргалка для облегчения понимания выполняемых действий, включая настройку масштаба частотных характеристик сигнала с помощью выше­описанного параметра n. Компоненты сигнала аналогичны первому примеру, за исключением параметра n, которы­й равен 100.

Применим к сигналу функции fft() и fftfreq() – и вычислим ДПФ и соответствующие её результатам частоты. Тем самым мы получим спектр рассматриваемого сигнала. Затем построим график, на котором выведем как исходный сигнал, так и полученный спектр.

y_fft = fft(x=xt)
freg_fft = fftfreq(n=n)

fig, axs = pyplot.subplots(1, 2)
axs[0].plot(time_vec, xt, label='Исходный сигнал')
axs[0].set_xlabel('Время (сек.)')
axs[0].set_ylabel('Амплитуда')
axs[0].legend(loc='best')
axs[1].stem((freg_fft.copy()[:n//2] / (1/n))[:50],
            (np.abs(y_fft.copy()[:n//2]) * (2/n))[:50],
            markerfmt='.')
axs[1].set_xlabel('Частота (Гц)')
pyplot.suptitle('Простой периодический синусоидальный сигнал и его спектр')
pyplot.show()

Результаты ДПФ удобно представить с помощью диаграммы со стеблями (stem plot) [ссылка]. Ранее была приведена шпаргалка, которая позволит наглядно понять, какие выполняются действия с «x и y» в качестве её аргументов, и как они изменяются. Тем не менее хотелось бы дать пояснения к представленному коду в функции stem().

Часть с freg_fft.

Сперва создаётся копия оригинального массива, чтобы избежать его изменений. Затем выбирается только первая половина значений частот, поскольку после выполнения ДПФ спектр имеет симметричную структуру, и положительные частоты содержат всю необходимую информацию о частотном содержании сигнала, в то время как отрицательные частоты являются их зеркальным отражением. Далее выполняется деление на «/(1/n)», чтобы частоты соответствовали диапазону фактического сигнала.

Часть с y_fft.

Сначала создаётся копия оригинального массива и берётся только первая половина значений комплексных чисел, полученных из ДПФ (из-за симметрии спектра); далее из этой половины значений извлекается модуль (=получаем амплитуды). Умножение на «*(2/n)» выполняется для масштабирования амплитуд к фактическим значениям сигнала. Это умножение на 2 учитывает вклад как положительных, так и отрицательных частот, поскольку амплитуды для действительных сигналов симметричны.

В конце применяется срез [:50] как к частотам, так и к амплитудам, чтобы отобразить только первые 50 значений для удобного восприятия информации.

В итоге мы легко можем оценить, какие частоты присутствуют в сигнале и с какой амплитудой: в данном сигнале имеется одна-единственная частота 10 Гц с амплитудой, равной 1. Зная это, мы легко может воспроизвести (=восстановить) исходный сигнал с помощью разных методов, которые будут рассмотрены далее. В этом примере давайте восстановим сигнал «вручную»:

A_inv, f_inv, phi = 1, 10, 0
x_inv = A_inv * np.sin(2*np.pi * f_inv * time_vec + phi)

pyplot.figure()
pyplot.plot(time_vec, xt, label='Исходный сигнал')
pyplot.plot(time_vec, x_inv, dashes=[10, 12],
            label='Восстановленный сигнал')
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
pyplot.legend(loc='best')
pyplot.title('Исходный и восстановленный "вручную" сигналы')
pyplot.grid(True)
pyplot.show()

Создадим сигнал посложнее:

n = 1000
time_vec = np.arange(0, 1, 1/n)

A1, f1 = 2, 2
xt1 = A1 * np.sin(2*np.pi * f1 * time_vec)
A2, f2 = 1.5, 6
xt2 = A2 * np.sin(2*np.pi * f2 * time_vec)
A3, f3 = 1, 20
xt3 = A3 * np.sin(2*np.pi * f3 * time_vec)

sum_signal = xt1 + xt2 + xt3

fig, axs = pyplot.subplots(2, 1, figsize=(14, 6))
axs[0].plot(time_vec, sum_signal)
axs[0].set_ylabel('Амплитуда')
axs[1].plot(time_vec, xt1, label='1-ая гармоника (2 Гц)')
axs[1].plot(time_vec, xt2, label='3-ая гармоника (6 Гц)')
axs[1].plot(time_vec, xt3, label='10-ая гармоника (20 Гц)')
axs[1].set_xlabel('Время (сек.)')
axs[1].set_ylabel('Амплитуда')
axs[1].legend(loc='best')
pyplot.suptitle('Суммированный сигнал и его гармоники')
pyplot.show()

Рассчитаем его спектр и выведем результат на графике:

y_fft = fft(x=sum_signal)
freg_fft = fftfreq(n=n)

fig, axs = pyplot.subplots(1, 2, figsize=(14, 3))
axs[0].plot(time_vec, sum_signal, label='Исходный сигнал')
axs[0].set_xlabel('Время (сек.)')
axs[0].set_ylabel('Амплитуда')
axs[0].legend(loc='best')
axs[1].stem((freg_fft.copy()[:n//2] / (1/n))[:50],
            (np.abs(y_fft.copy()[:n//2]) * (2/n))[:50],
            markerfmt='.')
axs[1].set_xlabel('Частота (Гц)')
pyplot.suptitle('Суммированный сигнал и его спектр')
pyplot.show()

С помощью функции ifft() восстановим сигнал по массиву комплексных чисел, полученных в результате выполнения ДПФ:

x_inv = ifft(y_fft)

pyplot.figure()
pyplot.plot(time_vec, sum_signal, label='Исходный сигнал')
pyplot.plot(time_vec, x_inv.real, dashes=[10, 12], label='Восстановленный сигнал')
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
pyplot.legend(loc='best')
pyplot.title('Исходный и восстановленный с помощью ifft() сигналы')
pyplot.grid(True)
pyplot.show()

Выполним фильтрацию сигнала и избавимся от самой высокой частоты, приняв её за шум.

Библиотека scipy.signal предоставляет множество функций для создания фильтров различных типов [ссылка], включая фильтры Баттерворта и Чебышёва. Но мы не станем использовать эти методы: вместо них мы рассмотрим два способа, которые хорошо подходят для наших целей, и просты для понимания. Оба способа имеют свои преимущества и недостатки; и позволяют работать как с фактическими, так и преобразованными значениями частот, как будет продемонстрировано далее.

Первый способ представляет собой фильтрацию по заданным частотам. Сперва определяется список fregs_for_stay с частотами, которые требуется оставить (в данном примере это 2 и 6 Гц и, соответственно, их зеркальные отражения). Далее создаётся маска, проверяющая вхождение полученных частот в фактическом масштабе в указанный список. Если маска истина, значение из y_fft сохраняется, в противном случае оно заменяется на 0.

fregs_for_stay = [2, 6]
mask = np.isin(np.abs(freg_fft) / (1/n), fregs_for_stay)
filtered_fft = np.where(mask, y_fft, 0)

x_inv = ifft(filtered_fft)

pyplot.figure()
pyplot.plot(time_vec, sum_signal, label='Исходный сигнал (xt1+xt2+xt3)')
pyplot.plot(time_vec, xt1+xt2, label='Сигнал с необходимыми частотами (xt1+xt2)')
pyplot.plot(time_vec, x_inv.real, dashes=[10, 12], label='Отфильтрованный сигнал')
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
pyplot.legend(loc='best')
text_1 = "Исходный и отфильтрованный восстановленный по ifft() "
text_2 = "сигналы (1 способ: прямое перечисление необходимых частот)"
pyplot.title(text_1 + text_2)
pyplot.grid(True)
pyplot.show()

Второй способ — частотная фильтрация сигнала по пороговым значениям. Сперва устанавливаются нижний (cutoff_low) и верхний (cutoff_high) пороги для фильтрации преобразованных частот. Затем создаётся маска, которая проверяет, находятся ли частоты в пределах заданных порогов. Значения из y_fft сохраняются, если маска истинна; в противном случае они заменяются на 0.

fig, axs = pyplot.subplots(1, 2, figsize=(14, 3))
axs[0].plot(time_vec, sum_signal, label='Исходный сигнал')
axs[0].set_xlabel('Время (сек.)')
axs[0].set_ylabel('Амплитуда')
axs[0].legend(loc='best')
axs[1].stem(freg_fft.copy()[:50],
            np.abs(y_fft.copy()[:50]) * (2/n),
            markerfmt='.')
axs[1].axvline(x=freg_fft[0], color='orange', linestyle='dashed')
axs[1].axvline(x=freg_fft[10], color='orange', linestyle='dashed')
axs[1].set_xlabel('Нормализованная (d=1) частота (Гц)')
pyplot.suptitle('Суммированный сигнал, его спектр и пороги фильтра')
pyplot.show()

cutoff_low = 0.0
cutoff_high = 0.01
mask = (np.abs(freg_fft) >= cutoff_low) & (np.abs(freg_fft) <= cutoff_high)
filtered_fft = np.where(mask, y_fft, 0)

x_inv = ifft(filtered_fft)

pyplot.figure()
pyplot.plot(time_vec, sum_signal, label='Исходный сигнал (xt1+xt2+xt3)')
pyplot.plot(time_vec, xt1+xt2, label='Сигнал с необходимыми частотами (xt1+xt2)')
pyplot.plot(time_vec, x_inv.real, dashes=[10, 12], label='Отфильтрованный сигнал')
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
pyplot.legend(loc='best')
text_1 = "Исходный и отфильтрованный восстановленный по ifft() "
text_2 = "сигналы (2 способ: с помощью фильтра с указанием нижнего и верхнего порогов)"
pyplot.title(text_1 + text_2)
pyplot.grid(True)
pyplot.show()

Следующий сигнал создадим с помощью функции np.cos(), и добавим равномерный шум. Для генерации равномерного шума применим функцию np.random.uniform(), которая создаст случайные числа, равномерно распределённые между двумя заданными границами: для примера возьмём -1.5 и 1.5.

n = 1000
time_vec = np.arange(0, 1, 1/n)
A, f = 1, 3
uniform_noise = np.random.uniform(-1.5, 1.5, size=time_vec.size)

xt = A * np.cos(2*np.pi * f * time_vec + uniform_noise)

pyplot.figure()
pyplot.plot(time_vec, xt)
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
pyplot.title('Зашумлённый периодический косинусоидальный сигнал')
pyplot.grid(True)
pyplot.show()

y_fft = fft(x=xt)
freg_fft = fftfreq(n=n)

fig, axs = pyplot.subplots(1, 2)
axs[0].plot(time_vec, xt, label='Исходный сигнал')
axs[0].set_xlabel('Время (сек.)')
axs[0].set_ylabel('Амплитуда')
axs[0].legend(loc='center right')
axs[1].stem(freg_fft.copy()[:30] / (1/n),
            np.abs(y_fft.copy()[:30]) * (2/n),
            markerfmt='.')
axs[1].set_xlabel('Частота (Гц)')
pyplot.suptitle('Зашумлённый периодический косинусоидальный сигнал и его спектр')
pyplot.show()

Обратите внимание, согласно спектру амплитуда сигнала на третьей частоте равна приблизительно 0.65, а не 1, как было указано при создании сигнала. Это произошло из-за добавления случайного шума к каждой точке сигнала, включая точку на третьей частоте, который может варьироваться от -1.5 до 1.5.

Создадим маску и отфильтруем сигнал, оставив только частоту 3 (и, следовательно, её зеркальное отражение n-3). Напоминаю, что частоты спектра распределены симметрично относительно его середины. У них одинаковая действительная часть, но мнимая часть имеет противоположные знаки. Это делает их комплексно-сопряжёнными числами, особенность которых заключается в том, что они позволяют упростить вычисления и представить сигналы в более компактной и удобной форме. В связи с этим важно отметить, что для получения исходного сигнала необходимо учитывать ещё и зеркальные частоты (с противоположной фазой). Без них восстановление сигнала может привести к искажению или неправильной интерпретации.

mask = (np.abs(freg_fft) / (1/n) == 3)
filtered_fft = np.where(mask, y_fft, 0)

x_inv = ifft(filtered_fft)

pyplot.figure()
pyplot.plot(time_vec, xt, label='Исходный сигнал')
pyplot.plot(time_vec, x_inv.real, linewidth=5, label='Отфильтрованный сигнал')
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
pyplot.legend(loc='best')
pyplot.title('Исходный и отфильтрованный восстановленный по ifft() сигналы')
pyplot.grid(True)
pyplot.show()

Перейдём к последнему и наиболее близкому к нашему случаю примеру сигнала. Сперва создадим для него несколько сезонных компонент — годовую, недельную и дневную, а также добавим гауссовский шум, который хорошо подходит для создания случайной изменчивости (помех в сигнале). Для генерации гауссовского (или нормального) шума применим функцию np.random.randn() с коэффициентов 2.5, которая создаст случайные числа, распределённые по нормальному закону со средним значением 0 и стандартным отклонением 2.5 соответственно коэффициенту.

Для их выявления будем использовать как спектрограмму, так и периодограмму из библиотеки scipy [ссылка]. С аргументом scaling=’density’ (по умолчанию) периодограмма представляет собой оценку спектральной плотности мощности сигнала [ссылка]. В отличие от спектрограммы, она менее чувствительна к шуму, что позволяет более наглядно выделить частоты, характеризующиеся значительными сезонными колебаниями (сезонные компоненты). В то же время высокочастотные гармоники становятся менее заметными, поэтому вместе с периодограммой будем рассматривать также и спектрограмму сигнала.

n = 26280
time_vec = np.arange(0, 1, 1/n)

A_year, A_week, A_day = 4, 3, 1.5
phi = -np.pi/2

xt_year = A_year * np.sin(2*np.pi * 3 * time_vec + phi)
xt_week = A_week * np.sin(2*np.pi * (3*52) * time_vec)
xt_day = A_day * np.sin(2*np.pi * (3*365) * time_vec )

norm_noise = 2.5 * np.random.randn(time_vec.size)

signal = xt_year + xt_week + xt_day + norm_noise

pyplot.figure()
pyplot.plot(time_vec, signal)
pyplot.xlabel('Время (сек.)')
pyplot.ylabel('Амплитуда')
text_1 = "Зашумлённый синусоидальный сигнал с годовой, "
text_2 = "недельной и дневной сезонностями"
pyplot.title(text_1 + text_2)
pyplot.grid()
pyplot.show()

y_fft = fft(x=signal)
freg_fft = fftfreq(n=n)

periodogramma = periodogram(x=signal)

fig, axs = pyplot.subplots(1, 2)
axs[0].stem(freg_fft.copy()[:2000] / (1/n),
            np.abs(y_fft.copy()[:2000]) * (2/n),
            markerfmt='.')
axs[0].set_xlabel('Частота (Гц)')
axs[0].set_ylabel('Амплитуда')
axs[0].set_title('Спектрограмма сигнала')
axs[1].stem((periodogramma[0][:2000] / (1/n)),
            (periodogramma[1][:2000] * (2/n)),
            markerfmt='.')
axs[1].set_xlabel('Частота (Гц)')
axs[1].set_ylabel('Мощность')
axs[1].set_title('Периодограмма сигнала')
pyplot.show()

По графикам видно, что в сигнале имеется три сезонные компоненты – годовая на частоте 3 Гц (3×1), недельная на частоте 156 Гц (3×52) и дневная на частоте 1095 Гц (3×365). Поскольку мы передали аргументу n функции fftfreq() общее количество точек сигнала n, то показаны совокупные результаты для всех 3-х периодов, отсюда и деление/умножение на 3 при их интерпретации.

Компоненту с самой низкой частотой называют основной. В нашем случае у основной частоты (годовая) самая большая амплитуда, поэтому это ещё и доминирующая частота. Периодограмма более сфокусировано выделила сезонные компоненты путём возведения в квадрат амплитуд ДПФ, в то время как шум стал менее заметным.

Для более удобной интерпретации полученных результатов применим такое понятие, как частота дискретизации (далее – ЧД). В контексте анализа временных рядов, ЧД определяет количество наблюдений (или выборок) за единицу времени. В этой связи важно точно определить, какой промежуток времени считать единицей измерения. В данном примере сигнал имеет 3 годовых периода, каждый из которых, напоминаю, насчитывает 8760 часовых наблюдений в год (3х24х365). В этой связи частота дискретизации 8760 выборок в год (=годовая единица времени) позволит корректно сопоставить результаты анализа с реальными данными.

Fs = 8760

fig, axs = pyplot.subplots(1, 2)
axs[0].stem(freg_fft.copy()[:2000] / (1/Fs),
            np.abs(y_fft.copy()[:2000]) * (2/n),
            markerfmt='.')
axs[0].set_xlabel('Частота (Гц)')
axs[0].set_ylabel('Амплитуда')
axs[0].set_title('Спектрограмма сигнала')
axs[1].stem((periodogramma[0][:2000] / (1/Fs)),
            (periodogramma[1][:2000] * (2/n)),
            markerfmt='.')
axs[1].set_xlabel('Частота (Гц)')
axs[1].set_ylabel('Мощность')
axs[1].set_title('Периодограмма сигнала')
pyplot.show()

Теперь, когда мы фактически уменьшили ЧД с 26280 Гц (3 года) до 8760 Гц (1 год), результаты стали более детальными в плане интерпретации и количественной оценки сезонных компонентов. Видно, что один год имеет годовую, недельную и дневную сезонности с частотами 1, 52 и 365 Гц соответственно. При этом частота дискретизации 8760 Гц гораздо более чем в два раза превышает максимальную частоту 365 Гц в сигнале, что удовлетворяет требованиям теоремы Котельникова, согласно которой ЧД должна быть не менее чем в два раза превышать максимальную частоту сигнала для полного восстановления непрерывного сигнала по его дискретным отсчётам. (Подробное описание многих понятий, которые используются при обработке сигналов, включая частоту заворота, являющуюся половиной ЧД, можно найти в специализированной литературе по данной теме.)

Обратите внимание на метки оси «х»: поскольку мы использовали ЧД Fs, то их кратность составила Fs/n. То есть в данном случае кратность меток снизилась с 1 до 0,3333. Это важно учитывать при использовании ЧД при фильтрации.

Дополнительно хотел бы добавить, что вначале хорошим решением будет использование всех точек сигнала. Это поможет избежать возможной путаницы при интерпретации результатов спектра.

Усложним пример, добавив к сигналу тенденцию, и посмотрим как изменится поведение частот сигнала.

trend = np.linspace(0.06, 2.5, n)
signal = xt_year + xt_week + xt_day + trend + norm_noise

y_fft = fft(x=signal)
freg_fft = fftfreq(n=n)

periodogramma = periodogram(x=signal)

fig, axs = pyplot.subplots(2, 2, figsize=(14, 6))
axs[0, 0].stem(freg_fft.copy()[:2000] / (1/Fs),
            np.abs(y_fft.copy()[:2000]) * (2/n),
            markerfmt='.')
axs[0, 0].set_xlabel('Частота (Гц)')
axs[0, 0].set_ylabel('Амплитуда')
axs[0, 0].set_title('Спектрограмма сигнала')

axs[1, 0].stem(freg_fft.copy()[:200] / (1/Fs),
            np.abs(y_fft.copy()[:200]) * (2/n),
            markerfmt='.')
axs[1, 0].set_xlabel('Частота (Гц)')
axs[1, 0].set_ylabel('Амплитуда')

axs[0, 1].stem((periodogramma[0][:2000] / (1/Fs)),
               (periodogramma[1][:2000] * (2/n)),
               markerfmt='.')
axs[0, 1].set_xlabel('Частота (Гц)')
axs[0, 1].set_ylabel('Мощность')
axs[0, 1].set_title('Периодограмма сигнала')

axs[1, 1].stem((periodogramma[0][:200] / (1/Fs)),
               (periodogramma[1][:200] * (2/n)),
               markerfmt='.')
axs[1, 1].set_xlabel('Частота (Гц)')
axs[1, 1].set_ylabel('Мощность')

pyplot.suptitle('Зашумлённый периодический сигнал с тенденцией и его спектр')
pyplot.show()

Согласно представленным графикам, наличие тенденции вносит значительные изменения в низкочастотный диапазон спектра сигнала: появились пики с различной амплитудой между нулевой и доминирующей частотами. Пик на нулевой частоте говорит о том, что в сигнале имеется постоянная составляющая, о ней мы более подробно поговорим далее. Также наблюдается искажение («дрожь») шума спектра сигнала. На периодограмме появился дополнительный пик рядом с доминирующим, который искажает интерпретацию.

На этом примере закончим рассмотрение примеров сигналов. Далее приступим к спектральному анализу временного ряда с температурными значениями.

df_temp = read_parquet('jena_climate_2009_2016_hour_grouped_data.parquet').asfreq('h')
temperature = df_temp['T (degC)'][:'2011'].to_numpy()

n = len(temperature)
time_vec = np.arange(0, 1, 1/n)
Fs = 8760

Выполним ДПФ и проанализируем результаты.

y_fft = fft(x=temperature)
freg_fft = fftfreq(n=n)

def plot_temperature_spectrum(data):
    """График со спектром температуры"""
    fig, axs = pyplot.subplots(3, 1, figsize=(14, 7))

    axs[0].plot(data['xt'], label='Исходный сигнал')
    axs[0].set_xlabel('Время (сек.)')
    axs[0].set_ylabel('Температура (град. C)')
    axs[1].stem(data['freg_fft'],
                data['y_fft'],
                markerfmt='.')
    axs[1].set_ylabel('Амплитуда')
    axs[1].grid(True)

    axs[2].stem(data['freg_fft'][:4000],
                data['y_fft'][:4000],
                markerfmt='.')
    axs[2].set_xlabel('Частота (Гц)')
    axs[2].set_ylabel('Амплитуда')
    axs[2].grid(True)

    pyplot.suptitle('Дискретное преобразование Фурье на исходном ВР')
    pyplot.show()

plot_temperature_spectrum(data={
    'xt': temperature,
    'freg_fft': freg_fft[:n//2] / (1/Fs),
    'y_fft': np.abs(y_fft[:n//2]) * (2/n)
    })

По спектру видно, что нулевая частота имеет самую большую амплитуду (в нашем случае равную 17). Это означает, что сигнал имеет так называемую постоянную составляющую. Наличие её в сигнале, равно как и тенденции, может сильно усложнить (или исказить) анализ сигнала и интерпретацию его спектральных характеристик. Поэтому для качественного проведения анализа и более точного выявления его компонентов их обычно извлекают.

На самом деле, ДПФ позволяет работать с данными и без их предварительной подготовки в виде удаления постоянной составляющей и тенденции, однако общее время обработки данных возрастёт. Искушение использовать все выявленные частоты, чтобы получить замечательные результаты, велико. Вопрос в том, целесообразно ли это? В этом исследовании предлагаю не поддаваться искушению — и попытаться смоделировать годовую сезонность, использовав только общие частоты для всех трёх периодов. Благодаря этому мы сэкономим время и вычислительные ресурсы. А затем проверим, как улучшиться предсказательная способность модели.

Для начала нам нужно решить в какой последовательности обработать (=извлечь) низкочастотные компоненты — тенденцию и постоянную составляющую. В литературе и на интернет-ресурсах описаны различные подходы к декомпозиции временных рядов. Например, в классических методах tsa.seasonal.seasonal_decompose [ссылка] и tsa.seasonal.STL [ссылка] библиотеки statsmodels определены процедуры, в которых итеративно выделяются компоненты временного ряда – тенденция, сезонность и остаток. При этом порядок извлечения компонентов может варьироваться в зависимости от конкретной задачи и используемого подхода, хотя в большинстве методов тенденция извлекается посредством сглаживания ряда, после чего уточняются сезонная и остаточная составляющие.

В ходе написания данного исследования мною были проведены эксперименты, направленные на корректное выделение сезонных колебаний. Результаты показали, что порядок извлечения компонентов имеет существенное значение: при последовательности «сначала тенденция, затем постоянная составляющая» на спектрограмме отсутствует пик на нулевой частоте, тогда как при обратном порядке наблюдается заметный ненулевой пик. Важно отметить, что извлечение тенденции осуществлялось несколькими способами. Сначала с использованием функции statsmodels.tsa.seasonal.seasonal_decompose, основанной на процедуре свёртки. Применение данной процедуры к большим временным рядам может привести к ряду проблем, особенно если не учитывать особенности граничных значений временного ряда в зависимости от выбранного режима алгоритма (применительно к функции seasonal_decompose). Затем использовался более гибкий метод под названием Локальное взвешенное сглаживание разброса (в англоязычной терминологии – Locally Weighted Scatterplot Smoothing, LOWESS), который применяется не только для сглаживания данных, но и для выявления тенденций и локальных изменений во временных рядах [ссылка]. Более того, с помощью функции ЛВСР выполняется разложение временного ряда на компоненты по методу «Сезонно-трендовое разложение с использованием ЛВСР» (в англоязычной терминологии – Seasonal-Trend decomposition using LOESS, STL) [ссылка].

Оба метода имеют свои плюсы и минусы, однако применительно к нашему временному ряду метод ЛВСР оказался предпочтительнее. На графике ниже представлено сравнение их результатов. Видно, что тенденция, извлечённая по методу свёртки, имеет загнутые концы. Это говорит о том, что для расчёта свёртки в точках, близких к краям ряда, часть оконного фильтра свёртки оказалась вне границ массива и заполнилась нулями, что исказило расчёт тенденции на границах. Напротив, тенденция, извлечённая по медоту ЛВСР лишена данного недостатка; отсюда экстраполяция её будущих значений не вызовет проблем.

Отдельно заметьте, что её динамика (при frac=0.5) отражает годовые колебания: подъёмы приходятся на середину лета, а спады — на середину зимы.

Рассмотрим оба способа.

Свёртка с использованием оконного фильтра

Сначала извлечём из данных тенденцию с помощью функции get_trend(). Подход, применённый в этой функции, аналогичен подходу извлечения тенденции в методе seasonal_decompose() [ссылка] из библиотеки statsmodels с тем лишь отличием, что мы сразу указываем, что наш период чётный (8760 точек).

Сперва создаётся фильтр в виде массива filt, который состоит из 0.5 в начале и в конце, а промежуток между началом и концом заполняется единицами. Затем этот массив делится на длину периода period, что нормализует его значения. Это важно для сглаживания, чтобы не изменять общий уровень сигнала. Функция convolve() из библиотеки numpy выполняет свёртку между исходными данными и созданным фильтром. Аргумент mode=’same’ означает, что результат свёртки будет иметь ту же длину, что и исходные данные. Это достигается за счёт дополнения данных нулями по краям перед свёрткой. Результатом функции get_trend() является массив чисел, представляющих собой извлечённую тенденцию. Поскольку в предыдущей части главы [ссылка] мы строили аддитивную модель, то далее вычтем тенденцию из временного ряда и тем самым получим обработанную температуру.

def get_trend(data, period):
    """Извлечение тенденции"""
    filt = np.array([0.5] + [1] * (period - 1) + [0.5]) / period
    return np.convolve(data, filt, mode='same')

trend = get_trend(data=temperature.copy(), period=8760)
processed_temp = temperature.copy() - trend

Отцентрируем обработанный сигнал вокруг нуля, убрав из него среднее значение. Затем заново выполним ДПФ и проанализируем результаты.

def get_centered_temperature(data):
    """Центрирование сигнала вокруг нуля"""
    mean_value = np.mean(data)
    return mean_value, data - mean_value

mean_value, processed_temp = get_centered_temperature(data=processed_temp.copy())

y_fft = fft(x=processed_temp)
freg_fft = fftfreq(n=n)

def plot_processed_temperature_spectrum(data):
    """График со спектром обработанной температуры"""
    fig, axs = pyplot.subplots(2, 1, figsize=(14, 7))

    axs[0].stem(data['freg_fft'],
                data['y_fft'],
                markerfmt='.')
    axs[0].set_ylabel('Амплитуда')
    axs[0].grid(True)

    axs[1].stem(data['freg_fft'][:data['axs1_slice']],
                data['y_fft'][:data['axs1_slice']],
                markerfmt='.')
    axs[1].set_xlabel('Частота (Гц)')
    axs[1].set_ylabel('Амплитуда')
    axs[1].grid(True)

    pyplot.suptitle('Дискретное преобразование Фурье на обработанном ВР')
    pyplot.show()

plot_processed_temperature_spectrum(data={
    'freg_fft': freg_fft[:6000] / (1/Fs),
    'y_fft': np.abs(y_fft[:6000]) * (2/n),
    'axs1_slice': 1200
    })

Итак, давайте проанализируем результаты. В начале отметим, что после выполнения вышеописанных действий основной годовой частотой сигнала стала 1-я частота, а её амплитуда хоть и уменьшилась до 9.4, но по-прежнему является самой большой из всех амплитуд. На нулевой частоте пик отсутствует.

Также имеется ещё одна основная частота — 365 (дневная). У этой частоты есть 4-е кратные ей гармоники. Обратите вниманием на их частотное представление. Наличие нескольких сопутствующих частот в гармонике говорит о её более сложной форме, чем простая синусоидальная структура. Подобное частотное представление П. Блюмфилд в своей книге охарактеризовал как «несинусоидальные колебания», «одним из аспектов которого могут быть более плоские впадины, чем пики». В данном случае при отдельном восстановлении сигнала по частотам 362-368 «гармоники 365» будет наблюдаться постепенное увеличение амплитуд к «шапке» годового периода с последующим зеркальным уменьшением, как показано на графике ниже. При этом видно, что у всех периодов поведение отличается.

Помимо этого сразу бросается в глаза «частокол» частот в диапазоне от 1 до – примерно – 80 с последующим плавным спадом в амплитудных значениях. Это означает, что годовые циклы имеют более сложную форму, чем кажущуюся. На графике ниже это наглядно продемонстрировано.

Локальное взвешенное сглаживание разброса

У функции lowess() библиотеки statsmodels есть важный параметр frac, который определяет долю данных (или ширину окна), в которой выполняется локальная регрессия. Это позволяет настроить уровень сглаживания в зависимости от структуры рассматриваемых данных. С помощью функции plot_lowess_compare() построим график, на котором отобразим исходный сигнал, извлечённую тенденцию по методу свёртки и тенденции по методу ЛВСР с различными значениями параметра frac, чтобы наглядно оценить, как разные значения frac влияют на сглаживание.

def plot_lowess_compare(data, frac_values):
    """График с исходным сигналом и тенденциями по свёртке и ЛВСР для разных значений frac"""
    pyplot.figure(figsize=(14, 6))

    pyplot.plot(data['time_vec'], data['temperature'], color='lightgray', label='Исходный сигнал')
    pyplot.plot(data['time_vec'], data['trend'], label='Свёртка (trend_convolve_same)')

    for frac in frac_values:
        lowess_res = lowess(data['temperature'], data['time_vec'], frac=frac)
        pyplot.plot(lowess_res[:, 0], lowess_res[:, 1], label=f'Локальная регрессия (lowess) (frac={frac})')

    pyplot.title('Сравнение тенденций по параметру frac')
    pyplot.xlabel('Время (сек.)')
    pyplot.ylabel('Температура (град. C)')
    pyplot.legend(loc='best')
    pyplot.grid(True)
    pyplot.show()

plot_lowess_compare(data={
    'time_vec': time_vec,
    'temperature': temperature.copy(),
    'trend': trend},
    frac_values=[0.1, 0.2, 0.3, 0.4, 0.5]
    )

По графику видно, что при малом значении frac ЛВСР более чувствительно к локальным изменениям, а при большом — более сглаженное. Значение frac=0.5 будет выбрано по причине того, что оно обеспечивает поведение, наиболее близкое к выполненной тенденции по процедуре свёртки; но в отличие от неё информация на границах временного ряда не теряется, что наглядно продемонстрировано на графике.

Итак, извлечём из данных тенденцию с помощью функции get_trend_lowess() и вычтем её из сигнала.

def get_trend_lowess(endog, exog, frac):
    """Извлечение тенденции по методу ЛВСР"""
    print("Выполнение расчёта тенденции по методу ЛВСР ...")
    trend = lowess(endog, exog, frac)
    return trend[:, 1]

trend_lowess = get_trend_lowess(endog=temperature.copy(),
                                exog=time_vec,
                                frac=0.5)
processed_temp = temperature.copy() - trend_lowess

Отобразим спектр обработанного сигнала.

y_fft = fft(x=processed_temp)
freg_fft = fftfreq(n=n)

plot_processed_temperature_spectrum(data={
    'freg_fft': freg_fft[:1200] / (1/Fs),
    'y_fft': np.abs(y_fft[:1200]) * (2/n),
    'axs1_slice': 80
    })

По спектру сигнала видно, что при извлечении тенденции по методу ЛВСР с уровнем сглаживания frac=0.5 нулевая частота практически отсутствует. Это говорит о том, что извлечённый тренд включит в себя постоянную составляющую и при её вычитании мы фактически вычли и среднее значение. Дополнительное (но не обязательное) вычитание среднего из полученного сигнала гарантирует отсутствие постоянной составляющей. На спектрограмме ниже это наглядно видно.

mean_value, processed_temp = get_centered_temperature(data=processed_temp.copy())
# mean = 0.017714

y_fft = fft(x=processed_temp)
freg_fft = fftfreq(n=n)

plot_processed_temperature_spectrum(data={
    'freg_fft': freg_fft[:6000] / (1/Fs),
    'y_fft': np.abs(y_fft[:6000]) * (2/n),
    'axs1_slice': 1200
    })

По результатам понятно, что анализ полученного спектра требует немало времени и значительных усилий для выявления необходимого количества частот с целью качественного прогнозирования годовой (да и дневной) сезонности. При этом необходимо учитывать объём обучающих данных и нерезиновые вычислительные ресурсы. В связи с этим я передаю задачу тщательной фильтрации и качественного восстановления сигнала с проведением множества интересных экспериментов на рассмотрение последующим исследователям и энтузиастам. Делитесь анализом в комментариях, с удовольствием обсудим. В дальнейшем будут продемонстрированы результаты, полученные с использованием только визуально выбранных частот, это позволит оценить эффективность данного подхода.

Внимательные читатели должно быть обратили внимание, что в примерах создания сигнала указывалась как функцию синуса, так и косинуса, и каждый раз функция ifft() верно восстанавливала сигнал, хотя мы не передавали ей информацию как он был построен. Давайте проверим, корректно ли мы вручную восстанавливаем сигнал. Отфильтруем временной ряд по 3-й частоте (n=n=len(temperature)), а затем сравним восстановление отфильтрованного сигнала по обоим способам — как с помощю ifft(), так и вручную.

def filter_signal(y_fft, n, fregs_for_stay=3):
    """Фильтрация сигнала по частоте №3"""
    freg_fft = fftfreq(n=n)
    mask = np.isin(np.abs(freg_fft) / (1/n), fregs_for_stay)
    filtered_fft = np.where(mask, y_fft, 0)
    return ifft(filtered_fft).real

x_inv = filter_signal(y_fft=y_fft,
                      n=n)

def show_ifft_result_plot(data):
    """График с обработанным исходным и отфильтрованным сигналами"""
    pyplot.figure()
    pyplot.plot(data['processed_signal'], label='Обработанный исходный сигнал')
    pyplot.plot(data['filtered_signal'], linewidth=3, label='Отфильтрованный сигнал')
    pyplot.xlabel('Время (сек.)')
    pyplot.ylabel('Амплитуда')
    pyplot.legend(loc='best')
    pyplot.title(data['title'])
    pyplot.grid(True)
    pyplot.show()

show_ifft_result_plot(data={
    'processed_signal': processed_temp,
    'filtered_signal': x_inv,
    'title': 'Обработанный исходный и восстановленный с помощью ifft() '
    'отфильтрованный сигналы'
    })

Восстановим сигнал вручную. Сначала с помощью функции fft_result_filter() отфильтруем данные по передаваемым в списке fregs_for_stay частотам. В основе функции уже известные вам действия. Для удобства предусмотрена фильтрация как по фактическим значениям частот (в этом случае явно указывается параметр d), так и преобразованными значениями частот по умолчанию (d=None). Затем выполняется ещё одна, результирующая, фильтрация по полученным частотам, чтобы оставить только те из них, у которых комплексные числа не равны нулю. Для этого будет использоваться dataframe (df) из библиотеки pandas, фильтрация по которому очень удобна и проста. Можно использовать и словарь python с выражением-генератором с фильтрацией значений, но команда выйдет не такой изящной и лёгкой для понимания.

Выполним фильтрацию по фактической частоте 1, определённой по спектрограмме с использованием частоты дискретизации, передав для этого параметру d значение Fs.

def fft_results_filter(fregs_for_stay, n, y_fft, d=None):
    """Фильтрация y_fft по частотам"""
    if d is None:
        freg_fft = fftfreq(n=n)
        mask = np.isin(np.abs(freg_fft) / (1/n), fregs_for_stay)

    else:
        if not isinstance(d, int):
            raise ValueError("d должно быть целым числом")
        freg_fft = fftfreq(n=n, d=1/d)
        mask = np.isin(np.abs(freg_fft), fregs_for_stay)

    filtered_fft = np.where(mask, y_fft, 0)

    df = DataFrame(data={
        'y_fft': filtered_fft,
        'fregs': freg_fft if d is None else freg_fft * 1/d
        })
    return df[df['y_fft'] != 0]

filtered_df = fft_results_filter(
    fregs_for_stay=[1],
    n=n,
    y_fft=y_fft,
    d=Fs)

Содержимое полученной переменной продемонстрировано на следующем рисунке.

Далее нам нужно восстановить сигнал по отфильтрованным данным с помощью функции get_signal_harmonics(). Эта функция использует векторные/матричные операции перемножения, что значительно ускоряет выполнение расчётов. Она имеет два аргумента: dataframe (отфильтрованные данные) и time (массив моментов времени).

Сначала из dataframe извлекаются столбцы y_fft и fregs, преобразуются в numpy-массивы и приводятся к фактическим значениям. Затем с помощью reshape(-1, 1) форма массивов изменяется с 1Д на 2Д (например, если в dataframe размерность столбца составляет (2,), то после reshape(-1, 1) мы получаем размерность (2, 1)). В итоге мы получаем две матрицы An и fn с размерностью (2, 1) в фактических масштабах значений.

Далее создаётся 1Д-массив комплексных чисел xt для хранения результата. Его начальная размерность составит (len(dataframe),), что соответствует количеству строк в dataframe. В данном примере размерность вектора составит (2,).

Затем выполняется перемножение массивов с использованием алгоритма автоматического расширения массивов numpy broadcasting. Поскольку time – это вектор с размерностью (26280,), то numpy «растянет» An и fn до размерности (2, 26280). Это означает, что каждое значение из An и fn будет применено ко всем значениям time для каждой гармоники, а time будет «размножен» по строкам, чтобы соответствовать размерности 2 (то есть у нас будет 2 строки и 26280 столбцов). Результатом будет массив xt с размерностью (2, 26280).

На последнем шаге создаётся словарь, в котором ключи – это индексы оригинального dataframe, а значения – строки массива xt. Каждая строка xt соответствует результату для каждого индекса и представляет собой одномерный массив длиной 26280 точек.

Итоговая структура данных – это словарь, где каждому индексу соответствует массив комплексных чисел.

def get_signal_harmonics(dataframe, time):
    """Получение гармоник отфильтрованного сигнала с помощью функции np.cos()"""
    An = (dataframe['y_fft'].to_numpy() * (1/n)).reshape(-1, 1)
    fn = (dataframe['fregs'].to_numpy() / (1/n)).reshape(-1, 1)
    xt = np.zeros(len(dataframe), dtype=complex)
    xt = An * np.cos(2 * np.pi * fn * time)
    return {str(index): xt_row for index, xt_row in zip(dataframe.index, xt)}

harmonics_dict = get_signal_harmonics(dataframe=filtered_df.copy(),
                                      time=time_vec)

def show_harmonics_processed_signal(data):
    """График сигнала, восстановленного вручную по гармоникам"""
    text_1 = 'Исходный обработанный сигнал'
    text_2 = 'Отфильтрованный сигнал x_inv'
    text_3 = 'Отфильтрованный сигнал df'

    pyplot.figure()
    pyplot.plot(data['time_vec'], data['processed_temp'], color='lightgray', label=text_1)
    pyplot.plot(data['time_vec'], data['x_inv'].real, label=text_2)
    pyplot.plot(data['time_vec'], data['sum_signal'].real, label=text_3, linestyle='dashed')
    pyplot.xlabel('Время (сек.)')
    pyplot.ylabel('Амплитуда')
    pyplot.legend(loc='best')
    pyplot.title('Сравнение восстановленных отфильтрованных сигналов')
    pyplot.grid(True)
    pyplot.show()

show_harmonics_processed_signal(data={
    'time_vec': time_vec,
    'processed_temp': processed_temp,
    'x_inv': x_inv,
    'sum_signal': np.sum(list(harmonics_dict.values()), axis=0)
    })

Результаты восстановления сигнала демонстрируют различия, которые объясняются использованием комплексной экспоненты в функции scipy.signal.ifft() [ссылка]. Подправим функцию и заново восстановим сигнал.

def get_signal_harmonics_exp(dataframe, time):
    """Получение гармоник отфильтрованного сигнала с помощью np.exp()"""
    An = (dataframe['y_fft'].to_numpy() * (1/n)).reshape(-1, 1)
    fn = (dataframe['fregs'].to_numpy() / (1/n)).reshape(-1, 1)
    xt = np.zeros(len(dataframe), dtype=complex)
    k = np.arange(len(time))
    xt = An * np.exp(2j * np.pi * fn * k / len(time))
    return {str(index): xt_row for index, xt_row in zip(dataframe.index, xt)}

harmonics_dict = get_signal_harmonics_exp(dataframe=filtered_df.copy(),
                                          time=time_vec)

show_harmonics_processed_signal(data={
    'time_vec': time_vec,
    'processed_temp': processed_temp,
    'x_inv': x_inv,
    'sum_signal': np.sum(list(harmonics_dict.values()), axis=0)
    })

Результаты восстановления сигнала стали идентичными. Мы устранили несоответствие в восстановлении и теперь функция работает как следует.

Давайте ещё раз посмотрим на спектрограмму.

plot_processed_temperature_spectrum(data={
    'freg_fft': freg_fft[:3000] / (1/n),
    'y_fft': np.abs(y_fft[:3000]) * (2/n),
    'axs1_slice': 100
    })

Для моделирования годовой сезонности я буду использовать диапазон частот от 1 до 19 включительно (при d=n). Сравнение исходного обработанного сигнала с восстановленным по этим частотам сигналом представлено на следующем графике.

fregs_for_stay = list(range(1, 20))  # d=n

filtered_df = fft_results_filter(
    fregs_for_stay=fregs_for_stay,
    n=n,
    y_fft=y_fft)

harmonics_dict = get_signal_harmonics_exp(dataframe=filtered_df.copy(),
                                          time=time_vec)

def show_harmonics_sum_signal(data):
    """График сигнала, восстановленного вручную по гармоникам"""
    text_1 = 'Обработанный исходный сигнал'
    text_2 = 'Отфильтрованный сигнал'

    pyplot.figure()
    pyplot.plot(data['time_vec'], data['processed_temp'], label=text_1)
    pyplot.plot(data['time_vec'], data['sum_signal'].real, label=text_2, linestyle='dashed')
    pyplot.xlabel('Время (сек.)')
    pyplot.ylabel('Амплитуда')
    pyplot.legend(loc='best')
    pyplot.title('Сравнение сигналов')
    pyplot.grid(True)
    pyplot.show()

show_harmonics_sum_signal(data={
    'time_vec': time_vec,
    'processed_temp': processed_temp,
    'sum_signal': np.sum(list(harmonics_dict.values()), axis=0)
    })

Отобразим на графике гармоники полученной годовой сезонности по отдельности и внимательно их рассмотрим.

def show_harmonics(harmonics_dict):
    """График с гармониками отфильтрованного сигнала"""
    df = DataFrame({key: value.real for key, value in harmonics_dict.items()})

    df.plot(subplots=True,
            layout=(7, 5),
            figsize=(14, 7),
            title='Гармоники отфильтрованного сигнала')

    pyplot.show()

show_harmonics(harmonics_dict=harmonics_dict)

Видно, что полученная годовая сезонность определяется гармониками с различными амплитудами, частотами и фазами. Поэтому – в данном случае – экстраполирование каждой гармоники отдельно от суммированного сигнала – это далеко не самая лучшая идея. Обычно при наличие одной сезонности экстраполируют «целостный», суммированный, сигнал, учитывающий совокупное поведение всех периодических составляющих. Такой подход сохраняет взаимосвязь между всеми гармониками, что, в свою очередь, позволяет моделировать сезонность более точно и с меньшими затратами на вычислительные ресурсы. Однако если имеется несколько сезонных компонент, то можно выбрать другой подход, при котором модель получает на вход сразу несколько экзогенных переменных, каждая из которых отвечает за свою сезонность – например, одну для годовой и другую для дневной, как в нашем случае. Такой метод позволяет явно разделить и смоделировать каждую циклическую компоненту, не сводя их к единому суммарному сигналу, и тем самым избавится от фактора «занижения» или искажения колебаний одной сезонной компоненты другой. В итоге модель может более точно учитывать особенности каждого временного диапазона, особенно если годовая и дневная сезонности существенно различаются по амплитуде и фазе, что должно повысить точность прогнозирования и интерпретируемость результатов.

Во избежания излишнего увеличения объёма этой части главы, ограничимся применением только первого подхода. (В конце данной части главы будет проведён опрос: если данное направление найдёт интерес у широкого читателя, то в будущих статьях рассмотрим и другие подходы.)

Давайте займёмся подготовкой входных данных.

def prepare_input_data(data):
    """Подготовка входных данных модели"""
    return {
        'endog_train': data['endog_train'],
        'endog_val': data['endog_val'],
        'exog_train': Series(data['exog_train'],
                             index=data['exog_train_index'],
                             name='Filtered_signal_train')
        }

input_data = prepare_input_data(data={
    'endog_train': df_temp['T (degC)'][:'2011'].copy(),
    'endog_val': df_temp['T (degC)']['2012-01-01':'2012-01-07'].copy(),
    'exog_train': np.sum(list(harmonics_dict.values()), axis=0) + mean_value + trend_lowess,
    'exog_train_index': df_temp[:'2011'].index
    })

Подготовка экзогенной валидационной части требует дополнительных действий. Сперва экстраполируем полученную суммированную годовую сезонность с помощью функции get_sum_signal_ext(). Используя функцию show_sum_signal_ext(), отобразим на графике обе сезонности – «тренировочную» и «валидационную».

def get_sum_signal_ext(dataframe, time):
    """Экстраполяция отфильтрованного суммированного сигнала годовой сезонности"""
    An = (dataframe['y_fft'].to_numpy() * (1/n)).reshape(-1, 1)
    fn = (dataframe['fregs'].to_numpy() / (1/n)).reshape(-1, 1)
    # xt = np.zeros(len(dataframe), dtype=complex)
    k = np.arange(len(time))
    return np.sum(list(An * np.exp(2j * np.pi * fn * k / len(time))), axis=0)

harmonics_values_ext = get_sum_signal_ext(dataframe=filtered_df.copy(),
                                          time=np.arange(1, 2, 1/n))

def show_sum_signal_ext(data):
    """График с исходной и экстраполированной сезонностями"""
    pyplot.figure()
    pyplot.plot(data['time_vec'], data['sum_signal'], label='Исходная')
    pyplot.plot(data['time_vec_ext'], data['sum_signal_ext'], label='Экстраполированная')
    pyplot.xlabel('Время (сек.)')
    pyplot.ylabel('Амплитуда')
    pyplot.legend(loc='best')
    pyplot.title('Сравнение сезонностей')
    pyplot.grid(True)
    pyplot.show()

show_sum_signal_ext(data={
    'time_vec': time_vec,
    'time_vec_ext': np.arange(1, 2, 1/n),
    'sum_signal': (np.sum(list(harmonics_dict.values()), axis=0)).real,
    'sum_signal_ext': harmonics_values_ext.real
    })

Следует понимать, что при использовании большого количества частот, определённых с помощью дискретного преобразования Фурье для восстановления сигнала, годовые сезонные колебания для всех рассматриваемых годов становятся менее унифицированными.

Далее экстраполируем тенденцию, извлечённую по методу ЛВСР. Не лишним будет напомнить, что Дж.Бокс и Г.Дженкинс писали про экстраполяцию: «Тренд можно подогнать полиномом, а сезонную компоненту — рядом Фурье. Прогноз производится путём экстраполяции этих подогнанных функций в будущее» [из книги «Анализ временных рядов: прогноз и управление /Дж.Бокс и Г.Дженкинс. Выпуск 1. – Издательство «Мир», Москва 1974. – 331 с.»]. В данном исследовании вместо полиномиальной регрессии мы будем использовать модель Хольта [ссылка], которая реализует двойное экспоненциальное сглаживание (=учитывает текущий уровень и тренд). Поскольку при параметре frac=0.5 мы получили чистую динамику тенденции без явной сезонной компоненты, но с захватом постоянной составляющей (будем рассматривать её как уровень), то применение данной модели оправдано. С определением параметра initialization_method=’estimated’ модель самостоятельно оценит начальный уровень и тенденцию, а с помощью метода optimized=True (по умолчанию) будут автоматически подобраны оптимальные настройки для их сглаживания.

def trend_extrapolation(data, steps=168, plot=True):
    """Выполнение прогноза тенденции с помощью модели Хольта и построение графика"""
    model = Holt(endog=data,
                 initialization_method='estimated')
    fitted_model = model.fit(optimized=True)
    trend_ext = fitted_model.forecast(steps=steps)

    if plot:
        pyplot.figure()
        pyplot.plot(np.arange(len(data)), data, label='Оригинальная')
        pyplot.plot(np.arange(len(data), len(data) + steps), trend_ext, label='Экстраполированная')
        pyplot.title('Сравнение тенденций')
        pyplot.legend()
        pyplot.show()

    return trend_ext

trend_ext = trend_extrapolation(data=trend_lowess)

Теперь мы имеем все необходимые компоненты для определения валидационной экзогенной переменной. Сложим компоненты переменной и добавим её в словарь с входными данными модели.

input_data['exog_val'] = Series(data=harmonics_values_ext[:168] + mean_value + trend_ext,
                                index=df_temp['2012-01-01':'2012-01-07'].index,
                                name='Filtered_signal_val')

Далее приступим к выполнению модели SARIMAX(5, 0, 1)(2, 1, 0, 24). Напоминаю, что при использовании экзогенных переменных в модели SARIMAX библиотеки statsmodels – statsmodels.tsa.statespace.sarimax – возникает ошибка при обучении (экзогенные параметры фактически игнорируются), поэтому применим модель ARIMA, у которой данная проблема не возникает. С учётом использования сезонного порядка и экзогенных переменных эта модель будет эквивалентна SARIMAX.

sarimax_model = ARIMA(endog=input_data['endog_train'],
                      exog=input_data['exog_train'].values.real,
                      order=(5, 0, 1),
                      seasonal_order=(2, 1, 0, 24),
                      trend=None,
                      enforce_stationarity=False,
                      enforce_invertibility=False)

sarimax_model = sarimax_model.fit()

preds = sarimax_model.get_prediction(start=input_data['endog_val'].index.min(),
                                     end=input_data['endog_val'].index.max(),
                                     exog=input_data['exog_val'].values.real)

mape = mean_absolute_percentage_error(y_true=input_data['endog_val'],
                                      y_pred=preds.predicted_mean)

Отобразим полученные результаты на графике.

def show_plot_sarimax_prediction(data, predictions, error_value):
    """График сравнения истинных значений с предсказанными"""
    pyplot.figure()

    pyplot.plot(data, label='Валидационные данные')

    label_mape = f" [mape: {error_value*100:.2f}%]"
    pyplot.plot(predictions, dashes=[4, 3],
                label='Предсказанные данные '+ label_mape)

    pyplot.xlabel('Дата наблюдения')
    pyplot.ylabel('Температура (град C)')
    pyplot.title('Предсказания на валидационных данных с использованием экзог.перем.')
    pyplot.grid(True)
    pyplot.legend(loc='best')
    pyplot.show()

show_plot_sarimax_prediction(data=input_data['endog_val'],
                             predictions=preds.predicted_mean,
                             error_value=mape)

# print(sarimax_model.summary())

Вероятно, полученные результаты могут сильно расстроить (или даже обидеть) читателя. Но не в нашем исследовании! График показывает, что нам удалось смоделировать (=учесть) годовую динамику: данная валидационная часть данных имеет тенденцию к плавному понижению, что фактически и отразилось на предсказании. С помощью следующего графика наглядно убедимся в этом.

def show_plot_input_data(data):
    """График входных данных"""
    fig, axs = pyplot.subplots(2, 1, figsize=(14, 8))

    axs[0].plot(data['endog_train'], label='Эндог. трен.')
    axs[0].plot(data['endog_val'], label='Эндог. вал.')
    axs[0].plot(data['exog_train'].index, np.real(data['exog_train']), label='Экзог. трен.')
    axs[0].plot(data['exog_val'].index, np.real(data['exog_val']), label='Экзог. вал.')
    axs[0].set_xlabel('Дата наблюдения')
    axs[0].set_ylabel('Температура (град. C)')
    axs[0].legend(loc='best')

    axs[1].plot(data['endog_val'].values, color='orange', label='Эндог. вал.')
    axs[1].plot(np.real(data['exog_val']), color='red', label='Экзог. вал.')
    axs[1].set_xlabel('Количество валидационных наблюдений')
    axs[1].set_ylabel('Температура (град. C)')
    axs[1].legend(loc='best')

    pyplot.suptitle('График входных данных')
    pyplot.show()

show_plot_input_data(data=input_data)

Более того, результаты исследования показывают, что использование данного подхода позволило улучшить точность предсказания по сравнению с моделью АРПСС(24, 1, 24), использованной в предыдущей части исследования [ссылка]: ошибка составила 67.37% против 63.52%. Однако предсказательная способность модели сильно уступает модели САРПСС, в которой в качестве экзогенной переменной используется целевой временной ряд, а для валидации — последние 168 точек обучающей выборки данных. Напомню: точность предсказаний той модели достигла 78.26%.

Другим немаловажным выводом является то, что дневная сезонность, выраженная четырьмя гармониками, оказывает значительное влияние на волатильность температуры (её изменчивость) как в течение суток, так и в течение недели. Наше первоначальное предположение о том, что мы можем отделаться моделированием лишь годовой сезонности, чтобы значительно улучшить предсказательную способность моделей, оказалось настолько потрясающим, что можно со смеху сложится напополам. Необходимо провести дополнительные интересные эксперименты, так как без учёта дневной сезонности дальнейшие действия не будут иметь большого смысла.

Среди возможных решений можно выделить следующие:

• Использовать суммированный сигнал в качестве экзогенной переменной, который будет учитывать как годовую, так и дневную сезонность.

• Моделировать сезонности в отдельных переменных: годовую, как это было сделано ранее, и дневную. Дневную сезонность можно дополнительно разделить по гармоникам и использовать их в качестве отдельных экзогенных переменных.

В следующей части исследования мы можем выполнить это вместе. Однако возникает вопрос: интересно ли вам развитие данного направление или предпочтительнее сразу перейти к рассмотрению модели глубокого обучения с использованием слоя ДКП (LSTM), которую я активно применяю в своей профессиональной деятельности? Чтобы это выяснить, ниже прикрепляется соответствующий опрос.

P.S.

Не могу удержаться, чтобы не поделиться ссылкой на очень интересную страницу «Интерактивное введение в преобразования Фурье».