авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 7 |

«Самарский научный центр Российской академии наук В.К. СЕМЁНЫЧЕВ, Е.В. СЕМЁНЫЧЕВ ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ РЯДОВ ДИНАМИКИ: СТРУКТУРЫ, МОДЕЛИ, ЭВОЛЮЦИЯ ...»

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

е) асимметричной экспоненциальной модели к данным добычи нефти в Албании а) б) Рис. 2.52. Применение моделей Хабберта (а) и кусочно-линейной Рис. 2.52. Применение модели Хабберта (а) и кусочно-линейной модели (б) модели (б) к данным добычи нефти в Венгрии к данным добычи нефти в Венгрии Нет оснований полагать, что для описания столь важной для эко номики мира динамики добычи нефти и прогнозирования ее «пика» ка кая-то из моделей (2.43)-(2.49) будет обладать универсальностью. Вряд ли одна модель будет применима ко всем временным рядам добычи вне зависимости от свойств конкретного месторождения, применяемых тех нологий, мировых цен на нефть, уровня агрегирования показателя и т.д.

Линейные модели (2.46) и (2.48), как и экспоненциальные модели (2.47) и (2.49) предполагают для известных методов идентификации по лучение длинных выборок после прохождения пика, позволяя обеспе чить только при этом условии высокую точность моделирования. Одна ко это же условие приводит к их малой прогностической ценности на наиболее интересном «падающем» отрезке ЖЦП. Недифференцируе мость моделей (2.46)-(2.49) затруднит нахождение пика импульсных моделей. Отсутствие колебательных компонент при моделировании ря да, изображенного на рисунке 2.52, невозможность моделирования мультилогистической динамики (например, повторных циклов в ЖЦП, на рисунке 2.51 а, г) не позволяют считать точность моделирования в 80% достаточной. О наиболее интересном для приложений параметре точности прогнозировании («привязанной» к положению точки начала прогнозирования и горизонту прогноза) сведения не обнаружены.

Сравним теперь по сложности идентификации, функциональным возможностям модель (2.45) и другие феноменологические импульсные модели для ЖЦП:

предложенную в [53] произведением полинома на экспоненци альную функцию T (t ) = T ( A0, A1,..., Am, t ) e 1t, (2.50) дробно-иррациональную функцию вида [34] A1t T (t ) = (2.51) A2 + A3t и экспоненциальную функцию с показателем в форме полинома второго порядка (по сути, другой формой записи модели Хабберта) [32, 107, 113, 116] T (t ) = e 0 +1t + 2t, (2.52) где 1, 2, A0, A1, A2, A3,..., Am – параметры модели (в общем случае нецелые числа), T ( A0, A1,..., Am, t ) = A0 + A1t +... + Amt m – алгебраический многочлен степени m.

Графики альтернативных моделей при наборе значений парамет ров представлены на рисунке 2.53. Видим, что модель (2.52) симмет рична относительно точки максимума, а модели (2.45), (2.50), (2.51) ас симетричны, причем у них рост быстрее спада.

T T 0 t (2.45) (2.50) t T T 0 t0 (2.52) (2.51) t Рис. 2.53. Общий вид известных импульсных моделей ЖЦП (2.45), (2.50)-(2.52) Все рассматриваемые импульсные модели являются существенно нелинейными по параметрам. Из рассмотренного комплекса моделей лишь (2.50) обладает возможностью идентификации с высокой точно стью на коротких выборках, но и то лишь при малых значениях степени многочлена ( m 4 ) на основе использования ARMA-моделей [53, 48].

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

Однако модель (2.51) допускает перспективное для приложений расширение в форме дробно-рациональной функции с четырьмя пара метрами ( P0, P1, Q1, Q2 ) и целыми значениями степеней аргумента:

P0 + Pt T (t ) =. (2.53) 1 + Q1t + Q2t При определенных значениях параметров модель (2.53) в первой координатной четверти принимает вид, соответствующий типичной форме импульсной модели ЖЦП. Значения параметров P0 и Q2 должны быть положительными. Кроме того, функция (2.53) должна быть нераз рывной, т.е. ее знаменатель не должен обращаться в ноль (1 + Q1k + Q2 ( k )2 0 ), что соответствует определенным ограничениям на параметры в виде неравенства Q12 4Q2.

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

Поэтому предложена альтернативная запись:

A(t B ) + C, (2.54) T (t ) = 1 + A(t B ) в которой параметры, A, B, C модели однозначно связаны с P0, P, Q1, Q2 следующими соотношениями:

C AB A 2 AB A P0 =,P=, Q1 =, Q2 =.

1 + AB 2 1 + AB 2 1 + AB 2 1 + AB Влияние параметров, A, B, C на форму кривой (2.54) иллюстри рует рисунок 2.54.

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

T T 0 C = 0 t t T T B A 0 t t Рис. 2.54. Зависимость формы тренда (2.54) от значений его параметров Тренд (2.54) имеет единственный максимум и, в случае асиммет рии, единственный минимум. Тем не менее, с его помощью можно мо делировать и мультилогистическую динамику с повторными циклами, используя взаимодействие нескольких трендов (рис. 2.55).

T T T1·(1 + T2) T1+ T2+ T T3 T T T T 0 t0 t (а) (б) Рис. 2.55. Получение повторных (мультилогистических) циклов путем взаимодействия нескольких трендов: (а) – аддитивного, (б) – прямо пропорционального мультипликативного Значения точек максимума и минимума могут быть рассчитаны по C C формуле: B ± +.

A A A Точки перегиба могут быть определены из решения уравнения третьей степени с помощью аналитических (или численных) методов:

( ) ( ) 3 A2 ( k ) + 3 AC A2 B ( k ) + 3 A2 B 2 2 ABC A k + +3 AB A2 B3 + 3 AB 2C C = 0.

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

Видим существенные преимущества по простоте использования предложенного расширения (2.54) в сравнении с другими известными моделями построения мультилогистической динамики: шестипара метрической би-логистической аддитивной моделью (2.27) Верхулста, би-логистической мультипликативной моделью (2.28) Верхулста и, тем более, с четырнадцатипараметрической би-логистической моде лью Капицы-Баренблата.

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

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

Можно считать, что в настоящее время уже возможна замена «качественного» обсуждения десятков используемых до настоящего времени графиков ЖЦП предлагаемым комплексом аналитических моделей, допускающих моделирование и краткосрочное прогнозиро вание траекторий.

Сформулируем задачи, которые для этого необходимо решить:

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

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

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

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

ГЛАВА 3. МОДЕЛИРОВАНИЕ И ПРОГНОЗИРОВАНИЕ РЯДОВ ДИНАМИКИ С ПОЛИНОМИАЛЬНЫМ ТРЕНДОМ И ЭВОЛЮЦИЕЙ ГАРМОНИК, МЕТОДИКА ИССЛЕДОВАНИЯ ТОЧНОСТИ 3.1. Модели роста в виде суммы линейного тренда и гармоник Как уже отмечалось выше, для моделей трендов в виде полиномов первого и второго порядков известны методы МНК-идентификации, ме тодика точечных и интервальных оценок точности при условии боль ших выборок, при аддитивной структуре стохастической компоненты, удовлетворяющей условиям Гаусса-Маркова [1, 2, 12, 71].

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

В случае линейного тренда рассмотрим модели:

Yk = A0 + A1k + A2 sin( k + ) + k, (3.1) Yk = A0 + A1k + Ai +1 sin(i k + i ) + k, (3.2) i = Yk = A0 + A1k + Ai +1 sin(i k + i ) + k, (3.3) i = где k – стохастическая компонента, которая здесь и далее отвечает принятым условиям Гаусса-Маркова (напомним, с упрощением в смыс ле требования лишь симметричности закона распределения).

Модель (3.1) является наиболее простой, и ее использование воз можно в относительно широком диапазоне значений объема выборки, в том числе на коротких выборках.

Линейный тренд в (3.1) может интерпретироваться и как первые два члена разложения в ряд Тейлора нелинейного тренда. Добавление гармоник в структуру колебательной компоненты модели (3.1) для мо делирования более сложных сезонных и/или циклических колебаний приводит к моделям (3.2) и (3.3).

При идентификации моделей (3.1)-(3.3) будем использовать двух этапную последовательную процедуру с применением на каждом этапе МНК или его расширения.

На первом этапе идентификации, как и для большинства рассмат риваемых в монографии моделей, будем конструировать ARMA-модели [17, 48, 50], используя Z -преобразование модели детерминированной компоненты ряда динамики.

Для конструирования ARMA-моделей можно применить анали тический подход, описанный в [48, 50, 57], но более целесообразно использовать систему компьютерной алгебры Maple, которая позво ляет проводить численные и символьные вычисления, упрощать гро моздкие математические выражения.

Автоматизированный вывод расчетных формул опирается на символьный процессор системы Maple 12, выполняющий необходи мое для конструирования ARMA-моделей Z -преобразование детер минированной компоненты модели, а также другие вычисления. Рас сматриваемые задачи Maple выполняет быстрее, чем известные про граммы Matlab и Mathematica.

С помощью Z -преобразования для (3.1) при k 4 получим сле дующую ARMA-модель четвертого порядка:

Yk = 1 ( Yk 1 2 Yk 2 + Yk 3 ) + 2 ( Yk 1 Yk 2 + Yk 3 ) Yk 4 + k, (3.4) где 1 = 2 cos, k – значения «новой» стохастической компоненты, образованной той же весовой суммой уровней стохастической компо ненты от k 1 до k 4, как и Yk в (3.4):

k = 1 ( k 1 2 k 2 + k 3 ) + 2 ( k 1 k 2 + k 3 ) k 4.

Стохастическая компонента k, так же, как и k, имеет нулевое математическое ожидание:

M [ k ] = M [ k ] (1 + 2) M [ k 1] + 2(1 + 1) M [ k 2 ] (1 + 2) M [ k 3 ] + (3.5) + M [ k 3 ] = 0 (1 + 2) 0 + 2(1 + 1) 0 (1 + 2) 0 + 0 = 0.

Ковариационная матрица для k имеет вид:

1 2 3 4 5 0... 1 2 3 4 5... 2 3 2 1 2 3 4... 3 2 1 2 3... cov( k ) = 4, (3.6) 5 4 3 2 1 2... 5 4 3 2 0............

...............

0... 0 0 0 0 1 = 14 2 + 51 2 + 12 2, 2 = 12 2 141 2 412 2, где 3 = 8 2 + 81 2 + 12 2, 4 = 4 2 21 2, 5 = 2.

При i j 4 M [ i j ] 0, т.е. имеет место автокорреляция уровней k. Принятое условие гомоскедастичности ошибки k обеспечивает и гомоскедастичность стохастической компоненты k :

M [ k ] = M [( i (1 + 2) i 1 + 2(1 + 1) i 2 (1 + 2) i3 + i 4 )2 ] = ( ) = M i2 + 4 i2 1 + 4 i2 2 + 4 i2 3 + i2 4 + 1 i2 1 + 41 i2 2 + 1 i2 3 = = 14 2 + 51 2 + 12 2.

Можно сказать, что модель (3.1) перепараметризована: вместо па раметра, нелинейно входящего в (3.1), мы получили параметр 1, ли нейно входящий в ARMA-модель (3.4), что значительно упрощает ее идентификацию и последующий расчет через принятые обозначения в (3.4).

Определим свойства МНК-оценки 1 из модели (3.4):

1N (Yk (1 + 2)Yk 1 + 2(1 + 1)Yk 1* = arg min N 4 k =5 (3.7) ) (1 + 2)Yk 3 + Yk 4 ) 2, где через N обозначен объем анализируемой выборки.

Примем, что:

a k = Yk 2(Yk 1 Yk 2 + Yk 3 ) + Yk 4, (3.8) bk = Yk 1 2Yk 2 + Yk 3, (3.9) ck = k 2( k 1 k 2 + k 3 ) + k 4, (3.10) d k = k 1 2 k 2 + k 3. (3.11) С учетом (3.8), (3.9) можно так записать выражение (3.7):

1N ( ak 1bk )2.

1 = arg min (3.12) N 4 k =5 Реализуем необходимое условие экстремума для (3.12), осущест вив его дифференцирование по параметру 1 :

2 N N ak bk 1 bk = 0, N 4 k =5 k = откуда МНК-оценка 1 будет равна:

N ak bk 1. (3.13) k = = N bk k = Для исследования смещенности оценки 1 подставим в (3.13) вы ражение ak в явном виде из исходной модели:

N ( ck 1d k + 1bk ) bk 1 = k =5. (3.14) N bk k = Математическое ожидание МНК-оценки 1 из (3.14) будет равно:

N N (ck 1d k + 1bk )bk (ck 1d k )bk k =5 = M k =5 + M [1 ] = M N N bk2 bk k =5 k = N 2 N bk (ck 1d k )bk + 1M kN5 = M k =5 N + 1.

= bk2 bk k =5 k = Видим, что оценка 1 имеет смещение, равное:

N ( ck 1d k ) bk M k =5 N. (3.15) bk k = Величина bk в (3.9) будет мала при относительно малом отличии друг от друга соседних четырех уровней ряда и, наоборот, будет велика при значительном их различии. В отношении стохастической компо ненты то же можно сказать о величинах c k в (3.10) и d k в (3.11) при со седних четырех и пяти уровнях ряда соответственно.

Смещение (3.15) можно уменьшать путем увеличения bk и уменьшения c k и d k, то есть следует достичь значительного различия уровней ряда и малого отличия стохастической компоненты.

Для этого можно применить прием прореживания выборки: уда ления из рассмотрения каждого i -того наблюдения, в результате чего получатся i новых прореженных выборок. Одновременно необходимо накладывать на полученную таким образом выборку условие равноуда ленности соседних наблюдений. На каждой из прореженных таким об разом выборок можно произвести идентификацию 1,i по формуле (3.13). Из идентифицированных на каждой из прореженных выборок оценок 1,i следует выбирать ту, которая обеспечивает минимум дис персии и рассчитать по ней МНК-оценку частоты:

= arccos 1,i. (3.16) При разных шагах прореживания (использования уровней ряда че рез два, через три и т.д. наблюдений) надо учесть это обстоятельство, рассчитывая коэффициент ARMA-модели по формуле:

1,i = 2 cos i.

Нетрудно показать, что прорежение выборки уменьшает и диспер сию оценки 1 (здесь через D обозначен оператор нахождения диспер сии):

N (ck 1d k + 1bk )bk D[1 ] = D k =5 = N bk k = N N (ck 1d k + 1bk )bk (ck 1d k + 1bk )bk k =5 = M k = =M N N 2 bk bk k =5 k = N N (ck 1d k )bk (ck 1d k )bk = M k =5 N 1 = + 1 M k =5 N 2 bk bk k =5 k = N (ck 1d k )bk = D k =5 N.

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

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

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

Идентификация любой ARMA-модели имеет ограничение по мультиколлинеарности – высокой взаимной коррелированности (линей ной связи) переменных в модели авторегрессии [12].

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

Строгая (perfect) мультиколлинеарность – наличие линейной функциональной связи между независимыми факторными переменными (иногда также между зависимой и определяемой переменной). Мульти коллинеарность (являясь, по сути, коррелированностью) нарушает одно из принятых условий Гаусса-Маркова и делает построение регрессии невозможным.

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

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

При этом однозначных критериев «силы» мультиколлинеарности не су ществует.

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

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

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

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

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

Иногда коэффициенты авторегрессии могут принять принципи ально неверные значения с точки зрения теории. Например, может ока заться, что определяемые через коэффициенты регрессии параметры должны рассчитываться при условиях cos 1 или exp() 0 и т.д.

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

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

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

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

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

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

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

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

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

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

После подстановки найденной на первом этапе идентификации оценки в модель (3.1) будем на втором этапе идентификации оцени вать оставшиеся параметры A0, A1, A2, модели (3.1) при помощи МНК, который формирует нормальную СЛАУ (с линейно входящими в нее параметрами A0, A1, A3, A4 ) из условия:

A0, A1, A3, A4 = N 2 (3.17) ( ) = arg min Yk A0 A1k A3 sin k A4 cos k, A0, A1, A3, A4 k =1 где A3 = A2 cos, A4 = A2 sin.

После получаемых из (3.17) МНК-оценок A0, A1, A3, A4 очевиден расчет:

A ( ) +( ) 2 A4, = = arctg.

A2 A A Для модели (3.2) в виде суммы линейного тренда и двух гармоник будет при k 6 справедлива ARMA-модель шестого порядка:

Yk (2 + 1 + 2 )Yk 1 + (31 + 2(1 + 2 ) + 12 )Yk (4 + 2(1 + 2 ) + 12 )Yk 3 + (3 + 2(1 + 2 ) + 12 )Yk 4 (3.18) (2 + 1 + 2 )Yk 5 + Yk 6 = k, где 1 = 2 cos 1, 2 = 2 cos 2, а «новая» стохастическая компонента определена соотношением того же вида, что и (3.18), но в уровнях сто хастической компоненты k :

k = k (2 + 1 + 2 ) k 1 + (31 + 2(1 + 2 ) + 12 ) k (4 + 2(1 + 2 ) + 12 ) k 3 + (3 + 2(1 + 2 ) + 12 ) k 4 (3.19) (2 + 1 + 2 ) k 5 + k 6.

Относительно k в (3.19) можно сделать выводы о ее центриро ванности, существовании автоковариации k для соседних шести на блюдений и гомоскедастичности.

Введем следующие обозначения:

(3.20) m = 1 + 2, (3.21) g = 12, ak = Yk 2(Yk 1 + Yk 5 ) + 3(Yk 2 + Yk 4 ) 4Yk 3 + Yk 6, bk = Yk 1 + 2(Yk 2 Yk 3 + Yk 4 ) Yk 5, ck = Yk 2 2Yk 3 + Yk 4, d k = k 2( k 1 + k 5 ) + 3( k 2 + k 4 ) 4 k 3 + k 6, ek = k 1 + 2( k 2 k 3 + k 4 ) k 5, f k = k 2 k 3 + k 4.

Применим МНК для ARMA-модели (3.18), реализуя условие:

N m, g = arg min ( ak + mbk + gck ).

(3.22) m, g k =7 Будем иметь в этом случае систему нелинейных уравнений отно сительно искомых параметров модели (в силу определения m в (3.20) и g в (3.21)), но, тем не менее, она может получить точное аналитическое решение.

Для этого полученные МНК-оценки m и g подставим (с учетом обозначений в (3.20) и (3.21)) в уравнение x 2 mx g = 0. Eсли корни вещественные и принадлежат интервалам значений 2 m 2, 2 g 2, то решение будет найдено из системы уравнений:

m = 1 + g = 12.

При этом необходимо рассмотреть два случая, так как 1 и 2 вхо дят в выражение симметрично.

N = arg min ( ak + (2 + 2 )bk + 22ck ).

Если 1 = 2, то 2 i = N = arg min ( ak + (2 + 2 )bk 22ck ).

Если 1 = 2, то 1 i = Дифференцируя (3.22) по переменным m и g, получим:

dN N N N (ak + mbk + gck ) 2 = 2 ak bk + m bk + g bk ck, dm k =7 k =7 k =7 k = dN N N N (ak + mbk + gck ) 2 = 2 ak ck + m bk ck + g ck.

dg k =7 k =7 k =7 k = Приравнивая выражения для частных производных к нулю и ре шая получаемую при этом нормальную СЛАУ, будем иметь:

akbk bk ck ak ck bk2, g= (3.23) ck2 bk2 ( bk ck ) akbk ck2 ak ck bk ck.

m= (3.24) ( bk ck ) bk2 ck Аналогично случаю с рассмотренной выше моделью (3.1), подста вим в (3.23) и в (3.24) a k в явном виде из исходной ARMA-модели и, находя математическое ожидание, выразим смещения оценок.

Смещение оценки g = 12 будет равно:

d + mek + gf k ) bk bk ck ( d k + mek + gf k ) ck bk ( k. (3.25) M ck bk ( bk ck ) 2 Смещение оценки m = 1 + 2 определит формула:

d + mek + gf k ) bk ck ( d k + mek + gf k ) ck bk ck ( k. (3.26) M ( bk ck ) ck bk 2 Как и в случае с моделью (3.1), величины bk и c k в (3.25) и (3.26) будут малы при относительно малом отличии друг от друга соседних пяти уровней ряда и, соответственно, велики при значительном их раз личии.

В отношении стохастической компоненты то же самое можно ска зать о величинах d k, ek, f k при соседних семи, пяти и трех наблюдениях соответственно.

Смещения (3.25) и (3.26) можно также уменьшить путем примене ния приема прореживания выборки, которое одновременно со смещени ем уменьшит и дисперсию оценок m и g :

( d k + mek + gf k mbk gck ) bk bk ck D[ g ] = D ck2 bk ( d k + mek + gf k mbk gck ) ck bk = (3.27) ( bk ck ) ( d + me + gf ) b k k k k bk ck ( d k + mek + gfk ) ck bk2.

= D ck bk ( bk ck ) 2 Аналогично выражается дисперсия и для m. Оценка оставшихся параметров модели (3.2) предполагает подстановку в нее обеих полу ченных оценок частот, составление и решение нормальных СЛАУ шес того порядка для оценки параметров линейного тренда, амплитуд и фаз гармоник, как это было сделано для (3.1).

ARMA-модель для наиболее сложной модели (3.3) будет иметь вид:

Yk = 2Yk 1 4Yk 2 + 6Yk 3 6Yk 4 + 6Yk 5 4Yk 6 + 2Yk Yk 8 + (1 + 2 + 3 )(Yk 1 2Yk 2 + 3Yk 3 4Yk 4 + 3Yk (3.28) 2Yk 6 + Yk 7 ) (12 + 13 + 23 )(Yk 2 2Yk 3 + 2Yk 2Yk 5 + Yk 6 ) + 123 (Yk 3 2Yk 4 + Yk 5 ) + k, где 1 = 2 cos 1, 2 = 2 cos 2, 3 = 2 cos 3, k – «новая» стохасти ческая компонента, по своей структуре аналогичная выражению (3.28), но в уровнях k.

Предыдущий пример показал трудность решения получающихся при применении МНК к (3.28) систем с зависимыми коэффициентами ( 1 2 + 13 + 2 3 ) при (Yk 2 2Yk 3 + 2Yk 4 2Yk 5 + Yk 6 ) и 1 2 3 при (Yk 3 2Yk 4 + Yk 5 ) в нормальных системах алгебраических уравнений.

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

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

Во-вторых, если параметрическая итерационная декомпозиция не обеспечивает решения, то можно использовать другие аналитические методы для решения получающихся в этих случаях уравнений: базисы Гребнера [13, 24, 74] или метода В.Н. Кублановской [27].

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

Понятие базиса Гребнера было введено австрийским математиком Бруно Бухбергером, который назвал конструкцию в честь своего науч ного руководителя – В. Гребнера, построившего первые примеры стан дартных базисов для свободного коммутативного случая в конце 40-х годов. Бухбергер предложил первую конструктивную процедуру вы числения базиса Гребнера, известную ныне как алгоритм Бухбергера.

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

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

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

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

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

Затем, в соответствии с классификацией, применяется алгоритм поис ка базиса Гребнера алгоритмом Бухбергера и его модификациями, появившимися позднее.

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

В системе Maple вычисление базиса Гребнера и сопутствующие таким вычислениям преобразования объединены в пакете «Groebner». Со временный пакет Maple позволяет работать с 500 000-значными чис лами (при соответствующей конфигурации машины).

Диапазон применимости базиса Гребнера для решения рассмат риваемых задач определился численным экспериментом: путем ис следования практики построения базиса Гребнера для систем N урав нений с N неизвестными полиномиальных многочленов степени M в системе Maple (таблица 3.1).

Таблица 3. Границы применения базиса Гребнера для решения систем полиномиальных уравнений с целыми коэффициентами в системе Maple Величина 1.. 10 50.. 100 1 000.. 10 000 10 000.. 100 коэффиц.

Степ\К-во* 2 3 4 5 5 7 2 3 4 5 6 2 3 4 5 6 2 3 4 2 1 1 1 1 2 3 1 1 1 1 3 1 1 1 1 3 1 1 1 3 1 1 1 3 1 1 3 1 1 3 1 2 4 1 1 3 1 2 3 1 3 1 5 1 2 3 1 3 1 3 2 6 1 3 1 3 1 3 2 7 1 3 1 3 1 3 3 8 2 3 2 3 3 3 3 9 2 3 2 3 3 3 3 10 2 3 3 3 3 3 3 * - по столбцам - количество переменных, по строкам - степень многочленов Обозначение результата вычислений:

1 Применение базиса Гребнера успешно Базис Гребнера вычисляется, но корни системы находятся 2 неверно из-за больших вчислительных погрешностей 3 Базис Гребнера не вычисляется Основная вычислительная проблема алгоритмов поиска базисов Гребнера – необходимость работы с большими целочисленными ко эффициентами, образующимися путем вычисления наименьшего об щего кратного коэффициентов исходной системы. Поэтому исследо вание границ практического применения проводится для систем урав нений, обладающих коэффициентами с различным количеством зна чащих цифр.

Результаты практического применения базиса Гребнера для ре шения систем полиномиальных уравнений в системе Maple показали, что базис Гребнера может применяться для решения систем не более, чем из 5-ти уравнений с 5-тью неизвестными (в случае многочленов 2-ой степени), и систем полиномов не более, чем 7-ой степени (в слу чае двух уравнений с двумя неизвестными).

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

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

Указанные методы использовались для нахождения параметров модели, состоящей из суммы линейного тренда и трех гармоник при помощи ARMA-моделей, использования формулы (3.28) для нахожде ния оценок трех частот 1, 2, 3, постановки полученных оценок час тот в (3.3) и последующего применения МНК для расчета оставшихся параметров модели.

Покажем пример количественного составления и решения поли номиального уравнения с помощью базиса Гребнера на примере моде лирования ВВП РФ за 2003-2006 гг. (рис. 3.1) моделью (3.3) с оценкой полученного прогноза на 2007-2008 гг. путем применения МНК и бази са Гребнера для модели (3.28).

2003 2004 2005 2006 2007 Рис. 3.1. Индекс реального ВВП РФ GDPEA_Q_DIRI, 2003.01= (данные Росстата URL: http: // www.gks.ru) Система нормальных уравнений при применении МНК будет иметь вид:

21247 + 105471 + 213592 + 213593 + 2 +53502 + 53503 + 1070112 + 209042 3 + 1070113 + 2 2 2 +10056123 + 274812 + 50282 3 + 502823 + 274813 + 22 2 2 +12232 3 + 244512 3 + 2445123 + 912212 3 = 0, 21247 + 213591 + 105472 + 213593 + +535012 + 53503 + 1070112 + 1070123 + 2090413 + +10056123 + 274812 2 + 27482 3 + 502812 3 + 2 +1223123 + 2445122 3 + 2445123 + 9122122 3 = 0, 2 2 21247 + 213591 + 213592 + 105473 + +535012 + 53502 + 2090412 + 1070123 + 1070113 + +1005612 3 + 502812 2 + 502812 + 27482 3 + 2748123 + 2 +1223122 + 24451223 + 244512 3 + 912212 2 3 = 0, 2 2 содержащий переменные 1, 2, 3 в степени не выше второй, что удовлетворяет требованиям таблицы 3.1 и позволяет рассчитывать на ее точное аналитическое решение.

Для вычисления базисов Гребнера использовали функция Groebner[Basis] системы Maple.

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

В рассматриваемом случае вычисления индекса ВВП России уравнение базиса Гребнера для идентификации параметров ARMA модели (3.28) примет вид:

+9578340658884866708628591791504712942773644600838160587689389422415452804786351251299231119+ +601366362885947812797806388505288862959176159336991738349738701721935736927869957767726901l8+ +226787925806526820724324513838610416662573989511170574297902324201130463923306725777818960117+ +603376634771468776851875018435706094686520423608544496505701596541055950186540487294839568116+ +1291513725644620084617569358839467554060254351176426842192292044842257849820814028662827760115+ +2374204470436624096492387908450287975214740561100063707059207576842588698010319670891359008114+ +3569958053597260799327741262330720166808437965965302946119304896036662501646697117768082304113+ +3921982688536177685312010965466206793664444799161315166529486782930776113982185162237415296112+ +2725966175638407370121903031026992120264494940872013324070279798404775234722338568033501952111+ +8120668252732964167156171125814452214241155066752412337631998632596763860489499648453519361l0+ -31990741587691847552702700937130186427105860565927826712989927403296929317966463578410700819 -36840988860839723855305646000586586923283812921764596090329582642208941627900516430157414418 -9926813124873649913013870101909013830741436219457892826503832884621149619875577836061900817 +90087016533516216421025565311447480924538949621673673786018138077759994102907845272371216+ 130991408290966855184594853744443847030781916840983573332140140199728429105894813199564815+ +18584375206546991692930439167126977903572474909453795011707193840716876474854323945472014+ +39708862311793447490793614751312126710423429079285085677452276753304768823983839641613+ -8543559285037587938145896375581760291823009617164461485189450987731127217537299251212 -12141776510175756264403948966041703917587277364414350191973696877265497044678082561 -5015986840051061219856165558985759717575367263985931934562811905593675963432960=0.

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

Корни уравнения относительно 1 будут равны:

2,01055;

2,00662;

1,62293;

0,75380;

0,01114;

0,02474;

0,13546;

0,45877.

1 = arccos, исключая не имею Проводя обратную замену щие экономического смысла комплексные значения, получим следую щие оценки частоты 1 :

1 = {2,517;

1,957;

1,576;

1,558;

1,503;

1,339}.

Вследствие симметрии системы такие же уравнения базиса Греб нера и такой же набор решений будет существовать и для 2, 3 ( 2, 3 ).

Истинные частоты системы определяются путем перебора из критерия максимального значения получаемого коэффициента детерминации R 2.

В рассматриваемом случае частоты будут равны:

1 = 1,503, 2 = 1,558, 3 = 2,517.

Коэффициенты C1, C 2, а также амплитуды колебаний A1, A2, A входят в модель динамики (3.3) линейно и находятся с помощью МНК из нормальной СЛАУ. В итоге, коэффициенты модели индекса реально го ВВП будут равны:

C1 = 108;

C2 = 2,132;

1 = 1, 4361;

A1 = 6, 406;

2 = 2, 297;

A2 = 13,87;

3 = 1,334.

A3 = 0,6295;

Полученный коэффициент детерминации модели равен R 2 = 0,99, а ошибки прогноза на 2007 и 2008 гг. оказались следующими: MAPE оценка 2,74%, а второй коэффициент Тейла T2 = 2, 41% (рис. 3.2).

2003 2004 2005 2006 2007 Рис. 3.2. Реальные данные индекса реального ВВП РФ и моделирование выражением (3.3) Из представленных результатов можно сделать вывод о высокой точности проведенного моделирования и прогнозирования довольно простой моделью, что позволяет рассчитывать на ее широкое исполь зование и в других приложениях.

3.2. Методика оценки точности и области применения методов идентификации Перейдем теперь к характеристике достигаемой предложенным инструментарием точности моделирования и прогнозирования моделя ми (3.1)-(3.3) на тестовых и реальных выборках.

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

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

Более общим подходом является методика Е.М. Четыркина [69] оценки точности метода идентификации логистической функции Вер хулста, использующая программную генерацию суммы анализируемо го (полезного) тренда с известными (задаваемыми) параметрами и ад дитивной помехи. Полезный сигнал имел при этом только один набор параметров, а соотношение мощностей генерируемой помехи и полез ного сигнала задавалось равным 5%, объем выборки был большим.

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

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

Объем детерминированной компоненты не превышает реальные выборки эволюционирующей динамики: принят менее 48 наблюдений.

На каждое сочетание параметров детерминированной компонен ты накладывается от 20 до 1000 соответсвующих выборок помехи. Под наложением понимается суммирование для ряда с аддитивной стохас тической компонентой, или умножение для структуры с мультиплика тивной стохастической компонентой. Результаты наложения образуют тестовые выборки. Дисперсия генерируемой детерминированной вы борки D [ Dk ] определяет мощность полезного сигнала.

Стохастическая помеха генерируется того же объема, что и детер минированная, причем с произвольными величинами математического ожидания и дисперсии. Дело в том, что обеспечить с высокой точно стью заданные числовые характеристики помехи с нормальным законом распределения на малых объемах выборки современными программны ми средствами практически невозможно. После генерации рассчитыва ется среднеарифметическое значение тестовой выборки для центриро вания помехи, рассчитывается стандартное отклонение выборки для нормирования. Умножая полученную нормированную и центрирован ную стохастическую тестовую выборку (её дисперсия равна единице, а математическое ожидание равно нулю) на коэффициент K n / s, будем иметь дисперсию стохастической компоненты D[ ], равную квадрату этого коэффициента.

При различных сочетаниях параметров генерируемой детермини рованной выборки и различных объемах выборки исследуется зависи мость точности моделирования (значений R 2 ) и точности прогнозиро вания (значений MAPE-оценки или K T 2 ) в функции от коэффициента «шум/сигнал» K n / s :

D[ ] K n/ s =, (3.29) D[ Dk ] где D[ ] – дисперсия компоненты k.

Из (3.29) очевидно, что коэффициент, на который следует умно жать нормированную и центрированную стохастическую компоненту, будет равен квадратному корню из произведения K n / s и D[ Dk ].

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

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

При проведении тестирования будем использовать часто встре чающиеся на практике объемы выборок: 24, 36 и 48 наблюдений.

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

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

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

Показ результатов исследований точности начнем с предложенно го метода параметризации модели (3.1).

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

Наложим ограничение на угловую частоту колебательной ком поненты. Для 24 наблюдений (2 года при помесячных измерениях) может принимать значения от 0,26;

для 36 – от 0,18;

для 48 – от 0,131.

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

В связи с этим может принимать значения до 0,9. Таким обра зом, определяются границы для параметра (таблица 3.2).

Таблица 3. Возможные значения параметра частоты при моделировании N (объем выборки) Границы изменения 24 0,26 0, 36 0,18 0, 48 0,131 0, Для определенности амплитуда колеблемости относительно трен да принята достаточно большой: в 30%.

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

Таблица 3. Параметры генерации выборки № A0 A1 A 1 0 0,6 5 0,9 2 1 0,8 7 0,8 0, 3 1 1,4 9 0,5 4 1 0,9 10 0,4 5 1 0,1 1 0,3 2, R 0, 0, 0, 0, 0,78 Kn/s 0 0,05 0,1 0,15 0,2 0,25 0, N=24 N=36 N= Рис. 3.3. Зависимость R 2 от коэффициента «шум/сигнал»

MAPE, % Kn/s 0 0,05 0,1 0,15 0,2 0,25 0, N=24 N=36 N= Рис. 3.4. Зависимость ошибки прогноза от коэффициента «шум/сигнал»

В диапазоне отношения мощностей «шум/сигнал» до 30% пред ложенный метод вызывает уменьшение значений R 2 не более, чем на 15%, причём устойчивость параметризации возрастает с увеличением числа наблюдений. Ошибка прогноза, представленная на рисунке 3.4, оказывается более чувствительной к действию помехи, но также укла дывается в категорию высокой точности.

Принятое предельным соотношение K n / s до 30% является во мно гих областях исследований (в измерительной технике, в радиотехнике) общепринятым при исследовании точности.

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

Числовые характеристики точности оценки параметров моделей и доверительных интервалов модели можно получить и в области реаль ных параметров модели. Для этого следует найти МНК-оценки пара метров модели и использовать их при генерации детерминированной модели (3.1).

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

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

Для осуществления моделирования с удовлетворительной точно стью временной ряд должен включать в себя хотя бы один период вы сокочастотной гармоники, т.е. той из двух гармоник, которая входит в модель с более высокой частотой. Для определенности колеблемость относительно тренда примем 30% для той из гармоник, которая модели рует сезонность. На амплитуду второй гармоники ограничения не на кладывали. Представляет интерес исследование влияния соотношения A и частот 1. Отношения амплитуд гармоник назначались амплитуд A A = [1,5;

2;

3;

5;

10], а соотношения значений частот – в диа в диапазоне A = [0,7;

0,8;

1,1;

1,5;

2]. Результаты проведенного моделирования пазоне и прогнозирования показали высокую точность: коэффициент детерми нации более 0,8;

а MAPE-оценка прогноза менее 15% при объеме вы борки в 36 и более наблюдений. Особо выделим случай, когда угловые частоты в модели (3.2) оказываются равными. Рисунок 3.5 демонстри рует график изменения коэффициента детерминации в функции коэф фициента «шум/сигнал» при принятых условиях. Видим, что оценки в проведенных экспериментах отражают высокую помехозащищенность предложенного метода идентификации.

R2 0, 0, 0, 0, 0,78 Kn/s 0 0,05 0,1 0,15 0,2 0,25 0, N=24 N=36 N= Рис. 3.5. Зависимость R 2 от коэффициента «шум/полезный сигнал»

при равенстве частот гармоник При коэффициенте «шум/сигнал» в 20% коэффициент детермина ции в среднем оказался равен 0,9, а ошибки прогноза менее 10%. Про ведены и исследования точности моделирования при различном соот ношении фаз гармоник: 1 = [0,5;

1;

2;

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

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


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

Исследование производились со следующими соотношениями па раметров модели: A2 = [1,5;

2;

3;

5;

10] ;

A2 = [9;

7;

4;

2,5;

0,5] ;

1 = [0,7;

0,8;

1,1;

1,5;

2] ;

A3 A = [2;

1,6;

1;

0,9;

0,85].

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

R2 0, 0, 0, 0, 0,78 Kn/s 0 0,05 0,1 0,15 0,2 0,25 0, N=24 N=36 N= Рис. 3.6. Зависимость R 2 от коэффициента «шум/сигнал» при соотношении амплитуд и частот гармоник 3;

4 и 1,1;

1 соответственно MAPE, % Kn/s 0 0,05 0,1 0,15 0,2 0,25 0, N=24 N=36 N= Рис. 3.7. Зависимость ошибки прогноза от коэффициента «шум/сигнал»

при соотношении амплитуд и частот гармоник 3;

4 и 1,1;

1 соответственно На рисунках 3.8, 3.9 и 3.10 показаны примеры использования мо делей (3.1), (3.2) и (3.3) предложенных методов их идентификации на основе ARMA-моделей и на реальных выборках на данных потреби тельского бюджета в среднем на душу населения, расходов на покупку товаров, расходов на покупку товаров и оплату услуг на душу населе ния Самарской области.

Рис. 3.8. Моделирование и прогнозирование потребительского бюджета в среднем на душу населения Самарской области суммой линейного тренда, одной гармоники и ошибки ( R 2 = 0,99 ) Рис. 3.9. Моделирование и прогнозирование расходов на покупку товаров населения Самарской области суммой линейного тренда, двух гармоник и ошибки ( R 2 = 0,95 ) Рис. 3.10. Моделирование и прогнозирование расходов на покупку товаров и оплату услуг населения Самарской области суммой линейного тренда, трех гармоник и ошибки ( R 2 = 0,9 ) На рисунках 3.11 и 3.12 показано сравнение точности моделирова ния и прогнозирования временных рядов общего объема денежных до ходов и расходов населения на основе предложенных ARMA-моделей и известной модели Хольта-Уинтерса. Прогноз производился на год впе ред: на 2009 г. Оказалось, что прогнозы, полученные на основе предло женных ARMA-моделей, требуют меньших объемов выборки (от 24 на блюдений или 2 лет) относительно известных моделей ARIMA, методов экспоненциального сглаживания и классической сезонной декомпози ции, предполагающих использование выборки объемом более 5-ти лет, или 60-ти ежемесячных наблюдений;

- применение ARMA-моделей дало увеличение точности модели рования на 10-20% коэффициента детерминации и на 1-3% - MAPE оценки прогноза.

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

В абсолютных показателях 1% ошибки этого показателя соответ ствует примерно 658 млн руб.

Таким образом, рынок потребления может быть переоценен или недооценен при допущении на этапе прогнозирования ошибки в один процент или примерно 0,66 млрд руб.

Рис. 3.11. Моделирование и прогнозирование общего объема денежных доходов населения Самарской области, млн руб. (на основе предложенной ARMA-модели и модели Хольта-Уинтерса) Рис. 3.12. Моделирование и прогнозирование общего объема денежных расходов населения Самарской области, млн руб. (на основе предложенной ARMA-модели и модели Хольта-Уинтерса) В таблице 3.4 приведено сравнение оценок точности прогнозиро вания (MAPE-оценок) для рядов общего объема денежных доходов в Самарской области при различном периоде наблюдений (год и месяц).

Таблица 3. Сравнение MAPE-оценок точности прогнозирования рядов общего объема денежных доходов населения с различным периодом наблюдений MAPE-оценка при прогнозировании общего объема денежных доходов населения:

а) с периодом б) с периодом наблюдений – месяц наблюдений – год Обобщенные параметрические 4,0% 0,64% ARMA ARIMA 4,8% 6,1% Модель экспоненциального сглаживания (в случае с еже 5,2% 2,6% месячными наблюдениями модель Хольта-Уинтерса) Декомпозиция Census II с экспоненциальным 7,9% сглаживанием Прогнозирование на основе ARMA-моделей, исходя из проведен ных исследований, позволило получить удовлетворительные по точно сти прогнозы (MAPE-оценка менее 5%), превосходящие, в том числе и на «длинных» выборках (объемом более 5 периодов), известные методы прогнозирования на основе моделей ARIMA, Хольта-Уинтерса, а также методов декомпозиции Census I и Census II.

Высокой точностью обладают и результаты моделирования под писного тиража газеты, условно названной «ИС», в г. Самара, прове денные помесячно с января 2002 г. по декабрь 2005 г. (рис. 3.13).

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

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

01. 06. 04. 09. 02. 07. 12. 05. 10. 11. Рис. 3.13. Подписной тираж газеты «ИС» (тыс. шт.) Затем наступает резкий сезонный спад (потери до 4 тысяч под писчиков) при летней «полугодовой» подписной кампании с 1 июля, последующий рост аналогичный росту 1-го полугодия и менее значи тельный (до 1000 подписчиков) спад при зимней подписной кампа нии.

1 2 3 4 5 6 7 8 9 10 11 - - - Рис. 3.14. Стационарная сезонная компонента тиража газеты «ИС»

Покажем, что можно предположить простую модель (3.1) дина мики подписного тиража и параметризировать ее на реальных данных с помощью ARMA-модели (таблица 3.5).

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

Таблица 3. Коэффициент детерминации, ошибка прогноза и аналитические выра жения моделей в функции объёма выборки в методе ARMA-моделей Ошибка Кол-во на- R прогноза, Аналитическое выражение модели блюдений % Т(t)=39765,07 12 0,9138 64, 2997,38t+7448,17cos(0,52t+3,68) 15 0,8661 17,65 T(t)=30907,15+84,62t+3522,06cos(0,52t-1,08) 18 0,8701 11,40 T(t)=32106,08-193,26t+2627,09cos(0,52t-1,24) 21 0,8758 5,69 T(t)=33473,48-466,08t+1440,58cos(0,52t-1,34) 24 0,8945 6,77 T(t)=33244,68-423,08t+1538,66cos(0,52t-1,5) 27 0,9011 3,78 T(t)=33560,55-488,85t+1725,11cos(0,52t-1,28) 30 0,9235 4,89 T(t)=33472,01-474,85t+1684,14cos(0,52t-1,23) 33 0,9362 4,60 T(t)=33747,24-505,91t+1474,71cos(0,52t-1,3) 36 0,9421 7,65 T(t)=33606,79-491,2t+1528,77cos(0,52t-1,39) 39 0,9531 6,56 T(t)=33595,62-489,81t+1465,35cos(0,52t-1,37) 42 0,9552 3,55 T(t)=33380,33-468,79t+1335,36cos(0,52t-1,24) 45 0,9624 7,27 T(t)=33399,96-470,36t+1323,24cos(0,52t-1,24) 48 0,9646 12,92 T(t)=33256,23-459,88t+1393,07cos(0,52t-1,33) Графическая иллюстрация сравнения точности моделирования с использованием предложенного ARMA-моделирования и метода клас сической сезонной декомпозиции по точности моделирования и прогно за приведена на рисунках 3.15 и 3.16.

Рис. 3.15. Сравнение значений R2 от числа наблюдений для методов классической сезонной декомпозиции и ARMA-моделирования MAPE, % 12 метод параметрической авторегрессии метод сезонной 4 декомпозиции 24 27 30 33 36 39 42 45 число наблюдений Рис. 3.16. Сравнение MAPE-оценок прогноза в функции количества наблюдений для методов классической сезонной декомпозиции и ARMA-моделирования Из рисунка 3.15 видно, что метод авторегрессии на короткой вы борке в 24 наблюдения уже обеспечивает высокую точность моделиро вания, причём нет необходимости в увеличении этого объёма. При оди наковой точности выигрыш по объёму выборки предложенного метода в сравнении с известным – более чем в два раза. Из рисунка 3.16 видим, что выигрыш в точности прогнозирования – полтора-два раза.

Еще одним примером может быть модель (3.2), которая визуально довольно близка к динамике годового сбора зерновых в Самарской об ласти в тыс. тонн (рис. 3.17) [60].

Рис. 3.17. Валовой сбор зерна по Самарской области в 1953-2005 гг.

(в хозяйствах всех категорий;

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

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

При моделировании использовались выборки годовых валовых сборов зерна разного объема: от 8 до 20 годовых наблюдений.

Напомним, что выборки, используемые в известных методах мно гомасштабной параметризации, требуют анализа от 4 до 12 периодов самой низкочастотной из двух гармоник, что с учётом теоремы Котель никова предполагает выборки от 24 до 100 наблюдений, т.е. от 24 до 100 лет.

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

В данном случае использовался, для сравнения результатов не ба зисы Гребнера, а другой метод идентификации (известный как метод вариации параметра Койка [71]): перебор значений одной из идентифи цируемых частот и выбор того из них, которое обеспечивает макси мальное значение коэффициента детерминации R 2.

Шаг перебора частоты обеспечивал 100 оценок параметров моде ли, из которых выбирались лучшие по принятому критерию. Набор па раметров модели, при котором было достигнуто максимальное значение величины R 2, использовался для прогноза на один шаг (на один год вперед), и при этом рассчитывалась MAPE-оценка прогноза.


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

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

Таблица 3. Результаты моделирования и прогнозирования валового сбора зерна в Самарской области на 2006 г.

В качестве критерия В качестве критерия Объем принимался принимался Max R выборки Max MF = Max( R 2 MAPE ) N MAPE, % MAPE, % R2 R 8 0,91 5,7 0,90 1, 10 0,88 11,5 0,83 8, 12 0,90 16, 14 0,68 27, 16 0,67 41, 18 0,75 42, 20 0,73 47, На рисунке 3.18 приведен результат моделирования сбора зерна при выборке в 12 наблюдений.

Реальная выборка Модель Прогноз 0 2 4 6 8 10 12 Рис. 3.18. Реальная выборка сбора зерна объёмом N =12, её моделирование выражением (3.2) и ошибка прогноза на один шаг При этом точность моделирования и прогнозирования оказалась выше на малых выборках (8-10 наблюдений), что обусловлено, видимо, высокой эволюцией ряда динамики.

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

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

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

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

Например, модели Yk = A0 + A1k + A2k 2 + A3 sin( k + ) + k будет соответствовать при k 5 ARMA-модель пятого порядка:

Yk 3Yk 1 + 4Yk 2 4Yk 3 + 3Yk 4 Yk 5 1(Yk 1 3Yk 2 + 3Yk 3 Yk 4 ) = = k 3 k 1 + 4 k 2 4 k 3 + 3 k 4 k 5 1( k 1 3 k 2 + 3 k 3 k 4 ), где 1 = 2 cos «перепараметризированный» в виде линейного коэф фициента ARMA-модели нелинейный параметр – частота колебаний гармонической компоненты.

Выше уже показан метод нахождения МНК-оценок 1, из нор мального уравнения, а МНК-оценок коэффициентов A0, A1, A2, A3 - из нормальной СЛАУ пятого порядка.

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

Yk = A0 + A1k + A2k + Ai+ 2 sin(i k + i ) + k i = при k 7 будем иметь ARMA-модель седьмого порядка:

Yk 3Yk 1 + 5Yk 2 7Yk 3 + 7Yk 4 5Yk 5 + 3Yk 6 Yk m(Yk 1 3Yk 2 + 4Yk 3 4Yk 4 + 3Yk 5 Yk 6 ) + + g (Yk 2 3Yk 3 + 3Yk 4 Yk 5 ) = = k 3 k 1 + 5 k 2 7 k 3 + 7 k 4 5 k 5 + 3 k 6 k m( k 1 3 k 2 + 4 k 3 4 k 4 + 3 k 5 k 6 ) + + g ( k 2 3 k 3 + 3 k 4 k 5 ), где m = 1 + 2, g = 12, 1 = 2 cos 1, 2 = 2 cos 2.

Для получения МНК-оценок 1, 2 (и МНК-оценок 1, 2 ) можно использовать базисы Гребнера, а для оценок A0, A1, A2, A3, A4, 1, 2 применить метод параметрической итерационной декомпозиции.

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

Программная реализации всех методов идентификации данного раздела приведена в [41].

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

( ) Yk = Tk 1 + S k + S kA + k = Tk + Tk S k + S kA + k, M M (3.30) M где Tk – полиномиальный тренд, Tk Sk – мультипликативная колеба тельная компонента, прямо пропорциональная тренду, S kM 1, SkA – не зависимая по отношению к тренду аддитивная колебательная компо нента, k – стохастическая компонента, которую будем считать отве чающей принятым условиям Гаусса-Маркова, k = 1,2,..., N, N – объем выборки [44, 54].

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

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

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

Например, они могут назначаться, исходя из экономической сущ ности моделируемого ряда динамики, или определяться из критериев качества моделирования.

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

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

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

Колебательные компоненты могут быть записаны рядом N M и N A гармоник с частотами q и r, амплитудами Aq, Bq, E r, Fr :

NM ( ) ( ) Sk = Aq sin q k + Bq cos q k, (3.31) M q = NA = Er sin ( r k ) + Fr cos ( r k ). (3.32) SkA r = Нелинейный тренд в модели (3.30) будем моделировать рядом Тейлора степени N T :

NT NT i Tk = Di +1 ( k ). (3.33) i = Известный классический метод сезонной декомпозиции [1, 4, 71] применяется только при сезонной колебательной компоненте, причем тре бует использования длительного интервала наблюдения.

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

Продемонстрируем идентификацию модели (3.30) с помощью пред ложенного метода параметрической итерационной декомпозиции тренд колебательной компоненты [59], который иллюстрирует рисунок 3.19.

Определение нулевого приближения параметров D10,..., DNT +1 тренда Tk при значении колебательных компонент равных нулю: S kM 0 = 0, S kA 0 = Сохранение да Удовлетворен критерий найденных точности?

параметров нет Уточнение параметров колебательных компонент:

уточнение значений колебательной компоненты YkS i = Yk Tki идентификация параметров модели YkS i = Tki 1 S kM i + S kA i + k Уточнение параметров тренда:

уточнение значений тренда Tki = Yk YkS i идентификация i приближения параметров D1i,..., DNT +1 тренда Tki i Увеличение индекса приближений i = i + Рис. 3.19. Алгоритм метода параметрической итерационной декомпозиции тренд-колебательных компонент Параметры D1i, D2,..., D NT +1 входят в модель тренда (3.33) линейно, i i поэтому их поиск легко может быть реализован с помощью МНК.

Результатом использования алгоритма параметрической итерацион ной декомпозиции тренд-сезонной компоненты является сведение задачи идентификации параметров модели (3.30) к идентификации параметров модели:

NM NT ( ) ( ) Di +1 ( k ) NT i + YkS = Aq sin q k + Bq cos q k q =1 i =0 (3.34) NA + Er sin ( r k ) + Fr cos ( r k ) + k.

r = Идентификацию параметров (3.34) проведем, как и ранее, с помо щью ARMA-моделирования. Для вывода общей формулы ARMA модели, соответствующей модели (3.34), вначале рассмотрим частный случай:

( ) 3 Yk = D1 ( k ) + D2 ( k ) + D3 ( k ) + D (3.35) ( A1 sin 1k + B1 cos 1k + A2 sin 2 k + B2 cos 2 k ) + k.

Получим при k 16 следующую ARMA-модель шестнадцатого по рядка:

Yk = 256Yk 8 µ14 µ 2 + 512 (Yk 7 + Yk 9 ) µ1 µ 2 + µ14 µ 4 34 (Yk 6 + 2Yk 8 + Yk 10 ) 384 µ14 µ 2 + 384 µ12 µ 2 + 1024 µ1 µ 2 + 2 4 + (Yk 5 + 3Yk 7 + 3Yk 9 + Yk 11 ) 128 µ14 µ 2 + 128 µ1µ 2 + 768 µ12 µ 2 + 768 µ1 µ 4 3 3 (Yk 4 + 4Yk 6 + 6Yk 8 + 4Yk 10 + Yk 12 ) 1 6 µ14 + 16 µ 2 + 256 µ1 µ 2 + 256 µ1µ 2 + 576 µ12 µ 2 + 4 3 3 + (Yk 3 + 5Yk 5 + 10Yk 7 + 10Yk 5 + 5Yk 11 + Yk 13 ) 32 µ1 + 32 µ 2 + 192 µ12 µ 2 + 192 µ1µ 3 3 (Yk 2 + 6Yk 4 + 15Yk 6 + 20Yk 8 + 15Yk 10 + 6Yk 12 + Yk 14 ) 24 µ12 + 24 µ 2 + 64 µ1µ 2 + (Yk 1 + 7Yk 3 + 2 1Yk 5 + + 35Yk 7 + 35Yk 9 + 21Yk 11 + 7Yk 13 + Yk 15 ) (3.36) [8 µ1 + 8 µ 2 ] ( 8Yk 2 + 28Yk 4 + 56Yk 6 + 70Yk 8 + 56Yk 10 + + + 28Yk 12 + 8Yk 14 + Yk 16 ) + k, где µ1 = C os 1, µ 2 = C os 2, k – остаток, описываемый той же моде лью (3.36), но в уровнях k.

Очевидна нецелесообразность использования модели (3.36) для идентификации: система нормальных уравнений будет нелинейной, имеет высокий порядок, требует больших выборок для обеспечения удовлетворительной точности: до 48 или 64 наблюдений (как уже отме чалось, в 3-4 раза больше порядка ARMA-модели).

Для упрощения записи ARMA-модели (3.36), можно отметить, что коэффициенты в сомножителях при однородных многочленах степеней µ1, µ 2 – биномиальные. Для выявления закономерности в коэффициен тах внутри квадратных скобок коэффициенты при слагаемых µ1m1 µ 2 2 за m пишем в таблицу 3.7 с разложением их на множители.

Таблица 3. Коэффициенты модели (3.36) при слагаемых µ1m1 µ 2 m µ1 µ12 µ13 µ 1 2 2 6 23 4 24 11 µ2 23 4 22 4 4 24 4 4 25 µ 2 2 6 23 6 4 2 4 6 6 25 6 4 2 6 µ 25 4 23 4 24 4 4 26 4 4 27 µ 2 6 1 24 1 25 1 4 27 1 4 28 Из таблицы 3.7 видим, что коэффициенты при степенях µ1, µ 2 мо гут быть записаны как произведения степеней двойки и биномиальных коэффициентов с чередованием знаков. После переноса Yk в правую часть уравнения (3.36) и учета выявленных закономерностей ARMA модель примет вид:

m,m µ1m2 µ 2 2 + k = 0, m (3.37) 12 m1,m2 = где m1 + m m1,m2 = ( 1) 2 m1 + m2 C4 1 C4 2 m1,m2, m m 42[ m1 + m2 ] C l m1,m2 = 42[ m1 + m2 ] Yk 2( 42)( m1 +m2 )2l, l = C ij – биномиальные коэффициенты.

Конструируя ARMA-модели с помощью программы Maple, можно отметить, что ARMA-модели, соответствующие динамическим моделям вида (3.35), сводимы к виду (3.37).

Проводя аналогичные рассуждения с включением колебательных компонент аддитивной структуры, получим общий вид ARMA-модели, соответствующей динамической модели колебательной компоненты YkS вида (3.34):

NT + 1 NM NA m1,m2,...,mN, p1, p2,..., pN µq rpr + k = 0, (3.38) mq m1,m2,...,mN =0 M A p1, p2,..., p N =0 q =1 r = A M где NM NA NM NA mq + pr +1 mq + pr + N M mq q =1 q =1 r =1 r = CN + = ( 1) m1,m2,...,mN, p1, p2,..., pN A (3.39) q = T M m1,m2,...,mN, p1, p2,..., p N A, M m,m,...,m =, p1, p2,..., pNA 1 2 NM NM NA ( NT +1) NM +NA mq + pr q=1 (3.40) l r = = Y S, C NM NA NM NA ( NT +1) NM +NAmq +pr k2(( NT +1) NM +NA ) mq +pr 2l l = q=1 r =1 q=1 r = – биномиальные коэффициенты, µ q = C os q, r = C os r, C ij k k, k 1,..., k N MaxRegr, – остаток, линейный относительно N MaxRegr = 2 ( N T + 1) N M + N A – порядок ARMA-модели.

Соответствие соотношений (3.38)-(3.40) динамической модели (3.34) проверено в практически важных случаях для всех сочетаний N T, N M, N A 3. Можно предполагать, что (3.38)-(3.40) будут соответст вовать модели (3.34) и при большей степени полинома тренда и при большем количестве гармоник.

Коэффициент при YkS в ARMA-модели (3.38)-(3.40) всегда равен 1, поэтому путем переноса слагаемого YkS в правую часть выражения (3.38) ARMA-модель запишем в явном виде. Идентификация парамет ров модели (3.38)-(3.40) возможна при длине выборок не менее N MaxRegr + 1.

ARMA-модель (3.38)-(3.40) позволяет найти переменные µq,r (косинусы частот гармоник колебательных компонент) с помощью МНК. Уравнение (3.38) – полином, относительно степеней µq,r, по этому соответствующий МНК представляет собой систему полиноми альных уравнений относительно µq,r. Метод решений таких систем в базисах Гребнера показан в разделе 3.1.

В рассматриваемой задаче можно привести и частный, но практи чески распространенный класс моделей с кратными частотами:

q = Kmq, r = Kar, (3.41) где Kmq, Kar – вектора натуральных чисел, задающие кратность частот q, r частоте.

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

Замена (3.41) позволяет перейти от µ q, r к одной неизвестной µ = cos ( ), а в МНК – от системы полиномиальных уравнений к более простому одному алгебраическому уравнению относительно µ, методы решения которых известны.

Подстановка (3.41) в функцию косинуса требует преобразования косинусов кратных частот. Использование тождества теории функций комплексного переменного:

cos n + i sin n = [ cos + i sin ], n (3.42) при раскрытии бинома с учетом Cos 2 = 1 Sin 2 позволяет получить уравнение замены µ n = cos n через µ :

floor ( n/2) n,l µ n2l, µn = (3.43) l = где floor ( n /2 ) 2 s l ( 1) Cn2sCsl, n,l = (3.44) s = C s – биномиальные коэффициенты, floor ( 1 n ) – целая часть числа l 1n.

Аналогично проводится замена для переменных r, определяю щих частоты аддитивно входящих гармоник. Замена первых шести кратных частот дает:

µ1 = µ, µ3 = 4µ 3 3µ, µ5 = 16µ 5 20µ 3 + 5µ, µ2 = 2µ 2 1, µ4 = 8µ 4 8µ + 1, µ6 = 32µ 6 48µ 4 + 18µ 2 1.

Замены (3.43) имеют вид полиномов. Следовательно, их подста новка в (3.38) приводит к вычислению степеней и произведений этих полиномов:

NS NF rpr, m µq q q=1 r = что, не умаляя общности результатов, сводит задачу к поиску произве дения полиномов:

(a µ )( ) + a2 µ m11 +... + a m1µ + am1+1 b1µ m 2 + b2 µ m 21 +... + bm 2 µ + a m 2+1, m которое, в свою очередь, вычисляется как полином c1µ m3 + c2 µ m31 +... + cm3µ + cm3+1, где вектор c – свертка векторов a и b.

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

Итак, идентификация модели (3.34) с кратными частотами (3.41) приводит к задаче идентификации параметра µ ARMA-моделей (3.38) (3.40), сводимой с помощью замен (3.43) и операции свертки для каж дой точки Yk, к виду k,1µ M + k,2 µ M 1 +... + k,M µ + k,M +1 = 0. (3.45) Параметр µ находится с помощью МНК из алгебраического урав нения:

N (k,1µ M +... + k,M µ + k,M +1 ) k =1. (3.46) ( ) + ( M 1) k,2 µ M 1 M M k,1µ +... + k,M = Решение уравнения (3.46) относительно µ – стандартная задача, реализуемая процедурами на многих языках программирования.

По корням данного уравнения можно определить µ, затем arccos µ и частоты q, на основе обратных замен (3.41).

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

Для повышения эффективности оценки частоты и повышения устойчивости ее МНК-идентификации целесообразно укрупнение шага дискретизации путем прореживания выборки, которое уже было рассмотрено в разделе 3.1, или, для сравнения, с помощью еще одного приема: консолидации значений (сглаживания нескольких значений) на определенном диапазоне по K diskr точкам с помощью усреднения абс цисс и ординат [53]:

K diskr Y( k 1)K = K diskr, Yk = j.

diskr + K diskr j = Данный прием по своим возможностям довольно близок к прореживанию: уменьшает автокорреляцию и остаточную дисперсию, что будет подробно показано ниже на примере одной из моделей ЖЦП.

Второй этап идентификации модели колебательной компоненты (3.34) заключается в идентификации амплитуд соответствующих гармоник Aq, Bq, Er, Fr. Зная частоты q,, идентификацию амплитуд r Aq, Bq, E r, Fr, входящих в модель (3.34) линейно, можно провести с помощью МНК.

Примером применения модели (3.30) может служить описание макроэкономических параметров–индекса реальных денежных доходов населения (HHI_Q_DIRI) и индекса реального ВВП (GDP_Q_DIRI). Рас сматривая поведение показателей на периоде между кризисами 1998 и 2008 годов, можно сравнить их динамику на первой и второй половинах периода – в 1999-2003 и в 2003-2007 годах.

Для описания индекса HHI_Q_DIRI лучшей моделью по значению коэффициента детерминации оказалась:

Yk = ( D1k + D2 ) (3.47) (1 + A1 sin k + B1 cos k + A2 sin 2 k + B2 cos 2 k ) + k.

Результаты описания HHI_Q_DIRI моделью (3.47) представлены на рисунке 3.20 и в таблице 3.8. Параметры модели (3.47), построенной по значениям индекса в 1999-2003 и в 2003-2007 гг., и амплитуды каж 2 дой из гармоник Vi = Ai + Bi представлены в таблице 3.9.

140 100 60 123412341234123412341234 а) б) 1999 2000 2001 2002 2003 2004 2003 2004 2005 2006 2007 Рис. 3.20. Индекс реальных денежных доходов населения HHI_Q_DIRI, 1993.04=100: а) 1999-2004 гг., б) 2003-2008 гг., маркеры – данные Росстата, сплошная линия – модель (3.47), точки – прогноз по модели (3.47) Таблица 3. Индекс HHI_Q_DIRI 1999 2000 2001 2002 2003 Год R2 MAPE Квартал 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 HHI_Q_DIRI, 64,8 71,4 72,2 87,8 74,8 84,3 85,9 96,3 80,6 90,5 95,5 104,8 88,5 99,1 104,1 120,9 104,6 114,3 117,6 138,7 118,3 124,0 128,8 153, 1992.04= Модель 1999-2003 + 62,9 71,0 77,3 81,6 71,6 82,4 86,2 96,3 81,1 93,2 95,2 110,1 92,1 103,0 104,9 122,5 105,1 111,6 115,5 133,0 120,2 119,0 127,1 141,6 0,9710 3,7% прогноз на 2004 год 2003 2004 2005 2006 2007 Год Квартал 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 HHI_Q_DIRI, 104,6 114,3 117,6 138,7 118,3 124,0 128,8 153,6 125,2 142,8 147,0 175,4 137,7 165,0 168,3 199,2 154,0 181,4 189,2 225,3 165,6 192,0 197,7 209, 1992.04= Модель 2003-2007 + 97,3 111,3 122,2 131,6 111,2 129,7 136,7 155,0 125,7 147,8 150,8 178,2 141,6 165,0 165,2 200,2 159,4 180,8 180,4 220,3 179,8 195,1 196,9 238,0 0,9747 6,0% прогноз на 2008 год Таблица 3. Параметры модели (3.47) D1 D2 omega A1 B1 A2 B2 V1 V 1999-2003 год 2,7445 68,5959 1,5082 -0,0193 -0,0709 0,0536 -0,0117 0,0734 0, 2003-2007 год 4,3991 108,2340 1,5254 -0,0349 -0,0804 0,0685 -0,0203 0,0876 0, Частота и амплитуды V1,V2 гармоник колебательной компонен ты близки между собой, а угловой коэффициент тренда возрастет в D120032007 / D1 2003 = 1,6 раз (таблица 3.9).

Индекс реального ВВП (GDP_Q_DIRI) лучше описывается моде лью:



Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 7 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.