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

Вычислительные технологии

Том 13, № 4, 2008

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

Б. Н. Луценко

Конструкт,орско-технологический институт вычислительной техники СО РАН,

Новосибирск, Россия e-mail: [email protected]

A new method of identification end estimation of parameters of the stochastic ARIMA model is considered. This method differs from the methods usually based on making use of correlation and partial correlation functions of the time series. It also makes possible to estimate greater dimension models with various season periods for subsequent forecasting of processes described using these models.

1. Формулировка задачи

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

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

2. Определения и теоретические предпосылки

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

количество параметров. Белый шум представляет собой последовательность независимых, одинаково распределенных случайных величин at с нулевым средним и постоянной дисперсией а2:

w ч w ч 2 / ч fa2 при т = 0,

Е(a) = 0; Е(at, ai+r) = а р(т) = j 0 при т

Обозначим сам случайный процесс, возможно подвергшийся предварительным преобразованиям, через Z(t). Отсчеты Z(t) эквидистантны, и в качестве временного дискрета принимается периодичность его регистрации. Если процесс стационарен и его среднее равно у, мы будем работать далее с центрированным процессом:

Zt = Zt - Е(Zt) = Zt - у.

При нестационарности из процесса может быть извлечена постоянная константа c и мы по-прежнему будем обозначать результат Zt:

ZZt = Zt - c.

Обозначим оператор сдвига на один такт назад через B (back), а вперед — через F (forward):

BZt = Zt—i; FZt = Zt+i-

Оператор взятия разности обозначим через V:

▽ Zt = (1 - B)Zt = Zt - Zt—1.

Оператором, обратным V, будет оператор суммирования (дискретный аналог интегрирования) .

Процесс авторегрессии (АР) порядка p описывается выражением вида

(1 - piB - ^2в2 - ... - wBp)Zt = at, или в более компактной записи:

%(B)Zt = at.

Процесс АР можно представить в иной форме, как

Zt = №P(B )] — 1at = Ф^ )at,

где Ф^) = £ фгBi; фо = 1.

i=0

Очевидно Ф(B)Ф(B) = 1, и, решая последовательно это уравнение, легко находятся все коэффициенты в разложении Ф^). В данной интерпретации процесс АР представляет собой выход линейного фильтра с передаточной функцией Ф^), на вход которого поступает белый шум.

Процесс АР будет стационарен, если корни его характеристического уравнения

Фр(х) = 1 - p1x - ... - ppxp = 0,

где x — комплексная переменная, будут лежать вне единичного круга в комплексной области. Выполнение условий стационарности обеспечивает и эргодичность процесса

(теорема Биркгофа—Хинчина [6]). Эргодичность чрезвычайно важна, поскольку реально приходится всегда иметь дело с ограниченным фрагментом единственной реализации процесса, по которому надлежит дать содержательные выводы о вероятностных параметрах получаемых оценок и прогноза.

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

^ = (1 - вг В - ... - вд В9 = ед (ВК

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

Процесс СС при конечной а2 будет всегда стационарен и обратим, т.е. представим сходящимся бесконечным рядом:

а = пгв

г=0

если корни его характеристического уравнения ед (ж) = 0 также лежат вне единичного круга в комплексной области.

И, наконец, процесс, объединяющий в себе качества АР и СС и именуемый процессом авторегрессии и скользящего среднего (АРСС):

фр (в)^ = ед (вк

И здесь Zt можно представить как выход линейного фильтра с передаточной функцией

[Фр(В )]-1ед (В ) = Ф(В ) = £ фгВг; ф = 1,

г=о

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

Фр(в )Ф(В) = ед (В).

Условия стационарности и обратимости процесса АРСС остаются прежними, т. е. присутствие оператора СС не отражается на условиях стационарности, а оператора АР — на условиях его обратимости.

Таким образом, однородный нестационарный процесс Zt, определяемый выражением

Фр(В) ▽ ^ = ед (B)at,

можно представить как результат кратного суммирования с соответствующими начальными условиями стационарного процесса:

= ▽

описываемого в свою очередь моделью АРСС:

Фp(B)wt = ед (В)аг.

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

3. Мультипликативная сезонная модель АРПСС

Модель задается своей структурой, в которую входят:

к — количество интервалов сезонности, включая для единообразия и единичный; для I = 1,... ,к;

Бг — значения интервалов сезонности;

рг — количество параметров АР для каждого интервала сезонности; ^ — количество параметров СС для каждого интервала сезонности; йг — количество сезонных разностей для каждого интервала сезонностей.

При этом оператор авторегрессии для периода сезонности Бг представим в виде

р'

ФРг (В8' ) = 1 ^ В*8', (1)

3 = 1

а оператор скользящего среднего —

<й

е« (в8' ) = 1 % в38'. (2)

3 = 1

В этих обозначениях сравнительно компактно мультипликативная модель АРПСС может быть представлена в виде

к к

Пфр- (В8' )▽% Ъ = Пе« (В8' К (3)

г=1 г=1

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

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

к

пв = т =£(р + Ф),

г=1

в = (^1,1, • • • , • • •, • • •, к,ркАд,... ,... ,9к,\,••• ,0кш) = (въ • • • (4)

А при использовании одношагового прогнозирования, основанного на разностной модели, удобнее ввести избыточные приведенные параметры модели, перемножив, соответственно, операторы в левой и правой частях выражения (3). Количество этих избыточных параметров для приведенного оператора АР будет равно

к

Пр = П(Рг + 1) - 1,

г=1

и для приведенного оператора СС соответственно —

к

П(* + !) - 1'

Пд =

г=1

Эти приведенные параметры ранжируют по степеням В. При этом, если обозначить приведенные параметры соответственно через и 0*, получим

к

ЦфРг(В3*) = £ №, ^ = 1, Го = 0, (5)

г=1 г=0

к П

ЦвЯг (В3 ) = £ 0*Вт, = 1, то = 0, (6)

и если обозначить

г=1 г=0

кк

^ = г е ^& х йг + 1,^), (7)

г=1 г=1

то исходная модель примет вид

и = £ 0**Вт*а*, (8)

г=0 г=0

похожий на модель АРСС, только теперь количество параметров избыточно, а степени операторов сдвига г и т^ следуют с неравномерным шагом.

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

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

Конечно, реально значимые р, ^ и ^ будут невелики.

Пример использования мультипликативной сезонной модели с тремя периодами сезонных изменений для прогнозирования описываемого ею процесса содержится в работе [7]. К сожалению, авторы не касаются вопросов идентификации модели, параметры которой они уже предполагают известными.

п

р

п

р

4. Схема решения задачи прогнозирования

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

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

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

Преобразованные таким образом данные подвергаются взятию сезонных разностей с переводом процесса АРПСС в АРСС. Следующий этап — идентификация модели по стационарному процессу АРСС.

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

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

5. Извлечение последовательности белого шума

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

Чтобы извлечь последовательность at из исходного процесса Wt на интервале t G (1, n), необходимо уже располагать значениями at и Wt при t < 0 , которых у нас изначально нет. Оценить их и позволяет предложенный М. Кенуем (M.H. Quenouille) так называемый метод складного ножа (jackknife). Суть операции состоит в том, что, используя заданный процесс Wt и задавшись значениями параметров модели, мы осуществляем последовательное одношаговое прогнозирование в прошлое процесса Wt с попутным вычислением последовательности белого шума et как разности между фактическими и прогнозируемыми на один шаг значениями Wt, поскольку реальные данные кроме инерционной компоненты, представляемой прогнозом, содержат и независимую от них случайную компоненту.

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

Ранее мы обращали внимание на то, что любой линейный процесс — АР, СС или АРСС — можно представить как выход линейного фильтра, на вход которого подается белый шум; когда это поступление прекращается, процесс постепенно (с учетом

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

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

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

При движении в отрицательном времени прогнозируемый процесс Ж* по мере угасания постепенно теряет признаки стохастичности (после прекращения действия сезонных компонент СС), перерождаясь в детерминированный процесс. При возвратном движении стохастичность процесса постепенно нарастает по мере приближения к £ = 0. Так что возможность извлечения а* и для £ < 0 вполне реальна.

Эта процедура позволяет погасить переходной процесс в а* до выхода а* в область положительных значений И для нас здесь это самое важное — обеспечить стационарный режим формирования а* на интервале прогнозной базы.

Конечно, для нас важно обеспечить стационарность а* для точки (параметров модели) в глобальном минимуме критерия или его окрестности и менее важно — в остальных.

Опишем теперь вкратце используемые формулы. Мы будем работать с сезонной моделью АРСС, представленной вариантом с избыточными приведенными параметрами (8) или в эквивалентном виде:

Пд Пр

= ^ 0гЧ_т - ^ ^+ а, £ > 0,

г=1 г=1

а при прогнозировании в обратном времени, т. е. после замены оператора В оператором

^, —

Пд Пр

^ = £ в*е*+тг - £ ^+ е*, £ < 0.

г=1 г=1

Прогноз для а* получается как условное математическое ожидание при заданных предшествующих моменту £ (или следующих за ним) значениях Ж* и а* (или е*), т.е., поскольку Е(а*) = Е(е*) = 0, полагается а* = 0 при £ > п и е* = 0 при £ < 1:

Пд Пр

= ^ в*а*_тг - ^ ^ , а* = Ж* - Ж* (9)

г=1 г=1

и при движении в обратном времени

Пд Пр

Ж = Е в*е*+тг - £ , е* = Ж* - Ж*. (10)

г=1 г=1

После завершения процедуры извлечения последовательности а* из Ж* вычисляется ее средний квадрат на интервале прогнозной базы:

в> = П£а?, (11)

г=1

а в последующем — и средняя дисперсия прогноза на интервале пролонгации.

6. Дихотомическое погружение в пространство параметров модели АРСС

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

Это метод выборочного зондирования всей области возможных значений параметров модели с постепенным отсевом неперспективных точек и связанных с ними областей и поэтапным отбором наиболее перспективных в плане минимизации критерия качества модели. Конечно, немыслимо детально исследовать пространства т-мерной размерности при т > 5. Поэтому, как и при планировании эксперимента, приходится экономно распоряжаться ресурсами. Из этих соображений и выбрано дихотомическое погружение, суть которого состоит в отборе для зондирования для каждой избранной точки из ее окрестности весьма ограниченного количества точек, названных нами дихотомическими по способу их размещения в этой окрестности. Количество этих точек при т-мерном пространстве параметров равно 2т. Они располагаются симметрично относительно центра области, являясь в свою очередь центрами 2т подобных же областей, но с уменьшенным в 2 раза по каждой из координат размером.

Для каждой из точек (а это вектор параметров модели АРСС) по заданному массиву

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

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

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

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

к

1 1 2 3

2 1 1 2 1 2 3

3 (—1, +1) (—2, +2) (—1, +1) (—3, +3) (—3, +1) (—1, +1)

4 2 4 2 6 4 2

5 0 0 0 0 —1 0

зондирования параметров модели АРСС. В строках 1-5 приводятся соответствующие данные:

1 — количество параметров АР или СС для одного периода сезонности;

2 — номера параметров по мере возрастания сдвига в прошлое;

3 — интервалы изменения параметров;

4 — размер стороны бруса для данного параметра;

5 — координата центра бруса по данному параметру.

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

Формирование координат дихотомических точек осуществляется по следующему алгоритму. После задания структуры модели определяется размерность вектора в ее параметров т = пв и вычисляется количество дихотомических точек т2 в окрестности любой точки т2 = 2т. Задается область возможных значений параметров модели в виде прямоугольного бруса с центром, выбираемым согласно строке 5 таблицы. Вводится диагональная матрица Вг с элементами-размерами сторон бруса Ь по каждому из т параметров модели:

Вг(г, г) = г = 1,..., т.

Если координаты некоторой ]-й точки, отобранной на 5-м этапе погружения для последующего обследования, обозначим во(5,^'), то координаты т2 дихотомических точек ее окрестности будут

в8О',к) = в0(5,^')+ X(к,о) ■ Вг ■ 2-1-*, к =1,...,т2. (12)

Здесь в0 и вэ — векторы строки размерности т, X(к, о) — к-я строка матрицы X (т2, т) относительных дихотомических координат с элементами +1 и — 1, получаемыми путем замены нулей на - 1 в бинарном представлении чисел (или номеров дихотомических точек) от 0 до т2 — 1.

По полученным координатам вэС?, к) точки с помощью метода складного ножа из последовательности ^ или извлекается последовательность аь, по которой вычисляется значение критерия. И эта операция производится для точек всех дихотомических окрестностей из отобранных на 5-м этапе погружения подмножества претендентов. Претенденты отбираются динамически заменой менее перспективных более перспективными, при этом участвуют и точки предшествующих этапов, если они конкурентоспособны. Количество оставляемых на достижение глобального минимума претендентов сокращается по мере погружения. Завершается процедура выбором единственной точки (вектора параметров АРСС), получаемой, возможно, усреднением координат точек, отобранных на последнем этапе.

7. Уточнение параметров модели методом модифицированного скорейшего спуска

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

Используя в качестве отправных точки, отобранные при дихотомическом погружении, на одном из ранних этапов мы попытались уточнить значения вектора параметров в, используя алгоритм скорейшего градиентного спуска [8]. Вскоре однако выяснилось, что этот алгоритм дает не очень значительное уточнение и зачастую останавливается, достигнув лишь локального минимума. Тогда мы обратились к алгоритму, использующему спуск по преобразованному градиенту [9]. Преобразование осуществлялось с помощью матрицы вторых частных производных (матрицы Гессе) от интегрального критерия качества модели по уточняемым параметрам. В качестве интегрального критерия была избрана средняя дисперсия прогноза на интервале пролонгации. Поначалу мы работали с оценкой дисперсии белого шума в качестве критерия.

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

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

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

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

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

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

^ = ^ фгО*-!, фо = 1,

г=0

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

V(1) = Б2; V(г) = V(г - 1) + Б2ф2_!, (13)

где Б2 представляет оценку дисперсии случайного белого шума или просто оценку его нормы, если он далеко не "белый", вычисляемую по формуле (11). Весовые коэффициенты фг определяются вектором 0. В качестве интегрального критерия можно взять некоторый линейный функционал от дисперсий прогноза ^(0) на интервале пролонгации Ь, в частности, просто среднюю величину дисперсии прогноза:

1 Ь

^(0) = ^ V(г)- (14)

г=!

Обозначим через

ft = ^, i = 1,...,m, (15)

dpi

компоненты вектора-строки градиента G, а через

h-- = д ^^ ij = 1 m (16)

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

Ge = G ■ H-1.

Если матрица H вырождена, т. е. rank H < m, операция обычного обращения заменяется операцией псевдообращения:

Ge = G ■ H+, (17)

где H + — псевдообратная матрица Мура—Пенроуза [10], удовлетворяющая условиям

HH+H = H; H+HH + = H+;

(18)

(HH+)T = HH+; (H+H )T = H+H, ^ J

обеспечивающим ее единственность. Здесь индекс "T" — символ транспонирования.

Если оценка вектора параметров модели, полученная на i-м шаге итерации, равна fli, то с помощью обычного скорейшего спуска уточненная оценка будет выглядеть следующим образом:

0i+1 = Pi — aiGi/ || Gi || •

Здесь Сг — вектор градиента, полученный на г-м шаге при в = вг; II || — его евклидова норма, а аг — некоторый коэффициент, смысл которого мы поясним несколько ниже.

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

вг+1 = вг - ТгСвг/ || Св; || •

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

Значение критерия И (в) вычисляется в равноудаленных на расстояние 8 точках вдоль градиента — на спуске и подъеме. Вначале на спуске в точках

вг - 3 • 8Свг/ || Св% ||, 3 = 1, 2,...,

вплоть до значения 3, при котором значение критерия превысит в п раз исходное значение ^(вО, где п = 1.5... 3 (и должно уменьшаться по мере приближения к глобальному минимуму). Затем такое же сканирование осуществляется на подъеме, т. е. для 3 = -1, -2,... Обычно в этом направлении расстояние до превышения ^(вО в п раз невелико.

Среди полученных точек находится точка с минимальным значением критерия. Поведение критерия в ее окрестности аппроксимируется параболой, и находится значение в, отвечающее минимуму параболы. Оно и является итоговым для данного спуска. Если минимальное значение критерия превышает исходное — аналогичное параболическое уточнение критерия ^(вО проводится для этой исходной точки, а соответствующее ему значение в принимается за вг+1. После этого осуществляется переход к скорейшему спуску по преобразованному матрицей Гессе градиенту, где операции производятся аналогично, и это завершает один цикл итерации. Итерации прекращаются, когда уменьшение интегрального критерия окажется меньше задаваемого порога е:

^(вг) - V(вг+1) <

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

8. Выбор модели и проверка ее адекватности

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

Адекватность модели оценивается по степени приближения к белому шуму последовательности аь [2, 6]. По-видимому, нельзя от реального процесса требовать идеальной "белизны", тем более при ограниченном объеме выборки.

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

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

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

9. Прогнозирование однородных нестационарных процессов АРПСС

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

Прогнозирование проводится с помощью линейной разностной модели последовательным продвижением на один шаг. На каждом шаге прогноз получается как условное математическое ожидание при известных значениях Zt и вычисленных заранее значениях последовательности а на интервале. Когда измеренные значения процесса и вычисленные значения а исчерпываются — вместо Zt используются их прогнозные значения а вместо а при Ь > п — их математические ожидания, т.е. нулевые значения. При этом в линейной разностной модели в качестве приведенных параметров АР будем использовать

к

ЦфРг (В* )▽% = £ , < = 1. (19)

г=1 г=0

То есть теперь в параметрах будут учтены и сезонные разности. Постоянная компонента сигнала вычисляется по оценке математического ожидания процесса после взятия сезонных разностей из Zt и по параметрам АР из (1):

во = С • ^, (20)

где

-. и к

^ = п-^л £ С = П(1 - £ ^).

г=ил г=1 ]=1

Прогнозирование осуществляется в соответствии с формулой

ид ир

Zt = во + £ в>4-т - £ , (21)

г=1 г=1

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

Точность прогноза Zt оценивается с использованием представления теперь уже нестационарного процесса процессом СС:

Zt = ^ Фг0,г-г, фо = 1.

г=1

Коэффициенты вычисляются по формулам

(1 для г = 0,

Е 3Ф— - в;, г > 1.

3 = 1

Дисперсия прогноза для г-го момента (г = 1,... ,Ь):

1-1

V (1) = # (22)

Квантили для доверительных интервалов будут зависеть от закона распределения белого шума о^. Если он оказывается нормальным, то 95%-й квантиль д будет 1.96, а 99%-й квантиль д — соответственно 2.58. И доверительные интервалы для прогноза в г-й точке будут

Zup(г) = Zn+1 + ^v/V(r); Zdn(l) = Zn+1 - дл/ПТ).

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

р(| ^„+11> AVV(Г)) > .

То есть для 95%-й вероятности пребывания в доверительном интервале потребуется квантиль д = 4.47, более чем в 2 раза превышающий соответствующий квантиль для нормального распределения. В последующем мы рассчитываем оценивать закон распределения ^ и по нему выбирать квантили, как и рекомендуется в [11, 12].

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

10. Иллюстративный пример

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

модели, составляет 5 недель (840 почасовых измерений). Использованная для описания процесса модель обладала следующей структурой: количество периодов сезонности k = 3; Si = {1, 24,168}, p = {2,1,1}, q = {3,1,1}, di = {0,0,1}, i = 1, 2, 3.

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

При данной структуре модели количество избыточных параметров АР почти в 3 раза превышает количество базовых (11 вместо 4), а избыточных параметров СС ровно в 3 раза больше базовых (15 вместо 5).

Для вектора параметров модели данного процесса была получена оценка:

в = (0.473, 0.099, 0.041, -0.798, -0.504, -0.077, -0.024, -0.057, -0.013)

с СКО белого шума Smi — 54.7.

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

Вычисления производились на PC AMD Sempron Processor 3000+, 1.60 ГГц, с ОЗУ 480 Мбайт и объемом дисковой памяти 74.5 Гбайта.

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

При задании на использование двух начальных приближений время счета составляло 2 с, пяти — 3 c, десяти — 6 c и для 25 начальных приближений — 13 с. То есть зависимость практически линейная. При этом величина критерия для данного ряда изменялась несущественно. Несколько изменялись оценки компонент вектора параметров модели, и это влекло изменение кумулятивной периодограммы и корреляционной функции для at. То есть варьируя в некоторых пределах количеством задаваемых начальных приближений, видимо, можно получить более адекватную прогнозируемому процессу модель. В данном примере был выбран вариант с 10 начальными приближениями. Он обеспечивал лучшее соответствие модели представляемому ряду в сравнении с вариантом с двумя начальными приближениями.

Адекватность полученной модели оценивалась по близости последовательности at, t £ (1.672) к реализации белого шума. Выборочная оценка плотности распределения вероятности at существенно отличается от нормальной. Она унимодальна, почти симметрична, в сравнении с нормальной плотностью более концентрированна в центре, но имеет затянутые хвосты. В связи с этим квантили, соответствующие 99%-му доверительному интервалу, равные -4.06 и 3.62, и соответствующие 95%-му доверительному

интервалу, равные -2.18 и 2.12, превышают соответствующие квантили для нормального распределения, а квантили для 90%-го доверительного интервала, равные -1.56 и 1.53, меньше, нежели у нормального распределения (-1.65,1.65).

Отклонение кумулятивной периодограммы а от кумулятивной периодограммы идеального белого шума не выходит за границы 75%-го доверительного интервала.

С автокорреляционной функцией а ситуация менее идеальная. На интервале в 330 точек в пяти точках она превышает уровень 3а, достигая в максимуме отклонения значения 4.3а, где а = еще в четырех точках выходит за границы 99%-го

доверительного интервала, а в остальных пребывает внутри него.

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

кВт/ч

9 янв 16 янв 23 янв 30 янв 6 фев 13 фев20 фев Часы

Рис. 1. Прогнозная база и две прогнозируемые недели: прогнозная база (жирная линия), две прогнозируемые недели (тонкая линия)

кВт/час

ПН ВТ ср чт пт об во ПН ВТ ор чт пт об во Дата- февраль 2006 з.

Рис. 2. Двухнедельный прогноз с коррекцией двух праздничных дней и доверительным интервалом на фоне фактических данных: фактические данные (средняя тонкая линия), скорректированный прогноз (жирная линия), верхняя и нижняя границы 95%-го доверительного интервала (тонкие линии)

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

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

На рис. 1 воспроизведен фрагмент ряда, по которому осуществляется идентификация модели и строится последующий прогноз. На рис. 2 приведен двухнедельный прогноз с 95%-м доверительным интервалом.

Заключение

Цель данной работы — развитие инструмента для прогнозирования процессов с сезонными колебаниями. Мультипликативную сезонную модель АРПСС, выбранную для описания этих процессов, характеризует повышенная устойчивость работы в сравнении с иными моделями с теми же периодами сезонности и тем же количеством параметров. Это обусловлено тем, что при пошаговом прогнозе используется текущее количество данных исходного процесса, равное количеству избыточных параметров модели, которое может в 2 раза или более превышать количество основных. Это преимущество модели достигается некоторым сужением класса описываемых ею процессов. Но использование избыточного количества параметров усложняет выбор начального приближения для базового, не избыточного вектора параметров. Традиционно выбор начального приближения основывается на использовании выборочной корреляционной функции стационарного процесса или При этом оценки параметров АР находятся решением линейных уравнений Юла—Уокера, а параметров СС — с применением итерационной процедуры Ньютона—Рафсона.

В случае мультипликативной модели связь с корреляционной функцией будет осуществляться именно через избыточные параметры. И чтобы перейти от них к базовым и оценить именно их, потребуется введение нелинейных уравнений связи между ними. И вместо линейной системы уравнений Юла—Уокера пришлось бы решать нелинейную систему уравнений. Аналогичная проблема, но в более сложном варианте возникла бы и с оценкой параметров СС. Кроме того, сезонность модели вынуждает работать с элементами корреляционной функции, существенно удаленными от начала координат. И это лишь усугубило бы низкую точность исходных выборочных данных.

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

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

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

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

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

В завершение два ключевых результата, которые следует выделить в работе:

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

— двухтактный алгоритм градиентного скорейшего спуска с улучшенной сходимостью.

Список литературы

[1] Андерсон X. Статистический анализ временных рядов. М.: Мир, 1976. 756 с.

[2] Бокс Дж., Дженкинс Г. Анализ временных рядов. Прогноз и управление. М.: Мир, 1974. Вып. 1. 406 с.

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

[4] Справочник по прикладной статистике / Под ред. Э. Ллойда и У. Ледермана. М.: Финансы и статистика, 1990. Т. 2. 528 с.

[5] Канторович Г.Г. Анализ временных рядов. Лекционные и методические материалы // Экономический журнал Высшей школы экономики. 2003. № 1. С. 134.

[6] Кендэл М. Временные ряды. М.: Финансы и статистика, 1981. 200 с.

[7] Бэкон Д.У., БроэкховЕн Л.Х. Предсказание временных рядов // Статистические методы для ЭВМ: Сб. науч. тр. / Под ред. К. Энслейна, Э. Релстона, Г.С. Уилфа. М.: Наука, 1986. С. 434-452.

[8] ВоЕводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984. 320 с.

[9] Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980. 280 с.

[10] Рао С.Р. Линейные статистические методы и их применения. М.: Наука, 1968. 548 с.

[11] Орлов А.И. Современная прикладная статистика // Заводская лаборатория. 1998. Т. 64. № 3. С. 52-60.

[12] Орлов А.И. Прикладная статистика XXI в. // Экономика XXI века. 2000. № 9. С. 3-27.

Поступила в редакцию 16 мая 2007 г., в переработанном виде —16 мая 2008 г.