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

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

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


Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |

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

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

Yk = D1k + D2 + (3.48) + E1 sin ( k ) + F1 cos ( k ) + E2 sin ( 2 k ) + F2 cos ( 2 k ) + k.

Результаты описания индекса реального ВВП моделью (3.48) представлены на рис. 3.21 и в таблице 3.10. Параметры построенных моделей и амплитуды каждой из гармоник Vi = Ei2 + Fi 2 приведены в таблице 3.11.

а) б) 2003 2004 2005 2006 2007 1999 2000 2001 2002 2003 Рис. 3.21. Индекс реального ВВП: а) 1999-2004 гг., GDP_Q_DIRI, 1994.01=100, б) 2003-2008 гг., GDPEA_Q_DIRI, 2003.01=100, маркеры – данные Росстата, сплошная линия – модель (3.48), точки – прогноз по модели (3.48) Частота и амплитуда второй гармоники V2 колебательной ком поненты остаются неизменными, амплитуда первой гармоники V1 в 2003-2007 гг. увеличивается.

Угловой коэффициент тренда возрастает при этом в D120032007 / D1 2003 = 1,32 раза (таблица 3.11).

Таблица 3. Индекс реального ВВП 1999 2000 2001 2002 2003 Год R2 MAPE Квартал 1 23 4 1 23 4 1 23 4 1 23 4 1 23 4 1 23 GDP_Q_DIRI, 88,5 94,0 109,2 105,0 98,6 103,6 120,7 113,6 103,2 108,8 127,9 118,7 107,1 113,4 133,6 125,9 115,2 122,4 142,2 135,4 123,6 131,6 152,3 144, 1994.01 = Модель 1999-2003 + 87,7 94,9 111,3 107,4 94,7 102,5 118,8 114,3 101,4 109,3 125,4 120,8 108,4 117,0 132,9 127,6 115,1 123,8 139,5 134,1 122,1 131,4 147,0 141,0 0,9805 1,8% прогноз на 2004 год 2003 2004 2005 2006 2007 Год Квартал GDPEA_Q_DIRI, 100,0 106,9 121,4 120,7 107,0 115,2 130,2 128,7 112,8 122,3 138,3 138,4 120,2 131,7 149,1 150,2 129,3 142,2 160,5 163,6 140,5 152,9 170,2 165, 2003.01= Модель 2003-2007 + 97,4 103,0 122,8 121,8 106,8 112,3 131,3 129,5 115,2 122,2 141,6 139,3 124,7 131,5 150,1 147,1 133,1 141,4 160,4 156,8 142,6 150,7 168,8 164,6 0,9756 1,1% прогноз на 2008 год Таблица 3. Параметры модели (3.48) описания индекса реального ВВП D1 D2 omega E1 F1 E2 F2 V1 V 1999-2003 год 1,7364 97,9032 0,7888 -0,2327 -0,0640 -4,6360 -10,1101 0,2413 11, 2003-2007 год 2,2951 107,5953 0,7902 0,4622 -0,2732 -7,0760 -9,9676 0,5369 12, Внешне кривые на рисунках 3.2, 3.20 и 3.21 похожи, однако ис пользуемый инструментарий различает их.

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

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

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

Модели (3.47), (3.48) с высокой точностью описывают многие макроэкономические показатели: например, индекс реального ВВП, индекс реальных денежных доходов населения: R 2 для них превыша ет 0,97.

MAPE-оценка прогноза на четыре квартала 2004 и 2008 годов не превышает 2% при описании индекса реального ВВП и не превышает 6% при описании индекса реальных денежных доходов населения.

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

Так, скорость роста индекса реальных денежных доходов насе ления в 1,58 раза превышала скорость роста индекса реального ВВП в 1999-2003 годах, и в 1,92 раза в 2003-2007 годах.

Известно утверждение, что зарплатоемкость российского ВВП в 2006 году достигла 33,3% процентов по сравнению с 23,6% в 2000 году. Опережающее в два раза повышение уровня оплаты труда по сравнению с ростом ВВП является тревожным сигналом для эко номики и указывает на необходимость ускоренной модернизации производства для обеспечения роста производительности труда, со поставимого с повышением уровня доходов населения.

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

( ) Yk = D1 ( k ) + D2k + D (1 + A1 sin ( k ) + B1 cos ( k ) + A2 sin ( 2 k ) + B2 cos ( 2 k ) ) + (3.48) + E1 sin ( k ) + F1 cos ( k ) + k.

Результаты идентификации модели (3.49) представлены на рисунке 3.22 и в таблице 3.12.

Параметры идентифицированных моделей и амплитуды каждой из гармоник V1 = A12 + B12, V2 = A2 + B2, V3 = A3 + B3 представлены в 2 2 2 таблице 3.13.

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

180 160 140 120 100 80 60 40 20 2000 2001 2002 2003 2004 2005 2006 2007 Рис. 3.22. Объем инвестиций в основной капитал в Самарской области, млн руб., маркеры – данные 2000-2008 гг., сплошная линия – модель (3.49) по данным 2000-2007 гг., точки – прогноз по модели (3.49) на 2008 год Таблица 3. Объем инвестиций в основной капитал в Самарской области 2000 2001 2002 2003 Год Квартал 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 Инвестиции в основной капитал 5 788 10 086 24 855 22 871 5 898 4 314 21 407 24 336 7 020 4 013 22 369 26 479 9 486 7 625 27 287 31 512 13 622 12 408 36 341 42 Самарской обл., млн. руб.

Модель 2000-2007 + 3 856 13 560 28 723 18 601 1 549 5 044 25 139 23 720 5 084 3 357 22 539 27 422 10 574 8 223 24 735 31 290 15 003 16 575 34 104 38 прогноз на 2008 год 2005 2006 2007 2008 Год Квартал 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 Инвестиции в основной капитал 18 468 16 188 51 061 59 461 22 807 20 108 71 071 80 768 29 570 25 989 91 039 106 618 44 885 30 340 106 568 142 Самарской обл., млн. руб.

Модель 2000-2007 + 18 027 23 671 50 405 53 516 22 532 25 313 70 132 78 101 33 658 20 410 86 871 112 662 56 869 12 626 93 118 154 прогноз на 2008 год 26% 0,9828 MAPE прогноза на R Таблица 3. Параметры модели (3.49) для описания объема инвестиций в основной капитал в Самарской области Тренд Адитивная компонета D1 D2 D3 E1 F1 V 83 -915 16 137 -14 055 -22 942 26 Мультипликативная компонента A1 B1 A2 B2 V1 V 0,9333 0,6217 -0,0136 0,0389 1,1214 0, Программа для идентификации модели (3.30) использована и для анализа динамики изменения средних потребительских цен на бензин на территории России (рис. 3.23). Применяя метод перебора (определяя последовательно степени полинома, количество аддитивных и мультип ликативных колебательных компонент, а также кратности их частот), выберем ту модель, для которой оценки точности моделирования (ко эффициента детерминации) и прогнозирования (MAPE-оценки) лучше.

Рис. 3.23. Полиномиальная модель тренда с аддитивно-мультипликативными колебательными компонентами динамики изменения цен на бензин на территории РФ Оказалось, что более точная модель состоит из полинома 1-й сте пени, 2-х гармоник аддитивной колебательной компоненты и 1-й гармо ники пропорционально мультипликативной колебательной компоненты с кратностью частот 3.

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

MAPE–оценка точности прогноза составила 0,9%, улучшая тем самым результат модели с полиномиальным трендом пятой степени более, чем в 20 раз, и для полинома второй степени – более, чем в 3 раза.

Рассмотрим возможность увеличения точности моделирования и прогнозирования за счет мониторинга эволюции модели.

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

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

Увеличилась и точность моделирования и прогнозирования: ко эффициент детерминации вырос – с 0,98305 до 0,9968 на первом отрез ке, до 0,9831 – на втором, а MAPE-оценка точности прогноза увеличи лась с 0,9% до 0,14% на первом отрезке, до 0,84% – на втором.

а) б) Рис. 3.24. Мониторинг эволюции полиномиальной модели тренда с аддитивно-мультипликативными колебательными компонентами динамики изменения цен на бензин на территории РФ: а) период с 15.06.2010 г.

по 25.10.2010 г.;

б) период с 1.11.2010 г. по 21.03.2011 г.

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

Наиболее точная модель динамики, описывающая весь интервал наблюдения, имеет вид полинома 4-ой степени с 2-мя аддитивными и 1-ой пропорционально мультипликативной гармониками (рис. 3.25).

Проводя визуальный анализ модели на рисунке 3.25, находим точку эволюции (рис. 3.26): для первого отрезка времени более точной модели соответствует полином 3-ей степени с 3-мя аддитивными гармониками и 1-ой пропорционально мультипликативной гармоникой, а на втором отрезке времени модель эволюционирует, принимая вид полинома 2-ой степени с 2-мя аддитивными и 1-ой пропорционально мультипликатив ной гармониками.

Рис. 3.25. Полиномиальная модель с аддитивно-мультипликативными колебательными компонентами динамики изменения цен на бензин на территории Самарской области (руб./л) Можно заключить, что разбиение интервала наблюдения и после довательная идентификация 2-х моделей оправдало себя в силу значи тельного увеличения точности моделирования и прогнозирования (таб лица 3.14).

а) б) Рис. 3.26. Мониторинг эволюции полиномиальной модели с аддитивно мультипликативными компонентами динамики изменения цен на бензин на территории Самарской области (руб./л): а) период с 15.06.2010 г. по 25.10.2010 г.;

б) период с 1.11.2010 г. по 21.03.2011 г.

Таблица. 3. Оценки моделирования и прогнозирования полиномиальными моделями с аддитивно-мультипликативными компонентами динамики изменения цен на бензин на территории Самарской области Период 15.06.2010 г.- 15.06.2010 г.- 1.11.2010 г. Оценка 21.03.2011 г. 25.10.2010 г. 21.03.2011 г.

Коэффициент 0,985 0,987 0, детерминации MAPE-оценка 1,66% 0,52% 0,93% Покажем, что модель (3.30) применима не только к российской статистике изменения цен на бензин, но и к аналогичным зарубежным показателям как макро-, так и мезоуровня.

Примером может служить идентификация модели, описываю щей динамику изменения средних потребительских цен на бензин на территории США за период с 15.06.2010 г. по 21.03.2011 г. (дол лар/галлон). Результат соответствующего моделирования цен на бен зин, а также мониторинга эволюции модели представлены на рисун ках 3.27 и 3.28.

Более точная модель на всем интервале наблюдения имеет вид полинома 2-ой степени, в то время как на отдельных участках она эволюционирует от полинома 1-ой степени до 3-ей степени. Также изменяется число аддитивных и мультипликативных гармоник. Точ ность моделирования и прогнозирования в обоих случаях увеличива ется. Аналогичные исследования относительно мезоуровня СЭС в США, округ Колумбия, дали схожие результаты.

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

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

Рис. 3.27. Полиномиальная модель с аддитивно-мультипликативными колебательными компонентами динамики изменения цен на бензин на территории США а) б) Рис. 3.28. Мониторинг эволюции полиномиальной модели с аддитивно мультипликативными компонентами динамики изменения цен на бензин на территории США: а) период с 15.06.2010 г. по 25.10.2010 г.;

б) период с 1.11.2010 г. по 21.03.2011 г.

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

Таблица. 3. Оценки моделирования и прогнозирования полиномиальными моделями с аддитивно-мультипликативными компонентами динамики изменения цен на бензин на территории округа Колумбия, США.

Период 15.06.2010 г.- 15.06.2010 г.- 1.11.2010 г. Оценка 21.03.2011 г. 25.10.2010 г. 21.03.2011 г.

Коэффициент детерминации 0,985 0,995 0, MAPE-оценка 1,99% 1,865% 1,889% 3.4. Модели роста в виде суммы полиномиального тренда и гармоник с независимо эволюционирующими моделями амплитуды колебательной компоненты Рассмотрим возможность идентификации еще одного возможного усложнения модели роста для моделирования эволюции: в виде суммы полиномиального тренда и гармоник с независимо эволюционирующи ми от тренда моделями амплитуды колебательной компоненты.

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

линейный (1.10):

C k = ( A0 + A1k ) sin ( k + ) ;

экспоненциальный (1.11):

Ck = Ae k sin ( k + ) ;

обобщенный экспоненциальный (1.12):

( ) C k = A0 + A1e k sin ( k + ) ;

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

Ck = A1e 1k sin ( k + ) + A2e 2k sin ( k + ).

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

NT N T i Yk = Di +1 ( k ) + Ck + k, i = не будет связан с уровнями тренда ( k удовлетворяет принятым услови ям Гаусса-Маркова).

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

При тренде в виде линейной функции и колебательной компонен те (1.10) будем иметь временной ряд:

Yk = A0 + A1k + ( A2 + A3 k ) sin ( k + ) + k, (3.50) которому при k 6 соответствует ARMA-модель:

Yk = 2Yk 1 3Yk 2 + 4Yk 3 3Yk 4 + 2Yk 5 Yk 6 + + 1 (2Yk 1 4Yk 2 + 4Yk 3 4Yk 4 + 2Yk 5 ) + (3.49) + 12 (Yk 2 2Yk 3 + Yk 4 ) + k, где 1 = 2 cos, а центрированную, автокоррелированную, гомоскеда стическую стохастическую компоненту k определит то же соотноше ние (3.51), но в уровнях k.

Компенсацию автокоррелированности k можно осуществить предложенными выше приемами, а оценку параметра 1 даст МНК идентификация (3.51).

Используя МНК-оценку 1, можно вычислить МНК-оценку по формуле = arccos.

Подставляя в (3.50), легко найти МНК-оценки остальных пара метров A0, A1, A2, A3, модели (3.50), решая нормальное СЛАУ шесто го порядка и проведя соответствующие расчеты.

Более сложной и более общей по приложениям в экономической практике может явиться модель:

Yk = A0 + A1k + A2 sin (1k + 1 ) + A3k sin ( 2 k + 2 ) + k. (3.50) Нетрудно получить при k 8 ARMA-модель восьмого порядка для поиска МНК–оценок ее коэффициентов 1 и 2 и последующего расче 1 1 1 и 2 = arccos.

= arccos та 2 Для обеспечения устойчивости решения в этом достаточно слож ном случае следует рекомендовать метод параметрической итеративной декомпозиции.

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

Достаточно выполнения нескольких (обычно, как показала прак тика, трех-четырех) итераций.

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

Yk = A0 + A1k + A2 sin (1k + 1 ) + A3e k sin ( 2 k + 2 ) + k, (3.51) где k – стохастическая компонента.

Целесообразно для (3.53) применить метод параметрической ите рационной декомпозиции.

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

Расчет оставшихся параметров модели (3.53) возможен из нор мальной СЛАУ шестого порядка при применении МНК.

Покажем возможность применения модели с эволюцией амплитуды одной колебательной компоненты и предложенный метод ее идентифика ции на микро-уровне СЭС: на реальных данных объема продаж лекарст венных средств «Pinkham Medicine Company» [114].

На рисунке 3.29 представлены графики объема продаж и модельного ряда:

Yk = 216105 251,761k + 1104517e 0,005k sin(0,0522k 0,798) + k.

Точностные характеристики модели: коэффициент детермина ции 0,77;

MAPE-оценка с горизонтом прогноза в одно наблюдение равна 0,2%, а MAPE-оценка прогноза на 1 3 выборки равна 10%. Можно пока зать, что данная модель дает в сравнении только с одной линейной моде лью тренда выигрыш MAPE-оценки прогноза в 4%.

Рис. 3.29. Моделирование объема продаж лекарственных средств «Pinkham Medicine Company»

Пример использования той же модели (3.53) (с добавлением еще одной гармоники) дадим на данных по ВНП США (в миллиардах дол ларов), взятых с 1951 по 1984 годы [114]. На рисунке 3.30 представлен график моделирования ВНП США выражением Yk = 633,77 + 36,4k + 8,5e0,0642k sin(0,9k 2,952) + + 16,1sin(0,6k 1,3) + k, для которого коэффициент детерминации оказался равен 0,997, а MAPE-оценка на горизонт наблюдения в три года равна 1%.

Примеры усложнения моделей можно продолжить: увеличивая порядок полиномов до второго, используя, наряду с моделью (1.10), (1.12) и модель (1.13) эволюции амплитуды гармоники.

Млрд. долларов Рис. 3.30. Моделирование ВНП США моделью с эволюцией амплитуды колебательной (циклической) компоненты Показанные приемы параметрической итеративной декомпози ции, использование базиса Гребнера демонстрируют во всех случаях достаточно высокую точность моделирования и прогнозирования предложенными моделями с предложенными моделями эволюции амплитуды колебательной компоненты.

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

Yk = (Tk + C k ) µ k = (Tk + C k )(1 + k ) = Tk + C k + (Tk + C k ) k, которая, в случае моделирования детерминированной компоненты сум мой линейной и гармонической колебательной компонентой, примет вид:

Yk = { A0 + A1k + A3 cos ( k + )}(1 + k ), (3.52) где k – стохастическая компонента удовлетворяющая принятым усло виям Гаусса-Маркова, причем 1 k 1.

Конструирование ARMA-модели в этом случае будет произво диться путем применения Z -преобразования к детерминированной Y компоненте Dk = k последующих действий по группировке и обрат 1+ k ному Z -преобразованию.

Получим при k 4 следующую ARMA-модель четвертого порядка для (3.54):

Yk = 1 (Yk 1 2Yk 2 + Yk 3 ) + 2 (Yk 1 Yk 2 + Yk 3 ) Yk 4 + k, (3.53) где:

k k k 1 k k = 1 Yk 1 2Yk 2 k + Yk 3 k + 1 + k 1 1 + k 2 1 + k k k 1 k 2 k k 4 – +2 Yk 1 k Yk 2 k + Yk 3 k Yk 1 + k 1 1 + k 2 1 + k 3 1 + k «новая» стохастическая компонента, 1 = 2 Cos.

Видим, что для модели (3.55) с мультипликативной стохастиче ской компонентой авторегрессионная часть модели имеет тот же самый вид, что и для модели (3.4) с аддитивной стохастической компонентой.

Различие заключается в стохастических компонентах k и k.

В (3.55) k является не только автокоррелированной (ее компенсация прореживанием показана), но и имеет явно выраженный гетероскеда стический характер: её значения (а, следовательно, и дисперсия) зависят от значений уровней Yk 1, Yk 2, Yk 3, Yk 4.

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

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

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

Для этого разделим (3.55) почленно на одно из значений, входя щих в ARMA-модель, например, на уровень Yk 1. В результате этого действия получим новую ARMA-модель в относительных единицах уровней:

2Y YY Y Y Yk = 1 1 k 2 + k 3 + 2 1 k 2 + k 3 k 4 + k, Yk 1 Yk 1 Yk 1 Yk 1 Yk 1 Yk где k 1 Y ( k 2 ) Yk 3 ( k k 3 ) k 2 k 2 k k = = 1 k + + 1 + k 1 Yk 1 (1 + k 2 ) Yk 1 (1 + k 3 ) Yk 1 k 1 Yk 2 ( k k 2 ) Yk 3 ( k k 3 ) Yk 4 ( k k 4 ) – +2 k + 1 + k 1 Yk 1 (1 + k 2 ) Yk 1 (1 + k 3 ) Yk 1 (1 + k 4 ) относительная невязка.

Тогда МНК-оценку 1 позволит определить условие Y Y Y 2Y Y Y N 1 = arg min k 1 1 k 2 + k 3 2 1 k 2 + k 3 + k 4, 1 k =4 Yk 1 Yk 1 Yk 1 Yk Yk Yk приводящее к линейному уравнению относительно 1.

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

Естественно, что коррекцию гетероскедастичности можно прово дить не только путём почленного деления выражения (3.55) на Yk 1, но и на другие уровни ряда, например, на Yk 2 или на Yk 3, которые, естест венно, не должны быть равными нулю.

После компенсации автокоррелированности приемом прорежива ния, очевидного расчёта по оценке 1 значения оценки частоты = arccos, представим (3.54) в виде:

{ )} (1 + ) = A k + A + ( Yk = A1k + A0 + A3 cos k + 1 k (3.54) { A k + A + A cos ( k + )}, + A4 cos k + A5 sin k + k 1 0 где использованы обозначения A4 = A3 sin, A5 = A3 cos.

{ )} ( k A1k + A0 + A3 cos k + Стохастическая компонента в (3.56) неавтокоррелирована, но вновь гетероскедастична. Для компен сации гетероскедастичности можно разделить левую и правую части этого выражения на величину уровня Yk.

Тогда будем иметь:

cos k sin k k +k, 1 = A1 + A0 + A4 + A5 (3.55) Yk Yk Yk Yk ( ) A1k + A0 + A3 cos k + где k = k – ряд уровней «новой» (относи Yk тельной) стохастической компоненты.

Параметры A1, A0, A4, A5, линейно входящие в (3.57), могут быть определены из нормальной СЛАУ четвёртого порядка, к которому при водит условие:

sin k N cos k k = arg min 1 A A0 A4 A A1, A2, A4, A5.

A1, A0, A4, A5 k =1 Yk Yk Yk Yk Что же касается МНК-оценок амплитуды и фазы гармоники, то они определятся по очевидным и уже знакомым формулам:

2 1/ ( ) +( ) = A4, = arctg ( A5 / A4 ).

A3 A Стохастическую компоненту k в (3.57) представим, после под становки Dk = (1 + k ){ A1k + A0 + A3 cos( k + )} из (3.54), в виде:

{ } k A1k + A0 + A3 cos( k + ) k =. (3.56) (1 + k ){ A1k + A0 + A3 cos(k + )} В силу того, что в (3.58) cos( k + ) 1 и cos(k + ) 1 можно сократить практически одинаковые сомножители в фигурных скобках числителя и знаменателя, прийти к выводу о гомоскедастичности и ма лости величины стохастической компоненты k при оценке параметров A1, A0, A3, :

k k.

1+ k Необходимый объём выборки для идентификации модели (3.54) с ARMA-моделью (3.55) четвертого порядка может составлять от 12 до наблюдений.

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

Yk = ( A0 A1k ) cos ( k + )(1 + k ) (3.57) + будем иметь при k 4 ARMA-модель четвертого порядка:

Yk = 21 (Yk 1 + Yk 3 ) 12Yk 2 Yk 4 2Yk 2 + k, (3.60) где k k 1 k k 3 k Yk 3 12Yk 2 k k = 21 Yk 1 + 1 + k 1 1 + k 3 1 + k k 4 k Yk 4 k 2Yk 2 k 1 + k 4 1 + k гетероскедастическая стохастическая компонента, 1 = 2 cos.

Для уменьшения гетероскедастичности перейдем к относительным величинам путем почленного деления (3.60), например, на Yk 2.

Тогда (3.60) примет вид:

Y Y Yk Y = 21 k 1 + k 3 12 k 4 2 + k, Yk 2 Yk 2 Yk 2 Yk где Y k 1 Yk 3 k k k 2 k k = 21 k 1 k k = + Yk 2 1 + k 1 Yk 2 1 + k 3 1 + k Yk 2 Y k 4 k k 4 k 2 k.

Yk 2 1 + k 4 1 + k В принципе, возможно деление (3.60) и на другие уровни опреде ляемого показателя в ARMA-модели: на Yk, Yk 1, Yk 3, Yk 4.

На первом этапе идентификации, после реализации приема проре живания, будем иметь нормальную СЛАУ второго порядка для очевид 1 и последующего расчета = arccos.

ной МНК-оценки На втором этапе МНК-идентификации, подставив в (3.59), вы полнив преобразование исходного выражения, получим:

( ) (1 + Yk = ( A1k + A0 ) cos k + )= k = A2 k cos k A3 k sin k + A4 cos k A5 sin k + (3.58) ( ) + k ( A1 k + A0 ) cos k +, где A2 = A1 cos, A3 = A1 sin, A4 = A0 cos, A5 = A0 sin, а стохастиче ( ) неавтокоррелирована, но ская компонента k ( A1k + A0 ) cos k + обладает свойством гетероскедастичности.

Для компенсации гетероскедастичности стохастической компо ненты разделим (3.61) почленно на Yk :

cos k sin k k k cos k A3 sin k + A4 + k, (3.59) 1 = A2 A Yk Yk Yk Yk ( ) k ( A1k + A0 ) cos k + где k = – гетероскедастическая стохастиче Yk ская компонента.

Из соотношения (3.62), решая соответствующую нормальную СЛАУ четвертого порядка, найдем оценки оставшихся параметров мо дели (3.59):

2 1/2 2 1/ ( ) +( ) ( ) +( ) 2 = A2 ( ), = arctg A5 / A4.

= A, A1 A3 A0 A Относительная невязка при параметризации модели (3.59) будет равна:

( )= ( ) k ( A1k + A0 ) cos k + k ( A1k + A0 ) cos k + k = = ( A1k + A0 ) cos ( k + ) (1 + k ) Yk ( ) k cos k + k =.

cos ( k + ) (1 + k ) 1 + k Видим, что стохастическая компонента k гомоскедастична. Од нако проведенные исследования на тестовых моделях (3.54) и (3.59) по уже показаннй выше методике оценки точности показали высокие пока затели точности моделирования на выборках от 12 до 16 наблюдений при соотношении Kn/s не более 10%. Характеристики точности прогно зирования на тестовых моделях несколько выше.

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

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

k 0, –1 –0,5 k 0 0, Рис. 3.31. Зависимость относительной среднеквадратической невязки k от стохастической компоненты k в моделях (3.54) и (3.59) Отметим, что при положительных значениях k соответствующие значения k остаются малыми (в пределах 0,5), а при k 0,5 абсолют ные значения k существенно возрастают. Таким образом, соотношение k k выполняется лишь при малых значениях k (менее 0,1) (или 1+ k 10% своего возможного диапазона).

Рассмотрим теперь статистические характеристики k в сравнении со статистическими характеристиками k более подробно (табл. 3.16).

Затем перейдем теперь к сравнению статистических характеристи ки k для частных случаев распределения k. Естественно начать с наи более часто рассматриваемого в научной литературе случая нормально го (центрированного и нормированного) симметричного распределения вероятностей стохастической компоненты N 0;

(табл. 3.17).

Таблица 3. Сравнение статистических характеристик k, k k k k ( 1;

1) k ( ;

0,5 ) Ограничения 1 x F ( x ) = F 1 = F Закон F ( x ) 1 x 1 x распределения 1 1 Плотность f ( x ) = f ( x ) f (1 x )2 1 x распределения Симметричность да нет распределения 0, M [ ] = 0 M [ ] = xf ( x )dx Матожидание D [ ] = const D [ ] = const D [ ] Дисперсия F ( 0 ) = 0, F ( 0 ) = 0, Медиана arg max ( F ( x ) ) = 0 ( ) arg max F ( x ) Мода x x Сравним теперь статистические характеристики распределения k со статистическими характеристиками (нормированного и центриро ванного) симметричного бета-распределения: 2 ( ;

) 1, = (примем = 4 ) (таблица 3.18).

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

Распределение k несимметрично: положительные значения огра ничены 0,5, а отрицательные не ограничены.

Таблица 3. Сравнение с нормальным распределением k k k ( 1;

1) лишь с вероятностью Ограни k ( ;

0,5 ) чения 0, ( ;

) ( ;

1) (1;

) F( ) F() Закон распределения 0, 0, 0 0. 0 0.5 –1 –0. –1 –0. F( ) F() Плотность распределения 0, 0, 0,5 – –1 M 1 –1 –0,5 0 0. Поскольку ограничение не выпол няется, существует разрыв в точке x = 1. Можно принять f (1) = M [ ] = xf ( x )dx + Матожидание + xf ( x )dx 0, M [ ] = 0 Значение M [ ] зависит от D [ ] :

чем больше D [ ], тем больше и аб солютное значение M [ ] Диспер По эмпирическим расчетам D [ ] D [ ] = 2 = сия превышает D [ ] в десятки раз Таблица 3. Сравнение с бета-распределением k k Ограни чения k ( 1;

1) k ( ;

0,5 ) F( ) F() Закон распределения 0,5 0, 0 0. –1 –0.5 0 0. –1 –0. f () f( ) Плотность распределения –1 1x –0,5 0 0, M –2 –1 0 0,5 f ( ) = 0 при 0, f ( ) = 0 при ( 1;

1) Матожи 1 0, дание M [ ] = xf ( x )dx = M [ ] = 0 = 1 + D [ ] = D [ ] ( 1) ( 1) D [ ] Дисперсия 1 D [ ] D [ ] = 4 ( 2 + 1) 0 1 2 3 4 5 6 7 8 Фактически это означает, что в сумме квадратов отклонений больший вес имеют те наблюдения, где реальные данные меньше мо дельных.

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

При относительной мощности мультипликативной стохастической помехи большей 10% целесообразно переходить от операции деления на уровни ряда динамики (по сути применения при МНК веса Wk = )к D(Yk ) взвешенному МНК. При этом полученные на «классическом» МНК оценки модели задают весовую функцию Wk =, {( )} A11 k A0 1 ) co s( k + + которую и будем использовать при следующем проходе для оценки взвешенных весов в МНК:

W k z +1 =.

{( )} A1 z k A0 z ) Cos( z k z + + Исследования показали возможность применения более сложного взвешенного МНК при 3-4 итерациях до соотношения Kn/s = 25-30%, ровно как показано в следующем параграфе для другой модели.

С помощью аддитивной модели (3.1) и мультипликативной моде ли (3.54) моделировалась динамика минимальных цен и максимальных цен на зерновые культуры в хозяйствах Самарской области.

Выяснилась определенная закономерность: например, что мульти пликативная модель для зерновых культур, обладая несколько худшей точностью моделирования (например, R2 = 0,89 для минимальных цен на рожь) по сравнению с аддитивной моделью ( R2 = 0,94 для минималь ных цен на рожь) обеспечила существенно лучшие показатели прогноза:

MAPE-оценка 0,18% для мультипликативной модели против 14,7% для аддитивной (рис. 3.32).

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

Уровень минимальной цены 04.2004 08.2004 11.2004 02.2005 05.2005 08.2005 09.2005 03. цена мультипликативная модель аддитивная модель Рис. 3.32. Динамика минимальных цен на рожь в Самарской области и результаты моделирования Визуальный анализ динамики цен на огурцы в РФ за 1998 2008 гг. (Росстат. URL: http://www.gks.ru.) показал, что в качестве модели могут быть рассмотрены или линейный тренд с гармониче ской колебательной компонентой аддитивной, или линейный тренд с мультипликативной гармонической компонентой:

Yk = (3702,6 k + 4519,8)(1 + 0,125Sin(0,6914 k + 0,17) + k, Yk = 3534, 2 k + 5102,6 + 1898Sin(0,9052 k 1,125) + k.

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

Таблица 3. Сравнение моделей мультипликативной и аддитивной структур для моделирования средних цен на огурцы в РФ Год 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 R2 MAPE Цены на огур 5066 10469 12436 16090 18628 22465 24321 26972 31383 36857 46893 цы, руб./тонна Мультипл.

модель, по 1998-2008 г. и 4615 9002 13415 17154 19824 21689 23650 26801 31795 38360 45274 50886 0,9920 1,9% прогноз на Аддитивная модель, по 1998-2008 г. и 3390 8223 13372 17602 20382 22287 24564 28175 33061 38187 42337 45039 0,9728 13,2% прогноз на Для описания цен за 1998-2009 гг. (использовано уже известное значение цены за 2009 г.) модель мультипликативной структуры будет иметь несколько другие параметры:

Yk = (3704, 6 k + 4459,9)(1 + 0,129Sin(0, 6914 k 0, 43) + k.

R 2 данной модели равен 0,992, а прогнозное значение средней це ны на 2010 год составляет 53011 руб./тонну.

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

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

Метод идентификации полиномов Dk первого и второго порядков, в которых стохастическая компонента аддитивна:

Yk = Dk + k, (3.60) известен [1, 12, 71] в форме линейной регрессии с использованием МНК, как и точечные, и интервальные оценки точности моделирования и прогнозирования.

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

Yk = D k (1 + k ) = Dk + Dk k = D k + k (3.61) где Dk – полиномы первого или второго порядков, k – стохастическая компонента, удовлетворяющая принятым условиям Гаусса-Маркова и неравенству 1 k 1, где k – «новая» гетероскедастическая стохас тическая компонента, у которой дисперсия зависит от номера наблюде ния:

2 22 2 2 D[ k ] = M [ k ] = M Dk µ k = Dk M µ k = Dk D[ µ k ].

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

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

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

n 1 D[ ] (Yk Dk ) min. (3.62) k =1 k Подставив (3.65) в условие МНК для (3.64), будем иметь:

n n 1 2 D[ ] (Yk Dk ) = D2 D[µ ] (Yk Dk ).

k =1 k =1 k k k Поскольку D[ k ] = const, то n D 2 (Yk Dk )2 min. (3.63) k =1 k Однако величина весов Dk заранее неизвестна, а оценить их вели чину можно одним из следующих способов.

Учитывая, что величина ошибки, как правило, мала по сравнению с уровнем ряда, в качестве оценок весов можно использовать сами уровни ряда Dk Yk2 [57]. Идентификацию можно при этом проводить в два этапа. На первом этапе модель идентифицируется с помощью МНК без весов. В результате получаются оценки Dk, которые затем использу ются в МВНК.

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

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

Yk = ( C 0 + C1k )(1 + k ), ( ) (1 + ).

Yk = C0 + C1k + C2 ( k ) k Для проведения компьютерного эксперимента зададим конкрет ные (истинные) значения параметров модели, длину выборки и величи ну случайной помехи, генерируемой с помощью датчика случайных чи сел.

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

Величину помехи будем задавать с помощью коэффициента «шум/сигнал», который задается для мультипликативных моделей в процентах (или в долях), отражая соотношение мощности помехи к мощности детерминированной части модели. Задавая конкретные зна чения K n / s, получим величину дисперсии помехи для конкретной моде ли. Будем генерировать выборки с уровнем шума 1%, 5%, 10%, 20%, 30%. Уровень шума до 10% будем считать малым.

Тот же показатель можно использовать для оценки точности мо делей, полученных в результате идентификации. Такой коэффициент «шум/сигнал» будем обозначать K n/ s. Логично предположить, что если модели идентифицируются верно, то в среднем K n/ s = K n/ s.

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

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

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

Разброс оценок параметров характеризуется коэффициентом ва риации, для оценки точности модели можно, как было указано выше, использовать K n/ s. Для характеристики точности прогноза будем ис пользовать второй коэффициент Тейла K T 2.

Прогноз с K T 2 20% считается достаточно точным, с K T 2 10% – точным. Очевидно, что при увеличении уровня шума точность полу чаемых оценок будет уменьшаться, поэтому будем рассматривать полу ченные оценки в зависимости от величины K n / s. Наборы параметров моделей, использованные в эксперименте, представлены в таблицах 3. и 3.21. Они отражают различные типы динамики, описываемые каждой из функцией.

Таблица 3. Значения параметров линейной модели Набор значений Параметр 1 2 3 50 100 700 C 10 10 -10 - C Таблица 3. Значения параметров параболической модели Набор значений Параметр 1 2 3 4 5 300 300 700 900 1500 C -15 5 -20 15 -5 C 0,7 0,7 0,3 -0,5 -0,3 -0, C На рисунке 3.33 а) и б) представлена зависимость коэффициентов вариации оценок параметров линейной модели от величины шума для каждого метода. Из графиков видно, что наименьший разброс оценок параметров обеспечивает итеративный подход.

Значения ряда можно использовать в качестве весов только при малом шуме.

а) б) Рис. 3.33. Зависимость коэффициентов вариации оценок а) C 0, б) C1 от K n / s Как показано на рисунке 3.34, наибольшую точность модели также обеспечивает третий метод. При использовании других подходов точ ность моделей существенно снижается с ростом зашумленности выбор ки. При малом шуме все методы дают приблизительно одинаковую точ ность.

50% 45% 40% 35% 30% 25% 20% 15% 10% 5% 0% 0% 5% 10% 15% 20% 25% 30% без весов вес Y[k] вес D[k] вес D[k] (итерат.) Рис. 3.34. Зависимость среднего K n/ s линейной модели от истинного K n / s На рисунке 3.35 представлена зависимость средней величины K T от K n / s. Все методы идентификации в среднем обеспечивают необходи мую точность прогноза.

0, 0, 0, 0, 0, 0, 0, 0, 0 0,05 0,1 0,15 0,2 0,25 0, без весов вес Y[k] вес D[k] вес D[k] (итерат.) истинное Рис. 3.35. Зависимость среднего K T 2 линейной модели от K n / s Зависимость коэффициентов вариации оценок параметров C0, C1, C2 параболической модели тренда от уровня шума представлены на рисунках 3.36-3.38. Методы, использующие оценки Dk в качестве ве сов, дают существенно меньший разброс оценок параметров.

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

10% 9% 8% 7% 6% 5% 4% 3% 2% 1% 0% 0% 5% 10% 15% 20% 25% 30% без весов вес Y[k] вес D[k] вес D[k] (итерат.) Рис. 3.36. Зависимость вариации С0 параболической модели от K n / s 100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0% 0% 5% 10% 15% 20% 25% 30% без весов вес Y[k] вес D[k] вес D[k] (итерат.) Рис. 3.37. Зависимость вариации С1 параболической модели от K n / s 45% 40% 35% 30% 25% 20% 15% 10% 5% 0% 0% 5% 10% 15% 20% 25% 30% без весов вес Y[k] вес D[k] вес D[k] (итерат.) Рис. 3.38. Зависимость вариации С2 параболической модели от K n / s На рисунке 3.39 представлена характеристика точности моделей для каждого из методов. Из графика видно, что три из четырех методов в среднем обеспечивают необходимую точность модели, даже если ко эффициенты вариации параметров велики.

Необходимую точность моделей обеспечивают все методы, кроме МВНК с Yk в качестве весов. Для других методов средний K n/ s меньше или равен истинному K n / s. Характеристика точности прогнозов пред ставлена на рисунке 3.40.

30% 25% 20% 15% 10% 5% 0% 0% 5% 10% 15% 20% 25% 30% без весов вес Y[k] вес D[k] вес D[k] (итерат.) Рис. 3.39. Зависимость среднего K n/ s параболической модели от K n / s 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0,05 0,1 0,15 0,2 0,25 0, без весов вес Y[k] вес D[k] вес D[k] (итерат.) истинное Рис. 3.40. Зависимость среднего K T 2 параболической модели от K n / s Таким образом, по совокупности критериев качества моделей можно сказать, что для идентификации моделей с мультипликативной структурой стохастической компоненты лучше использовать МВНК с оценками Dk в качестве весов. При этом нет необходимости организо вывать итеративную процедуру, поскольку она не дает существенного улучшения качества идентификации.

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

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

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

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

ГЛАВА 4. МОДЕЛИРОВАНИЕ И ПРОГНОЗИРОВАНИЕ РЯДОВ ДИНАМИКИ, ОПИСЫВАЕМЫХ ЭКСПОНЕНТАМИ И ИХ СОЧЕТАНИЯМИ С ГАРМОНИКОЙ 4.1. Использование модели в виде обобщенной экспоненциальной функции Экспонента в виде модели (1.18) Tk = A1 exp( 1k ) относится к числу едва ли не самых распространенных математических моделей тренда, входит в состав значительного числа других более сложных мо делей трендов.

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

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

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

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

Если 0, то имеем тренд с возрастающими значениями (уров нями), причём это возрастание не просто ускоренное, а с возрастающим ускорением и с возрастающими производными более высоких порядков.

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

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

Во второй главе была представлена и обобщенная экспоненциаль ная функция (2.5):

T (t ) = (1 exp ( t ) ), которая при очевидных переобозначениях примет вид:

T (t ) = A0 + A1 exp( 1t ).

Модель (2.5) может интерпретироваться как простейшая балансо вая модель Кейнса, включающая в себя доход, расход, потребление и инвестиции, а её квадрат описывает неоклассическую модель динамики фондовооруженности в зависимости от нормы амортизации, прироста во времени трудовых ресурсов [10, 14, 23, 25].

Воспроизводство при непрерывном темпе потребления, динамиче скую межотраслевую модель роста В. Леонтьева, экономический рост СЭС при различных траекториях потребления описывает уже сумма экспонент (обычно от двух до четырёх) с различными темпами роста и/или убывания [14].

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

Рассмотрим ряд динамики в виде суммы тренда (2.5) и стохасти ческой компоненты k, отвечающей принятым условиям Гаусса Маркова:

Yk = A0 + A1 exp( 1k ) + k. (4.1) В [78] предложен метод МНК-идентификации, основанный на конструировании ARMA-модели второго порядка, которую при k представим в виде:

Yk = 1 (Yk 1 Yk 2 ) + Yk 1 + k 1 ( k 1 k 2 ) + k 1 = (4.2) = 1 (Yk 1 Yk 2 ) + Yk 1 + k, где k = k 1 ( k 1 k 2 ) + k 1 – «новая» гомоскедастическая, центри рованная и автокоррелированная стохастическая компонента, 1 = exp( 1 ).

Заметим, что использование первых разностей в (4.2) представит ARMA-модель в абсолютных цепных приростах П ц.

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

Подстановка 1 в (4.1) позволит на втором этапе идентификации найти МНК-оценки А0 и A1 из решения соответствующей нормальной СЛАУ второго порядка.

При идентификации модели (1.18) повторяем первый этап иден тификации, а на втором этапе МНК-оценка A1 будет осуществляться из нормального линейного уравнения.

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

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

Динамические свойства производственных и экономических сис тем проявляются в сдвиге времени (временном лаге) от осуществления затрат (от принятия управленческого решения) до получения результа тов. Выделяют следующие типы лагов [12, 71]:

технологический лаг – от начала проектирования объекта до его ввода в действие;

строительный лаг – от начала работ по сооружению объекта до их завершения;

лаг отдачи – от реализации капиталовложений до освоения про ектной мощности;

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

психологический лаг – инерция в поведении людей;

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

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


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

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

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

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

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

Запаздывание проявляется и в процессах потребления – в отстава нии моментов изменения интенсивности потребления или спроса на те или иные товары и услуги в зависимости от момента изменения их цен и доходов потребителей. Это отставание приписывается, в частности, и «психологической инерции» потребителей.

Для моделирования всего многообразия запаздываний используют отдельно или в сочетании друг с другом звено идеальной задержки и динамические преобразователи первого и/или второго рода [23].

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

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

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

q exp ( p з ) W ( p) =, 1 + T1 p где q R, p – комплексная переменная.

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

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

Yk = q + q1 exp( k ) + k, где q1 = q exp( з ), =.

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

На первом этапе идентификации находим, на втором – q, q1, а затем по результатам обоих этапов рассчитывается и параметр запазды вания:

ln q1 + ln q з =.

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

Таким образом, необходимо в (4.2) анализировать не менее, чем 6-8 наблюдений.

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

Таблица 4. Истинные значения параметров модели Минимальное Максимальное Шаг Параметр значение значение варьирования 10 50 A 10 20 A 1 -0,1 0,2 0, Исследование проводилось на выборках объемом 24, 36 и 48 на блюдений. Для каждого сочетания параметров модели генерировалось 20 выборок. Коэффициент «шум/сигнал» изменялся от 0 до 0,35. В об щей сложности для исследования было сгенерировано 43 200 выборок.

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

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

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

наблюдений, присваиваемые средним значениям аргумента интервала сглаживания, что позволяет также уменьшить автокорреляцию k. Есте ственно, что при обоих приемах компенсации автокорреляции объем выборки уменьшается, что требует поиска компромисса между обосно ванием увеличения объема используемой статистической выборки и ус ловием стационарности модели на этой выборке. При сглаживании вы бираются МНК-оценки параметров на том шаге сглаживания, при кото ром получим лучшие показатели критерия точности моделирования и/или прогнозирования. Результаты исследования качества идентифи кации показаны на рисунке 4.1 а) и б), сведены в таблицу 4.2.

а) б) Рис. 4.1. Зависимость R 2 от K n / s при использовании прореживания (а) и сглаживания (б) выборки Видим, что качество идентификации остается достаточно высоким даже при мощности помехи равной 35%.

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

модель хорошо описывает тестовую выборку.

Таблица 4. Значение R 2 при различной величине шума и объеме исходной выборки R 2 при прореживании R 2 при сглаживании Истинный Kn/s R2 n = 24 n = 36 n = 48 n = 24 n = 36 n = 0,00 1,00000 1,00000 1,00000 1,00000 1,00000 1,00000 1, 0,05 0,95238 0,95614 0,95484 0,95402 0,95556 0,95474 0, 0,10 0,90909 0,91539 0,91311 0,91229 0,91510 0,91273 0, 0,15 0,86957 0,87779 0,87524 0,87338 0,87721 0,87529 0, 0,20 0,83333 0,84314 0,84047 0,83905 0,84247 0,83954 0, 0,25 0,80000 0,81157 0,80869 0,80591 0,81033 0,80703 0, 0,30 0,76923 0,78280 0,77652 0,77713 0,78013 0,77600 0, 0,35 0,74074 0,75586 0,74934 0,74680 0,75200 0,74850 0, При этом объем выборки незначительно влияет на качество иден тификации. И прореживание, и сглаживание выборки в данном случае дали приблизительно одинаковый результат.

Результаты исследования качества прогнозирования показаны на рисунке 4.2 а) и б). Из рисунка 4.2 видно, что средняя ошибка прогнози рования не превышает 15% даже при значительной зашумленности вы борки, а при величине шуме до 20% составляет менее 10%.

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

а) б) Рис. 4.2. Зависимость MAPE-оценки от K n / s при использовании прореживания (а) и сглаживания (б) выборки Для проведения такого исследования генерировались 10 000 выбо рок объемом 36 наблюдений, а глубину прогноза назначали в 12 наблю дений, мощность шума в 10%. Результаты численного эксперимента при ведены в таблицах 4.3 и 4.4. Средние значения оценок параметров оказа лись близки к их истинным значениям. Наибольшим разбросом обладают оценки параметра 1, коэффициент вариации которого составляет около 15%. Вариация R 2 составила менее 1%.

При этом ни на одной из 10 000 выборок коэффициент детермина ции не составил менее 0,89, а ошибка прогноза не превысила 20%.

Таблица 4. Оценки параметров модели и показателей качества идентификации при использовании прореживания выборки R Параметр MAPE A0 A Истинное значение 100 90 0,1 0,90836 0, Математическое 99,489 90,67 0,1004 0,91353 0, ожидание Среднеквадратическое 3,7621 4,8125 0,015204 0,0091789 0, отклонение Коэффициент 0,037814 0,053077 0,15143 0,010048 0, вариации Минимальное 74,419 75,101 0,049061 0,89932 0, значение Максимальное 109,1 108,5 0,17087 0,96199 0, значение Таблица 4. Оценки параметров модели и показателей качества идентификации при использовании сглаживании выборки объемом 36 наблюдений, глубиной прогноза 12 наблюдений и мощностью шума 10% R Параметр MAPE A0 A Истинное значение 100 90 0,1 0,9082 0, Математическое 99,464 90,626 0,10036 0,91322 0, ожидание Среднеквадратическое 3,8796 4,9008 0,015782 0,0093105 0, отклонение Коэффициент 0,039005 0,054077 0,15725 0,010195 0, вариации Минимальное 75,786 73,398 0,047603 0,89473 0, значение Максимальное 109,49 108,32 0,16709 0,95968 0, значение Близость графиков точностных характеристик при выборках в 24, 36 и 48 наблюдений побудила провести исследования и на меньших вы борках: в 12 и 6 наблюдений. При этом оказалось, что качество модели рования на выборках объемом 12 и 6 наблюдений практически такое же, что и показано на рисунках 4.1 и 4.2. Однако погрешность прогнозиро вания на выборке в 6 наблюдений с горизонтом прогноза в 2 наблюде ния практически в два раза больше, чем при 12 наблюдениях, что не по зволяет рекомендовать использование таких выборок при соотношении «шум/сигнал» более 20%.


Значительный интерес представляет количественная оценка влия ния шага прореживания на точность моделирования и прогнозирования, которая была выполнена на тестовых выборках объемом 36 наблюдений в широком диапазоне значений параметров модели. Шаг прореживания изменялся от 1 (без прореживания) до 12 (максимально допустимый шаг при заданном объеме выборки). Зависимость коэффициента детермина ции от коэффициента шум/сигнал для выборки приведена на рисунке 4.3. При шагах прореживания 4-5 качество идентификации становится приемлемым, а при 6-12 шагах прореживания результаты практически не различаются.

без прореживания 0.9 шаг = 0.8 шаг = 0.7 шаг = 0.6 шаг = шаг = 0. шаг = 0. шаг = 0. шаг = 0. шаг = 0. шаг = 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 шаг = Рис. 4.3. Зависимость R 2 от K n / s при различных шагах прореживания Зависимость коэффициента детерминации от шага прореживания при различных мощностях шума показана на рисунке 4.4, из которого видно, что качество идентификации не уменьшается при увеличении шага прореживания от одного шага до предельно допустимого.

0. 0. 0. 0. 01 2 3 4 5 6 7 8 9 10 11 шум 30% шум 10% шум 20% Рис. 4.4. Зависимость R2 от шага прореживания при K n / s = 0,1;

0, 2;

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

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

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

Yk = ( A0 + A1 exp( 1k ) ) µ k, где уровни µ k заданы в относительных единицах (в долях, в процентах).

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

Подобное утверждение приведено, например, и в [12].

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

Yk = ( A0 + A1 exp(1k ))(1 + k ) = A0 + A1 exp(1k ) +, (4.3) + ( A0 + A1 exp(1k )) k где стохастическая компонента k удовлетворяет принятым условиям Гаусса-Маркова.

В (4.3) новая стохастическая компонента ( A0 + A1 exp( 1k )) k пропорциональна уровням тренда, т.е. свойство гетероскедастичности сохранено. В соответствии с методикой, уже изложенной выше, сконст руируем при k 2 для (4.3) ARMA-модель второго порядка:

(4.4) Yk = Yk 1 + 1 (Yk 1 Yk 2 ) + k, где k = Yk ( k 1 k 2 + k 1 + k 2 ) + Yk 1( k k 2 + k + k 2 ) + – + 1 {Yk 1( k k 2 + k + k 2 ) Yk 2 ( k k 1 + k + k 1)} «новая» стохастическая компонента (центрирована, обладает автокор реляцией, гетероскедастична), а 1 = exp( 1 ).

Автокорреляцию k в (4.4) можно компенсировать приемом про реживания. Сложнее вопрос компенсации гетероскедастичности k.

Учтем, что условие Yk k в выражении для k позволяет пренеб речь малыми величинами второго порядка, содержащими произведения k 1 k 2, k k 2, k k 1.

Тогда для компенсации гетероскедастичности стохастической ком поненты k можно предложить прием почленного деления (4.4) на Yk 1, что дает:

Y Yk = 1 + 1 1 k 2 + k, (4.5) Yk 1 Yk Yk Y где k = (k + k 2 )(1 + 1 ) (k 1 + k 2 ) (k + k 1 ) 1 k 2.

Yk 1 Yk Заметим, что гетероскедастичность k определяется теперь лишь Y Y вариацией отношений k и k 2, т.е. зависит от цепных темпов роста Yk 1 Yk уровней ряда наблюдений Yk. В силу того, что для экспоненциального тренда они постоянны [4], то и для ряда Yk можно ожидать малой ва риации дисперсии k, т.е. ее практической гомоскедастичности.

Теперь можно к (4.5) применить МНК для нахождения оценки и рассчитать 1 с учетом обозначений в (4.4). Подставляя 1 в (4.4) и деля его на Yk, получим линейную регрессию с гомоскедастической стохастической компонентой A0 A1 exp(1 k ) ( A0 + A1 exp(1 k ) k 1= + + = Yk Yk Yk A0 A1 exp(1 k ) ( A0 + A1 exp(1 k ) k = + + ( A0 + A1 exp(1 k )(1 + k ) Yk Yk A0 A1 exp(1 k ) k + +, (1 + k ) Yk Yk из которой, решая соответствующую нормальную СЛАУ второго по рядка для минимизации относительной среднеквадратической невязки, определим МНК-оценки A0, A1.

Можно ожидать, что с ростом мощности стохастической компо ненты k (увеличения отношения мощностей помеха/полезный сигнал) более 10% компенсация гетероскедастичности путем деления на Yk (т.е. по сути применения взвешенного МНК с весовой функции ) Wk Yk будет недостаточной.

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

На первой итерации параметры модели будем оценивать «обыч ным» МНК. Полученные оценки параметров A0 1, A11,11 формируют взвешенный МНК с весовой функцией:

Wk =, ) ( A11 e 1 k A0 + который будет использоваться при следующем проходе для оценки МНК с весовой функцией:

W k z +1 = ) ( z A0 z + A1 z e 1 k до достижения требуемой точности.

На рисунке 4.5 показана зависимость величин второго коэффици ента Тейла от K n / s для различных весовых функций в идентификации обобщенной показательной функции с мультипликативной стохастиче ской компонентой, а также K T 2 для истинных значений параметров, за даваемых программой генерации. Видим, что все методы дают сущест венно менее точные прогнозы по сравнению с исходными моделями, хотя все они обеспечивают приемлемую точность прогнозирования при значениях K n / s до 15%.

18% 16% 14% 12% 10% 8% 6% 4% 2% 0% 0% 5% 10% 15% 20% 25% 30% без весов вес Y[k] вес D[k] вес D[k] (итерат.) истинное Рис. 4.5. Зависимость второго коэффициента Тейла от K n / s Таким образом, использование итерационной процедуры и ис пользование весовой функции Wk = 2 не дали в дан ) ( 1 1 1k A0 + A1 e ном случае существенного увеличения точности.

Целесообразность применения взвешенного МНК для нескольких рассмотренных моделей (суммы и произведения линейного тренда с мультипликативной колебательной компонентой, полиномиального тренда с мультипликативной стохастической компонентой) зависит от модели и от конкретного набора ее параметров.

Вывод аналитических выражений для оценки смещения и дисперсии оценки параметра 1 при реализации приемов компенсации автокорреля ции для модели с трендом (1.18) и для модели (4.2) оказался в принципе возможен.

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

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

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

Yk = C0 + A0 sh (0k + 0 ) + k, (4.6) Yk = C0 + A0 ch ( 0 k + 0 ) + k. (4.7) Гиперболические синус и косинус представляют собой соответст венно, полуразность и полусумму экспонент с одинаковым по модулю, но противоположным по знаку показателем степени. На рисунке 4. представлены графики моделей (4.6) и (4.7) в сравнении с обобщенной экспонентой.

Y (t ) C0 + A A exp ( 0t + 0 ) C0 + C C0 + A0 sh ( 0t + 0 ) C0 + A0 sh ( 0t + 0 ) 0 t Рис. 4.6. Графики обобщенной экспоненты, гиперболического синуса и гиперболического косинуса При С 0 0, 0 0 точка перегиба и минимум функций (4.6) и (4.7) соответственно смещаются в область первого координатного угла, а мо дели в области малых значений аргумента приобретают характер, суще ственно отличающийся от динамики обобщенной экспоненты.

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

Осуществим ARMA-идентификацию моделей (4.6) и (4.7), выпол нив предварительно замену переменных:

A0 sh (0k + 0 ) = A1 sh 0k + A2 ch 0k, (4.8) A0 ch (0k + 0 ) = A1 ch 0k + A2 sh 0k, (4.9) где A1 = A0 ch 0, A2 = A0 sh 0.

В обоих случаях приходим к одинаковой ARMA-модели:

Yk = (1 + ) (Yk 1 Yk 2 ) + Yk 3 + k, (4.10) где = 2ch (0 ), k = k k 3 (1 + ) ( k 1 k 2 ) – «новая» сто хастическая компонента.

Осуществим идентификацию параметра, используя МНК (ком поненту k считаем удовлетворяющей принятым условиям Гаусса Маркова) и приемы прореживания выборки для уменьшения автокорре лированности k. Из линейного нормального уравнения определим, затем с учетом обозначений в (4.10) можем рассчитать оценку 0 = arch, подставляя которую в (4.6) и (4.7), определим и МНК оценки параметров C0, A01, A02, линейно входящих в (4.8) и (4.9).

A ()() 2 2 = Arth.

= Далее очевиден расчет:, A0 A01 A A Было проведено исследование качества идентификации моделей (4.6) и (4.7) на тестовых выборках, в ходе которого было сгенерировано более двух сотен тысяч рядов динамики с различными значениями па раметров моделей (диапазон изменения параметров в 10-20 раз), раз личными соотношениями мощностей помехи и полезного сигнала в диапазоне от 0 до 30%.

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

Необходимая точность идентификации обеспечивается только при нахождении данных точек внутри выборки, на некотором удалении от ее краев (на 20-30% длины выборки). Кроме того, для модели (4.6) точ ность прогноза оказывается выше, когда точка перегиба находится в первой половине выборки, а для модели (4.7) – когда точка минимума находится во второй половине выборки.

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

Влияние других параметров на точность прогнозирования значи тельно меньше, но можно отметить, что второй коэффициент Тейла K T увеличивается при больших величинах C 0 и уменьшается при больших величинах A0. Покажем, что при выборе значений параметров, соответ ствующих описанным условиям, и при достаточном объеме выборки модели (4.6) и (4.7) позволяют получить высокую точность прогноза.

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

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

а) 1 б) 0, 0,95 0, 0,9 0, 0,85 0, 0,8 0, 0, 0, 0, 0% 5% 10% 15% 20% 25% 30% 0% 5% 10% 15% 20% 25% 30% 24 36 48 24 36 Рис. 4.7. Зависимость средних R 2 (а) и K T 2 (б) от K n / s для модели (4.6) а) 1 б) 0, 0, 0, 0, 0,85 0, 0, 0, 0, 0,7 0% 5% 10% 15% 20% 25% 30% 0% 5% 10% 15% 20% 25% 30% 24 36 24 36 Рис. 4.8. Зависимость средних значений R 2 и K T 2 (а) от K n / s (б) для модели (4.7) Из графиков видна высокая точность моделирования, но заметим, что при объеме выборки равном 24 значениям и менее существенно па дает точность прогноза. В целом, точность прогнозирования выше для модели (4.7).

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

4.2. Модели рядов динамики в виде квазиполиномов Квазиполиномы (определение А.Г. Гранберга [14]) описывают большое количество широко употребляемых в практике трендов вида:

m T (t ) = Ai t Bi exp( i t ) Cos(i t + i ), (4.11) i = где i – частоты, i – фазы гармонических компонент, Ai, Bi, i – пара метры, m – порядок квазиполинома.

Частными случаями (4.11) (при m = 1, i = i = Bi = 0 ) являются рассмотренная выше обобщенная экспоненциальная функция (2.5) T (t ) = A0 + A1 exp( 1t ) и экспоненциальная функция (1.18) (при A0 = 0 ) T (t ) = A1 exp( 1t ).

Применение в экономической практике нашли, например, сле дующие модели:

T (t ) = A1 exp( 1t ) + A2 exp( 2t ) ;

T (t ) = ( A1 + A2t ) exp( 1t ) ;

T (t ) = A1 Cos(1t + 1 ) + A2 Cos( 2t + 2 ) ;

T (t ) = A0 + A1 Cos(1t + 1 ) + A2 Cos( 2 t + 2 ) ;

T (t ) = A0 t + A1 Cos(1t + 1 ) + A2 Cos( 2t + 2 ) ;

T ( t ) = ( A0 + A0 t ) Cos (1t + 1 ) ;

T (t ) = A0 + A2t + A3Cos ( 3t + 3 ) ;

m T (t ) = ( Ai + Bit ) exp( it ), где обычно m i = и др.

Для каждой из приведенных моделей при постулировании адди тивной стохастической компоненты, удовлетворяющей принятым усло виям Гаусса-Маркова, можно сконструировать соответствующую ARMA-модель и получить МНК-оценки ее параметров [48].

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

Из предложенного выше комплекса моделей обращение к исполь зованию базиса Гребнера необходимо, например, для моделей тренда m T (t ) = A1cos (1t + 1 ) + A2 cos ( 2t + 2 ) и T (t ) = ( Ai + Bit )exp( it ).

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

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

T (t ) = A1 exp( 1t ) + A2 exp( 2t ), T (t ) = ( A1 + A2t ) exp( 1t ), T (t ) = A1 exp( 1t ) + A2 exp( 2t ) + ( A3 + A4t ) exp( 3t ) они могут представлять весьма разнообразные по форме графики трен дов, например, такие, как показаны на рисунках 4.9-4.11.

10 T(t) 10 T(t) 10 T(t) A A2 8 6 4 2 t t t 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 Рис. 4.9. Вид функции T (t ) = ( A1 + A2t ) exp( 1t ) при A1, A2 0, 0 1 100 T(t) 20 T(t) 5 t t 01 2 3 4 56 7 8 9 01 2 3 4 5 6 7 8 9 а) б) Рис. 4.10. Вид функции T (t ) = A1 exp( 1t ) + A2 exp( 2t ) при:

а) A1 = 9;

A2 = 11;

1 = 0,5;

2 = 0, 4 ;

б) A1 = 9;

A2 = 11;

1 = 0,1;

2 = 0, T(t) 100 T(t) 2 t t 01 2 3 4 56 7 8 9 01 2 3 4 5 6 7 8 9 б) а) Рис. 4.11. Вид функции T (t ) = A1 exp( 1t ) + A2 exp( 2t ) + ( A3 + A4t ) exp( 3t ) при: а) A1 = 9;

A2 = 11;

A3 = 15;

A4 = 1;

1 = 0,5;

2 = 0,4, 3 = 1;

б) A1 = 9;

A2 = 11;

A3 = 10;

A4 = 50;

1 = 0,1;

2 = 0, 2, 3 = 0, Известно, что модель:

T (t ) = A1 exp( 1t ) + A2 exp( 2t ) используется при 1 0 и 2 0 для описания связи между валовым продуктом и капитальными вложениями в экономической системе с учетом лага задержки при реализации капитальных вложений [28, 38].

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

Еще одним приложением рассматриваемых моделей может быть описание воспроизводства национального дохода, характеристика ди намики потребления и накопления в экономических системах: при этом используется модель (1.18), если потребление отсутствует, или модель T (t ) = A1 exp( 1t ) + A2 exp( 2t ) для двух случаев: если 2 = (т.е. потребление постоянно), или если 1 2 (т.е. потребление увели чивается с увеличением темпа прироста дохода).

Обратим также внимание на отличие модели T (t ) = A1 exp( 1t ) + A2 exp( 2t ) (имеющей более быструю динамику спа да уровней) на начальном участке аргумента (рис. 4.10 а)) от экспонен циальной функции. В предыдущем параграфе, напомним, отличие от экспоненциальной функции на малых значениях аргумента искали с помощью гиперболического синуса и гиперболического косинуса.

В модели T (t ) = ( A1 + A2t ) exp( 1t ) множитель ( A1 + A2t ) может, в зависимости от сочетания знаков и значений параметров A1 и A2, уси лить или уменьшить скорость изменения экспоненциальной функции (рис. 4.9), сделать функцию немонотонной (рис. 4.9, рис. 4.11).

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

Дело в том, что модели:

T (t ) = A1 exp( 1t ) + A2 exp( 2t ), (4.12) T (t ) = ( A1 + A2t ) exp( 3t ) (4.13) в сумме со стохастической компонентой дают одну и ту же ARMA модель (при k 2 ):

Yk = 1Yk 1 2Yk 2 + k, (4.14) где коэффициенты ARMA-модели определены как:

1 = exp( 1 ) + exp( 2 ), 2 = exp ( (1 + 2 ) ) для (4.12), но:

1 = 2 exp( 3 ), 2 = exp ( 2 3 ) для (4.13).

Тогда условием отнесения тренда ряда к модели (4.12) может быть выполнение следующей системы неравенств для коэффициентов ARMA-модели (4.12) (и для их оценок) на плоскости коэффициентов ARMA-моделей:

0 1 2, 2, 0 2 0, а условием принятия модели (4.13) являются соотношения:

0 1 2,.

2 = 0, Несколько более сложным, но и более тонким описанием, чем мо дель (4.1) временного лага, является модель запаздывания второго по рядка с трендом:

(4.15) Y (t ) = A0 + A1 exp( 1t ) + A2 exp( 2t ) + k, A0 2 exp(1 з ) A exp(2 з ) 1, A1 = 0 1, 1 =, 2 =, где A1 = 1 2 2 1 T1 T T1, T2 – постоянные времени.



Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |
 





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

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