Книги, научные публикации Pages: | 1 | ... | 4 | 5 | 6 | 7 | 8 | ... | 10 |

Оглавление Введение................................... 11 I Введение в социально-экономическую статистику 15 1. Основные понятия 17 1.1. Краткая историческая справка.. ...

-- [ Страница 6 ] --

Ряд d-х разностей (разностей d-го порядка) получается как степень разност ного оператора, то есть применением разностного оператора d раз.

При d = 2 получается 2 = (1 - L)2 = 1 - 2L + L2, поэтому 2xt = =(1 - 2L + L2)xt = xt - 2xt-1 + xt-2.

Для произвольного порядка d следует использовать формулу бинома Ньютона:

d d!

k k d =(1 - L)d = (-1)kCd Lk, гд е Cd =, k!(d - k)!

k= d k так что dxt =(1 - L)dxt = (-1)kCd xt-k.

k= 11.8. Модели регрессии с распределенным лагом Часто при моделировании экономических процессов на изучаемую переменную xt влияют не только текущие значения объясняющего фактора zt, но и его ла ги. Типичным примером являются капиталовложения: они всегда дают результат с некоторым лагом.

Модель распределенного лага можно записать следующим образом:

q xt = + jzt-j + t = + (L)zt + t. (11.18) j= 376 Глава 11. Основные понятия в анализе временных рядов q j где q Ч величина наибольшего лага, (B) = jL Ч лаговый многочлен, j=0 j t Ч случайное возмущение, ошибка. Коэффициенты j задают структуру лага q и называются весами. Конструкцию jzt-j часто называют скользящим j= средним переменной zt 10.

Рассмотрим практические проблемы получения оценок коэффициентов j в модели (11.18). Модель распределенного лага можно оценивать обычным ме тодом наименьших квадратов, если выполнены стандартные предположения ре грессионного анализа. В частности, количество лагов не должно быть слишком большим, чтобы количество регрессоров не превышало количество наблюдений, и все лаги переменной zt, т.е. zt-j (j =0,..., q), не должны быть коррелированы с ошибкой t.

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

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

1) Для каждого конкретного q оценивается модель (11.18), и из нее берется t-статистика для последнего коэффициента, т.е. q. Эти t-статистики рассматри ваются в обратном порядке, начиная с q = Q (и заканчивая q =0). Как только t-статистика оказывается значимой при некотором наперед заданном уровне, то следует остановиться и выбрать соответствующую величину q.

2) Следует оценить модель (11.18) при q = Q. Из этой регрессии берутся F -статистики для проверки нулевой гипотезы о том, что коэффициенты при по следних Q - q +1 лагах, т.е. q,..., Q, одновременно равны нулю:

H0 : j =0, j = q,..., Q.

Соответствующие F -статистики рассчитываются по формулам:

(RSSQ - RSSq-1)/(Q - q +1) Fq =, RSSQ/(T - Q - 2) где RSSrЧ сумма квадратов остатков из модели распределенного лага при q = r, T Ч количество наблюдений. При этом при проведении расчетов для сопостави мости во всех моделях надо использовать одни и те же наблюдения Ч те, которые использовались при q = Q (следовательно, при всех q используется одно и то же T ). Эти F -статистики рассматриваются в обратном порядке от q = Q до q = (в последнем случае в модели переменная z отсутствует). Как только F -статистика Другое часто используемое название Ч линейный фильтр.

11.9. Условные распределения оказывается значимой при некотором наперед заданном уровне, то следует оста новиться и выбрать соответствующую величину q.

3) Для всех q от q =0 до q = Q рассчитывается величина информационно го критерия, а затем выбирается модель с наименьшим значением этого инфор мационного критерия. Приведем наиболее часто используемые информационные критерии.

Информационный критерий Акаике:

RSS 2(n +1) AIC =ln( ) +, T T где RSS сумма квадратов остатков в модели, T Ч фактически использовавшееся количество наблюдений, n Ч количество факторов в регрессии (не считая кон станту). В рассматриваемом случае n = q +1, а T = T0 - q, гд е T0 Ч количество наблюдений при q =0.

Байесовский информационный критерий (информационный критерий Шварца):

RSS (n +1) lnT BIC =ln( ) +.

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

11.9. Условные распределения Условные распределения играют важную роль в анализе временных рядов, осо бенно при прогнозировании. Мы не будем вдаваться в теорию условных распре делений, это предмет теории вероятностей (определения и свойства условных рас пределений см. в Приложении A.3.1). Здесь мы рассмотрим лишь основные пра вила, по которым можно проводить преобразования. При этом будем использовать следующее стандартное обозначение: если речь идет о распределении случайной величины X, условном по случайной величине Y (условном относительно Y ), то это записывается в виде X|Y.

Основное правило работы с условными распределениями, которое следует за помнить, состоит в том, что если рассматривается распределение, условное отно сительно случайной величины Y, тос Y и ее функциями следует поступать так же, как с детерминированными величинами. Например, для условных математических ожиданий и дисперсий выполняется E ((Y ) +(Y )X|Y ) =(Y ) +(Y )E(X|Y ), var ((Y ) +(Y )X|Y ) =2(Y )var(X|Y ).

378 Глава 11. Основные понятия в анализе временных рядов Как и обычное безусловное математическое ожидание, условное ожидание представляет собой линейный оператор. В частности, ожидание суммы есть сумма ожиданий:

E (X1 + X2|Y ) =E(X1|Y ) +E(X2|Y ).

Условное математическое ожидание E(X|Y ) в общем случае не является де терминированной величиной, т.е. оно является случайной величиной, которая мо жет иметь свое математическое ожидание, характеризоваться положительной дис персией и т.п.

Если от условного математического ожидания случайной величины X еще раз взять обычное (безусловное) математическое ожидание, то получится обычное (безусловное) математическое ожидание случайной величины X. Таким образом, действует следующее правило повторного взятия ожидания:

E (E(X|Y )) = E(X).

В более общей форме это правило имеет следующий вид:

E (E(X|Y, Z)|Y ) =E(X|Y ), что позволяет применять его и тогда, когда второй раз ожидание берется не полно стью, т.е. не безусловное, а лишь условное относительно информации, являющейся частью информации, относительно которой ожидание бралось первый раз.

Если случайные величины X и Y статистически независимы, то распределе ние X, условное по Y, совпадает с безусловным распределением X. След ова тельно, для независимых случайных величин X и Y выполнено, в частности, E(X|Y ) =E(X), var(X|Y ) =var(X).

11.10. Оптимальное в среднеквадратическом смысле прогнозирование: общая теория 11.10.1. Условное математическое ожидание как оптимальный прогноз Докажем в абстрактном виде, безотносительно к моделям временных рядов, общее свойство условного математического ожидания, заключающееся в том, что оно минимизирует средний квадрат ошибки прогноза.

Предположим, что строится прогноз некоторой случайной величины x на ос нове другой случайной величины, z, и что точность прогноза при этом оценивается 11.10 Оптимальное в среднеквадратическом смысле прогнозирование на основе среднего квадрата ошибки прогноза = x - xp(z), гд е xp(z) Чпро гнозная функция. Таким образом, требуется получить прогноз, который бы мини мизировал E 2 = E (x - xp(z))2.

Оказывается, что наилучший в указанном смысле прогноз дает математическое ожидание x, условное относительно z, т.е. E (x|z), которое мы будем обозначать x(z). Докажем это. Возьмем произвольный прогноз xp(z) и представим ошибку прогноза в виде:

x - xp(z) = =(x - x(z)) + (x(z) - xp(z)).

Найдем сначала математическое ожидание квадрата ошибки, условное относи тельно z:

E 2|z = E (x - x(z))2|z + +2E [(x - x(z))(x(z) - xp(z))|z] +E ( - xp(z)t)2|z.

x(z) При взятии условного математического ожидания с функциями z можно обра щаться как с константами. Поэтому x(z) E ( - xp(z))2|z =(x(z) - xp(z)) и E [(x - x(z))(x(z) - xp(z))|z] =E (x - x(z)|z)(x(z) - xp(z)) = =(x(z) - x(z))(x(z) - xp(z)) = 0.

Используя эти соотношения, получим E 2|z = E (x - x(z))2 |z +(x(z) - xp(z))2.

Если теперь взять от обеих частей безусловное математическое ожидание, то (по правилу повторного взятия ожидания) получится E 2 = E (x - x(z))2 + E ( - xp(z))2.

x(z) Поскольку второе слагаемое неотрицательно, то E (x - xp(z))2 = E 2 E (x - x(z))2.

380 Глава 11. Основные понятия в анализе временных рядов Другими словами, средний квадрат ошибки прогноза достигает минимума при xp(z) = =E (x|z).

x(z) Оптимальный прогноз xp(z) = =E (x|z) является несмещенным. Дей x(z) ствительно, по правилу повторного взятия ожидания E (E (x|z)) = E (x).

Поэтому E = E (x - xp(z)) = E (x) - E (E (x|z)) = 0.

11.10.2. Оптимальное линейное прогнозирование Получим теперь формулу для оптимального (в смысле минимума среднего квад рата ошибки) линейного прогноза. Пусть случайная переменная z, на основе кото рой делается прогноз x, представляет собой n-мерный вектор: z =(z1,..., zn).

Без потери общности можно предположить, что x и z имеют нулевое математи ческое ожидание. Будем искать прогноз x в виде линейной комбинации zj:

xp(z) =1z1 +... + nzn = z, где =(1,..., n) Ч вектор коэффициентов. Любой прогноз такого вида яв ляется несмещенным, поскольку, как мы предположили, Ex =0 и Ez =0.

Требуется решить задачу минимизации среднего квадрата ошибки (в данном случае это эквивалентно минимизации дисперсии ошибки):

E (x - xp(z))2 min!

Средний квадрат ошибки можно представить в следующем виде:

E (x - xp(z))2 = E x2 - 2 zx + zz = x - 2 Mzx + Mzz, где x = Ex2 Ч дисперсия x, Mzx = E [zx] Ч вектор, состоящий из ковариаций zj и x, а Mzz = E [zz ] Ч ковариационная матрица z. (Напомним, что мы рас сматриваем процессы с нулевым математическим ожиданием.) Дифференцируя по, получим следующие нормальные уравнения:

-2Mzx +2Mzz =0, откуда - = Mzz Mzx.

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

Таким образом, оптимальный линейный прогноз имеет вид:

- xp(z) =z Mzz Mzx. (11.19) Ошибка оптимального линейного прогноза равна - = x - xp(z) =x - z Mzz Mzx.

Эта ошибка некоррелирована с z, то есть с теми переменными, по которым делается прогноз. Действительно, умножая на z и беря математическое ожидание, получим -1 - E (z) =E zx - zz Mzz Mzx = Mzx - MzzMzz Mzx, т.е.

E (z) =0.

Средний квадрат ошибки оптимального прогноза равен 2 -1 -1 - E 2 = E (x - xp(z))2 = x - 2MxzMzz Mzx + MxzMzz MzzMzz Mzx.

После преобразований получаем 2 - E 2 = x - MxzMzz Mzx. (11.20) Несложно увидеть аналогии между приведенными формулами и формулами МНК. Таким образом, данные рассуждения можно считать одним из возможных теоретических обоснований линейного МНК.

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

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

E 2 = x - Mxz.

382 Глава 11. Основные понятия в анализе временных рядов В общем случае оптимальный линейный прогноз (11.19) не совпадает с услов ным математическим ожиданием E (x|z). Другими словами, он не является опти мальным среди всех возможных прогнозов. Пусть, например, z имеет стандартное нормальное распределение: z N(0, 1), а x связан с z формулой x = z2 - 1.

Тогда, поскольку x и z некоррелированы, то = 0, и оптимальный линейный прогноз имеет вид xp(z) = 0 при среднем квадрате ошибки прогноза равном E (z2 - 1)2 =2. В то же время прогноз по нелинейной формуле xp(z) =z2 - будет безошибочным (средний квадрат ошибки прогноза равен 0).

В частном случае, когда совместное распределение x и z является многомер ным нормальным распределением:

x 0 x Mxz N,, z 0n Mzx Mzz оптимальный линейный прогноз является просто оптимальным. Это связано с тем, что по свойствам многомерного нормального распределения (см. Приложение A.3.2) условное распределение x относительно z будет иметь следующий вид:

-1 2 - x|z N z Mzz Mzx, x - MxzMzz Mzx.

- Таким образом, E (x|z) = z Mzz Mzx, что совпадает с формулой оптимального линейного прогноза (11.19).

11.10.3. Линейное прогнозирование стационарного временного ряда Пусть xt Ч слабо стационарный процесс с нулевым математическим ожида нием. Рассмотрим проблему построения оптимального линейного прогноза этого процесса, если в момент t известны значения ряда, начиная с момента 1, т.е. толь ко конечный ряд x =(x1,..., xt). Предположим, что делается прогноз на шагов вперед, т.е. прогноз величины xt+. Для получения оптимального линейного (по x) прогноза можно воспользоваться формулой (11.19). В случае стационарного временного ряда ее можно переписать в виде:

xt() =x -1t,, (11.21) t 11.10 Оптимальное в среднеквадратическом смысле прогнозирование где 0 1 t- 1 0 t- t =...

....

...

t-1 t-2 Ч автоковариационная матрица ряда (x1,..., xt), а вектор t, составлен из ковариаций xt+ с (x1,..., xt), т.е.

t, =(t+-1,..., ).

Можно заметить, что автоковариации здесь нужно знать только с точностью до множителя. Например, их можно заменить автокорреляциями.

Рассмотрим теперь прогнозирование на один шаг вперед. Обозначим через t вектор, составленный из ковариаций xt+1 с (x1,..., xt), т.е. t =(t,..., 1) = = t,1. Прогноз задается формулой:

t xt(1) = x -1t = x t = txt-i.

t i i= Прогноз по этой формуле можно построить только если матрица t неособенная.

Коэффициенты t, минимизирующие средний квадрат ошибки прогноза, задаются i нормальными уравнениями t = t или, в развернутом виде, t t|k-i| = k, k =1,..., t.

i i= Ошибка прогноза равна = xt+1 - xt(1).

Применив (11.20), получим, что средний квадрат этой ошибки равен E 2 = 0 - t -1t.

t Заметим, что 0 - t -1t = |t+1| / |t|, т.е предыдущую формулу можно пере t писать как E 2 = |t+1| / |t|. (11.22) 384 Глава 11. Основные понятия в анализе временных рядов Действительно, матрицу t+1 можно представить в следующей блочной форме:

t t t+1 =.

t По правилу вычисления определителя блочной матрицы имеем:

|t+1| = 0 - t -1t |t|.

t Если |t+1| =0, т.е. если матрица t+1 вырождена, то средний квадрат ошибки прогноза окажется равным нулю, т.е. оптимальный линейный прогноз будет без ошибочным. Процесс, для которого существует такой безошибочный линейный прогноз, называют линейно детерминированным.

Укажем без доказательства следующее свойство автоковариационных матриц:

если матрица t является вырожденной, то матрица t+1 также будет вырожден ной.

Отсюда следует, что на основе конечного отрезка стационарного ряда (x1,..., xt) можно сделать безошибочный линейный прогноз на один шаг впе редв том и только в том случае, если автоковариационная матрица t+1 является вырожденной ( |t+1| =0).

Действительно, пусть существует безошибочный линейный прогноз. Возможны два случая: |t| = 0 и |t| = 0. Если |t| = 0, то средний квадрат ошибки прогноза равен |t+1| / |t|, откуд а |t+1| =0, еслиже |t| =0, то из этого также следует |t+1| =0.

Наоборот, если |t+1| =0, тонайд ется s ( s t)такое, что |s+1| =0, но |s| =0.

Тогда можно сделать безошибочный прогноз для xt+1 на основе (x1+t-s,..., xt), а, значит, и на основе (x1,..., xt).

При использовании приведенных формул на практике возникает трудность, связанная с тем, что обычно теоретические автоковариации k неизвестны. Тре буется каким-то образом получить оценки автоковариаций. Обычные выборочные автоковариации ck здесь не подойдут, поскольку при больших k (сопоставимых с длиной ряда) они являются очень неточными оценками k. Можно предложить следующий подход11:

1) Взять за основу некоторую параметрическую модель временного ряда. Пусть Ч соответствующий вектор параметров. Рассчитать теоретические автоковари ации для данной модели в зависимости от параметров: k = k().

Этот подход, в частности, годится для стационарных процессов ARMA. В пункте 14.8 дается альтернативный способ прогнозирования в рамках модели ARMA.

11.10 Оптимальное в среднеквадратическом смысле прогнозирование 2) Оценить параметры на основе имеющихся данных. Пусть b Ч соответству ющие оценки.

3) Получить оценки автоковариаций, подставив b в формулы теоретических автоковариаций: k k(b).

4) Использовать для прогнозирования формулу (11.21), заменяя теоретические автоковариации полученными оценками автоковариаций.

11.10.4. Прогнозирование по полной предыстории.

Разложение Вольда Можно распространить представленную выше теорию на прогнозирование ря д а в случае, когд а в момент t известна полная предыстория t =(xt, xt-1,... ).

Можно определить соответствующий прогноз как предел прогнозов, полученных на основе конечных рядов (xt, xt-1,..., xt-j), j = t, t - 1,..., -. Без доказа тельства отметим, что этот прогноз будет оптимальным в среднеквадратическом смысле. Если рассматривается процесс, для которого |t| =0 t, топоаналогиис (11.22) средний квадрат ошибки такого прогноза равен |t+1| E 2 = lim.

t |t| |t+1| |t| Заметим, что всегда выполнено 0 <, т.е. средний квадрат ошиб |t| |t-1| ки не увеличивается с увеличением длины ряда, на основе которого делается про гноз, и ограничен снизу нулем, поэтому указанный предел существует всегда.

Существуют процессы, для которых |t| = 0 t, т.е. для них нельзя сделать безошибочный прогноз, имея только конечный отрезок ряда, однако |t+1| lim =0.

t |t| Такой процесс по аналогии можно назвать линейно детерминированным. Его фак тически можно безошибочно предсказать, если имеется полная предыстория про цесса t =(xt, xt-1,... ).

Если же данный предел положителен, то линейный прогноз связан с ошибкой:

E 2 > 0. Такой процесс можно назвать регулярным.

Выполнены следующие свойства стационарных рядов.

A. Пусть xt Ч слабо стационарный временной ряд, и пусть t Чошибкиод ношагового оптимального линейного прогноза по полной предыстории процесса (xt-1, xt-2,...). Тогд аошибки t являются белым шумом, т.е. имеют нулевое ма 386 Глава 11. Основные понятия в анализе временных рядов тематическое ожидание, не автокоррелированы и имеют одинаковую дисперсию:

E (t) =0, t, E (st) =0, при s = t, E 2 =, t.

B. Пусть, кроме того, xt является регулярным, т.е. E 2 = > 0. Тогд а он представим в следующем виде:

xt = it-i + vt, (11.23) i= где 0 = 1, i < ;

процесс vt здесь является слабо стационарным, i= линейно детерминированным и не коррелирован с ошибками t: E (svt) =0 при s, t. Такое представление единственно.

Утверждения A и B составляют теорему Вольда. Эта теорема является одним из самых фундаментальных результатов в теории временных рядов. Утверждение B говорит о том, что любой стационарный процесс можно представить в виде так называемого линейного фильтра от белого шума12 плюс линейно детерминиро ванная компонента. Это так называемое разложение Вольда.

Доказательство теоремы Вольда достаточно громоздко. Мы не делаем попытки его излагать и даже обсуждать;

отсылаем заинтересованных читателей к гл. 7 книги Т. Андерсона [2].

Заметим, что коэффициенты разложения i удовлетворяют соотношению E (t-ixt) E (t-ixt) i = =.

E t-i Для того чтобы это показать, достаточно умножить (11.23) на t-i и взять мате матическое ожидание от обеих частей.

Разложение Вольда имеет в своей основе прогнозирование на один шаг впе ред. С другой стороны, если мы знаем разложение Вольда для процесса, то с по мощью него можно делать прогнозы. Предположим, что в момент t делается прогноз на шагов вперед, т.е. прогноз величины xt+ на основе предыстории t =(xt, xt-1,... ). Сдвинем формулу разложения Вольда (11.23) на периодов вперед:

xt+ = it+ -i + t+.

i= Можно назвать первое слагаемое в 11.23 также процессом скользящего среднего бесконечного порядка MA(). Процессы скользящего среднего обсуждаются в пункте 14.4.

11.10 Оптимальное в среднеквадратическом смысле прогнозирование Второе слагаемое, vt+, можно предсказать без ошибки, зная t. Из первой суммы первые слагаемых не предсказуемы на основе t. При прогнозирова нии их можно заменить ожидаемыми значениями Ч нулями. Из этих рассуждений следует следующая формула прогноза:

xt() = it+-i + t+.

i= Без доказательства отметим, что xt() является оптимальным линейным прогно зом. Ошибка прогноза при этом будет равна - it+-i.

i= Поскольку t Ч белый шум с дисперсией, то средний квадрат ошибки про гноза равен - 2 i.

i= Напоследок обсудим природу компоненты vt. Простейший пример линейно детерминированного ряда Ч это, говоря неформально, случайная константа:

vt =, где Ч наперед заданная случайная величина, E =0.

Типичный случай линейно детерминированного ряда Ч это, говоря неформаль но, случайная синусоида:

vt = cos(t + ), где Ч фиксированная частота, и Ч независимые случайные величины, причем имеет равномерное распределение на отрезке [0;

2].

Это два примера случайных слабо стационарных рядов, которые можно без ошибочно предсказывать на основе предыстории. Первый процесс можно модели ровать с помощью константы, а второй Ч с помощью линейной комбинации синуса и косинуса:

cos(t) + sin(t).

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

388 Глава 11. Основные понятия в анализе временных рядов 11.11. Упражнения и задачи Упражнение 1.1. Дан временной ряд x =(5, 1, 1, -3, 2, 9, 6, 2, 5, 2).

Вычислите среднее, дисперсию (смещенную), автоковариационную и авто корреляционную матрицы.

1.2. Для временного ряда y =(6, 6, 1, 6, 0, 6, 6, 4, 3, 2) повторите упраж нение 1.1.

1.3. Вычислите кросс-ковариации и кросс-корреляции для рядов x и y из преды дущих упражнений для сдвигов -9,..., 0,..., 9.

1.4. Для временного ряда x =(7, -9, 10, -2, 21, 13, 40, 36, 67, 67) оцените параметры полиномиального тренда второго порядка. Постройте точечный и интервальный прогнозы по тренду на 2 шага вперед.

1.5. Сгенерируйте 20 рядов, задаваемых полиномиальным трендом третьего по рядка t =5+4t-0.07t2 +0.0005t3 длиной 100 наблюдений, с добавлением белого шума с нормальным распределением и дисперсией 20.

Допустим, истинные значения параметров тренда неизвестны.

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

б) Для 20 рядов оцените полиномиальный тренд третьего порядка по пер вым 50 наблюдениям. Вычислите оценки параметров тренда и их ошиб ки. Сравните оценки с истинными значениями параметров.

в) Проведите те же вычисления, что и в пункте (б), для 20 рядов, используя 100 наблюдений. Результаты сравните.

г) Используя предшествующие расчеты, найдите точечные и интерваль ные прогнозы на три шага вперед с уровнем доверия 95%.

1.6. Найдите данные о динамике денежного агрегата M0 в России за 10 последо вательных лет и оцените параметры экспоненциального тренда.

1.7. Ряд x =(0.02, 0.05, 0.06, 0.13, 0.15, 0.2, 0.31, 0.46, 0.58, 0.69, 0.78, 0.81, 0.95, 0.97, 0.98) характеризует долю семей, имеющих телевизор. Оцените параметры логиcтического тренда.

1.8. По ряду x из упражнения 3 рассчитайте ранговый коэффициент корреляции Спирмена и сделайте вывод о наличии тенденции.

11.11 Упражнения и задачи Таблица 11. Расходы на рекламу 10 100 50 200 20 70 100 50 300 Объем продаж 1011 1030 1193 1149 1398 1148 1141 1223 1151 1.9. Дан ряд:

x =(10, 9, 12, 11, 14, 12, 17, 14, 19, 16, 18, 21, 20, 23, 22, 26, 23, 28, 25, 30).

а) Оцените модель линейного тренда. Остатки, полученные после исклю чения тренда, проверьте на стационарность с использованием рангового коэффициента корреляции Спирмена.

б) Рассчитайте для остатков статистику Бартлетта, разбив ряд на 4 интер вала по 5 наблюдений. Проверьте однородность выборки по дисперсии.

в) Рассчитайте для остатков статистику ГолдфельдаЧКвандта, исключив 6 наблюдений из середины ряда. Проверьте однородность выборки по дисперсии. Сравните с выводами, полученными на основе критерия Бартлетта.

1.10. По данным таблицы 11.1 оцените модель распределенного лага зависимости объема продаж от расходов на рекламу с лагом 2.

Определите величину максимального лага в модели распределенного лага, используя различные критерии ( t-статистики, F -статистики, информацион ные критерии).

Задачи 1. Перечислить статистики, использующиеся в расчете коэффициента автокор реляции, и записать их формулы.

2. Чем различается расчет коэффициента автокорреляции для стационарных и нестационарных процессов? Записать формулы.

3. Вычислить значение коэффициента корреляции для двух рядов:

x =(1, 2, 3, 4, 5, 6, 7,... ) и y =(2, 4, 6, 8, 10, 12, 14,... ).

4. Посчитать коэффициент автокорреляции первого порядка для ряда x =(2, 4, 6, 8).

5. Есть ли разница между автокорреляционной функцией и трендом автокорре ляции?

390 Глава 11. Основные понятия в анализе временных рядов 6. Записать уравнения экспоненциального и полиномиального трендов и при вести формулы для оценивания их параметров.

7. Записать формулу для оценки темпа прироста экспоненциального тренда.

8. Привести формулу логистической кривой и указать особенности оценивания ее параметров.

9. Оценить параметры линейного тренда для временного ряда x =(1, 2, 5, 6) и записать формулу доверительного интервала для прогноза на 1 шаг вперед.

10. Дан временной ряд: x =(1, 0.5, 2, 5, 1.5). Проверить его на наличие тренда среднего.

11. Пусть L Ч лаговый оператор. Представьте в виде степенного ряда следую щие выражения:

2 -1, 5 2.8 - а) ;

б) ;

в) ;

г).

1 - 0.8L 1 - 0.9L 1+0.4L 1+0.5L Рекомендуемая литература 1. Айвазян С.А. Основы эконометрики. Т.2. Ч М.: Юнити, 2001.

2. Андерсон Т. Статистический анализ временных рядов. Ч М.: Мир, 1976.

(Гл. 1, 3, 7).

3. Бокс Дж., Дженкинс Г. Анализ временных рядов. Прогноз и управление.

Вып. 1. Ч М.: Мир, 1974. (Гл. 1).

4. Кендалл М. Дж., Стьюарт А. Многомерный статистический анализ и вре менные ряды. Ч М.: Наука, 1976. (Гл. 45Ц47).

Маленво Э. Статистические методы эконометрии. Вып. 2. Ч М.: Статисти ка, 1976. (Гл. 12).

5. Магнус Я.Р., Катышев П.К., Пересецкий А.А. Эконометрика Ч начальный курс. Ч М.: Дело, 2000. (Гл. 12).

6. Enders Walter. Applied Econometric Time Series. Ч Iowa State University, 1995. (Ch. 1).

7. Mills Terence C. Time Series Techniques for Economists, Cambridge University Press, 1990 (Ch. 5).

8. Wooldridge Jeffrey M. Introductory Econometrics: A Modern Approach, 2nd ed., Thomson, 2003 (Ch. 10).

Глава Сглаживание временного ряда 12.1. Метод скользящих средних Одним из альтернативных по отношению к функциональному описанию тренда вариантов сглаживания временного ряда является метод скользящих или, как еще говорят, подвижных средних.

Суть метода заключается в замене исходного временного ряда последователь ностью средних, вычисляемых на отрезке, который перемещается вдоль времен ного ряда, как бы скользит по нему. Задается длина отрезка скольжения (2m +1) по временной оси, т.е. берется нечетное число наблюдений. Подбирается полином p t = aktk (12.1) k= к группе первых (2m +1) членов ряда, и этот полином используется для опре деления значения тренда в средней (m +1)-й точке группы. Затем производится сдвиг на один уровень ряда вперед и подбирается полином того же порядка к группе точек, состоящей из 2-го, 3-го,..., (2m +2)-го наблюдения. Находится значение тренда в (m +2)-й точке и т.д. тем же способом вдоль всего ряда до последней группы из (2m +1) наблюдения. В действительности нет необходимости строить полином для каждого отрезка. Как будет показано, эта процедура эквивалентна нахождению линейной комбинации уровней временного ряда с коэффициентами, 392 Глава 12. Сглаживание временного ряда которые могут быть определены раз и навсегда и зависят только от длины отрезка скольжения и степени полинома.

Для определения коэффициентов a0, a1,..., ap полинома (12.1) с помощью МНК по первым (2m +1) точкам минимизируется функционал:

m = (xt - a0 - a1t -... - aptp)2. (12.2) t=-m Заметим, что t принимает условные значения от -m до m. Это весьма удоб ный прием, существенно упрощающий расчеты. Дифференцирование функционала по a0, a1,..., ap дает систему из p +1 уравнения типа:

m m m m m a0 tj + a1 tj+1 + a2 tj+2 + + ap tj+p = xttj, t=-m t=-m t=-m t=-m t=-m j =0, 1,..., p. (12.3) Решение этой системы уравнений относительно неизвестных параметров m a0, a1,..., ap (i =0, 1,..., 2p) облегчается тем, что все суммы ti при t=-m нечетных i равны нулю. Кроме того, т. к. полином, подобранный по 2m +1 точ кам, используется для определения значения тренда в средней точке, а в этой точке t =0, то, положив в уравнении (12.1) t =0, получаем значение тренда, равное a0.

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

Система нормальных уравнений (12.3), которую нужно разрешить относитель но a0, разбивается на две подсистемы: одну Ч содержащую коэффициенты с чет ными индексами a0, a2, a4,..., другую Ч включающую коэффициенты с нечет ными индексами a1, a3, a5,.... Решение системы относительно a0 зависит от чис m m ленных значений ti и линейных функций от x типа xttj.

t=-m t=-m В итоге, значением тренда в центральной точке отрезка будет средняя ариф метическая, взвешенная из значений временного ряда от x-m до xm c весовыми коэффициентами t, которые зависят от значений m и p:

m a0 = txt.

t=-m Указанная формула применяется для всех последующих отрезков скольжения, с вы числением значений тренда в их средних точках.

Продемонстрируем рассматриваемый метод на примере полинома второй степени и длины отрезка скольжения, равной пяти точкам. Здесь надо свести к минимуму сумму:

= (xt - a0 - a1t - a2t2)2.

t=- 12.1. Метод скользящих средних Получается система уравнений:

2 2 2 a0 + a1 t + a2 t2 = xt, t=-2 t=-2 t=-2 t=- 2 2 2 a0 t + a1 t2 + a2 t3 = xtt, t=-2 t=-2 t=-2 t=- 2 2 2 a0 t2 + a1 t3 + a2 t4 = xtt2.

t=-2 t=-2 t=-2 t=- Для конкретных значений сумм при ap система уравнений приобретает вид:

5a0 +10a2 = xt, t=- 10a1 = xtt, t=- 10a0 +34a2 = xtt2.

t=- Решение этой системы относительно a0 дает следующий результат:

2 1 a0 = 17 xt - 5 xtt2 = (-3x-2 +12x-1 +17x0 +12x1 - 3x2).

35 t=-2 t=- Весовые коэффициенты для полиномов 2Ц5 степени и длины отрезка скольжения от 5 до 9 представлены в таблице 12.1.

Таблица 12.1. Фрагмент таблицы Каудена для весов t Длина Степени полинома отрезка скольже ния 2m +1 m p =2, p =3 p =4, p = 5 2 (-3, 12, 17, 12, -3) 7 3 (-2, 3, 6, 7, 6, 3, -2) (5, -30, 75, 131, 75, -30, 5) 9 4 (-21, 14, 39, 54, 59, 54, (15, -55, 30, 135, 179, 135, 231 30, -55, 15) 39, 14, -21) 394 Глава 12. Сглаживание временного ряда Метод скользящих средних в матричной форме Введем следующие обозначения:

m 1. cj = xttj.

t=-m Так как xt и tj известны, то cj также известно для каждого j =0, 1,..., p.

m 2. i = ti, i =0, 1,..., 2p. Тогда t=-m 0, если i Ч нечетно, 2m + i =, если i =0, 1i +2i +... + mi, если i Ч четно.

В таких обозначениях система (12.3) принимает вид:

0 1 p a0 c 1 2 p+1 a1 c =.

.....

......

.....

p p+1 2p ap cp В краткой записи эта система выглядит как Ma = c, где матрица M Ч известна, кроме того, ее элементы с нечетными индексами рав ны нулю, вектор c также известен.

Из полученной системы следует a = M-1c.

Теперь можно использовать формулу Крамера для нахождения ak:

det Mk+ ak =, det M 12.1. Метод скользящих средних где матрица Mk+1 получается из матрицы M заменой (k +1)-го столбца векто ром c.

Таким образом, det M1 det M2 det Mp+ a =,,...,.

det M det M det M Рассмотрим частный случай, когда m =2 и p =2, т.е. временной ряд аппроксими руется полиномом второй степени:

t = a0 + a1t + a2t2.

Система уравнений, которую нужно решить относительно ak, имеет вид:

2 2 2 a0 tj + a1 tj+1 + a2 tj+2 = xttj, t=-2 t=-2 t=-2 t=- где x-2, x-1, x0, x1, x2 Чизвестны, j =0,..., p. Находим i:

0, если i Ч нечетно, i =, если i =0, 1i +2i, если i Ччетно.

Тогда 0 1 2 5/2 0 M = =.

1 2 3 0 5 2 3 4 5 0 Находим определители:

25 det M = c0 0 5 2 det M1 = =5(17c0 - 5c2) = 17 xt - 5 xtt2, c1 5 t=-2 t=- c2 0 5/2 c0 35 det M2 = = c1 = xtt, 0 c1 2 t=- 5 c2 396 Глава 12. Сглаживание временного ряда 5/2 0 c0 2 c 25 det M3 = =25 - c0 = xtt2 - xt.

0 5 c 2 2 t=-2 t=- 5 0 c Отсюда:

2 det M1 a0 = = 17 xt - 5 xtt2, det M t=-2 t=- det M2 a1 = = xtt, det M t=- 2 det M3 1 a2 = = xtt2 - xt.

det M 14 t=-2 t=- Таким образом, 3 12 17 12 a0 = - x-2 + x-1 + x0 + x1 - x2, 35 35 35 35 a1 = -0, 2x-2 - 0, 1x-1 +0, 1x1 +0, 2x2, 1 1 1 1 a2 = x-2 - x-1 - x0 - x1 + x2, 7 14 7 14 и каждый из этих коэффициентов получается как взвешенная средняя из уровней временного ряда, входящих в отрезок.

Оценки параметров a1, a2,..., ap необходимы для вычисления значений трен д авпервых m и последних m точках временного ряда, поскольку рассмотренный способ сглаживания ряда через a0 сделать это не позволяет.

Размерность матрицы M определяется степенью полинома: (p +1) (p +1), пределы суммирования во всех формулах задаются длиной отрезка скольжения.

Следовательно, для выбранных значений p и m можно получить общее решение в виде вектора (a0, a1,..., ap).

Свойства скользящих средних m 1. Сумма весов t в формуле a0 = txt равна единице.

t=-m Действительно, пусть все значения временного ряда равны одной и той же m m константе c. Тогд а txt = c t должна быть равна этой константе t=-m t=-m m c, а это возможно только в том случае, если t =1.

t=-m 2. Веса симметричны относительно нулевого значения t, т.е. t = -t 12.1. Метод скользящих средних Это следует из того, что весовые коэффициенты при каждом xt зависят от tj, а j принимает только четные значения.

3. Для полиномов четного порядка p =2k формулы расчета a0 будут теми же самыми, что и для полиномов нечетного порядка p =2k +1.

Пусть p =2k +1, тогда матрица коэффициентов системы (12.3) при неизвест ных параметрах a0, a1,..., ap будет выглядеть следующим образом:

m m m m t0 t t2k t2k+ t=-m t=-m t=-m t=-m m m m m t t2 t2k+1 t2k+ t=-m t=-m t=-m t=-m....

.....

....

m m m m t2k t2k+1 t4k t4k+ t=-m t=-m t=-m t=-m m m m m t2k+1 t2k+2 t4k+1 t4k+ t=-m t=-m t=-m t=-m Для нахождения a0 используются уравнения с четными степенями t при a0, следовательно, половина строк матрицы, включая последнюю, в расчетах участво вать не будет.

В этом блоке матрицы, содержащем коэффициенты при a0, a2, a4,..., по следний столбец состоит из нулей, так как его элементы Ч суммы нечетных сте пеней t. Таким образом, уравнения для нахождения a0 при нечетном значении p =2k +1 в точности совпадают с уравнениями, которые надо решить для нахож дения a0 при меньшем на единицу четном значении p =2k:

m m m t0 t2 t2k t=-m t=-m t=-m m m m t2 t4 t2k+ t=-m t=-m t=-m.

...

....

...

m m m t2k t2k+2 t4k t=-m t=-m t=-m 4. Оценки параметров a1,..., ap тоже выражены в виде линейной комбинации уровней временного ряда, входящих в отрезок, но весовые коэффициенты в этих формулах в сумме равны нулю и не симметричны.

Естественным образом возникает вопрос, какой степени полином следует вы бирать и какой должна быть длина отрезка скольжения. Закономерность такова:

чем выше степень полинома и короче отрезок скольжения, тем ближе расчетные 398 Глава 12. Сглаживание временного ряда значения к первоначальным данным. При этом, помимо тенденции могут воспро изводиться и случайные колебания, нарушающие ее смысл. И наоборот, чем ниже степень полинома и чем длиннее отрезок скольжения, тем более гладкой является сглаживающая кривая, тем в большей мере она отвечает свойствам тенденции, хотя ошибка аппроксимации будет при этом выше.

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

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

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

st = xt + st-1, (12.4) где st Ч значение экспоненциальной средней в момент t, Ч параметр сглаживания (вес последнего наблюдения), 0 <1, =1 -.

Экспоненциальную среднюю, используя рекуррентность формулы (12.4), мож но выразить через значения временного ряда:

st = xt + (xt-1 + st-2) =xt + xt-1 + 2st-2 =... = = xt + xt-1 + 2xt-2 +... + jxt-j +... + t-1x1 + ts0 = t- = jxt-j + ts0, (12.5) j= t Ч количество уровней ряда, s0 Ч некоторая величина, характеризующая на чальные условия для первого применения формулы (12.4) при t =1. Вкачестве s0 можно использовать первое значение временного ряда, т.е. x1.

Так как < 1, то при t величина t 0, а сумма коэффициентов t- j 1.

j= Действительно, 1 j = =(1 - ) =1.

1 - 1 - j= 12.2. Экспоненциальное сглаживание Тогда последним слагаемым в формуле (12.5) можно пренебречь и st = jxt-j = (1 - )jxt-j.

j=0 j= Таким образом, величина st оказывается взвешенной суммой всех уровней ря да, причем веса уменьшаются экспоненциально, по мере углубления в историю процесса, отсюда название Ч экспоненциальная средняя.

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

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

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

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

Выделяют два вида моделей, которые можно изобразить схематично:

1. Модель с аддитивным сезонным эффектом, предложенная Тейлом и Вей джем (Theil H., Wage S.):

xt = ft + gt + t, (12.6) где ft отражает тенденцию развития процесса, gt, gt-1,..., gt-k+1 Ч аддитив ные коэффициенты сезонности;

k Ч количество опорных временных интервалов (фаз) в полном сезонном цикле;

t Чбелыйшум.

2. Модель с мультипликативным сезонным эффектом, разработанная Уин терсом (Winters P.R.):

xt = ft mt t, (12.7) где mt, mt-1,..., mt-k+1 Ч мультипликативные коэффициенты сезонности.

400 Глава 12. Сглаживание временного ряда В принципе, эта модель после логарифмирования может быть преобразована в модель с аддитивным сезонным эффектом.

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

И в том, и в другом случае обновление параметров модели производится по схе ме экспоненциального сглаживания. Оба варианта допускают как наличие тенден ции (линейной или экспоненциальной), так и ее отсутствие.

Множество комбинаций различных типов тенденций с циклическими эффекта ми аддитивного и мультипликативного характера можно представить в виде обоб щенной формулы:

ft = fd1 +(1- f )d2, где ft Ч некоторый усредненный уровень временного ряда в момент t после устранения сезонного эффекта, f Ч параметр сглаживания, 0

xt, Ч если сезонный эффект отсутствует, xt - gt-k, Ч в случае аддитивного сезонного эффекта, d1 = xt, Ч в случае мультипликативного сезонного эффекта.

mt-k Таким образом, d1 представляет собой текущую оценку процесса xt, очищенную от сезонных колебаний с помощью коэффициентов сезонности gt-k или mt-k, рассчитанных для аналогичной фазы предшествующего цикла.

ft-1, Ч при отсутствии тенденции, d2 = ft-1 + ct-1, Ч в случае аддитивного роста, ft-1 rt-1, Ч в случае экспоненциального роста.

В этой формуле ct-1 Ч абсолютный прирост, характеризующий изменение сред него уровня процесса, или аддитивный коэффициент роста, rt-1 Ч коэффициент экспоненциального роста.

Например, для модели с аддитивным ростом и мультипликативным сезонным эффектом подойдет график, изображенный на рисунке 12.1а, а для модели с экс поненциальным ростом и аддитивным сезонным эффектом Ч график на рисунке 12.1б.

12.2. Экспоненциальное сглаживание Примеры графиков для некоторых типов адаптивных сезонных моделей a) б) xt xt t t Модель с аддитивным ростом Модель с экспоненциальным ростом и мультипликативным сезонным эффектом и аддитивным сезонным эффектом Рис. 12.1. Графики некоторых типов временных рядов Адаптация всех перечисленных параметров осуществляется с помощью экспо ненциального сглаживания:

gt = g(xt - ft) +(1- g)gt-k, xt mt = m +(1- m)mt-k, ft ct = c(ft - ft-1) +(1- c)ct-1, ft rt = r +(1- r)rt-1, ft- где 0

Первые две формулы представляют собой линейную комбинацию текущей оценки коэффициента сезонности, полученной путем устранения из исходного уровня процесса значения тренда ( xt - ft и xt/ft), и оценки этого параметра на аналогичной фазе предшествующего цикла ( gt-k и mt-k). Аналогично, две последние формулы являются взвешенной суммой текущей оценки коэффициен та роста (соответственно, аддитивного ft - ft-1 и экспоненциального ft/ft-1) и предыдущей его оценки ( ct-1 и rt-1).

Очевидно, что в случае отсутствия тенденции и сезонного эффекта получается простая экспоненциальная средняя:

ft = f xt +(1- f )ft-1.

402 Глава 12. Сглаживание временного ряда Рассмотрим для иллюстрации модель Уинтерса с аддитивным ростом и мульти пликативным сезонным эффектом:

xt ft = f +(1- f )(ft-1 + ct-1), mt-k xt mt = m +(1- m)mt-k, (12.8) ft ct = c(ft - ft-1) +(1- c)ct-1.

Расчетные значения исследуемого показателя на каждом шаге, после обновле ния параметров ft, mt и ct, получаются как произведение ft mt.

Прежде чем воспользоваться полной схемой экспоненциального сглаживания (12.8), а сделать это можно начиная с момента t = k +1, необходимо получить начальные, отправные значения перечисленных параметров.

Для этого с помощью МНК можно оценить коэффициенты f1 и c1 регрессии:

xt = f1 + c1t + t, и на первом сезонном цикле (для t =1,..., k) адаптацию параметров произвести по усеченному варианту:

ft = f xt +(1- f )ft-1, xt mt =, t =1,..., k, ft ct = c(ft - ft-1) +(1- c)ct-1, gt = xt - ft.

Задача оптимизации модели сводится к поиску наилучших значений парамет ров f, m, c, выбор которых определяется целями исследования и характером моделируемого процесса. Уинтерс предлагает находить оптимальные уровни этих коэффициентов экспериментальным путем, с помощью сетки значений f, m, c (например, (0, 1;

0, 1;

0, 1), (0, 1;

0, 1;

0, 2),... ). В качестве критерия сравне ния вариантов рекомендуется стандартное отклонение ошибки.

12.3. Упражнения и задачи Упражнение 1.1. Сгенерируйте 20 рядов по 100 наблюдений на основе полиномиального тренда t =5 + 4t - 0, 07t2 +0.0005t3 с добавлением белого шума с нор мальным распределением и дисперсией 20.

Таблица 12.2. Производство природного газа в СССР (миллиардов кубических футов) январь февраль март апрель май июнь июль август сентябрь октябрь ноябрь декабрь 1971 653.1 589.5 653.1 610.7 610.7 583.2 600.1 614.2 600.1 642.5 642.5 670. 1972 670.8 649.5 695.4 664.5 638.9 621.3 620.7 619.4 624.8 653.1 663.6 706. 1973 720.1 656.6 734.2 691.9 688.3 688.4 691.2 701.2 653.1 673.2 720.1 673. 1974 720.1 709.5 776.6 737.7 741.3 723.7 724.6 758.9 760.0 808.4 811.2 882. 1975 864.9 871.9 868.4 861.2 864.8 833.3 833.1 829.6 829.6 840.2 900.1 953. 1976 953.1 914.3 967.2 921.3 917.8 916.8 924.9 924.9 917.8 988.4 974.3 1009. 1977 1048.4 960.2 960.2 1048.4 998.9 956.6 984.9 995.5 999.0 1175.5 1180.0 1190. 1978 1129.6 1129.4 1126.1 1076.7 1080.2 1034.3 1062.5 1064.7 1023.7 1147.2 1136.7 1196. 1979 1230.0 1220.0 1220.0 1175.5 1182.5 1140.2 1157.8 1161.4 1164.9 1249.6 1250.6 1306. 1980 1309.6 1232.0 1306.1 1246.1 1256.7 1200.2 1246.1 1260.2 1270.8 1270.0 1323.8 1376. 1981 1419.1 1299.0 1420.0 1345.0 1313.0 1271.0 1270.0 1334.0 1334.0 1430.0 1430.0 1460. 1982 1504.0 1380.0 1528.5 1436.7 1457.9 1412.0 1419.1 1436.7 1447.3 1546.1 1528.5 1623. 1983 1627.4 1486.1 1652.0 1528.5 1570.8 1517.9 1514.4 1539.1 1482.6 1648.5 1648.5 1747. 1984 1747.4 1648.5 1757.9 1680.3 1697.9 1623.8 1669.7 1697.9 1694.4 1821.5 1803.8 1870. 1985 1930.9 1775.6 1941.5 1853.2 1892.1 1765.0 1825.0 1846.2 1870.9 1990.9 1962.7 2047. 1986 2075.6 1895.6 2118.0 1983.9 2005.0 1906.2 1959.2 1969.7 1976.8 2103.9 2089.8 2188. 1987 2221.3 2030.6 2210.7 2083.6 2118.9 2012.9 2048.3 2048.3 2083.6 2223.9 2259.2 2330. 1988 2369.6 2221.3 2366.1 2224.8 2275.0 2146.6 2118.9 2189.5 2189.5 2357.0 2394.8 2447. 1989 2510.0 2300.0 2391.0 2333.0 2336.0 2187.7 2208.0 2279.0 2200.0 2500.0 2484.0 2495. 1990 2630.0 2400.0 2420.0 2391.0 2430.0 2250.0 2340.0 2340.0 2250.0 2500.0 2450.0 2460. 12.3. Упражнения и задачи 404 Глава 12. Сглаживание временного ряда а) Проведите сглаживание сгенерированных рядов с помощью полинома первой степени с длиной отрезка скольжения 5 и 9.

б) Выполните то же задание, используя полином третьей степени.

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

1.2. Имеются данные о производстве природного газа в СССР (табл. 12.2).

а) Постройте графики ряда и логарифмов этого ряда. Чем они различают ся? Выделите основные компоненты временного ряда. Какой характер носит сезонность: аддитивный или мультипликативный? Сделайте вывод о целесообразности перехода к логарифмам.

б) Примените к исходному ряду метод экспоненциального сглаживания, подобрав параметр сглаживания.

в) Проведите сглаживание временного ряда с использованием адаптивной сезонной модели.

Задачи 1. Сгладить временной ряд x =(3, 4, 5, 6, 7, 11), используя полином первого порядка с длиной отрезка скольжения, равной трем.

2. Записать формулу расчета вектора коэффициентов для полинома третьей степени с помощью метода скользящей средней в матричной форме с рас шифровкой обозначений.

3. В чем специфика аппроксимации первых m и последних m точек временного ряда при использовании метода скользящих средних?

4. Найти параметры адаптивной сезонной модели для временного ряда x =(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4,... ).

5. Изобразить график временного ряда с аддитивным ростом и мультиплика тивным сезонным эффектом.

6. Изобразить график временного ряда с экспоненциальным ростом и аддитив ным сезонным эффектом.

7. Записать модель с экспоненциальным ростом и мультипликативным сезон ным эффектом, а также формулу прогноза на 5 шагов вперед.

12.3. Упражнения и задачи Рекомендуемая литература 1. Андерсон Т. Статистический анализ временных рядов. Ч М.: Мир, 1976.

(Гл. 3).

2. Лукашин Ю.П. Адаптивные методы краткосрочного прогнозирования. Ч М.: Статистика, 1979. (Гл. 1, 2).

3. Кендалл М. Дж., Стьюарт А. Многомерный статистический анализ и вре менные ряды. Ч М.: Наука, 1976. (Гл. 46).

4. Маленво Э. Статистические методы эконометрии. Вып. 2. Ч М.: Статисти ка, 1976. (Гл. 11, 12).

5. Mills Terence C. Time Series Techniques for Economists, Cambridge University Press, 1990 (Ch. 9).

Глава Спектральный и гармонический анализ 13.1. Ортогональность тригонометрических функций и преобразование Фурье временного ряда Как известно, тригонометрические функции cos t и sin t являются периодиче скими с периодом 2:

cos(t +2) =cos t, sin(t +2) =sin t.

Функции cos(t - ) и sin(t - ) периодичны с периодом 2/. Действи тельно, cos(t - ) =cos(t +2 - ) =cos ((t +2/) - ), sin(t - ) =sin(t +2 - ) =sin((t +2/) - ).

Величина /2, обратная периоду, называется линейной частотой, назы вают угловой частотой. Линейная частота равна числу периодов (не обязательно целому), содержащемуся в единичном интервале, то есть именно такое число раз функция повторяет свои значения в промежутке [0, 1].

Рассмотрим функцию:

R cos(t - ) =R(cos t cos +sint sin ) = cos(t) + sin(t), 13.1 Ортогональность тригонометрических функций.

где = R cos, = R sin или, что эквивалентно, R = 2 + 2, tg = Коэффициент R, являющийся максимумом функции R cos(t - ) называется амплитудой этой функции, а угол называется фазой.

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

Две функции (t) и (t), определенные на конечном множестве {1,..., T}, называются ортогональными, если их скалярное произведение, определенное как сумма произведений значений (t) и (t) в этих точках, равно нулю:

T (t) (t) =0.

t= Система T тригонометрических функций в точках t {1,..., T} 2j T cjt =cos t, j =0, 1,...,, T (13.1) 2j T - sjt =sin t, j =1,..., T ортогональна, т.е. скалярное произведение векторов T T (cj, ck) = cjtckt =0, j = k, 0 j, k, (13.2) t= T T - (sj, sk) = sjtskt =0, j = k, 0

Для доказательства этого утверждения полезны следующие равенства T 0, при j =0, 2j cos t = (13.5) T T, при j =0, T, t= T 2j sin t =0, (13.6) T t= 408 Глава 13. Спектральный и гармонический анализ истинность которых легко установить, выразив тригонометрические функции через показательные с использованием формул Эйлера:

ei =cos i sin, (13.7) cos = (ei + e-i), (13.8) sin = (ei - e-i). (13.9) 2i Итак, при j = T T 2j 2j 2j cos t = ei T t + e-i T t = T t=1 t= 2j 2j 1 1 - ei2j 1 1 - e-i2j = ei T + e-i T =0, 2j 2j 1 - ei T 2 1 - e-i T где предпоследнее равенство получено из формулы суммы геометрической прогрес сии, а последнее Ч из формулы (13.7), т.к.

ei2j =cos(2j) i sin(2j) =1.

T 2j Очевидно, что при j =0, T cos t = T.

t=1 T Равенство (13.6) доказывается аналогично. При доказательстве соотношений (13.2Ц13.4) используются утверждения (13.5, 13.6).

Таким образом, T T T 2j 2k 1 2(j - k) 1 2(j + k) (cj, ck) = cos tcos t = cos t+ cos t = T T 2 T 2 T t=1 t=1 t= T 0, j = k, 0 j, k, T T = (13.10), j = k, 0

T T T 2j 2k 1 2(j - k) 1 2(j + k) (sj, sk) = sin t sin t = cos t - cos t = T T 2 T 2 T t=1 t=1 t= T - 0, j = k, 0

2 13.1 Ортогональность тригонометрических функций T 2j 2k (cj, sk) = cos t sin t = T T t= T T 1 2(j + k) 1 2(j - k) = sin t + sin t =0. (13.12) 2 T 2 T t=1 t= Мы доказали выполнение (13.2Ц13.4) для указанного набора функций, полу чив одновременно некоторые количественные их характеристики. Таким образом, 2j 2j функции cos t и sin t образуют ортогональный базис и всякую функцию, T T в том числе и временной ряд {xt}, определенный на множестве {1,..., T}, можно разложить по этому базису, т.е. представить в виде конечного ряда Фурье:

[T/2] 2j 2j xt = j cos t + j sin t, (13.13) T T j= или, вспоминая (13.1), кратко [T/2] xt = (jcjt + jsjt), j= где 0 и [T/2] при четном T отсутствуют (т.к. sin 0 = 0, sin t =0).

Величину 2j/T = j называют частотой Фурье, а набор скаляров j и j ( j =0, 1,..., [T/2]) Чкоэффициентами Фурье.

Если cjt и sjt Ч элементы векторов cj и sj, стоящие на t-ом месте, то, пе реходя к векторным обозначениям, (13.13) можно переписать в матричном виде:

x = C S, (13.14) где x =(x1,..., xT ), =(0,..., [T/2]), =(1,..., [(T -1)/2]), C = {cjt}, j =0, 1,..., [T/2], t =1,..., T, S = {sjt}, j =1,..., [(T - 1)/2], t =1,..., T.

410 Глава 13. Спектральный и гармонический анализ Перепишем в матричной форме свойства ортогональности тригонометрических функций, которые потребуются при вычислении коэффициентов Фурье:

c sk =0, j, k, j c 1T =0, k =0, k s 1T =0, k, k (13.15) c ck = s sk =0, j = k, j j c ck = s sk = T/2, k =0, T/2, k k c c0 = T, c cT/2 = T, для четных T, T/ где 1T =(1,..., 1) Ч T -компонентный вектор.

Для нахождения коэффициентов Фурье скалярно умножим c на вектор x и, j воспользовавшись изложенными свойствами ортогональности (13.15), получим:

c x = c C S =(c c0,..., c c[T/2], c s1,..., c s[(T -1)/2]) = j j j j j j T T = jc cj = j, для j =0,.

j 2 Таким образом, T 2 2 2j T j = c x = xt cos t, для j =0,, j T T T t= T 1 0 = c x = xt, (13.16) T T t= T 1 T/2 = c x = (-1)txt, для четных T.

T/ T T t= Аналогично находим коэффициенты j:

T 2 2 2j j = s x = xt sin t. (13.17) j T T T t= 13.2. Теорема Парсеваля 13.2. Теорема Парсеваля Суть теоремы Парсеваля состоит в том, что дисперсия процесса xt разлагается по частотам соответствующих гармоник следующим образом:

T/2- 2 var(xt) = Rj + RT/2, для четных T, (13.18) j= (T -1)/ var(xt) = Rj, для нечетных T. (13.19) j= Покажем, что это действительно так. Из (13.14) мы имеем:

C x x = = C S S C C C S = = S C S S C = = 0 S [T/2] [(T -1)/2] = C + S = 21 1T + 2c cj + j s sj, 0 T j j j j=1 j= где c и s Ч диагональные матрицы. Таким образом, если T Ччетно, то T/2 T/2- x x = 21 1T + 2c cj + j s sj = 0 T j j j j=1 j= T/2-1 T/2-1 T/2- T T T 2 = 2T + 2 + j + 2 T = 2T + (2 + j )+2 T = 0 j T/2 0 j T/ 2 2 j=1 j=1 j= T/2-1 T/2- T T 2 2 2 = 2T + Rj + 2 T = R0T + Rj + RT/2T. (13.20) 0 T/ 2 j=1 j= Аналогично для нечетных T :

(T -1)/2 (T -1)/ x x = 21 1T + 2c cj + j s sj = 0 T j j j j=1 j= 412 Глава 13. Спектральный и гармонический анализ (T -1)/2 (T -1)/ T T = 2T + 2 + j = 0 j 2 j=1 j= (T -1)/2 (T -1)/ T T 2 2 = 2T + (2 + j ) =R0T + Rj. (13.21) 0 j 2 j=1 j= Разделим уравнения (13.20) и (13.21) на T и перенесем в левые части R0. Сучетом того, что R0 = 2 =2, получаем выражения для дисперсии процесса xt.

x T/2- x x 2 2 var(xt) = - R0 = Rj + RT/2, для четных T, (13.22) T j= (T -1)/ x x 2 var(xt) = - R0 = Rj, для нечетных T. (13.23) T j= Таким образом, вклад в дисперсию процесса для T/2-й гармоники равен RT/2, ад ля k-й гармоники, k = T/2, равен Rk.

Следовательно, наряду с определением коэффициентов Фурье для k-й гармо ники, можно определить долю этой же гармоники в дисперсии процесса.

13.3. Спектральный анализ Введем понятия периодограммы и спектра.

Периодограммой называют последовательность значений {Ij}:

T T Ij = (2 + j ), j =0, 1,...,, j 2 T T т.е. Ij равно квадрату амплитуды j-ой гармоники, умноженному на, Ij = Rj.

Величина Ij называется интенсивностью на j-ой частоте.

На практике естественнее при вычислении периодограммы использовать цен трированный ряд xt = xt - x. При этом меняется только I0. Для центрированного ряда 0 = 0, поэтому I0 = 2 = 0. Все остальные значения периодограммы не меняются, что следует из (13.5) и (13.6) Ч влияние константы на остальные значения обнуляется. В оставшейся части главы мы будем использовать только центрированный ряд.

В определении периодограммы принципиальным является то, что гармониче ские частоты fj = j/T (j = 0, 1,..., [T/2]) изменяются дискретно, причем наиболее высокая частота составляет 0, 5 цикла за временной интервал.

13.3. Спектральный анализ Вводя понятие спектра, мы ослабляем это предположение и позволяем частоте изменяться непрерывно в диапазоне 0 - 0.5 Гц ( 0.5 цикла в единицу времени).

Обозначим линейную частоту через f и введем следующие обозначения:

T T 2 f = xt cos(2ft), f = xt sin(2ft) T T t=1 t= и 2 Rf = 2 + f.

f Функция T T 2 p(f) = Rf = 2 + f = f 2 2 T T = xt cos(2ft) + xt sin(2ft), (13.24) T t=1 t= где 0 f, называется выборочным спектром. Очевидно, что значения пе риодограммы совпадают со значениями выборочного спектра в точках fj, то есть p(fj) =Ij.

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

И периодограмму, и спектр представляют для наглядности в виде графика, на оси ординат которого Ч интенсивность Ij или p(f), на оси абсцисс Ч частота j fj = или f, соответственно. График выборочного спектра часто называют T спектрограммой.

Спектрограмма нужна для более наглядного изображения распределения дис k персии между отдельными частотами. Если частоте f = соответствует пик на T спектрограмме, то в исследуемом ряду есть существенная гармоническая состав T ляющая с периодом 1/f =.

k Целью спектрального анализа является определение основных существенных гармонических составляющих случайного процесса путем разложения дисперсии процесса по различным частотам. Спектральный анализ позволяет исследовать смесь регулярных и нерегулярных спадов и подъемов, выделять существенные гармоники, получать оценку их периода и по значению спектра на соответствующих частотах судить о вкладе этих гармоник в дисперсию процесса.

414 Глава 13. Спектральный и гармонический анализ Исследования показывают, что наличие непериодического тренда (тренда с бес конечным периодом) дает скачок на нулевой частоте, т.е. в начале координат спек тральной функции. При наличии циклических составляющих в соответствующих частотах имеется всплеск;

если рядслишком зазубрен, мощность спектра пере мещается в высокие частоты.

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

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

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

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

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

Теорема ВинераЧХинчина:

T -1 T - p(f) =2 c0 +2 ck cos 2fk =2s2 1+2 rk cos 2fk, (13.25) k=1 k= где rk = ck/c0 = ck/s2 Ч выборочные автокорреляции.

Доказательство.

Объединим коэффициенты Фурье f, f в комплексное число df = f - if, где i Ч мнимая единица. Тогда T T T p(f) = 2 + f = (f - if )(f + if) = df df, (13.26) f 2 2 где df комплексно сопряжено с df.

13.4. Связь выборочного спектра с автоковариационной функцией Используя формулы для f и f, получим:

T T 2 df = xt cos 2ft - i sin 2ft = xte-i2ft.

T T t=1 t= Точно так же T df = xtei2ft.

T t= Подставляя полученные значения df и df в выражение (13.26), получаем:

T T T 2 p(f) = xte-i2ft xtei2ft = 2 T T t=1 t = T T = xtxt e-i2f(t-t ). (13.27) T t=1 t = Произведем замену переменных: пусть k = t - t. Так как автоковариация равна T ck = xtxt-k, T t=k+ что тождественно T -k ck = xtxt+k, T t= то выражение (13.27) преобразуется следующим образом:

T T T -1 T 2 xtxt e-i2f(t-t ) = e-i2fk xtxt-k = T T t=1 t =1 k=-(T -1) t=k+ T -1 T T - =2 e-i2fk xtxt-k =2 e-i2fkck.

T k=-(T -1) t=k+1 k=-(T -1) Тогда T -1 T - p(f) =2 e-i2fkck =2 ck cos 2fk - i sin 2fk = k=-(T -1) k=-(T -1) T - =2 c0 +2 ck cos 2fk, 0 f.

k= 416 Глава 13. Спектральный и гармонический анализ Теперь допустим, что выборочный спектр, характеризующийся эмпирическими значениями частоты, амплитуды и фазы, вычислен для ряда из T наблюдений и мы можем неоднократно повторить этот эксперимент, соответственно собрав множество значений j, j и p(f) по повторным реализациям. Тогда среднее значение p(f) будет равно:

T - E(p(f)) = 2 E(c0) +2 E(ck)cos2fk, (13.28) k= где c0 и ck Ч эмпирические значения автоковариации. C учетом того, что E(ck) при больших T стремится к теоретической автоковариации k, получим, переходя к пределу, так называемую теоретическую спектральную плотность, или спектр мощности:

p(f) = lim E(p(f)), 0 f T или p(f) =2 0 +2 k cos(2fk) = k= =22 1+2 k cos(2fk), (13.29) k= где k = k/0 = k/2 Ч теоретические автокорреляции.

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

Иногда более удобно использовать автокорреляции: разделим обе части p(f) на дисперсию процесса, 2, и получим нормированный спектр g(f):

g(f) =2 1+2 k cos 2fk, 0 f.

k= Если процесс представляет собой белый шум, то, согласно приведенным фор мулам, p(f) = 22, т.е. спектральная плотность белого шума постоянна. Этим объясняется термин белый шум. Подобно тому, как в белом цвете смешаны в одинаковых объемах все цвета, так и белый шум содержит все частоты, и ни одна из них не выделяется.

13.5. Оценка функции спектральной плотности 13.5. Оценка функции спектральной плотности На первый взгляд, выборочный спектр, определенный как T - p(f) = 2 cke-i2fk = (13.30) k=-(T -1) T - = 2 c0 +2 ck cos 2fk, 0 f, (13.31) k= является естественной и правильной оценкой функции спектральной плотности:

+ p(f) =2 ke-i2fk =2 0 +2 k cos 2fk, 0 f. (13.32) k=- k= Известно, что выборочная автоковариация ck Ч это асимптотически несме щенная оценка параметра k, таккак lim E(ck) =k, T и поэтому выборочный спектр есть также асимптотически несмещенная оценка функции спектральной плотности.

Однако дисперсия оценки (выборочного спектра) не уменьшается по мере роста размера выборки. Это означает, что рассматриваемая оценка несостоятельна. В то время как график функции теоретической спектральной плотности стационарного стохастического процесса гладкий, Ч график выборочного спектра, построен ный на основе эмпирических данных, неровный. Использование данной оценки в качестве оценки функции спектральной плотности может привести к ложным выводам.

Поведение выборочного спектра иллюстрируют спектрограммы на рис. 13.1 а), б), в). Гладкая жирная линия соответствует теоретической спектральной плотно сти случайного процесса, а неровная тонкая линия Ч оценке по формуле (13.24).

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

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

418 Глава 13. Спектральный и гармонический анализ а) б) T = T = в) г) T = 1000 T = 1000, сглаженная оценка Рис. 13. 1) Взвешивание ординат периодограммы.

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

Сглаживающая оценка определяется в форме M M j j - k ps = k p = kIj-k. (13.33) T T k=-M k=-M Здесь {-M,..., -1, 0, 1,..., M } Ч 2M +1 коэффициентов скользящего среднего, которые в сумме должны давать единицу, а также должны быть симмет ричными, в смысле -k = k. Как правило, коэффициент 0 максимальный, а остальные коэффициенты снижаются в обе стороны. Таким образом, наиболь ший вес в этой оценке имеют значения выборочного спектра, ближайшие к данной j частоте, в связи с чем данный набор коэффициентов называют спектральным T окном. Через это окно мы как бы смотрим на функцию спектральной плотности.

Формула сглаженной спектральной оценки определяется только для значений j = M,..., [T/2] - M. Для гармонических частот с номерами j =0,..., M - и j = [T/2] - M +1,..., [T/2] оценка не определена, поскольку при дан ных значениях j величина (j - k) может принимать значения -M,..., -1;

[T/2] + 1,..., [T/2] + M. Однако проблема краевых эффектов легко решается, 13.5. Оценка функции спектральной плотности если доопределить функцию выборочного спектра (и, соответственно, периодо грамму), сделав ее периодической. Формула (13.24) дает такую возможность, по скольку основана на синусоидах и косинусоидах. Таким образом, будем считать, что выборочный спектр определен формулой (13.24) для всех частот f (-, +).

При этом ясно, что выборочный спектр будет периодической функцией с перио дом 1, зеркально-симметричной относительно нуля (и относительно любой часто ты вида i/2, гд е i Ч целое число):

p(f - i) =p(f), i =..., -1, 0, 1,...

и p(-f) =p(f).

Формулу (13.33) несложно обобщить так, чтобы можно было производить сгла живание не только в гармонических частотах. Кроме того, в качестве расстояния между лусредняемыми частотами можно брать не 1/T, а произвольное > 0.

Таким образом приходим к следующей более общей формуле:

M ps(f) = k p(f - k). (13.34) k=-M Просматривая поочередно значения выборочного спектра и придавая наиболь ший вес текущему значению, можно сгладить спектр в каждой интересующей нас точке.

2) Взвешивание автоковариационных функций.

Известно, что при увеличении лага k выборочные автоковариации ck явля ются все более неточными оценками. Это связано с тем, что k-я автоковариация вычисляется как среднее по T -k наблюдениям. Отсюда возникает идея в формуле (13.25) ослабить влияние дальних автоковариаций за счет применения понижаю щих весов так, чтобы с ростом k происходило уменьшение весового коэффициента при ck.

Если ряд весов, связанных с автоковариациями c0, c1,..., cT -1, обозначить как m0, m1,..., mT -1, оценка спектра будет иметь вид:

T - pc(f) =2 m0c0 +2 mkck cos 2fk, где 0 f. (13.35) k= Набор весов mk, используемый для получения такой оценки, получил название корреляционное, илилаговое окно.

420 Глава 13. Спектральный и гармонический анализ При использовании корреляционного окна для уменьшения объема вычислений удобно брать такую систему весов, что mk = 0 при k K, гд е K < T. Тогд а формула (13.35) приобретает вид K- pc(f) =2 m0c0 +2 mkck cos 2fk. (13.36) k= Корреляционное окно удобно задавать с помощью функции m(), заданной на интервале [0;

1], такой, что m(0) = 1, m(1) = 0. Обычно функцию выбирают так, чтобы между нулем и единицей эта функция плавно убывала. Тогда понижающие веса mk при k =0,..., K вычисляются по формуле mk = m(k/K).

Ясно, что при этом m0 = 1 (это величина, с помощью которой мы взвешиваем дисперсию в формуле (13.35)) и mK =0.

Наиболее распространенные корреляционные окна, удовлетворяющие пере численным свойствам Ч это окна Парзена и ТьюкиЧХэннинга (см. рис. 13.2).

Окно ТьюкиЧХэннинга:

m(z) = (1 + cos(z)).

Окно Парзена:

1 - 6z2 +6z3, если z [0;

];

m(z) = 2(1 - z)3, если z [1;

1].

Можно доказать, что сглаживание в частотной области эквивалентно пони жающему взвешиванию автоковариаций. Чтобы это сделать, подставим в (13.34) выборочный спектр, выраженный через комплексные экспоненты (13.30):

M T - ps(f) =2 j cke-i2(f-j)k.

j=-M k=-(T -1) Меняя здесь порядок суммирования, получим T -1 M ps(f) =2 cke-i2fk jei2jk.

k=-(T -1) j=-M 13.5. Оценка функции спектральной плотности m(z) окно Парзена 0. окно ТьюкиЧХэннинга 0. 0. 0. z 0.2 0.4 0.6 0.8 Рис. 13.2. Наиболее популярные корреляционные окна Введя обозначение M M mk = jei2jk = 0 + j ei2jk + e-i2jk j=-M j= M = 0 +2 j cos(2jk), j= придем к формуле T - ps(f) = 2 mkcke-i2fk = k=-(T -1) T - = 2 m0c0 +2 mkck cos 2fk, k= которая, очевидно, совпадает с (13.35).

Кроме того, без доказательства отметим, что подбором шага ивесов j мы всегда можем сымитировать применение корреляционного окна (13.35) с помощью (13.34). Существует бесконечно много способов это сделать.

В частности, как несложно проверить, частотное окно ТьюкиЧХэннинга полу 1 чим, если возьмем M =1, 0 =, 1 = -1 = и =.

2 2K 422 Глава 13. Спектральный и гармонический анализ Особого внимания требует вопрос о том, насколько сильно нужно сглаживать спектральную плотность. Для корреляционных окон степень гладкости зависит от того, насколько быстро убывают понижающие веса. При фиксированной функ ции m() все будет зависеть от параметра K Ччемменьше K, тем более гладкой будет оценка.

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

При этом принято говорить о ширине окна, илиширине полосы.

Если ширина окна слишком большая, то произойдет пересглаживание, оцен ка будет сильно смещенной, и пики спектральной плотности станут незаметными.

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

Таким образом, ширина окна выбирается на основе компромисса между смещением и дисперсией.

13.6. Упражнения и задачи Упражнение Сгенерируйте ряд длиной 200 по модели:

xt =10 +0.1t +4sin(t/2) - 3cos(t/2) + t, где t Ч нормально распределенный белый шум с дисперсией 3. Предположим, что параметры модели неизвестны, а имеется только сгенерированный ряд xt.

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

1.2. Выделите тренд и гармоническую составляющую, сравните их параметры с истинными значениями.

Упражнение Для временного ряда, представленного в таблице 12.2 (с. 403), выполните следующие задания.

2.1. Исключите из временного ряда тренд.

13.6. Упражнения и задачи 2.2. Остатки ряда, получившиеся после исключения тренда, разложите в ряд Фу рье.

2.3. Найдите коэффициенты j и j разложения этого ряда по гармоникам.

2.4. Постройте периодограмму ряда остатков и выделите наиболее существенные гармоники.

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

2.6. Вычислите выборочные коэффициенты автоковариации и автокорреляции для ряда остатков после исключения тренда.

2.7. Найдите значение периодограммы для частоты 0.5 разными способами (в том числе через автоковариационную функцию).

Упражнение Используя данные таблицы 12.2, выполните следующие задания.

j 3.1. С помощью периодограммы вычислите оценку спектра на частоте fj =, 2K K = 4. Получите сглаженную оценку спектра с помощью спектрального окна ТьюкиЧХэннинга.

3.2. Рассчитайте автокорреляционную функцию rj, j = 1,..., 4. Оцените спектр с помощью корреляционного окна ТьюкиЧХэннинга при K = 4.

Сравните с предыдущим результатом.

3.3. Оцените спектр с помощью корреляционного окна Парзена при K =4.

3.4. Постройте график оценки спектра для корреляционного окна ТьюкиЧ j Хэннинга в точках fj =, j =0,..., 20.

Задачи 1. Записать гармонику для ряда x =(1, -1, 1, -1, 1, -1,... ).

2. Пусть временной ряд xt имеет гармонический тренд: t = 3 cos(t) + +4sin(t). Найти значения амплитуды, фазы и периода.

3. Записать ортогональный базис, по которому разлагается исследуемый про цесс xt в ряд Фурье для T =6, T =7.

424 Глава 13. Спектральный и гармонический анализ 4. Записать ковариационную матрицу для гармонических переменных, состав ляющих ортогональный базис, если T =6, T =7.

5. Привести формулу расчета коэффициентов гармонических составляющих временного ряда.

T 6. Что описывает формула: I(f) = (2 + f ), гд е 0 f ? Почему f f 2 меняется в указанном диапазоне значений?

7. Как соотносятся понятия интенсивности и амплитуды, периодограммы и спек тра?

8. Вывести формулу для определения периодограммы на нулевой частоте.

9. Как связаны выборочный спектр и автокорреляционная функция для чисто случайного процесса? Записать формулу с расшифровкой обозначений.

10. Пусть для ряда из 4-х наблюдений выборочная автокорреляционная функция 1 1 равна: r1 = 2, r2 = 2, r3 = 2, дисперсия равна 1. Вычислить значение выборочного спектра на частоте 4.

11. Как соотносится выборочный спектр с автоковариационной функцией, спек тральными и корреляционными окнами?

12. Пусть для ряда из 4-х наблюдений выборочная автокорреляционная функция равна: r1 = r2 = r3 = - 2, 1 4, 4, дисперсия равна 1. Вычислить значение сглаженной оценки выборочного спектра на частоте с помощью окна Парзена с весами mk, k =1, 2.

13. По некоторому временному ряду рассчитана периодограмма:

1 1 I(0) = 2, I =6, I =1, I =4.

6 3 Найти оценки спектральной плотности для тех же частот с использованием окна ТьюкиЧХэннинга.

14. Записать уравнение процесса с одной периодической составляющей для ча стоты 0.33, амплитуд ы 2 ифазы 0.

15. Изобразить графики спектра для стационарных и нестационарных процессов.

13.6. Упражнения и задачи Рекомендуемая литература 1. Андерсон Т. Статистический анализ временных рядов. Ч М.: Мир, 1976.

(Гл. 4, 9).

2. Бокс Дж., Дженкинс Г. Анализ временных рядов. Прогноз и управление.

Вып. 1. Ч М.: Мир, 1974. (Гл. 2).

3. Гренджер К., Хатанака М. Спектральный анализ временных рядов в эконо мике. Ч М.: Статистика, 1972.

4. Дженкинс Г., Ваттс Д. Спектральный анализ и его приложения. Ч М.:

Мир, 1971.

5. Кендалл М. Дж., Стьюарт А. Многомерный статистический анализ и вре менные ряды. Ч М.: Наука, 1976. (Гл. 49).

6. Маленво Э. Статистические методы эконометрии. Вып. 2. Ч М.: Стати стика, 1976. (Гл. 11, 12).

7. Бриллинджер Д. Временные ряды. Обработка данных и теория. Ч М.: Мир, 1980. (Гл. 5).

8. Hamilton James D., Time Series Analysis. Ч Princeton University Press, 1994.

(Ch. 6).

9. Judge G.G., Griffiths W.E., Hill R.C., Luthepohl H., Lee T. Theory and Practice of Econometrics. Ч New York: John Wiley & Sons, 1985. (Ch. 7).

Глава Линейные стохастические модели ARIMA 14.1. Модель линейного фильтра Стационарный стохастический процесс {xt} с нулевым математическим ожи данием иногда полезно представлять в виде линейной комбинации последователь ности возмущений t, t-1, t-2,..., т.е.

xt = t + 1t-1 + 2t-2 +... = it-i, (14.1) i= или с использованием лагового оператора:

xt =(1 +1L + 2L2 + )t, где 0 =1 и выполняется |i| <, (14.2) i= т.е. ряд абсолютных значений коэффициентов сходится.

Уравнение (14.1) называется моделью линейного фильтра, а линейный опе ратор:

(L) =1 +1L + 2L2 + = iLi, i= 14.1. Модель линейного фильтра преобразующий t в xt, Чоператором линейного фильтра.

Компактная запись модели линейного фильтра выглядит следующим образом:

xt = (L)t.

Предполагается, что последовательность { t } представляет собой чисто слу чайный процесс или, другими словами, белый шум (см. стр. 353). Напомним, что автоковариационная и автокорреляционная функции белого шума имеют очень простую форму:

, k =0, 1, k =0, k = = k 0, k =0, 0, k =0, а его спектральная плотность имеет вид p(f) =20 =2 = const.

Таким образом, белый шум легко идентифицируется с помощью графиков автокор реляционной функции и спектра. Часто предполагается, что последовательность {t} состоит из независимых одинаково распределенных величин. Упростить ана лиз помогает дополнительное предположение о том, что {t} имеет нормальное распределение, т.е. представляет собой гауссовский белый шум.

Данная модель не является произвольной. Фактически, согласно теореме Воль да, любой слабо стационарный ряд допускает представление в виде модели линей ного фильтра, а именно: разложение Вольда ряда xt1. Следует помнить, однако, что разложение Вольда единственно, в то время как представление (14.1), вообще говоря, неоднозначно2. Таким образом, разложение Вольда представляет процесс в виде модели линейного фильтра, в то время как модель линейного фильтра не обя зательно задает разложение Вольда.

Как мы увидим в дальнейшем, модель линейного фильтра (14.1) применима не только к стационарным процессам, таким что выполняется (14.2), Ч с соответ ствующими оговорками она упрощает анализ и многих нестационарных процессов.

Если процесс {xt} подчинен модели (14.1), то при выполнении условия (14.2) он имеет математическое ожидание, равное нулю:

E(xt) = iE(t-i) =0.

i= В разложении Вольда произвольного стационарного процесса может присутствовать также пол ностью предсказуемая (линейно детерминированная) компонента. Однако такая компонента, если ее свойства известны, не создает больших дополнительных сложностей для анализа.

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

428 Глава 14. Линейные стохастические модели ARIMA Если требуется, чтобы математическое ожидание xt не было равно нулю, то урав нение модели линейного фильтра должно включать константу:

xt = + it-i = + (L)t.

i= Выведем формулы для автоковариаций рассматриваемой модели:

k = E(xtxt+k) =E it-i jt+k-j = i=0 j= = ijE(t-it+k-j) = ii+k. (14.3) i=0 j=0 i= Здесь учитывается, что для белого шума, j = i + k, E(t-it+k-j) = 0, j = i + k.

Заметим, что из (14.2) следует сходимость возникающих здесь рядов. Это го ворит о том, что данное условие подразумевает стационарность.

Действительно, пусть (14.2) выполнено. Тогда существует индекс I, такой что |i| 1 при i >I (иначе бы ряд не сошелся). Тогда I |ii+k| |i||i+k| |i||i+k| + |i+k|.

i=0 i=0 i=0 i=I+ Поскольку оба слагаемых здесь конечны, то |ii+k| <.

i= Ясно, что модель линейного фильтра (14.1) в общем виде представляет в основ ном теоретический интерес, поскольку содержит бесконечное число параметров.

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

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

Частными случаями применения рациональной аппроксимации являются мо дели авторегрессии AR(p) и скользящего среднего MA(q). В общем случае полу чаем смешанные процессы авторегрессии Ч скользящего среднего ARMA(p, q).

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

14.2. Влияние линейной фильтрации на автоковариации и спектральную плотность Пусть два стационарных процесса c нулевым математическим ожиданием {zt} и {yt} связаны между собой соотношением:

zt = jyt-j, j т.е. zt получается применением к yt линейного фильтра (L) = jLj.

j Пределы суммирования не указываем, поскольку они могут быть произвольными3.

y Пусть k Ч автоковариации процесса yt, а py(f) Ч его спектральная плот ность. Найдем те же величины для zt. Автоковариации zt равны z k = E(ztzt-k) =E jyt-j syt-k-s = s j y = jsE(yt-jyt-k-s) = jsk+s-j.

s s j j Спектральную плотность zt можно записать в виде z pz(f) =2 kei2fk.

k= В том числе, при соответствующих предположениях, пределы суммирования могут быть беско нечными. Кроме того, zt здесь может зависеть от опережающих значений yt.

430 Глава 14. Линейные стохастические модели ARIMA Подставим сюда формулы для ковариаций:

y pz(f) = 2 jsk+s-jei2fk = s k=- j y = 2 js k+s-jei2fk.

s j k= Произведем здесь замену k = k + s - j:

y pz(f) = 2 js k ei2f(k +j-s) = s j k = y = 2 k ei2fk jei2fjse-i2fs = s k =- j = 2 k ei2fk jei2fj se-i2fs.

s k =- j Первый множитель Ч это спектральная плотность yt. Два последних множи теля представляют собой сопряженные комплексные числа, поэтому их произве дение равно квадрату их модуля. Окончательно получим pz(f) =py(f) je-i2fj (14.4) j или 2 pz(f) =py(f) j cos 2fj + j sin 2fj. (14.5) j j Данную теорию несложно применить к модели линейного фильтра (14.1). Для этого заменяем zt на xt, а yt на t. Автоковариации xt равны k = ijk+j-i.

i=0 j= Поскольку 0 =, и k =0 при k =0, то k = ii+k, i= 14.3. Процессы авторегрессии что совпадает с полученной ранее формулой (14.3).

Для спектральной плотности из (14.4) получаем p(f) =2 je-i2fj 2. (14.6) j= Поскольку |e-i2fj| =1, то ряд здесь сходится и | p(f)| <.

14.3. Процессы авторегрессии В модели авторегрессии текущее значение процесса xt представляется в виде линейной комбинации конечного числа предыдущих значений процесса и белого шума t:

xt = 1xt-1 + 2xt-2 + + pxt-p + t, (14.7) при этом предполагается, что текущее значение t не коррелировано с лагами xt.

Такая модель называется авторегрессией p-го порядка и обозначается AR(p) (от английского autoregression).

Используя лаговый оператор L, представим уравнение авторегрессии в виде:

(1 - 1L - 2L2 -... - pLp)xt = t, или кратко, через лаговый многочлен (L) =1 - 1L - 2L2 -... - pLp:

(L)xt = t.

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

xt = (L)t, где (L) =-1(L), т.е. (L) Ч оператор, обратный оператору (L).

Удобным и полезным инструментом для изучения процессов авторегрессии яв ляется характеристический многочлен (характеристический полином) p (z) =1 - 1z - 2z2 -... - pzp =1 - jzj j= и связанное с ним характеристическое уравнение 1 - 1z - 2z2 -... - pzp =0.

432 Глава 14. Линейные стохастические модели ARIMA Как мы увидим в дальнейшем, от того, какие корни имеет характеристическое уравнение, зависят свойства процесса авторегрессии, в частности, является ли процесс стационарным или нет.

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

Процесс Маркова Процессом Маркова (марковским процессом) называется авторегрессионный процесс первого порядка, AR(1):

xt = xt-1 + t, (14.8) где t представляет собой белый шум, который не коррелирует с xt-1. Зд есь мы упростили обозначения, обозначив = 1.

Найдем необходимые условия стационарности марковского процесса. Предпо ложим, что процесс {xt} слабо стационарен. Тогда его первые и вторые момен ты неизменны. Находя дисперсии от обеих частей (14.8), получим, учитывая, что cov(xt-1, t) =0:

var(xt) =2var(xt-1) +var(t) или 2 2 x = 2x +.

Ясно, что при || 1, с учетом > 0, правая часть этого равенства должна быть больше левой, что невозможно. Получаем, что у стационарного марковского процесса || < 1.

Пусть, с другой стороны, || < 1. Представим xt через белый шум { t }. Это можно осуществить с помощью последовательных подстановок по формуле (14.8):

xt = xt-1 + t = (xt-2 + t-1) +t = 2xt-2 + t-1 + t, потом xt = 2(xt-3 + t-3) +t-1 + t = 3xt-3 + 2t-3 + t-1 + t, ит.д.

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

xt = t + t-1 + 2t-2 +....

14.3. Процессы авторегрессии Это же представление можно получить с использованием оператора лага. Урав нение (14.8) запишется в виде (1 - L)xt = t.

Применив к обеим частям уравнения (1 - L)-1, получим xt =(1 - L)-1t =(1 +L + 2L2 +... )t = t + t-1 + 2t-2 +....

В терминах модели линейного фильтра (14.1) для марковского процесса i = i. Поэтому |i| = ||i = <, (14.9) 1 -|| i=0 i= т.е. условие стационарности модели линейного фильтра (14.2) выполняется при || < 1.

Можно сделать вывод, что условие стационарности процесса Маркова имеет следующий вид:

|| < 1.

Свойства стационарного процесса AR(1):

1) Если процесс xt слабо стационарен, то его математическое ожидание неиз менно, поэтому, беря математическое ожидание от обеих частей (14.8), получим E(xt) =E(xt), откуда E(xt) =0.

Если добавить в уравнение (14.8) константу:

xt = + xt-1 + t, то E(xt) = + E(xt) и E(xt) =.

1 - 2) Найдем дисперсию процесса Маркова, используя полученное выше уравне 2 2 ние x = 2x + :

var(xt) =0 = x =. (14.10) 1 - 434 Глава 14. Линейные стохастические модели ARIMA Можно также применить общую формулу для автоковариаций в модели линей ного фильтра (14.3) с k =0:

2 2 2 0 = x = i = (1 + 2 + 4 +... ) =.

1 - i= При || 1 дисперсия процесса {xt}, вычисляемая по этой формуле, неограни ченно растет.

3) Коэффициент автоковариации k-го порядка по формуле (14.3) равен 2 k = ii+k = ii+k = i=0 i= 2 = 2i+k = k 2i = k.

1 - i=0 i= Можно вывести ту же формулу иным способом. Ошибка t некоррелирована не только с xt-1, но и с xt-2, xt-3 и т.д. Поэтому, умножая уравнение процесса (14.8) на xt-k при k > 0 и беря математическое ожидание от обеих частей, получим E(xtxt-k) =E(xt-1xt-k) или k = k-1, k > 0.

Таким образом, учитывая (14.10), получим k = 0k = k. (14.11) 1 - 4) Коэффициент автокорреляции, исходя из (14.11), равен:

k k = = k.

При 0 <1 автокорреляционная функция имеет форму затухающей экс поненты (рис. 14.1а), при -1 <0 Ч форму затухающей знакопеременной экспоненты (рис. 14.1б).

Если >1, процесс Маркова превращается во взрывной процесс.

В случае =1 имеет место так называемый процесс случайного блуждания, который относится к разряду нестационарных.

14.3. Процессы авторегрессии а) б) k k 1 k k 0 < 1 < -1 < 1 < - Рис. 14. Процесс Юла Процессом Юла называют авторегрессию второго порядка AR(2):

xt = 1xt-1 + 2xt-2 + t, (14.12) или, через лаговый оператор:

(1 - 1L - 2L2)xt = t.

Для стационарности процесса авторегрессии AR(2) необходимо, чтобы корни 1, 2 характеристического уравнения 1 - 1z - 2z2 =0, которые, вообще говоря, могут быть комплексными, находились вне единичного круга на комплексной плоскости, т.е. |1| > 1, |2| > 1.

Неформально обоснуем условия стационарности AR(2), разложив характеристиче ский полином (z) на множители (по теореме Виета 12 = - ):

z z (z) =-2(1 - z)(2 - z) = 1 - 1 - =(1 - G1z)(1 - G2z), 1 1 где мы ввели обозначения G1 = и G2 =.

1 Такое разложение позволяет представить уравнение AR(2) ввид е:

(1 - G1L)(1 - G2L)xt = t. (14.13) 436 Глава 14. Линейные стохастические модели ARIMA Введем новую переменную:

vt =(1 - G2L)xt, (14.14) тогда уравнение (14.13) примет вид:

(1 - G1L)vt = t.

Видим, что ряд vt является процессом Маркова:

vt = G1vt-1 + t.

Для того чтобы процесс vt был стационарным, коэффициент G1 по модулю должен быть меньше единицы, т.е. |1| > 1.

С другой стороны, из (14.14) следует, что xt имеет вид процесса Маркова с ошиб кой vt:

xt = G2xt-1 + vt.

Условие стационарности такого процесса имеет аналогичный вид: |2| > 1.

Приведенные рассуждения не вполне корректны, поскольку vt не является белым шумом. (Укажем без доказательства, что из стационарности vt следует стационар ность xt.) Кроме того, корни 1, 2 могут быть, вообще говоря, комплексными, что требует изучения свойств комплексного марковского процесса. Более коррект ное обоснование приведенного условия стационарности будет получено ниже для общего случая AR(p).

Условия стационарности процесса AR(2), |1| > 1, |2| > 1, можно переписать в эквивалентном виде как ограничения на параметры уравнения авторегрессии:

2 + 1 < 1, 2 - 1 < 1, (14.15) - 1 <2 < 1.

Проверим эти условия. Для этого рассмотрим два случая.

1) Пусть корни характеристического уравнения вещественные, то есть 2 +42 0. Тогда для выполнения условий |1| > 1, |2| > 1 необходимым - требованием является |12| > 1 или > 1. В таком случае один из корней обязательно лежит вне отрезка [-1, 1]. Для того чтобы и второй корень не попал в этот отрезок, необходимо и достаточно, чтобы значения характеристического по линома (L) =1 - 1L - 2L2 вточках -1 и 1 были одного знака. Это условие можно описать неравенством:

(-1) (1) > 0, или (1 + 1 - 2)(1 - 1 - 2) > 0.

14.3. Процессы авторегрессии a) 2 б) 1 Ц1 - 1 1 Ц1 - Рис. 14. Таким образом, случай вещественных корней описывается системой:

2 +42 0, |2| < 1, (1 + 1 - 2)(1 - 1 - 2) > 0.

Если связать 1 и 2 с координатными осями, то область, соответствующую данной системе, можно изобразить на рисунке (см. рис. 14.2а).

2) Если корни комплексные, то они имеют одинаковую абсолютную величину:

|1| = |2| = -.

И тогда вместе с отрицательностью дискриминанта достаточно условия:

|2| < 1 (область решений см. на рис. 14.2б).

Если объединить случаи 1) и 2), то общее решение как раз описывается системой неравенств (14.15) и соответствующая область на координатной плоскости пред ставляет собой треугольник, ограниченный прямыми:

2 = 1 +1, 2 = -1 +1, 2 = -1.

Автокорреляционная и автоковариационная функция AR(p) Для стационарного процесса авторегрессии:

xt = 1xt-1 + 2xt-2 + + pxt-p + t (14.16) 438 Глава 14. Линейные стохастические модели ARIMA можно вывести формулу автокорреляционной функции. Умножив обе части урав нения на xt-k:

xt-kxt = 1xt-kxt-1 + 2xt-kxt-2 +... + pxt-kxt-p + xt-kt, и перейдя к математическим ожиданиям, получим уравнение, связывающее коэф фициенты автоковариации различного порядка:

k = 1k-1 + 2k-2 +... + pk-p, k > 0. (14.17) Это выражение является следствием того, что соответствующие кросс ковариации между процессом и ошибкой равны нулю: E(xt-kt) =0 при k >0, т.к. xt-k может включать лишь ошибки j для j t - k.

Делением уравнения (14.17) на 0 получаем важное рекуррентное соотноше ние для автокорреляционной функции:

k = 1k-1 + 2k-2 +... + pk-p, k > 0. (14.18) Подставляя в выражение (14.18) k =1,..., p, получаем, с учетом симметрич ности автокорреляционной функции, так называемые уравнения ЮлаЧУокера (YuleЧWalker) для AR(p):

1 = 1 + 21 +... + pp-1, 2 = 11 + 2 +... + pp-2, (14.19) p = 1p-1 + 2p-2 +... + p.

Мы имеем здесь p линейных уравнений, связывающих p автокорреляций, 1,..., p. Из этой системы при данных параметрах можно найти автокорреляции.

С другой стороны, при данных автокорреляциях из уравнений ЮлаЧУокера можно найти параметры 1,..., p. Замена теоретических автокорреляций выборочны ми дает метод оценивания параметров процесса AR(p).

В частности, для процесса Юла получим из (14.19) 1 1 =, 2 = + 2.

1 - 2 1 - Зная 1,..., p, все последующие автокорреляции k ( k >p) можем найти по рекуррентной формуле (14.18).

Для нахождения автоковариационной функции требуется знать 0, дисперсию процесса xt. Если умножить обе части (14.16) на t и взять математические ожи дания, получим, что E(txt) = E(2) =. Далее, умножая обе части (14.16) t на xt и беря математические ожидания, получим 0 = 11 + 22 +... + pp +.

14.3. Процессы авторегрессии Значит, если известны автокорреляции, то дисперсию xt можно вычислять по формуле:

0 = x =. (14.20) 1 - 11 - 22 -... - pp Автоковариации затем можно вычислить как j = jx.

С учетом полученного дополнительного уравнения можно записать вариант уравнений ЮлаЧУокера для автоковариаций:

0 = 11 + 22 +... + pp +, 1 = 10 + 21 +... + pp-1, 2 = 11 + 20 +... + pp-2, (14.21) p = 1p-1 + 2p-2 +... + p0.

В этой системе имеется p+1 уравнений, связывающих p+1 автоковариацию, что позволяет непосредственно вычислять автоковариации при данных параметрах.

Заметим, что 14.17 и 14.18 имеют вид линейных однородных конечно разностных уравнений, а для подобных уравнений существует общий метод на хождения решения. Решив уравнение 14.18, можно получить общий вид автокор реляционной функции процесса авторегрессии. Проведем это рассуждение более подробно для процесса Юла, а затем рассмотрим общий случай AR(p).

Вывод формулы автокорреляционной функции процесса Юла Из соотношения (14.18) для p =2 получаем:

k = 1k-1 + 2k-2, k > 0. (14.22) или, используя лаговый оператор, который в данном случае действует на k, (L)k =0, гд е (z) = 1 - 1z - 2z2 Ч характеристический полином про цесса Юла. Как мы видели, этот характеристический полином можно представить ввид е:

(z) =(1 - G1z)(1 - G2z), 1 где G1 = и G2 =, а 1, 2 Ч корни характеристического уравнения.

1 Таким образом, k удовлетворяет уравнению:

(1 - G1L)(1 - G2L)k =0. (14.23) 440 Глава 14. Линейные стохастические модели ARIMA Найдем общее решение этого уравнения. Введем обозначение:

k =(1 - G2L)k. (14.24) Для k имеем (1 - G1L)k =0 или k = G1k-1. Поэтому k = G1k-1 = = Gk-11.

В свою очередь, из (14.24), с учетом того, что 0 =1 и 1 = 1(1 - 2), следует, что 1 = 1 - G20 = - G2.

1 - Поскольку, по теореме Виета, G1 + G2 = 1, G1G2 = -2, получаем выра жение для 1 через корни характеристического уравнения:

G1 + G2 G1(1 - G2) 1 = - G2 =. (14.25) 1+G1G2 1+G1G Возвращаясь к формуле (14.24), имеем, исходя из рекуррентности соотноше ния:

k = G2k-1 + k = G2(G2k-2 + k-1) +k =... = k- = Gk + Gk-11 + Gk-22 +... + G2k-1 + k = Gk + Gk-1-ss+1.

2 2 2 2 s= Подставляя сюда k = Gk-11, получим s k- G1 (G1/G2)k - k = Gk + 1Gk-1 = Gk + 1Gk-1 = 2 2 2 G2 (G1/G2) - s= Gk - Gk 1 = Gk + 1 1 2 = Gk + 1 - Gk.

G1 - G2 G1 - G2 1 G1 - G2 Таким образом, общее решение уравнения (14.22) имеет вид:

k = A1Gk + A2Gk, (14.26) 1 где коэффициенты A1 и A2 вычисляются по формулам:

G1(1 - G2) G2(1 - G2) 2 A1 =, A2 = -, (14.27) (G1 - G2)(1 + G1G2) (G1 - G2)(1 + G1G2) причем A1 + A2 =1.

14.3. Процессы авторегрессии 2 k k k k 3 k k Рис. 14. В стационарных процессах корни характеристического уравнения лежат вне единичного круга. Следовательно, |G1| < 1 и |G2| < 1, и автокорреляцион ная функция состоит из совокупности затухающих экспонент, что на рисунке 14. соответствует областям 1, 2, 3 и 4, лежащим выше параболической границы 2 +42 0.

При этом, если оба корня положительны или доминирует по модулю от рицательный корень (соответственно, положительное G), автокорреляционная функция затухает, асимптотически приближаясь к экспоненте (области 1 и 4 на рис. 14.3). Когда же оба корня отрицательны или доминирует по модулю поло жительный корень (или отрицательное G), автокорреляционная функция затухает по экспоненте знакопеременно (области 2 и 3 на рис. 14.3).

Если корни разного знака и совпадают по модулю, то затухание k происходит в области положительных значений, но имеет колебательный характер:

G1 = -G2, следовательно A1 = A2 =0.5 и 0, если k Ч нечетное, k = Gk, если k Ч четное.

Рассмотрим случай, когда корни 1 = G-1 и 2 = G-1 Чкомплексные.

1 Тогда автокорреляционная функция процесса авторегрессии второго порядка будет представлять собой затухающую синусоиду:

dk sin(k + ) k =, sin 442 Глава 14. Линейные стохастические модели ARIMA 1 1+d где d = -2, =arccos, =arctg tg.

2 -2 1 - d Подтвердим справедливость формулы. Поскольку коэффициенты характеристи ческого многочлена Ч действительные числа, то 1 и 2 Ч комплексно сопряженные, поэтому G1 и G2 тоже являются комплексно-сопряженными. Далее, G1, как и любое комплексное число, можно представить в виде:

G1 = dei = d cos + i d sin, где d = |G1| ( = |G2|) Чмод уль, а Ч аргумент комплексного числа G1. Соот ветственно, для G2 как для сопряженного числа верно представление:

G2 = de-i = d cos - i d sin.

Для нахождения d воспользуемся тем, что -2 = G1 G2 = dei de-i = d2.

ei + e-i G1 + G Значит, d = -2. Кроме того, cos = =, поэтому, учиты 2 2d вая, что G1 + G2 = 1, получаем:

cos =.

2 - Подставим теперь G1 = dei и G2 = de-i в выражения (14.27) для A1 и A2:

dei(1 - d2e-2i) ei - d2e-i A1 = = = (dei - de-i)(1 + d2eie-i) (ei - e-i)(1 + d2) 1+d 1+i tg (1 - d2)cos + i(1 + d2)sin 1 - d = =.

2i(1 + d2)sin 2i Пусть Ч такой угол, что 1+d tg = tg.

1 - d Тогда 1+i tg cos + i sin ei A1 = = =.

2i 2i sin ei - e-i Отсюда найдем A2:

e-i A2 =1 - A1 = -.

ei - e-i 14.3. Процессы авторегрессии Несложно увидеть, что A1 и A2 являются комплексно-сопряженными.

Подставим найденные A1 и A2 в(14.26):

ei e-i k = A1Gk + A2Gk = dkeik - dke-ik = 1 ei - e-i ei - e-i ei(k+) - e-i(k+) = dk.

ei - e-i Таким образом, подтверждается, что автокорреляционную функцию можно записать вформе dk sin(k + ) k =.

sin Нахождение автокорреляционной функции AR(p) с помощью решения конечно-разностного уравнения Аналогичным образом можно изучать автокорреляционную функцию процес са AR(p) при произвольном p. Запишем уравнение (14.18) с помощью лагового оператора, действующего на k:

(L)k =0, k > 0, (14.28) где (L) =1 - 1L - 2L2 -... - pLp.

Рассмотрим характеристическое уравнение:

(z) =1 - 1z - 2z2 -... - pzp =0.

Пусть i (i =1,..., p) Ч корни этого уравнения. Мы будем предполагать, что все они различны. Характеристический многочлен (z) можно разложить следующим образом:

p (z) =-p (i - z).

i= Обозначим через Gi значения, обратные корням характеристического урав нения: Gi =. Тогда, учитывая, что 1 2... p = - и соответственно p i G1 G2... Gp = -p, имеем:

p (1 - Giz) p p (z) =-p - z = -p i=1 = (1 - Giz).

p Gi i= Gi i= i= 444 Глава 14. Линейные стохастические модели ARIMA Исходя из этого, перепишем уравнение (14.28) в виде:

(1 - G1L)(1 - G2L)... (1 - GpL)k =0. (14.29) Из теории конечно-разностных уравнений известно, что если все корни i различны, то общее решение уравнения (14.18) имеет вид:

k = A1Gk + A2Gk +... + ApGk, k > -p, (14.30) 1 2 p где Ai Ч некоторые константы, в общем случае комплексные. (Обсуждение ре шений линейных конечно-разностных уравнений см. в Приложении A.4.) Проверим, что это действительно решение. Подставим k в (14.29):

p p (1 - GiL)k = (1 - GiL)(A1Gk + A2Gk +... + ApGk) = 1 2 p i=1 i= p p p = (1 - GiL)A1Gk + (1 - GiL)A2Gk +... + (1 - GiL)ApGk =0.

1 2 p i=1 i=1 i= Для доказательства использовался тот факт, что LGk = Gk-1 и поэтому j j p (1 - GiL)Gk = (1 - GiL)(1 - GjL)Gk =0.

j j i=1 i =j Формулы для коэффициентов A1,..., Ap можно получить из условий:

0 =1, k = -k, откуда p p p Ai Ai =1 и AiGk =, k =1,..., p- 1.

i Gk i i=1 i=1 i= Другой способ состоит в том, чтобы вычислить 1,..., p-1 из уравнений ЮлаЧ Уокера (14.19), а затем составить на основе (14.30) при k =0,..., p- 1 систему линейных уравнений, откуда и найти A1,..., Ap.

Если все корни характеристического уравнения удовлетворяют условию |i| > 1, то |Gi| < 1 i и все слагаемые в (14.30) затухают с ростом k. Если же для какого-то корня выполнено |i| < 1, то (при условии, что Ai =0) соответ ствующее слагаемое луходит на бесконечность. Если |i| =1, то соответствую щее слагаемое не затухает. Из этих рассуждений следует условие стационарности AR(p) Ч все корни соответствующего характеристического уравнения по модулю должны быть больше единицы.

14.3. Процессы авторегрессии Если корень i = G-1 действителен, элемент AiGk в (14.30) убывает с ро i i стом k экспоненциально, коль скоро |i| > 1. Если же есть пара комплексно сопряженных корней i = G-1, j = G-1, то соответствующие коэффициен i j ты Ai, Aj также будут сопряженными и в составе автокорреляционной функции появится экспоненциально затухающая синусоида (см. вывод автокорреляционной функции процесса Юла).

Таким образом, из соотношения (14.30) следует, что в общем случае автокорре ляционная функция стационарного процесса авторегрессии является комбинацией затухающих экспонент и затухающих синусоид.

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

На практике же для коротких рядов различительная сила автокорреляционной функции не очень высока. Однако часто изучение автокорреляционной функции является хорошим заделом исследования системы.

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

Спектр стационарного процесса авторегрессии В главе 13 мы определили спектральную плотность стационарного процесса как косинус-преобразование Фурье автоковариационной функции (13.29):

p(f) =2 0 +2 k cos 2fk. (14.31) k= Из этой общей формулы найдем спектральную плотность для стационарного процесса Маркова ( |1| < 1). Автоковариационная функция этого процесса имеет вид:

k k =.

1 - Подставляя эти автоковариации в формулу (14.31), получим p(f) = 1+2 k cos 2fk.

1 - k= Воспользовавшись представлением косинуса через комплексную экспоненту (формулами Эйлера) Ч 2cos2fk = ei2fk + e-i2fk, после несложных преобра 446 Глава 14. Линейные стохастические модели ARIMA зований получим:

2 k 2 k p(f) = 1ei2f + 1e-i2f - 1, 1 - k=0 k= т.е.

p(f) = (1z)k + (1z-1)k - 1, 1 - k=0 k= где мы ввели обозначение z = ei2f.

По формуле бесконечной геометрической прогрессии 2 1 p(f) = + - 1.

1 - 2 1 - 1z 1 - 1z- Произведение знаменателей двух дробей можно записать в разных формах:

(1 - 1z)(1 - 1z-1) =1 - 1(z + z-1) +2 = =1 - 21 cos 2f + 2 = 1 - 1e-i2f.

Приведя к общему знаменателю, получим 2 1 - 1z-1 +1- 1z - 1+1(z + z-1) - p(f) = = 1 - 2 (1 - 1z)(1 - 1z-1) 2 1 - =.

1 - 2 (1 - 1z)(1 - 1z-1) Таким образом, спектральная плотность марковского процесса равна p(f) =, 1 - 21 cos 2f + или, в другой форме, p(f) =.

1 - 1e-i2f Ошибку t авторегрессии произвольного порядка AR(p) можно выразить в виде линейного фильтра от yt:

p t = xt - jxt-j, j= 14.3. Процессы авторегрессии поэтому для вычисления спектра авторегрессионного процесса можно воспользо ваться общей формулой (14.4), характеризующей изменение спектра при приме нении линейного фильтра.

Спектральная плотность белого шума t равна 2. Применение формулы (14.4) дает p 2 = p(f) 1 - je-i2fj, j= откуда p(f) =.

1 - 1e-i2f - 2e-i4f - -pe-i2pf В частном случае процесса Юла формула спектральной плотности имеет вид:

p(f) = = 1 - 1e-i2f - 2e-i4f =.

1+2 + 2 - 21(1 - 2)cos2f - 22 cos 4f 1 Разложение Вольда и условия стационарности процессов авторегрессии Как уже говорилось, модель AR(p) можно записать в виде модели линейного фильтра:

xt = it-i = (L)t, i= где (L) =-1(L). Если процесс авторегрессии стационарен, то это разложение Вольда такого процесса.

Найдем коэффициенты модели линейного фильтра i процесса авторегрессии.

Для этого в уравнении p (L)(L) = iLi 1 - jLj = i=0 j= 448 Глава 14. Линейные стохастические модели ARIMA приравняем коэффициенты при одинаковых степенях L. Получим следующие урав нения:

0 = 1, 1 - 01 = 0, 2 - 11 - 02 = 0,...

p - p-11 -... - 1p-1 - 0p = 0, p+1 - p1 -... - 2p-1 - 1p = 0,...

Общая рекуррентная формула имеет следующий вид:

p i = ji-j, i > 0, (14.32) j= где 0 =1 и i =0 при i <0.

Это разностное уравнение, которое фактически совпадает с уравнением для автокорреляций (14.18). Соответственно, если все корни характеристического уравнения различны, то общее решение такого уравнения такое же, как указа но в (14.30), т.е.

i = B1Gi + B2Gi +... + BpGi, i > -p, 1 2 p где Bi Ч некоторые константы. Коэффициенты Bi можно вычислить, исходя из известных значений i при i 0.

Очевидно, что если |Gj| < 1 j, то все слагаемые здесь экспоненциально затухают, и поэтому ряд, составленный из коэффициентов i сходится абсолютно:

|i| <.

i= Таким образом, указанное условие гарантирует, что процесс авторегрессии явля ется стационарным. Это дополняет вывод, полученный при анализе автокорреля ционной функции.

Итак, условием стационарности процесса AR(p) является то, что корни i характеристического уравнения лежат вне единичного круга на комплексной плос кости.

Для процесса AR(2) имеем два уравнения для коэффициентов B1, B2:

B1 B + = -1 =0, G1 G B1 + B2 = 0 =1, 14.3. Процессы авторегрессии G1 G откуда B1 = и B2 = -. Таким образом, в случае процесса Юла G1-G2 G1-G имеем Gi+1 - Gi+ 1 i =, i > -1, G1 - G и xt = (Gi+1 - Gi+1)t-i.

G1 - G2 i=0 1 Оценивание авторегрессий Термин авторегрессия для обозначения модели (14.7) используется потому, что она фактически представляет собой модель регрессии, в которой регрессорами служат лаги изучаемого ряда xt. По определению авторегрессии ошибки t явля ются белым шумом и некоррелированы с лагами xt. Таким образом, выполнены все основные предположения регрессионного анализа: ошибки имеют нулевое ма тематическое ожидание, некоррелированы с регрессорами, не автокоррелированы и гомоскедастичны. Следовательно, модель (14.7) можно оценивать с помощью обычного метода наименьших квадратов.

Отметим, что при таком оценивании p начальных наблюдений теряются. Пусть имеется ряд x1,..., xT. Тогда регрессия в матричной записи будет иметь следую щий вид:

xp+1 xp x1 p+ xp+2 xp+1 x2 p+.

= +.

....

p xT xT -1 xT -p T Как видим, здесь используется T - p наблюдений.

Оценки МНК параметров авторегрессии равны = M-1m, где = =(1, 2,..., p), аматрицы M и m, как нетрудно увидеть, фактически состоят из выборочных автоковариаций ряда xt. Отличие от стандартных выборочных ав токовариаций состоит в том, что используются не все наблюдения.

Можно рассматривать данную регрессию как решение уравнений ЮлаЧ Уокера для автоковариаций (14.21) (или, что эквивалентно, уравнений для автокор реляций (14.19)), где теоретические автоковариации заменяются выборочными.

450 Глава 14. Линейные стохастические модели ARIMA Действительно, уравнения ЮлаЧУокера (14.21) без первой строки записываются ввид е =, где 1 1 1 2 p- 2 1 1 1 p- =, =,.....

......

.....

p p-1 p-2 p-3 откуда =-1.

Замена теоретических значений автоковариаций k выборочными автоковари ациями ck позволяет найти параметры процесса авторегрессии. Ясно, что при этом можно использовать и стандартные формулы выборочных автоковариаций. Тем са мым, мы получаем еще один из возможных методов оценивания авторегрессий Ч метод моментов.

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

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

Пусть {xt} Ч произвольный стационарный ряд (не обязательно авторегрес сия) и j Ч его автокорреляции. Применим к нему уравнения ЮлаЧУокера (14.19), как если бы процесс представлял собой авторегрессию k-го порядка, и найдем по автокорреляциям коэффициенты. Если обозначить j-й коэффици ент уравнения авторегрессии порядка k через kj, то уравнения ЮлаЧУокера 14.3. Процессы авторегрессии принимают вид:

1 1 2... k-1 k1 1 1 1... k-2 k2 2 1 1... k-3 k3 = 3.

......

.......

......

k-1 k-2 k-3... 1 kk k Частная автокорреляция k-го порядка определяется как величина kk, полученная из этих уравнений.

Решение этих уравнений соответственно для k = 1, 2, 3 дает следующие результаты (здесь используется правило Крамера):

1 1 1 1 1 1 2 2 - 2 1 11 = 1, 22 = =, 33 =.

- 1 1 1 1 1 1 2 1 Частная автокорреляционная функция рассматривается как функция частной автокорреляции от задержки k, гд е k =1, 2,....

Для процесса авторегрессии порядка p частная автокорреляционная функция {kk} будет ненулевой для k p и равна нулю для k >p, то есть обрывается на задержке p.

Значение выборочного частного коэффициента автокорреляции kk вычисля ется как МНК-оценка последнего коэффициента в уравнении авторегрессии AR(k).

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

452 Глава 14. Линейные стохастические модели ARIMA 14.4. Процессы скользящего среднего Другой частный случай модели линейного фильтра, широко распространенный в анализе временных рядов, Ч модель скользящего среднего, когд а xt линейно зависит от конечного числа q предыдущих значений :

xt = t - 1t-1 - 2t-2 -... - qt-q. (14.33) Модель скользящего среднего q-го порядка обозначают MA(q) (от английского moving average).

Данную модель можно записать и более сжато:

xt = (L)t, через оператор скользящего среднего:

(L) =1 - 1L - 2L2 -... - qLq. (14.34) Легко видеть, что процесс MA(q) является стационарным без каких-либо огра ничений на параметры j.

Действительно, математическое ожидание процесса E(xt) =0, а дисперсия 2 2 2 0 =(1 +1 + 2 +... + q), т.е. равна дисперсии белого шума, умноженной на конечную величину (1 + 1 + 2 + 2 +... + q).

Остальные моменты второго порядка (k, k) также от времени не зависят.

Автоковариационная функция и спектр процесса MA(q) Автоковариационная функция MA(q) (-k + 1k+1 +... + q-kq), k =1, 2,..., q, k = (14.35) 0,k > q.

В частном случае для MA(1) имеем:

2 0 =(1 +1), 1 = -1, k =0, k > 1, 14.4. Процессы скользящего среднего и автоковариационная матрица, соответствующая последовательности x1, x2,..., xT, будет иметь следующий трехдиагональный вид:

1+1 -1 0 -1 1+1 -1 = 0 -1 1+1 0.

....

.....

....

0 0 0 1+ В общем случае автоковариационная матрица процесса скользящего среднего порядка q имеет q ненулевых поддиагоналей и q ненулевых наддиагоналей, все же остальные элементы матрицы равны нулю.

Автокорреляционная функция имеет вид:

-k + 1k+1 +... + q-kq, k =1, 2,..., q, 1+1 +... + q k = (14.36) 0,k > q.

Таким образом, автокорреляционная функция процесса MA(q) обрывается на задержке q, и в этом отличительная особенность процессов скользящего среднего.

С другой стороны, частная автокорреляционная функция, в отличие от авторе грессий, не обрывается и затухает экспоненциально. Например, для MA(1) частная автокорреляционная функция имеет вид 1 - p -1 2p+2.

1 - Ясно, что модель скользящего среднего является частным случаем модели ли нейного фильтра (14.1), где j = -j при j = 1,..., q и j = 0 при j > q.

Фактически модель линейного фильтра является моделью MA().

Формула спектра для процесса скользящего среднего следует из общей форму лы для модели линейного фильтра (14.6):

p(f) =2 - 1e-i2f - 2e-i4f -... - qe-i2qf.

454 Глава 14. Линейные стохастические модели ARIMA Соответственно, для MA(1):

2 2 p(f) =2 1 - 1e-i2f =2(1 + 1 - 21 cos 2f);

для MA(2):

p(f) =2 - 1e-i2f - 2e-i4f = 2 2 =2 1+1 + 2 - 21(1 - 2)cos2f - 22 cos 4f.

Обратимость процесса MA(q) Авторегрессию, как мы видели выше, можно представить как MA(). С другой стороны, процесс скользящего среднего можно представить в виде AR().

Рассмотрим, например, MA(1) (будем для упрощения писать вместо 1):

xt = t - t-1, (14.37) Сдвигом на один период назад получим t-1 = xt-1 + t-2 и подставим в (14.37):

xt = t - xt-1 - 2t-2.

Далее, t-2 = xt-2 + t-3, поэтому xt = t - xt-1 - 2xt-2 - 3t-3.

Продолжая, получим на k-м шаге xt = t - xt-1 - 2xt-2 - -kxt-k - k+1t-k-1.

Если || < 1, то последнее слагаемое стремится к нулю при k. Переходя к пределу, получаем представление AR() для MA(1):

xt = - jxt-j + t. (14.38) j= С помощью лагового оператора можем записать это как (L)xt = t, где (L) =(1 - L)-1 = jLj.

j= 14.4. Процессы скользящего среднего В то время как процесс (14.37) стационарен при любом, процесс (14.38) стационарен только при || < 1. При || 1 веса -j в разложении (14.38) растут (при || =1 не меняются) по абсолютной величине по мере увеличения j.

Тем самым, нарушается разумная связь текущих событий с событиями в прошлом.

Говорят, что при || < 1 процесс MA(1) является обратимым, а при || 1 Ч необратимым.

В общем случае уравнение процесса MA(q) в обращенной форме можно запи сать как t = -1(L)xt = (L)xt = jxt-j.

j= Процесс MA(q) называется обратимым, если абсолютные значения весов j воб ращенном разложении образуют сходящийся ряд. Стационарным процесс MA(q) является всегда, но для того, чтобы он обладал свойством обратимости, параметры процесса должны удовлетворять определенным ограничениям.

Выведем условия, которым должны удовлетворять параметры 1, 2,..., q процесса MA(q), чтобы этот процесс был обратимым.

- Пусть Hi, i =1,..., q Ч корни характеристического уравнения (L) = (будем предполагать, что они различны). Оператор скользящего среднего (L) через обратные корни характеристического уравнения можно разложить на мно жители:

q (L) = (1 - HiL).

i= Тогда обратный к (L) оператор (L) можно представить в следующем виде:

q 1 Mi (L) =-1(L) = =. (14.39) q 1 - HiL i= (1 - HiL) i= Каждое слагаемое (14.39) можно, по аналогии с MA(1), представить в виде беско нечного ряда:

Mi j = Mi Hi Lj, i =1,..., q, 1 - HiL j= который сходится, если |Hi| < 1.

456 Глава 14. Линейные стохастические модели ARIMA Тогда процесс MA(q) в обращенном представлении выглядит как q j t = Mi Hi Ljxt, i=1 j= и он стационарен, если корни характеристического уравнения (L) =0 лежат вне единичного круга. Иными словами, MA(q) обладает свойством обратимости, если - для всех корней выполнено |Hi | > 1, т.е. |Hi| < 1 i. Если же для одного из корней |Hi| 1, то ряд не будет сходиться, и процесс MA(q) будет необратимым.

Для каждого необратимого процесса MA(q), у которого корни характеристи ческого уравнения не равны по модулю единице, существует неотличимый от него обратимый процесс того же порядка. Например, процесс MA(1) (14.37) с || > можно записать в виде xt = t - t-1, 1 - L где t = t является белым шумом. Мы не будем доказывать, что t 1 - 1/ L является белым шумом, поскольку это технически сложно. Вместо этого мы укажем на простой факт: пусть t Ч некоторый белый шум. Тогда процесс t- t-1 имеет такую же автоковариационную функцию, как и процесс xt, заданный уравнением 2 (14.37), если дисперсии связаны соотношением = 2. Для того чтобы в этом убедиться, достаточно проверить совпадение дисперсий и автоковариаций первого порядка (остальные автоковариации равны нулю).

В общем случае процесса MA(q), чтобы сделать его обратимым, требуется обратить все корни характеристического уравнения, которые по модулю меньше q единицы. А именно, пусть (L) = (1-HiL), Ч характеристический многочлен i= где |Hi| < 1 при i =1,..., m и |Hi| > 1 при i = m +1,..., q. Тогд а q m - (L) = (1 - HiL) (1 - Hi L) (14.40) i=1 i=m+ Ч характеристический многочлен эквивалентного обратимого процесса.

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

14.5. Смешанные процессы авторегрессии Ч скользящего среднего 14.5. Смешанные процессы авторегрессии Ч скользящего среднего ARMA (модель БоксаЧДженкинса) На практике иногда бывает целесообразно ввести в модель как элементы ав торегрессии, так и элементы скользящего среднего. Это делается для того, чтобы с использованием как можно меньшего числа параметров уловить характеристики исследуемого эмпирического ряда. Такой процесс называется смешанным процес сом авторегрессии Ч скользящего среднего и обозначается ARMA(p, q):

xt = 1xt-1 +... + pxt-p + t - 1t-1 -... - qt-q, (14.41) или, с использованием оператора лага, (1 - 1L - 2L2 -... - pLp)xt =(1 - 1L - 2L2 -... - qLq)t.

В операторной форме смешанная модель выглядит так:

(L)xt = (L)t, где (L) Ч оператор авторегрессии, (L) Ч оператор скользящего среднего.

Модель (14.41) получила название модели БоксаЧДженкинса, поскольку бы ла популяризирована Дж. Боксом и Г. Дженкинсом в их известной книге Анализ временных рядов [3]. Методология моделирования с помощью (14.41) получила название методологии БоксаЧДженкинса.

Автокорреляционная функция и спектр процесса ARMA(p, q) Рассмотрим, как можно получить автоковариационную и автокорреляционную функции стационарного процесса ARMA(p, q), зная параметры этого процесса.

Для этого умножим обе части уравнения (14.41) на xt-k, гд е k 0, и перейд ем к математическим ожиданиям:

E(xt-kxt) =1E(xt-kxt-1) +2E(xt-kxt-2) +... + pE(xt-kxt-p) + + E(xt-kt) - 1E(xt-kt-1) - 2E(xt-kt-2) -... - qE(xt-kt-q).

Обозначим через s кросс-ковариацию изучаемого ряда xt иошибки t с за держкой s, т.е.

s = E(xtt-s).

458 Глава 14. Линейные стохастические модели ARIMA Поскольку процесс стационарен, то эта кросс-ковариационная функция не зависит от момента времени t. В этих обозначениях E(xt-kt-j) =j-k.

Получаем выражение для автоковариационной функции:

k = 1k-1 +... + pk-p + -k - 11-k -... - qq-k. (14.42) Так как xt-k зависит только от импульсов, которые произошли до момента t - k, то j-k = E(xt-kt-j) =0 при j

Для того чтобы найти остальные нужные нам кросс-ковариации, 0,..., q, необходимо поочередно умножить все члены выражения (14.41) на t, t-1,..., t-q и перейти к математическим ожиданиям. В итоге получится следующая систе ма уравнений:

0 =, 1 = 10 - 1, 2 = 11 + 20 - 2,...

Общая формула для всех 1 s p имеет вид:

s = 1s-1 + + s0 - s.

При s >p (такой случай может встретиться, если p

Отсюда рекуррентно, предполагая и параметры и известными, найдем s.

Далее, зная s, по аналогии с уравнениями ЮлаЧУокера (14.21) по формуле (14.42) при k = 0,..., p с учетом того, что -k = k найдем автоковариа ции 0,..., p. Остальные автоковариации вычисляются рекуррентно по формуле (14.42).

Автокорреляции рассчитываются как k = k/0. Заметим, что если требуется найти только автокорреляции, то без потери общности можно взять ошибку t с единичной дисперсией: =1.

Если в уравнении (14.42) k >q, то все кросс-корреляции равны нулю, поэтому k = 1k-1 +... + pk-p, k > q. (14.43) 14.5. Смешанные процессы авторегрессии Ч скользящего среднего Поделив это выражение на 0, выводим уравнение автокорреляционной функции для k >q:

k = 1k-1 +... + pk-p, (14.44) или (L)k =0, k > q.

Таким образом, начиная с некоторой величины задержки, а точнее, когда q

По аналогии с AR(p) условия стационарности ARMA(p, q) определяются кор нями характеристического уравнения (L) =0: если эти корни лежат вне единич ного круга, то процесс стационарен.

В качестве примера рассмотрим процесс ARMA(1, 1):

xt = xt-1 + t - t-1, (14.45) или через лаговый оператор:

(1 - L)xt =(1 - L)t.

Процесс стационарен, если -1 <1, и обратим, если -1 <1.

Для вывода формулы автокорреляционной функции умножим (14.45) на xt-k и перейдем к математическим ожиданиям:

E(xt-kxt) =E(xt-kxt-1) +E(xt-kt) - E(xt-kt-1), или k = k-1 + -k - 1-k. (14.46) Исследуем поведение автоковариационной функции при различных значениях па раметра k.

Заметим, что граничные условия у AR(p) другие, поэтому автокорреляционные функции не будут совпадать.

460 Глава 14. Линейные стохастические модели ARIMA При k = 0 = -1 + 0 - 1. (14.47) Чтобы найти второе слагаемое, умножим уравнение процесса (14.45) на t и возьмем математическое ожидание:

0 = E(xtt) =E(xt-1t) +E(tt) - E(tt-1) =. (14.48) Аналогичным способом распишем E(xtt-1):

1 = E(xtt-1) =E(xt-1t-1) +E(t-1t) - E(t-1t-1) =( - ).

Равенство E(xt-1t-1) = подтверждается так же, как (14.48).

Итак, принимая во внимание, что k = -k, выражение для дисперсии запи сывается как 2 0 = 1 + - ( - ). (14.49) При k =1 равенство (14.46) преобразуется в 1 = 0 + -1 - 0.

Используя ранее приведенные доводы относительно математических ожиданий, стоящих в этом уравнении, имеем:

1 = 0 -. (14.50) При k k = k-1.

Выразим автоковариации в (14.49) и (14.50) через параметры модели и.

Получим систему уравнений 2 0 = 1 + - ( - );

1 = 0 - ;

и решим ее относительно 0 и 1.

Решение имеет вид:

1 - 2 + 0 =, 1 - ( - )(1 - ) 1 =.

1 - 14.5. Смешанные процессы авторегрессии Ч скользящего среднего k k k k k - k -1 1 k k k k - k k - - Рис. 14.4. График автокорреляционной функции процесса ARMA(1, 1) С учетом того, что k = k/0, получаем выражения для автокорреляционной функции процесса ARMA(1, 1):

( - )(1 - ) 1 =, 1 - 2 + k = k-11, k 2.

На рисунке 14.4 изображены графики автокорреляционной функции процесса ARMA(1, 1) при различных сочетаниях значений параметров и.

По аналогии с процессами AR(p) и MA(q) выводится формула спектра процесса ARMA(p, q). Пусть {t} Ч такой процесс, что t = t - 1t-1 -... - qt-q.

Тогда xt, описываемый уравнением (14.41), можно записать в виде xt = 1xt-1 +... + pxt-p + t.

462 Глава 14. Линейные стохастические модели ARIMA По формуле (14.4), с одной стороны, p(f) =px(f) =p(f) 1 - 1e-i2f - 2e-i4f -... - pe-i2pf, а с другой стороны, для процесса скользящего среднего {t} выполняется p(f) =.

1 - 1e-i2f - 2e-i4f -... - qe-i2qf Таким образом, получаем 1 - 1e-i2f - 2e-i4f -... - qe-i2qf p(f) =2. (14.51) - 1e-i2f - 2e-i4f -... - pe-i2pf Представление процесса ARMA в виде MA() и функция реакции на импульсы Так же, как и в случае авторегрессии, стационарный процесс ARMA можно записать в виде модели линейного фильтра, или, другими словами, скользящего среднего бесконечного порядка MA():

(L) xt = t = t + 1t-1 + 2t-2 +... = it-i = (L)t, (14.52) (L) i= где 0 =1.

Коэффициенты i представляют собой функцию реакции на импульсы для процесса ARMA, т.е. i является количественным измерителем того, как неболь шое изменение (лимпульс) в t влияет на x через i периодов, т.е. на xt+i, что можно символически записать как dxt+i i =.

dt Один из способов вычисления функции реакции на импульсы сводится к ис пользованию уравнения (z)(z) =(z), или (1 - 1z - 2z2 -... - pzp)(1 + 1z + 2z2 +... ) =(1 - 1z -... - qzq), 14.6. Модель ARIMA из которого, приравнивая коэффициенты в левой и правой частях при одинаковых степенях z, можно получить выражения для i.

Более простой способ состоит в том, чтобы продифференцировать по t урав нение ARMA-процесса, сдвинутое на i периодов вперед, p q xt+i = jxt+i-j + t+i - jt+i-j, j=1 j= p dxt+i dxt+i-j = j - i, dt j=1 dt где 0 = -1 и j =0 при j >q. Таким образом, получим рекуррентную формулу для i = dxt+i/dt:

p i = ji-j - i.

j= При расчетах по этой формуле следует положить 0 =1 и i =0 при i <0.

Если процесс ARMA является обратимым5, то полученное представление в виде MA() является разложением Вольда этого процесса.

14.6. Модель авторегрессии Ч проинтегрированного скользящего среднего ARIMA Характерной особенностью стационарных процессов типа ARMA(p, q) являет ся то, что корни i характеристического уравнения (L) =0 находятся вне еди ничного круга. Если один или несколько корней лежат на единичной окружности или внутри нее, то процесс нестационарен.

Теоретически можно предложить много различных типов нестационарных мо делей ARMA(p, q), однако, как показывает практика, наиболее распространен ным типом нестационарных стохастических процессов являются интегрированные процессы или, как их еще называют, процессы с единичным корнем. Единичным называют корень характеристического уравнения, равный действительной едини це: i =1.

Разложение Вольда необратимого процесса, у которого некоторые корни характеристического уравнения по модулю больше единицы, такое же, как у эквивалентного обратимого процесса. Ошибки однопериодных прогнозов, лежащие в основе разложения Вольда, при этом не будут совпадать сошибкамимод ели t.

464 Глава 14. Линейные стохастические модели ARIMA мнимая часть действительная часть - - Рис. 14.5. Корни характеристического уравнения процесса xt =2.8xt-1 - 3.1xt-2 +1.7xt-3 - 0.4xt-4 + t +0.5t-1 - 0.4t- на комплексной плоскости Рассмотрим в качестве примера следующий процесс ARMA(4, 2):

xt =2.8xt-1 - 3.1xt-2 +1.7xt-3 - 0.4xt-4 + t +0.5t-1 - 0.4t-2. (14.53) Характеристическое уравнение этого процесса имеет следующие корни: 1 =1 + + i, 2 = 1 - i, 3 = 1.25, 4 = 1. Все корни лежат за пределами единичного круга, кроме последнего, который является единичным. Эти корни изображены на рисунке 14.5.

Оператор авторегрессии этого процесса можно представить в следующем виде:

1 - 2.8L +3.1L2 - 1.7L3 +0.4L4 =(1 - 1.8L +1.3L2 - 0.4L3)(1 - L) = =(1 - 1.8L +1.3L2 - 0.4L3), где =1 - L Ч оператор первой разности.

Введем обозначение wt = xt = xt - xt-1. Полученный процесс {wt} явля ется стационарным процессом ARMA(3, 2), задаваемым уравнением:

wt =1.8wt-1 - 1.3wt-2 +0.4wt-3 + t +0.5t-1 - 0.4t-2. (14.54) В общем случае, если характеристическое уравнение процесса ARMA(p + d, q) содержит d единичных корней, а все остальные корни по модулю больше единицы, то d-я разность этого временного ряда wt = dxt = (L)(1 - L)dxt может быть представлена как стационарный процесс ARMA(p, q):

Pages: | 1 | ... | 4 | 5 | 6 | 7 | 8 | ... | 10 |

Книги, научные публикации

Blog

Оглавление Введение................................... 11 I Введение в социально-экономическую статистику 15 1. Основные понятия 17 1.1. Краткая историческая справка.. ...