Выделение периодической компоненты временного ряда (пример)

Материал из MachineLearning.

(Различия между версиями)

Версия 13:53, 22 мая 2011

Содержание

1 Введение
2 Постановка задачи
3 Пути решения задачи
4 Вычислительный эксперимент
- 4.1 Исследование модельного зашумленного синуса.
- 4.2 Исследование величины коэффициента корреляции в зависимости от накладываемого шума.
5 Исследование величины коэффициента корреляции при уменьшении числа отсчетов за период.
6 Исследование величины коэффициента корреляции при сокращении временного ряда.
7 Исследование реального временного ряда электрокардиограммы.
8 Исследование применения МНК для прогнозирования реальных временных рядов.
9 Заключение
10 См. также
11 Литература

Введение

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

Автокорреляционная функция исследует временной ряд на наличие периодической компоненты, сдвигая ряд на несколько временных отсчетов и сравнивая с самим собой.

Метод наименьших квадратов (МНК) оценивает параметры для тригонометрической аппроксимации данного ряда. Так как любая последовательность, обладающая периодичностью может быть разложена в ряд Фурье, необходимо принять коэффициенты перед синусами и косинусами за коэффициенты регрессии и оценить их величину. Если найденная корреляция (коэффициент при определенном синусе или косинусе) велика, то можно заключить, что существует строгая периодичность на соответствующей частоте в данных.

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

Постановка задачи

Дан временной ряд $X_t$ , где $n$ - длина временного ряда, $t\in\{1,...,n\}$ - номер отсчета. Предполагаем, что в рассматриваемом временном ряде нет пропущенных значений, и он имеет периодические составляющие с периодом $T=\{\tau_1,..\tau_p\}$ . Работа состоит из трех следующих ступеней.

Во-первых, тестирование на модельной задаче. Дан зашумлённый синус с известным периодом. Необходимо исследовать изменение коэффициента корреляции в следующих ситуациях:

       1. при увеличении шума;
       2. при уменьшении числа отсчётов на период;
       3. при сокращении длины временного ряда.

Во-вторых, тестирование на реальном временном ряде. Дан временной ряд электрокардиограммы, включающий периодическую компоненту со сложным строением. Необходимо исследовать его на наличие временных периодичностей, используя алгоритмы автокорреляционной функции и МНК.

В-третьих, необходимо выяснить пригодность метода наименьших квадратов для прогнозирования временных рядов.

Для контроля качества алгоритма прогноза будем выделять во временном ряде $l$ последовательных значений (контрольную выборку), которые алгоритм будет прогнозировать по предыдущим значениям. В качестве критерия качества прогноза будем минимизировать следующий функционал:

$Q = \sum\limits_{t=1}^{l}|\hat{X_t}-X_t|,$

где $\hat{X_t}$ - прогнозируемое значение в $t$ -ый момент времени, $X_t$ - фактическое значение.

Пути решения задачи

Опишем используемые в работе алгоритмы.

Автокорреляционная функция.

Автокорреляционная функция - это характеристика временного ряда, которая помогает находить его повторяющиеся участки, скрытые из-за наложений шума или других помех.

Для дискретного временного ряда $X_1, X_2, ..., X_n$ с известными матожиданием $\mu$ и дисперсией $\sigma$ автокорреляционную функцию можно рассчитать по следующей формуле:

$R(w)=\frac{1}{(n-w)\sigma ^2}\sum_{t=1}^{n-w} (X_t-\mu)\cdot(X_{t+w}-\mu),$

где $n$ - длина временного ряда, $w$ - текущая задержка во времени. Таким образом получим функцию $R(w)$ , зависящую от лагов (задержек во времени). Исследуя ее на экстремальные значения, получим искомые значения периодов $T=\{\tau_1,..\tau_p\}$ .

Оценка периода осложняется тем, что в некоторой окрестности оцениваемого периода, наблюдаются локальные максимумы коэффициентов корреляции. Следовательно, необходимо усреднение коэффициентов корреляции, а также удаление <<близких>> и кратных периодов (<<близкими>> в работе считаются периоды, отличающиеся друг от друга менее, чем на величину $\delta$ ).

Ряд Фурье.

Сделаем предположение о наличии периодичности в предлагаемом ряде и обратимся к теореме.

Теорема 1.

Если некоторая периодическая функция с периодом $2j$ на интервале $[-j, j]$ удовлетворяет условиям Дирихле (имеет конечное число экстремумов и точек разрыва I рода), то она может быть представлена в виде суммы ряда Фурье (разложена в ряд Фурье).

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

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

$X_t=\frac{a_0}{2}+\sum_{k=1}^{n} {a_k}\cos({\lambda_k}t) + {b_k}\sin({\lambda_k}t),$

где $\lambda_k=2\pi\eta_k$ , $\eta_k=\frac{k}{n}$ , $k=1,2,.., n$ .

Коэффициенты $a_k$ и $b_k$ определяются следующими рядами:

$a_k=\frac{2}{n} \sum_{i=1}^n {X_i}\cos(\lambda_k t), k=0,1,2.., n;$ $b_k=\frac{2}{n} \sum_{i=1}^n {X_i}\sin(\lambda_k t), k=1,2.., n.$

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

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

Тригонометрическая интерполяция методом наименьших квадратов.

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

Введем непрерывную функцию $\varphi(x)$ для аппроксимации дискретной зависимости $g(x_i)$ , $i =1,...,n$ . Будем считать, что $\varphi(x)$ построена при условии наилучшего квадратичного приближения, если:

(1) $M=\sum_{i=1}^n (\varphi(x_i)-g(x_i))^2=\min.$

Рассмотрим случай линейной аппроксимации:

$\varphi(x)=c_{0}\varphi_{0}(x)+c_{1}\varphi_{1}(x)+...+c_{m}\varphi_{m}(x),$

где $\varphi_{0},...,\varphi_{m}$ - произвольные базисные функции, $c_0,...,c_m$ - неизвестные коэффициенты. Количество базисных функций должно быть меньше количества заданных точек для того, чтобы их суперпозиция определялась единственным образом.

Для решения задачи линейной аппроксимации в общем случае следует найти условия минимума суммы квадратов отклонений. Это можно свести к задаче поиска корня системы уравнений $\frac{\partial M}{\partial c_k}=0$ , $k=1,...,m$ . Вычисление данных производных, при учёте равенства (1) приведёт к следующей системе алгебраических уравнений: $\left\{ \begin{array}{lcl} \sum_{i=1}^n (c_{0}\varphi_{0}(x)+c_{1}\varphi_{1})(x)+\ldots+c_{m}\varphi_{m}(x)-g_i)\varphi_{0}(x)=0;\\ &&\\ \sum_{i=1}^n (c_{0}\varphi_{0}(x)+c_{1}\varphi_{1})(x)+\ldots+c_{m}\varphi_{m}(x)-g_i)\varphi_{1}(x)=0; &&\\ \ldots \\ &&\\ \sum_{i=1}^n (c_{0}\varphi_{0}(x)+c_{1}\varphi_{1})(x)+\ldots+c_{m}\varphi_{m}(x)-g_i)\varphi_{m}(x)=0.\\ \end{array} \right.$

Далее следует решить полученную СЛАУ относительно коэффициентов $c_0,...,c_m$ . Для решения СЛАУ обычно составляется расширенная матрица коэффициентов, которую называют матрицей Грамма, элементами которой являются скалярные произведения базисных функций и столбец свободных коэффициентов:

$\begin{pmatrix} (\varphi_0,\varphi_0)&(\varphi_0,\varphi_1)&\ldots&(\varphi_0,\varphi_m)&(\varphi_0,g)\\ (\varphi_1,\varphi_0)&(\varphi_1,\varphi_1)&\ldots&(\varphi_1,\varphi_m)&(\varphi_1,g)\\ \ldots&\ldots&\ddots&\ldots&\ldots\\ (\varphi_m,\varphi_0)&(\varphi_m,\varphi_1)&\ldots&(\varphi_m,\varphi_m)&(\varphi_m,g)\\ \end{pmatrix} </p>$ $(\varphi_j,\varphi_k)=\sum_{i=1}^n \varphi_{j}(x_i)\varphi_{k}(x_i)(\varphi_{j}g)=\sum_{i=1}^n \varphi_{j}(x_i)g(x_i),$

где $j=0,...,m$ , $k=0,...,m$ .

После того как с помощью, например, метода Гаусса найдены коэффициенты $c_0,...,c_m$ , можно построить аппроксимирующую кривую или вычислить координаты заданной точки. Таким образом, задача аппроксимации решена.

Для того чтобы сформировать ортонормированный базис, коэффициенты нормировки будут следующими: для $a_0$ - $(4/n)^{0.5};$ для $a_k$ и $b_k$ - $(2/n)^{0.5},$ где $a_k$ - коэффициенты нормировки для косинусов, $k=0,...,m$ ; $b_k$ - коэффициенты нормировки для синусов, $k=1,...,m$ ; $n$ - количество отсчетов. В качестве функции $g(x)$ рассмотрим дискретный ряд $X_t$ .

При работе с временными рядами необходимо учесть возможное наличие тренда или присутствие постоянного слагаемого. Обе эти составляющие исключим из данных, поскольку они могут привести к большим погрешностям при подсчёте функционала $Q = \sum\limits_{t=1}^{l}|\hat{X_t}-X_t|$ . Пользуясь методом наименьших квадратов мы учтём и тренд и наличие постоянного слагаемого. Необходимо заметить, что тригонометрическая интерполяция основана на разложении в ряд Фурье, поэтому она также не годится для выявления периодической компоненты ряда. В данной работе она используется для нахождения тренда, содержащего периодическую компоненту. Оставшиеся точки исследуются при помощи автокорреляционной функцией.

Вычислительный эксперимент

В данной части работы происходит сравнение алгоритмов автокорреляционной функции и метода наименьших квадратов. Также определяется пригодность МНК для оценки параметров при построении тригонометрической интерполяции ряда.

Исследование модельного зашумленного синуса.

В качестве модельных данных будем использовать функцию $X(t)=0.3\sin(2 \pi t/7)$ . Количество отсчетов $n=100$ . При исследовании временных рядов автокорреляционной функцией для поиска <<близких>> периодов будем считать $\delta=1\%$ .

Исследование величины коэффициента корреляции в зависимости от накладываемого шума.

Была определена зависимость коэффициента корреляции в зависимости от шума (величина шума определялась в процентах от максимального значения функции). Распределение шума - равномерное. Результаты представлены на графике.

Максимальное значение функции $f_{max}=0.2925$ . Максимальный коэффициент корреляции $R_{max}=0.9250$ соответствует 10\%-ому зашумлению. Минимальный коэффициент корреляции $R_{min}=0.6378$ соответствует 100\%-ому зашумлению.

Исследование величины коэффициента корреляции при уменьшении числа отсчетов за период.

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

Максимальный коэффициент корреляции $R_{max}=0.9295$ соответствует 100-та отсчетам за период. Минимальный коэффициент корреляции $R_{min}=0.4355$ соответствует 7-ми отсчетам за период.

Исследование величины коэффициента корреляции при сокращении временного ряда.

Расчет был произведен начиная с длины ряда $n=100$ до $n=10$ . Величина шага составляла 10 временных отсчетов.

Максимальный коэффициент корреляции $R_{max}=0.9295$ соответствует $n=100$ . Минимальный коэффициент корреляции $R_{min}=0.3333$ соответствует $n=10$ , где $n$ - длина ряда.

Выводы:

1. при увеличении шума коэффициент корреляции уменьшается;

2. при уменьшении количества отсчётов за период коэффициент корреляции уменьшается;

3. при уменьшении длины временного ряда коэффициент корреляции уменьшается.

Данные выводы напрямую следуют из устройства корреляционной функции.

Исследование реального временного ряда электрокардиограммы.

Рассмотрим реальный временной ряд со сложной периодической составляющей. Количество отсчетов $n=4170$ .

После применения корреляционной функции к первоначальному ряду получим следующие коэффициенты корреляции:

Корреляция	0.6	0.6	0.4	0.4	0.4	0.3	0.2	0.1	0.1
Период	43	121	453	163	205	288	2371	2089	2811

При анализе таблицы и графика получим, что корреляционная функция определила <<локальный>> период временного ряда, тогда как глобальная периодичность осталась не выявленной. Стоит заметить, что локальный период также определён не точно и коэффициенты корреляции очень малы.

Применим метод наименьших квадратов для получения тригонометрического тренда данного временного ряда. Аппроксимация происходит с помощью двух гармоник.

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

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

Корреляция	0.5	0.2	0.2	0.1	0.1
Период	459	2700	3160	1320	868

Заметим, что хотя коэффициенты корреляции малы, периодичность ряда определена верно. Малость коэффициентов объясняется сложным строением ряда.

Исследование применения МНК для прогнозирования реальных временных рядов.

Метод наименьших квадратов при тригонометрической интерполяции основывается на добавлении новых гармоник для лучшего совпадения реальной функции и её аппроксимации. Критерием схожести служил функционал: $Q = \sum\limits_{t=1}^{l}|\hat{X_t}-X_t|,$

где $\hat{X_t}$ - прогнозируемое значение в $t$ -ый момент времени, $X_t$ - фактическое значение.

В данной работе расчет велся начиная с двух гармоник. Остановка алгоритма происходила в двух случаях:

1. номер гармоники больше половины длины ряда;

2. относительная ошибка на один отсчёт составляет менее 1\%.

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

Длина входного ряда 3000. Горизонт прогноза - 1000 точек. Остановка метода наименьших квадратов произошла на 148 гармонике. Функционал качества составил $Q=163.6817$ на 1000 временных отсчётов, что составляет $16.3\%$ . Результат прогнозирования ряда приведён на графике.

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

Заключение

В работе рассмотрена зависимость коэффициента корреляции от различных входных параметров. Исследованы способы получения периодичности временных рядов, а также метод нахождения скрытых периодичностей. Рассмотрена возможность применения метода наименьших квадратов при оценке параметров для тригонометрической интерполяции и прогнозирования реального временного ряда. Анализ проводился на модельных данных и на реальном временном ряде электрокардиограммы.Необходимый для повторения вычислительного эксперимента код можно найти по [1].

См. также

Литература

Gujarati, D. Basic Econometrics, 4th ed / D. Gujarati. - The McGraw-Hill Companies, 2004.

А.М. Тер-Крикоров, М.И. Шабунин. Курс математического анализа. Ed. by . Ходан. - ФИЗМАТ-ЛИТ, 2001.

В.В. Стрижов, Г.О. Пташко. Алгоритмы поиска суперпозиций при выборе оптимальных регрессионных моделей. - Вычислительный центр им. А.А. Дородницына РАН, 2007.

В.В. Витязев. Спектрально-корреляционный анализ равномерных временных рядов. - Издательство С.-Петербургского университета, 2001.

А.И. Орлов. Прикладная статистика. - Издательство "Экзамен", 2004.

Ю.В. Попов. О выделении периодической компоненты из временного ряда показателя количества катастроф.//Проблемы безопасности полетов. - 2008.

В.В. Стрижов. Методы индуктивного порождения регрессионных моделей. - Вычислительный центр им. А.А. Дородницына Российской Академии наук, 2008. - P. 37.

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%92%D1%8B%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D0%BF%D0%B5%D1%80%D0%B8%D0%BE%D0%B4%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%BE%D0%B9_%D0%BA%D0%BE%D0%BC%D0%BF%D0%BE%D0%BD%D0%B5%D0%BD%D1%82%D1%8B_%D0%B2%D1%80%D0%B5%D0%BC%D0%B5%D0%BD%D0%BD%D0%BE%D0%B3%D0%BE_%D1%80%D1%8F%D0%B4%D0%B0_%28%D0%BF%D1%80%D0%B8%D0%BC%D0%B5%D1%80%29»

@@ Строка 126: / Строка 126: @@
 [[Изображение:PeriodicComponents1.png‎|Модельный временной ряд X(t)=0.3*sin(2*pi*t/7)]]
-==Исследование величины коэффициента корреляции в зависимости от накладываемого шума.==
+===Исследование величины коэффициента корреляции в зависимости от накладываемого шума.===
 Была определена зависимость коэффициента корреляции в зависимости от шума (величина шума определялась в процентах от максимального значения функции). Распределение шума - равномерное. Результаты представлены на графике.