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

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

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


Pages:     | 1 |   ...   | 4 | 5 || 7 |

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

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

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

1Yk = 11Yk 1 11Yk 1 + 1 k, (4.16) где 1 = ln( 1 ), 2 = ln( 2 ), 1 – первая разность, 1 k – первая раз ность уровней k, аналогичная по форме модели (4.16) для наблюдений определяемой переменной Yk, но состоящая из уровней k.

Из (4.16) легко получить МНК-оценки 1, 2, решая соответст вующую СЛАУ второго порядка, и, подставив их в (4.15), рассчитать МНК-оценки A0, A1, A2, а затем найти величину запаздывания по фор муле A1 (1 2 ) ln A0 з =.

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

Рассмотрим теперь идентификацию модели колебательной компо ненты ряда динамики из двух гармоник, например, получаемую после реализации метода детрендирования, который в методе параметриче ской итерационной декомпозиции из более сложной модели дает Y (t ) = A1 cos(1t + 1 ) + A2 cos(2t + 2 ) + k = A3 cos 1t + A4 sin 1t + (4.17) + A5 cos 2t + A6 sin 2t + k, где A3 = A1 cos 1, A4 = A1 sin 1, A5 = A2 cos 2, A6 = A2 sin 2.

Модели (4.17) (она также относится к квазиполиномам) соответст вует при k 4 следующая ARMA-модель:

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

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

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

С учетом обозначений в (4.17) также очевиден и МНК-расчет из нормальной СЛАУ четвертого порядка оставшихся для идентификации параметров A1, A2, 1, 2 модели в (4.17).

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

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

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

Для другой модели:

Y (t ) = A0 + A1cos (1t + 1 ) + A2 cos ( 2t + 2 ) + (t ) будет справедлива ARMA-модель (4.18), но уже в первых разностях. На втором этапе идентификации порядок СЛАУ увеличится до пятого.

Требуемая длительность реализации 15-20 наблюдений.

Модель:

Y (t ) = A0t + A1cos (1t + 1 ) + A2 cos ( 2t + 2 ) + (t ) уже была рассмотрена в третьей главе.

Теперь обратим внимание на то, что простейшая модель Y (t ) = A1cos (1t + 1 ) + (t ) (4.19) допускает помехозащищенную идентификацию на 6-8 наблюдениях, что укладывается в длительность одного периода или даже его части, в зависимости от величины шага дискретизации. Известные методы оцен ки периода (полупериода) модели (4.19) гармоники требуют обычно фиксации минимально двух последовательных переходов динамической траектории через ноль (или, что сложнее для регистрации в силу мень шего значения производной и, соответственно, худшего различия со седних уровней ряда) – через максимум.

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

4.3. Квазиполиномы, сочетающие экспоненту с гармонической компонентой В [48] рассмотрены шесть квазиполиномов, сочетающих в своей детерминированной части экспоненту с гармонической компонентой:

T (t ) = A1 exp( 1t ) Cos(1t + 1 ) + A2 exp( 2t ) + A3 exp( 3t ), T (t ) = A1 exp( 1t ) Cos(1t + 1 ) + A2 exp( 2t )( A3 + A4 k ), T (t ) = A1 exp( 1t ) Cos(1t + 1 ) + A2 exp( 2t )Cos ( 2t + 2 ), T (t ) = A1 exp( 1t ) + A2 exp( 2t ) + ( A3 + A4t ) exp( 3t ), T (t ) = A1 exp( 1t ) + A2 exp( 2t ) + ( A3 + A4t ) exp( 3t ), T (t ) = A1 exp( 1t ) + A2 exp( 2t ) + A3 exp( 3t ) + A4 exp( 4t ).

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

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

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

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

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

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

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

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

(4.20) Yk = A + B exp( k ) cos ( k + ) + k.

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

Если последняя – один год, то период цикла будет более 6 лет.

Пример функции (4.20) представлен на рисунке 4.12.

40 Dk k 0 10 20 30 40 50 60 70 80 90 Рис. 4.12. Вид детерминированной компоненты модели (4.20) Dk = A + B exp( k ) cos ( k + ) при A = 20;

B = 1;

= 0,03 ;

= 0, 27;

= 0, Можно показать, модели (4.20) соответствует при k 3 ARMA модель третьего порядка:

(4.21) Yk = 1Yk 1 12Yk 2 + 2Yk 3 + k, где 1 = 2 exp( )Cos + 1, 2 = exp(2 ), k – «новая» стохастическая компонента, образованная из уровней k так же, как и Yk в (4.21).

Использование базиса Гребнера, приема прореживания выборки для компенсации автокорреляции k позволят определить МНК-оценки 1, 2, а затем, с учетом обозначений в (4.21), найти и,.

Подставляя найденные оценки, в (4.20), из нормальной СЛАУ третьего порядка, соответствующей МНК-идентификации моде ли:

Yk = A + B Cos exp( k )Cos k B Sin exp( k )Sin k + k найдем A, B,.

Еще одной широко используемой в практике экономических ис следований моделью является зависимость спроса Yk от дохода и цены потребительского блага с учетом инерции потребителей [23]:

Yk = D exp( 1k ) cos (1k + 1 ) + k. (4.22) Цена предполагается в (4.22) неизменной, и исследуется инерция потребителей, проявляющаяся в запаздывании реакции на изменение дохода.

Вид функции (4.22) представлен на рисунке 4.13.

16 Dk k 0 10 20 30 40 50 60 70 80 90 - - Рис. 4.13. Вид детерминированной компоненты модели (4.22) при D = 18;

1 = 0,18;

1 = 5,9;

1 = 1, С учетом уже представленных выше результатов для модели (4.20) идентификация (4.22) очевидна, она даже проще.

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

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

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

При этом, как правило, совмещаются нормативный и дескриптив ный подходы [14], осуществляется эконометрическое моделирование для оперативной идентификации и коррекции макроэкономических мо делей.

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

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

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

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

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

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

Yk = C1 exp(1k ) + C 2 exp( 2 k ) cos( k + ) + k. (4.23) Выражению (4.23), как можно показать, соответствует при k ARMA-модель авторегрессии третьего порядка наблюдений:

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

Дальнейший алгоритм идентификации очевиден: использование базиса Гребнера для получения МНК-оценок 1, 2, 3, расчет МНК оценок 1, 2,, затем расчет из нормальной СЛАУ третьего порядка МНК-оценок A1, A2,.

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

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

Можно показать, что модели (4.25) будет соответствовать при k 4 следующая ARMA-модель шестого порядка:

Yk = Yk 1 Yk 2 + Yk 3 + µ1Yk 1 µ 2Yk 2 + µ3Yk, (4.26) µ 4Yk 4 + µ5Yk 5 µ6Yk 6 + k где µ1 = 1 + 2 + 12, µ2 = 2 + 312 + 12 + 122, µ3 = 1 + 12 + 13 + + 312 + 122, µ4 = 1 + 12 + 13 + 3122 + 132, µ5 = 12 + 2 + + 132, µ6 = 13, 1 = exp(1), 2 = 2cos2, k – «новая» стохастическая компонента, удовлетворяющая принятым ус ловиям Гаусса-Маркова, имеющая тот же вид, что и Yk, но в уровнях k.

Из анализа таблицы 3.1 можно сделать вывод о возможности приме нения базиса Гребнера для МНК-оценок 1, 2 модели (4.26) и, соответст венно, 1,, а также последующего определения МНК-оценок A0, A1, A2, 2 из нормальной СЛАУ шестого порядка. Чтобы уменьшить порядок используемой для идентификации нормальной СЛАУ с шестого до четвертого, тем самым уменьшив вычислительные погрешности, можно рекомендовать прием параметрической итерационной декомпозиции для (4.25) на тренд и колебательную компоненту и применение ARMA моделирования для каждой из них.

Количественный пример моделирования ряда динамики с использо ванием тренда в виде обобщенной экспоненциальной функции и колеба тельной компоненты проведем на примере данных официальной статисти ки населения г.о. Самара. Численность населения в интервале с 1970 г. по 1989 г. росла (рис. 4.14), а во втором (с 1990 по 2009 гг.) она демонстрирует тенденцию к уменьшению (рис. 4.15) [66]. Выбор модели производился из набора функций, отвечающих визуальным тенденциям динамического ря да. Некоторые из них выходили за рамки рассмотренных до этого парагра фа монографии моделей. Для первого диапазона с 1970 г. по 1989 г. выбор модели осуществлялся из следующего семейства пяти функций:

1) Yk = C + A0e0k + k ;

( ) 2) Yk = C + A0e0k (1 + k ) ;

= ( C + A e ) (1 + A sin ( k + ) ) + 0k 3) Yk k;

4) Yk = C + A0e 0k + A1e1k sin (1k + 1 ) + k ;

5) Yk = C + ( B0 + B1k ) e 0k + k.

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

1) Yk = 1291,5 268,0e0,110k + k, R 2 = 0,994 ;

2) Yk = (1296,8 271,5e0,105k ) (1 + k ), R 2 = 0,994 ;

3) Yk = (1307,0 276,7e0,097 k ) (1 + 0,0063sin ( 0, 473k 3,13) ) + k, R 2 = 0,998 ;

4) Yk = 1292,9 265,0e0,108k + 12,5e0,067 k sin ( 0,495k + 2,83) + k, R 2 = 0,999 ;

5) Yk = 1260,1 + ( 224,5 47,8k ) e0,259k + k, R 2 = 0,999.

Таблица 4. Результаты моделирования и прогнозирования населения г.о. Самара с 1970 по 1989 гг.

Модель Годы Исходный 1 2 3 4 1970 1035,3 1023,5 1025,3 1030,2 1031,7 1035, 1971 1048,4 1051,4 1052,3 1052,7 1052,8 1049, 1972 1071,7 1076,4 1076,6 1073,4 1072,5 1069, 1973 1092 1098,8 1098,5 1093,2 1091,8 1974 1112,6 1118,9 1118,2 1112,4 1111,4 1112, 1975 1133,5 1136,9 1136 1131,5 1131,2 1133, 1976 1151,3 1153 1152 1150,1 1150,5 1977 1166 1167,4 1166,4 1167,8 1168,7 1168, 1978 1178,8 1180,3 1179,4 1183,9 1184,9 1183, 1979 1203,3 1191,9 1191,1 1198 1198,5 1196, 1980 1211,9 1202,3 1201,6 1209,4 1209,4 1207, 1981 1217,5 1211,6 1211,1 1218,2 1217,7 1216, 1982 1220,2 1219,9 1219,6 1224,7 1223,8 1224, 1983 1229,8 1227,4 1227,3 1229,3 1228,7 1230, 1984 1233,6 1234,1 1234,2 1232,9 1232,9 1236, 1985 1236,7 1240 1240,4 1236,4 1237,1 1240, 1986 1241,9 1245,4 1246 1240,6 1241,7 1244, 1987 1247,3 1250,2 1251,1 1245,8 1246,8 1247, 1988 1253,2 1254,5 1255,6 1252,3 1252,3 1249, 1989 1257,3 1258,4 1259,7 1259,7 1257,8 1251, R 0,99416 0,99406 0,99847 0,99866 0, MAPE 0,15% 0,18% 0,39% 0,38% 0,56% Прогноз на оценка 1 шаг KT2 0,00103 0,0013 0,00276 0,00264 0, MAPE 0,22% 0,25% 0,37% 0,47% 0,56% Прогноз на оценка 2 шага KT2 0,00152 0,00178 0,00292 0,0035 0, MAPE 0,35% 0,44% 0,44% 0,73% 0,48% Прогноз на оценка 3 шага KT2 0,00247 0,00311 0,00355 0,00524 0, MAPE 0,65% 0,70% 0,51% 1,00% 0,36% Прогноз на оценка 4 шага KT2 0,00456 0,00496 0,00421 0,00725 0, MAPE 0,99% 1,07% 0,82% 0,63% 0,28% Прогноз на оценка 5 шагов KT2 0,007 0,00755 0,00658 0,00446 0, Объем используемой выборки в данном моделировании был равен двадцати наблюдениям. Если исходить из точности моделирования, то, казалось бы, следует выбрать четвертую модель.

Однако лучший прогноз на три шага (года) дала первая модель, а лучший прогноз на четвертый и пятый шаги дала пятая модель. Исходя из задачи обеспечения в исследовании высоких показателей точности моделирования и прогнозирования, можно воспользоваться комплекс ным критерием для различных i-моделей MFi = max ( Ri2 MAPEi ), при чем в MF можно задавать «веса» (важность) точности моделирования и прогнозирования. Обычно принимают, что сумма весов равна единице.

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

Рис. 4.14. Исходный ряд динамики населения г.о. Самара в 1970-1989 гг. и его моделирование Погрешность прогнозирования (характеризуемая MAPE-оценкой) менее 1,07%. Точность моделирования практически одинакова: более 99% статистических данных объясняются каждой из рассматриваемых моделей.

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

Тенденции уменьшения населения г.о. Самара в диапазоне с 1990 г.

по 2009 г. моделировались следующими шестью функциями:

1) Yk = C + A0e0k + k ;

2) Yk = C + A0 e 0 k + A1 sin (1k + 1 ) + k ;

3) Yk = C + A0 e 0 k + A1e1k sin (1k + 1 ) + k ;

4) Yk = C + ( B0 + B1k ) e0k + k ;

5) Yk = C + ( B0 + B1k ) e 0k + A1 sin (1k + 1 ) + k ;

6) Yk = C + ( B0 + B1k ) e 0 k + A1 sin (1k + 1 ) + A2 sin ( 2 k + 2 ) + k.

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

0,060 k + k, R 2 = 0,971 ;

1) Yk = 1081,2 + 153,3e 2) Yk = 1599, 2 381,3e 0,012 k + 10,1sin ( 0,342 k + 0,88 ) + k, R 2 = 0,983 ;

3) Yk = 1090, 4 + 143, 42e 0,066 k + 11,1e 0,092 k sin ( 0, 440 k 0,63 ) + k, R 2 = 0,983 ;

4) Yk = 1119,9 + (107,0 + 19,7 k ) e 0,193k + k, R 2 = 0,982 ;

5) Yk = 1128,6 + ( 94,1 + 27,8k ) e 0,240 k + 3, 4sin (1,596k + 0, 20 ) + k, R 2 = 0,988 ;

6) Yk = 1041,6 + (191,5 0,043k ) e 0,043k + 4,6sin ( 0,508k 1,00 ) + +3,5sin (1, 272 k 1,98 ) + k, R 2 = 0,987.

Таблица 4. Результаты моделирования и прогнозировании населения г.о. Самара в диапазоне с 1990 по 2014 гг.

Модель Исходный ряд 1 2 3 4 5 1990 1225,5 1234,5 1225,7 1227,2 1226,8 1223,3 1991 1225,1 1225,6 1223 1222,7 1224,3 1227,8 1220, 1992 1220,2 1217,2 1219,1 1218,3 1219,4 1220,4 1219, 1993 1212,9 1209,3 1213,9 1213,4 1213 1211,7 1215, 1994 1203,8 1201,8 1207,7 1207,4 1205,8 1208,3 1206, 1995 1200,7 1194,8 1200,4 1200,4 1198,2 1202,1 1197, 1996 1191,7 1188,2 1192,5 1192,6 1190,7 1189,3 1191, 1997 1187,2 1182 1184,4 1184,4 1183,4 1179,3 1187, 1998 1177,8 1176,1 1176,3 1176,2 1176,5 1176,4 1180, 1999 1170,8 1170,6 1168,6 1168,6 1170 1171,5 1169, 2000 1158,9 1165,4 1161,8 1161,8 1164,1 1160,9 1158, 2001 1146,4 1160,5 1156 1156,1 1158,7 1154,1 2002 1157,9 1155,9 1151,3 1151,6 1153,8 1154,3 1153, 2003 1155,4 1151,5 1147,6 1148 1149,5 1151,7 1150, 2004 1144,2 1147,4 1144,8 1145,1 1145,6 1143,6 1144, 2005 1133,4 1143,6 1142,7 1142,7 1142,2 1139,7 1139, 2006 1143,3 1140 1140,8 1140,5 1139,2 1142,1 1139, 2007 1139 1136,5 1138,7 1138,2 1136,5 1141 1140, 2008 1135,4 1133,3 1136,2 1135,8 1134,2 1134,4 1136, 2009 1134,7 1130,3 1132,8 1133,1 1132,2 1132,4 1128, 2010 - 1127,4 1128,3 1130,3 1130,5 1136,2 1119, 2011 - 1124,7 1122,5 1127,3 1129 1135,5 1116, 2012 - 1122,2 1115,6 1124,4 1127,7 1129,8 1114, 2013 - 1119,8 1107,6 1121,6 1126,5 1129,1 1110, 2014 - 1117,6 1098,8 1119,1 1125,6 1133,4 1103, R2 0,97053 0,98327 0,98332 0,98179 0,98779 0, MAPE 0,63% 0,71% 0,52% 0,66% 0,27% 0,27% Прогноз на 1 шаг KT2 0,0045 0,00189 0,00503 0,00188 0,00369 0, MAPE 0,66% 0,99% 0,80% 0,52% 0,39% 0,24% Прогноз на 2 шага KT2 0,00478 0,00703 0,00598 0,00177 0,00371 0, MAPE 1,07% 1,02% 0,60% 3,41% 0,28% 0,18% Прогноз на 3 шага KT2 0,00769 0,00753 0,00491 0,00135 0,02621 0, MAPE 1,42% 1,40% 1,84% 0,70% 1,64% 0,36% Прогноз на 4 шага KT2 0,01025 0,01006 0,01459 0,00266 0,00553 0, MAPE 1,02% 0,93% 1,45% 0,45% 1,56% 0,27% Прогноз на 5 шагов KT2 0,00772 0,00712 0,0129 0,00275 0,00373 0, Объем используемой выборки в этом случае равен 19-ти наблюде ниям. Точность моделирования оказалась высокой для всех моделей:

коэффициент детерминации моделей более 0,97. Если и в этом случае исходить только из точности моделирования, то следует предпочесть пятую модель в виде уменьшающейся логистической функции Рамсея и колебательной гармонической компоненты. Обратим внимание, однако, на то, что четвертая модель имеет практически постоянную лучшую точность прогнозирования. Эта «устойчивость» может быть основанием для предпочтения модели в виде произведения обобщенной экспонен циальной функции и линейной функции.

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

Для сравнения отметим, что использование стандартной програм мы Excel позволяет идентифицировать ряд динамики известным спосо бом (при предположении мультипликативной структуры стохастиче ской компоненты µ k с логнормальным законом распределения, выпол нении операции логарифмирования, применении затем МНК для полу чаемой линейной регрессии), т.е. экспоненциальную модель Yk = A1 exp( 1k ) µ k.

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

R 2 = 0,78.

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

Проведенное моделирование показало, что при краткосрочном прогнозе динамика численности жителей г.о. Самара отрицательна.

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

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

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

Моделирование в рядах динамики колебательной компоненты, пропорциональной обобщенному экспоненциальному тренду ( ) Yk = C + A0e 0k (1 + A sin ( k + ) ) + k, дало следующий результат:

( ) Yk = 1097,8 + 121,8e 0,129 (1 + 0, 003sin ( 0,886 k 1,57 ) ) + k.

При этом коэффициент детерминации оказался равным 0,99921, а MAPE-оценка прогноза на один шаг вперед - 0,05%, второй коэффици ент Тейла 0,03%. MAPE-оценка прогноза на два шага вперед оказалась равной 0,06%, второй коэффициент Тейла - 0,04%. MAPE-оценка про гноза на три шага вперед равна 0,05%, а второй коэффициент Тейла – 0,03%. Для четырех шагов вперед получим оценку точности прогнози рования 0,04% и 0,03%, на пять – 0,038% и 0,031% соответственно.

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

Yk = C + A0 e 0k + A1e1k sin (1k + 1 ) + k дало следующий результат:

Yk = 1130,5 + 129,7 e 0,326 k + 37,0e 0,348 k sin ( 0,665k 1,63 ) + k с коэффициентом детерминации 0,99961.

MAPE-оценка прогноза на один шаг вперед оказалась равной 0,04%, а второй коэффициент Тейла – 0,038%;

MAPE-оценка прогноза на два шага вперед равна 0,034%;

второй коэффициент Тейла – 0,025%.

MAPE-оценка прогноза на три шага вперед - 0,032%, а второй коэффи циент Тейла - 0,023%, на четыре шага вперед получим 0,037% и 0,028%, на пять – 0,034% и 0,026%.

На рисунке 4.16 представлены абсолютные ошибки трех предло женных моделей для исходного ряда.

Ошибка 1 Ошибка 2 Ошибка Рис. 4.16. Абсолютные ошибки моделирования динамики населения г.о. Самара На рисунке 4.17 представлены графики MAPE-ошибок прогноза численности населения в абсолютных единицах, а на рисунке 4.18 – графики ошибки прогноза по второму коэффициенту Тейла в абсолют ных единицах.

Рис. 4.17. Графики MAPE-ошибок прогноза численности населения Рис. 4.18. Графики ошибки прогноза численности населения по второму коэффициенту Тейла Очевидны преимущества моделирования и прогнозирования численности населения г.о. Самара моделями с учетом эволюциони рующей колебательной компоненты в анализируемом ряде динамики (второй и третьей моделями).

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

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

ГЛАВА 5. МОДЕЛИРОВАНИЕ И ПРОГНОЗИРОВАНИЕ ЖИЗНЕННОГО ЦИКЛА ПРОДУКТА И СОЦИАЛЬНОЙ ДИНАМИКИ 5.1. Выбор метода идентификации логистической динамики моделью Верхулста Из почти пятидесяти приведенных во второй главе монографии известных моделей логистической динамики показателя Yk наиболее распространенной до настоящего времени можно считать симметрич ную модель Верхулста (ее называют часто классической):

A + k.

Yk = (5.1) 1 + A1e k Для стохастической компоненты k в (5.1) будем считать спра ведливыми принятые условия Гаусса-Маркова.

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

До настоящего времени не решена задача сравнения точности мо делирования и прогнозирования методов идентификации модели (5.1) в динамическом диапазоне сочетаний параметров модели, при различных соотношениях мощности помехи k и полезного сигнала.

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

A Yk =, 1 + A1e k k A A 1 = A1e k k ln 0 1 = ln A1 k + ln k.

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

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

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

1 1 1 = ln, ln m 1Yk 2 Yk Yk Yk 2 где m – длина отрезка;

1 Yk Yk 1, A0 = m 1Yk 1 + Yk Yk 1 1 ) ( 1 1 e A Yk Yk A1 = 1.

1 ) ( 1 me + Yk Yk 1 Методы Фишера и Готеллинга [69] предполагают использование Y dY = Y дифференциального уравнения логистической функции:.

dk A В обоих методах используются приближения расчета производной.

Метод Фишера использует приближенное вычисление темпов прироста:

N dY 1 1 Yk +1 1 Yk +,, B = arg min ln + BYk 1, A0 = ;

ln dk Yk 2 Yk 1,B k =1 2 Yk 1 B где N – объем выборки.

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

dY Y = Yk Yk 1.

dk k В частности, для метода Готеллинга будем иметь:

Yk Yk = Yk + k, где k – «новая» (не равная k в (5.1)) стохастиче A ская компонента. Оценки параметров, B находятся с помощью МНК:

N ( ), 2, B = arg min Yk Yk 1 Yk 1 + BYk 12 =.

A B,B k = Yk +1 Yk Метод Юла использует регрессию темпов прироста на Yk Yk +1 :

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

Нелинейно входящие в уравнение параметры заменяют линейны e 1 = Q, = P и применяют МНК:

ми e A N 1 Y Y Q, P = arg min k +1 k Q + PYk +1.

Yk Q, P k =1 Метод Родса [69] использует взятие разности между соседними обратными значениями ряда, а идентификацию осуществляет относи (1 e ) и e тельно параметров :

A (1 e ) + e 1 + k, = Yk +1 A0 Yk где k – «новая» стохастическая компонента.

Метод Нейра основан на идентификации параметров регрессии разности соседних обратных значений ряда на их сумму, а решение e 1 e и :

осуществляет относительно параметров e + 1 A + e 1 e 1 2 e 1 1 1 + k, = + + 1 A0 e + 1 Yk +1 Yk Yk +1 Yk e где k – «новая» стохастическая компонента.

Все перечисленные методы предполагают последовательное вы числение оценок параметров модели: на первом этапе производится расчет оценок параметров, A0, а на втором этапе находят оценки па раметра A1.

В работах [64, 69] предлагается оценивать значение параметра A A A 1 = A1e k, ln A1 = k + ln 0 1 ), а затем ис методом моментов ( Yk Yk кать средние значения:

A N N N ln Y0 ln A1 k N ( N + 1) N k, что, с учетом k = k = k =1 k =, дает = + N N N k = оценку:

1 N ( N + 1) N A + ln 0 1.

= exp A N 2 k =1 Yk Однако при наличии в выборке значений, превышающих найден ную оценку уровня насыщения логистической кривой A0, метод момен тов становится уже неработоспособным, поскольку возникает отрица тельное число под знаком логарифма.

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

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

1 e k N A1 = arg min A.

Yk A0 A k = A Однако нужно учитывать, что в этом случае получим гетероскеда стическую стохастическую компоненту в модели:

( ) A0 + k 1 + A1e k 1 + A1e k A + k = = Yk =,.

( ) 1 + A1e k 1 + A1e k Yk A0 + k 1 + A1e k ( ) k 1 + A1e k 1 1 + A1e k + k, k = =.

Yk A0 A0Yk Для компенсации гетероскедастичности и, соответственно, уменьшения неэффективности оценок параметров модели можно рас смотреть возможности использования метода ARMA-моделирования.

Для его реализации осуществим замену переменных модели 1 A, A1 = 10, сконструируем разностную схему (при k 2 ) детер A0 = A00 A минированной части модели:

Dk = ( Dk 1 Dk 2 ) + Dk 1, где = e, Dk – детерминированная часть модели ( Dk = A00 + A10e k ).

Первый метод реализации этого метода (назовем его ARMA I) учитывает соотношение Dk = k и приводит к следующей модели Yk авторегрессии:

= Gk = ( Gk 1 Gk 2 ) + Gk 1 + k, Yk где k = k k 1 + ( k 2 k 1 ) - гетероскедастическая стохастическая компонента.

Оценка параметра (или параметра ) находится с помощью взвешенного МНК для компенсации гетероскедастичности.

В качестве оценок весов wk можно, во-первых, использовать сами уровни ряда динамики Gk, если величины помехи малы по сравнению с уровнями ряда. Или, во-вторых, реализовывать взвешенный МНК в два этапа: на первом модель идентифицируется с помощью МНК, а на вто ром этапе полученные оценки Dk используются в качестве весов. Тре тий способ задания весов включает использование какого-либо метода непараметрического сглаживания: вначале оцениваются значения k = Yk Yk, где Yk – полученные сглаженные значения ряда. Затем по соседним m точкам строится ряд оценок дисперсии k. Заметим, в каче стве недостатка, что при этом теряются ( m 1) значения ряда:

1 m1 1 m1 k i m k i.

Sk = m 1 i =0 i = После этого полученные оценки дисперсии используются в каче стве весов wk :

n 2( k G ( Gk 1 Gk 2 ) Gk 1 ).

= arg min k =2 wk Второй метод (его назовем ARMA II) учитывает соотношение Dk = :

Yk k (Yk 1 Yk )Yk 2 (Yk 2 Yk 1 )Yk = Yk 1 k 2 + Yk 2 k k 1 k 2 k ( k 1 k 2 ), а затем через МНК позволяет найти:

(Y Y ) Y n = arg min 2 (Yk 1 Yk ) k 2 k 1 k.

Yk k =2 wk Второй этап идентификации для обоих ARMA-методов одинаков:

находятся оценки параметров A00, A10 с помощью взвешенного МНК, а затем вычисляются МНК-оценки параметров A0, A1 :

).

( N 1 = arg min 2 Dk A00 A10e k A00, A A00, A10 k =1 wk Отметим, что модель (5.1) также может быть также идентифици рована с помощью численного решения МНК, например, методом Ле венберга-Марквардта [9].

Постановка задачи идентификации в этом случае выглядит так:

найти такое значение вектора параметров модели, которое бы опреде лило локальный минимум функции ошибки E :

N E = (Yk f ( ) ), k = где f ( ) – регрессионная модель.

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

( )) ( J T (Y f ( ) ), = J T J + diag J T J где J – якобиан функции f ( ) в точке, 0 – параметр регуляриза ции, назначаемый на каждой итерации алгоритма.

Проведем теперь количественное исследование точности модели рования и прогнозирования для модели Верхулста (5.1) девятью по следними методами: Фишера, Готеллинга, Юла, Родса, Нейра, трех сумм, ARMA I, ARMA II и методом Левенберга-Марквардта.

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

Соотношение мощностей помехи и полезного сигнала K n / s варьи ровалось от 0 до 0,3 с шагом 0,05.

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

Была рассмотрена и точность оценивания отдельных параметров модели: смещение, среднеквадратическое отклонение, коэффициент ва риации оценок параметров.

Для этого было сгенерировано 10 000 выборок объемом 36 наблю дений, глубина (горизонт) прогноза назначалась в 12 наблюдений, а ко эффициент «шум/сигнал» был принят в 10%.

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

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

Применение взвешенного МНК в методах ARMA I и ARMA II не улучшило существенно качество идентификации.

Результаты идентификации для методов Родса, Юла и Нейра не приводятся, поскольку они не дали удовлетворительного результата да же в случае добавления шума, мощность которого не превышает 5% от мощности полезного сигнала: полученные с помощью данных методов модели объясняют менее 10% исходных данных.

Таблица 5. Значения параметров модели, использованные при генерации тестовых выборок Параметр Минимальное значение Максимальное значение Шаг 50 100 A 50 200 A 0,1 0,8 0, Таблица 5. Результаты исследования качества оценивания отдельных параметров модели R Параметр MAPE A0 A Метод Готеллинга Истинное значение 50 50 0,3 0,9082 0, Математическое ожи 47,1300 134,4400 0,4654 0,6717 0, дание Среднеквадратическое 1,3398 889,8900 0,1125 0,1793 0, отклонение Коэффициент вариа 0,0284 6,6193 0,2417 0,2670 0, ции Метод трех сумм Истинное значение 50 50 0,3 0,90842 0, Математическое ожи 49,4240 267,5400 0,3228 0,9020 0, дание Среднеквадратическое 1,7138 4972,0000 0,0788 0,0159 0, отклонение Коэффициент вариа 0,0347 18,5840 0,2441 0,0177 0, ции Метод ARMA II Истинное значение 50 50 0,3 0,9082 0, Математическое ожи 56,6130 109,5400 0,2094 0,5635 0, дание Среднеквадратическое 576,0500 1356,3000 0,2326 0,3114 0, отклонение Коэффициент вариа 10,1750 12,3830 1,1106 0,5527 1, ции Метод Левенберга-Марквардта Истинное значение 50 50 0,3 0,908 0, Математическое ожи 50,4380 47,7970 0,2831 0,9146 0, дание Среднеквадратическое 1,4700 42,0710 0,0461 0,0097 0, отклонение Коэффициент вариа 0,0291 0,8802 0,1630 0,0106 0, ции Анализ полученных результатов показывает, что только метод Ле венберга-Марквардта в рассмотренном достаточно широком динамиче ском диапазоне значений параметров до соотношения мощности помехи в 30% от мощности полезного сигнала дает приемлемый результат по точности моделирования и прогнозирования.

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

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

В качестве примера рассмотрен жизненный цикл продуктов ком пании Electronic Arts (EA), которая является разработчиком компьютер ных видеоигр, в частности, популярной серии «Need For Speed». Ис пользована статистика индекса поисковых запросов (SVI, Search Volume Index) сервиса Google [94].

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

A + C0 + C1k + k.

Yk = 1 + A1e k Однако для идентификации этой модели приведенные выше мето ды напрямую не могут быть применены. Этот факт потребовал обраще ния к методу параметрической итерационной декомпозиции [59].

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

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

Таблица 5. Результаты моделирования и прогнозирования жизненного цикла продукции компании EA NFS Pro Street NFS Undercover 1,04 1, Yk = Yk = Модель 1 + 83437e 1,6 k 1 + 14534e 1,9 k 0,06 + 0,08k + k 0,11 + 0,13k + k R2 0,99 0, MAPE 2,71% 1,68% а) б) Рис. 5.1. Данные статистики поисковых запросов на продукцию компании EA и результаты моделирования и прогнозирования ее жизненного цикла:

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

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

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

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

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

Для сравнения альтернативных вариантов с различной доходно стью принято использовать относительный показатель риска: отноше ние среднеквадратического отклонения показателя экономической эф фективности инвестиций к ожидаемой средней величины показателя (математического ожидания) – коэффициент вариации (CV ) [1, 2]:

, (5.2) CV = Pcp где Pcp – среднее ожидаемое значение показателя эффективности проек та (математическое ожидание).

Модели уровней Yk рядов наблюдений добычи нефти описывают, например, моделями:

YkStepenExp = A1e 1k ( k ) (5.3) + k,, A1 ( k ) + k, YkDrobnRatz = (5.4) A2 + A3 ( k ) где 1, 2, A1, A2, A3 – параметры модели (в общем случае нецелочислен ные).

Напомним, что во второй главе монографии функции (5.3) и (5.4) использовались в качестве моделей ЖЦП.

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

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

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

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

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

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

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

Yk n stepen = Tn ( A1, A2,..., An +1, k ) e 1k + An+ 2 + k, (5.5) где Tn ( A1, A2,..., An +1, k ) = A1 ( k )n + A2 ( k )n 1 +... + An k + An +1 – алгебраи ческий многочлен степени n, k – стохастическая компонента.

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

Идентификацию будем как всегда реализовывать в два этапа: на первом этапе сконструируем ARMA-модель и применим МНК для оценки нелинейного параметра модели 1, а на втором этапе, с исполь зованием полученной МНК-оценки 1, осуществим МНК-оценку дру гих оставшихся линейных параметров Ai, i = 1, 2,..., n + 2.

Продемонстрируем предложенный метод моделирования на при мере шести пластов Турнейского яруса месторождения Самарской об ласти. Для построения модели (5.5) рассмотрим возможность использо вания многочленов Tn от первой до пятой степени.

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

Yk1stepen = (135,187 k 250,367 ) e 0,22241k + 134,52, R12stepen = 0,70498 ;

( ) Yk2stepen = 28,378(k )2 20,594k 113,777 e0,25106k + 101,0136, R2 stepen = 0,90932 ;

( ) Yk3stepen = 20,13168(k )3 53,633(k )2 1,033k 123,093 e0,43127k + + 125,643, R3 stepen = 0,94692 ;

( ) Yk4stepen = 5,942(k )4 20,399(k )3 + 42,933(k )2 94,725k 125, e0,53691ti + 129,402, R4 stepen = 0,94315 ;

( Yk5stepen = 15,633(k )5 165,437(k ) 4 + 682,619(k )3 1144,757(k ) 2 + + 487,313k 134,560 ) e0,78122 k + 138,122.

R5 stepen = 0,92206.

На рисунке 5.2 представлены графические результаты идентифи кации истории добычи нефти рассмотренным семейством моделей (5.5).

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

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

Yi Yi Y 5 stepen i 350 Y2 stepen i 300 Y1 stepen i Y 3 stepen i Y 4 stepen i a) б) 0 5 10 15 20 25 30 35 t 0 5 10 15 20 25 30 35 t i i Рис. 5.2. Описание уровня годовой добычи нефти из пласта моделями на базе многочленов: а) 1 и 2 степени, б) 3 (пунктир), 4 и 5 степеней Покажем решение для случая использования в (5.5) полинома третьего порядка с параметрами 1, A1, A2, A3, A4, A5 :

( ) 3 Yk = A1 ( k ) + A2 ( k ) + A3 ( k ) + A4 e 1k + A5 + k. (5.6) Динамической модели (5.6) соответствует при k 4 ARMA-модель четвертого порядка:

Yk = (Yk 5 Yk 4 ) µ14 + 4 (Yk 3 Yk 4 ) µ1 + (5.7) + 6 (Yk 3 Yk 2 ) µ12 + 4 (Yi 1 Yi 2 ) µ1 + Yk 1 + k, где µ1 = e1, k – «новая» автокоррелированная стохастическая компо нента того же вида, что и (5.7), но в уровнях k.

Для (5.7) нормальное уравнение реализации МНК для поиска µ примет вид полиномиального уравнения с одним неизвестным:


l1µ1 + l2 µ1 + l3µ1 + l4 µ14 + l5 µ1 + l6 µ12 + l7 µ1 + l8 = 0, 7 6 5 (5.8) где N l1 = 4 Yk25 2Yk 4Yk 5 + Yk24, k = N l2 = 28 Yk 5Yk 3 Yk 5Yk 4 Yk 4Yk 3 + Yk24, k = N l3 = 12 11Yk 4Yk 3 + 4Yk24 + 4Yk23 3Yk 5Yk 2 + 3Yk 5Yk 3 + 3Yk 4Yk 2, k = N l4 = 20 6Yk23 + Yk 1Yk 5 + 7Yk 4Yk 2 Yk 5Yk 2 6Yk 3Yk k = Yk 1Yk 4 6Yk 4Yk 3 ], N l5 = 4 52Yk 3Yk 2 + 18Yk23 + Yk 1Yk 5 17Yk 1Yk 4 + 16Yk 4Yk 2 + k = +16Yk 1Yk 3 YkYk 5 + Yk Yk 4 + 18Yk22, N l6 = 12 6Yk 3Yk 2 + 6Yk22 + Yk Yk 4 6Yk 1Yk 4 + 7Yk 1Yk k = 6Yk 1Yk 2 YkYk 3 ], N l7 = 4 11Yk 1Yk 2 + 4Yk22 + 3Yk 1Yk 3 + 4Yk22 3Yk Yk 3 + 3YkYk 2, k = N l8 = 4 YkYk 2 + Yk21 Yk 1Yk 2 YkYk 1.

k = После компенсации автокоррелированности и решения уравнения (5.8) с относительно МНК-оценки µ1 можно осуществить расчет 1 = ln µ1.

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

Второй этап МНК-оценки параметров A1, A2,..., An + 2 очевиден и по требует решения нормальной СЛАУ пятого порядка.

Рассчитанный предложенным методом геологический риск паде ния добычи нефти для пластов Турнейского яруса шести месторожде ний составил 14% со среднеквадратичным отклонением 6%. Отдельные показатели CV оказались равными 11,43%;

15,18%;

26,09%;

12,41%;

9,58%;

10,46%.

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

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

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

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

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

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

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

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

Yk = C1e1k + C2e 2k + k, (5.9) где стохастическая компонента k удовлетворяет принятым условиям Гаусса-Маркова.

При предложении модели (5.9) исходили из предпосылки, что од на из экспонент модели уровней ряда Yk, например, Yk1 = C1e1k, форми рует в большей мере тренд кривой жизненного цикла на этапе роста, а другая ( Yk2 = C2e 2k ) – на этапе спада (рис. 5.3).

Sales Y2=C2e2k k Yk=Y1+Y k k - Y1=C1e1k k - - 2 4 6 8 t Рис. 5.3. Иллюстрация моделирования жизненного цикла продаж автомобилей Мустанг (тыс. шт.) моделью (5.9) ЖЦП описывает и более сложный комплекс моделей:

Yk = ( C1 + D1k ) e1k + ( C2 + D2 k ) e 2 k + k, (5.10) Yk = С1e1k + С2e2k + С1e1k ( A1 sin (1k ) + B1 cos (1k ) ) + (5.11) + С2e2k ( A2 sin (2k ) + B2 cos (2k ) ) + k, Yk = С1e1k + С1e1k ( A1 sin (1k ) + B1 cos (1k ) ) + (5.12) + A2 sin (2k ) + B2 cos (2k ) + k, Yk = D1 + D2k + e1k ( A1 sin (1k ) + B1 cos (1k ) ) + (5.13) + A2 sin (2k ) + B2 cos (2k ) + k, Yk = D1 + D2k + D3 ( k ) + e1k ( A1 sin (1k ) + B1 cos (1k ) ) + (5.14) + A2 sin (2k ) + B2 cos (2k ) + k.

Отметим, что (5.11), (5.12), (5.13) и (5.14) могут, в зависимости от конкретных приложений, моделировать колебательные компоненты временного ряда жизненного цикла сезонными гармониками с некрат ными частотами, а также применяться к задаче моделирования эконо мических циклов, рассмотренных, например, в работе [98].

Модели (5.9) и (5.10) можно рассматривать как развитие извест ных моделей [82, 100], а (5.5) – как расширение модели Рамсея (2.33).

Идентификация модели (5.5) показана в предыдущем параграфе, а моделям (5.9) (при k 2 ) и (5.10) (при k 4 ) соответствуют ARMA модели второго и четвертого порядков соответственно:

Yk = (1 + 2 ) Yk 1 1 2Yk 2 + k, (5.15) ( ) Yk = 2 (1 + 2 ) Yk 1 12 + 412 + 2 Yk 2 + 212 (1 + 2 ) Yk (5.16) 122 Yk + k где 1 = e1, 2 = e 2.

Нетрудно видеть, что в (5.15) и (5.16) стохастические компоненты k и k также удовлетворяют принятым условиям Гаусса-Маркова, за исключением наличия автокорреляции, компенсация которой может быть обеспечена, например, прореживанием выборки.

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

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

Для модели (5.11) тренд описывается моделью (5.9), а колебатель ная компонента – выражением:

YkS = C1e1k ( A1 sin (1k ) + B1 cos (1k ) ) + (5.17) + C2e1k ( A2 sin (2k ) + B2 cos (2k ) ) + k.

Коэффициенты 1, 2 при параметрической итерационной деком позиции будут известны по результатам идентификации тренда ARMA моделью (5.9) (через 1, 2 очевиден расчет параметров 1, 2 ), а линей ные параметры модели (5.9) C1, С2 определит МНК.

ARMA-модель колебательной компоненты (5.17) будет иметь при k 4 вид:

( ) Yk = 2 ( µ11 + µ22 ) Yk 1 12 + 2 + 4µ1µ212 Yk 2 + (5.18) ( ) 122 µ2 Yk 3 122 Yk + 12 µ1 + k + где µ1 = cos (1 ), µ 2 = cos ( 2 ).

В модели (5.18) коэффициенты 1, 2 находят итерациями после идентификации тренда, а неизвестными являются только коэффициенты µ1, µ 2, которые входят в уравнение нелинейно.

Модели (5.12)-(5.14) при применении параметрической итераци онной декомпозиции состоят, соответственно, из трендов:

YkT = C1e 2t + k, (5.19) YkT = D1 + D2 k + k, (5.20) YkT = D1 + D2 k + D3 ( k ) + k (5.21) и колебательной компоненты YkS = C1e1k ( A1 sin (1k ) + B1 cos (1k ) ) + A2 sin (2k ) + (5.22) + B2 cos (2k ) + k ARMA-модель тренда (5.19), как уже было показано в третьей гла ве, имеет при k 1 вид Yk = 1Yk 1 + k, (5.23) а тренды (5.20), (5.21) линейны относительно своих параметров, кото рые легко находятся с помощью классического МНК. Напомним, что выше показана возможность идентификации таких трендов и при структуре мультипликативной стохастической компоненты.

ARMA-модель колебательной компоненты (5.22) при k 4 имеет вид:

( ) ( ) Yk = 2 ( µ2 + µ11 ) Yk 1 1 + 12 + 4µ1µ21 Yk 2 + 2 µ11 + 2µ212 Yk (5.24) 12Yk 4 + k.

В случае модели (5.12) параметр 1 в (5.23) может быть найден из идентификации тренда, а в случае моделей (5.13), (5.14) в модели (5.24) будут неизвестны все три параметра первого этапа идентификации 1, µ1, µ 2.

Заменой 1 = µ11, 2 = 12 степень многочлена модели (5.23) может быть понижена:

Yk = (Yk 2 + Yk 4 ) 2 + 2Yk 3µ2 2 + 2 (Yk 1 + Yk 3 ) 1 4Yk 2 1µ2 + (5.25) + 2Yk 1µ2 Yk 2 + k.

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

Выше уже была показана данными таблицы 3.1 возможность ис пользования для (5.15), (5.16), (5.18) и (5.25) решения систем нормаль ных нелинейных полиномиальных уравнений в базисах Гребнера.

Рассмотрим теперь устойчивость предложенных методов иденти фикации моделей (5.9)-(5.14) в диапазонах значений параметров моде лей и соотношений дисперсий «шум/полезный сигнал».

Устойчивость методов моделирования и прогнозирования иссле дована на тестовых выборках в диапазоне параметров, указанных в таб лице 5.4 для моделей (5.9)-(5.14), обозначенных M1-M3, а также в табли це 5.5 – для моделей (5.12)-(5.16), обозначенных M4-M6.

Соотношения дисперсий «шум/сигнал» K n / s задавались в диапазо не от 0 до 30%.

На рисунках 5.4 и 5.5 представлены графики коэффициента детер минации в функции соотношения мощностей помехи и полезного сигнала R 2, а также второго коэффициента Тейла в функции соотношения мощно стей помехи и полезного сигнала K T 2, являющихся результатом усредне ния 40 случайных реализаций шума для данных рисунка 5.4 и 20 реализа ций шума для данных рисунка 5.5.


Значения коэффициентов R 2, K T 2 говорят о достаточно высокой точности моделирования и прогнозирования моделей (5.9)-(5.11), (5.12) (5.14) при уровнях шума K n / s 10% ( R 2 0,93, а K T 2 20% ) и о доволь но высоком качестве моделирования и прогнозирования моделей (5.9) (5.11), (5.12)-(5.14) при уровнях шума до K n / s 30% ( R 2 0,85, а K T 2 30% ).

Точность идентификации модели (5.12) оказалась несколько ниже, чем точность идентификации других моделей.

Таблица 5. Диапазон параметров модельных задач исследования устойчивости алгоритмов идентификации моделей (5.9)-(5.11) к шуму 1 1 2 С1 D1 A1 B1 С2 D2 A2 B M Min -15 -0,3 3,00 -0, Max -3 -0,1 15,00 -0, M Min -15 -1,0 -0,3 3,00 0,20 -0, Max -3 -0,2 -0,1 15,00 1,00 -0, M Min -15 -0,3 0,01 0,01 0,3 3,00 -0,10 0,01 0,01 0, Max -3 -0,1 0,10 0,10 1,0 15,00 -0,03 0,10 0,10 1, 1 R2 T2,% 0.95 M3 M1 M1 M 0.9 M M 0.85 K N/S,% K N/S,% 0. 0 10 20 30 0 10 20 Рис. 5.4. Коэффициент детерминации R2 и второй коэффициент Тейла в T случае идентификации моделей (5.9) – M1, (5.10) – M2, (5.11) – M при различных сочетаниях дисперсий шума и сигнала K n / s Таблица 5. Диапазон параметров модельных задач исследования устойчивости алгоритмов идентификации моделей (5.12)-(5.14) к шуму 1 1 D1 D2 D3 С1 A1 B1 A2 B M Min 20 -0,10 1,2 0,1 0,1 0,6 10 Max 50 -0,01 1,5 0,7 0,7 0,9 30 M Min 20 5 -0,10 1,2 10 10 0,6 10 Max 50 10 -0,01 1,5 30 30 0,9 30 M Min 20 5 0,1 -0,10 1,2 10 10 0,6 10 Max 50 10 1,0 -0,01 1,5 30 30 0,9 30 1 R2 T2,% M 0. M 0. M6 M M 0. M KN/S,% KN/S,% 0. 0 10 20 30 0 10 20 Рис. 5.5. Коэффициент детерминации R2 и второй коэффициент Тейла T2 в случае идентификации моделей (5.12) – M4, (5.13) – M5, (5.14) – M при различных сочетаниях дисперсий шума и сигнала K n / s Исследуем теперь точностные предложенных методов идентифи кации моделей (5.9)-(5.14) на реальных данных.

Коэффициенты моделей (5.9)-(5.11) идентификации модели жиз ненного цикла автомобилей Мустанг производства «Форд», а также ко эффициенты детерминации R 2 и второго коэффициента Тейла K T 2 све дены в таблицу 5.6, а результаты построения моделей показаны на ри сунке 5.6.

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

Таблица 5. Идентификация жизненного цикла автомобилей Мустанг D 1 1 1 С2 D2 2 2 R Модель С1 A1 B1 A2 B2 T2, % Mustang M1 -1341 -0,83411 1604 -0,37337 0,979 22, M2 -688 -97 -1,11334 953 280 -0,46185 0,979 22, M3 -1130 -0,89715 0,085 -0,126 1,5006 1394 -0,34818 -0,020 -0,100 1,9792 1,000 10, Sales 600 Sales 600 Sales Mustang Mustang Mustang 400 200 M1 M2 M t0 t0 t 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 Рис. 5.6. Идентификация жизненного цикла автомобилей Мустанг.

Серый цвет – исходные данные, черная толстая линия – идентифицированная модель, черная тонкая линия – прогноз Исследуем интерес пользователей системы Google к сотовым те лефонам марок Nokia N70, N73, N82 и N95, используя данные о за просах пользователей, которые еженедельно агрегируются в стати стике Google Trends [94]. Смоделируем с помощью (5.9)-(5.11) ежеме сячную динамику индекса запросов (таблица 5.7, рисунок 5.7). Усред няя коэффициенты детерминации и Тейла по всем маркам телефонов, получим, что первая модель описывает индекс запросов с R 2 = 0,88, K T 2 = 12,3% ;

вторая с R 2 = 0,95, K T 2 = 14,6% ;

а третья – с R 2 = 0,97, K T 2 = 5,3%. Наилучшее приближение, как видим из приведенного ри сунка 5.7, дала модель (5.11).

Таблица 5. Идентификация интереса пользователей Google к моделям телефонов компании Nokia R 1 1 2 Модель С1 D1 A1 B1 С2 D2 A2 B2 T2, % N M1 -12,68 -0,06428 13,11 -0,05172 0,920 8, M2 1,25 0,59 -0,36200 -0,91 0,25 -0,06383 0,961 12, M3 -3,59 -0,08129 -0,089 -0,078 0,2736 4,06 -0,04022 0,012 -0,080 0,1166 0,975 7, N M1 -6,81 -0,08963 6,97 -0,04631 0,945 2, M2 -3,11 0,04 0,01212 3,39 0,35 -0,03977 0,966 16, M3 -5,76 -0,09569 -0,006 -0,101 0,1531 5,92 -0,04392 -0,006 -0,082 0,1342 0,965 3, N M1 -2,93 -0,11783 2,82 -0,07348 0,877 7, M2 -0,80 -0,35 -0,32285 0,82 0,09 -0,09630 0,950 12, M3 -4,25 -0,11295 0,006 -0,040 0,9326 4,15 -0,08170 -0,005 -0,015 0,2379 0,985 4, N M1 -6,55 -0,09093 6,71 -0,03445 0,787 30, M2 -12,69 -1,09 -0,13646 13,18 -0,20 -0,03170 0,922 16, M3 -2,73 -0,13739 -0,033 -0,295 0,7671 3,07 -0,01952 0,270 -0,194 0,7671 0,950 5, Search index 1.5 Search index 1.5 Search index 1. N70 N N 1 0. 0.5 0. M1 M2 M 0 t t t 1 2 3 1 2 3 4 1 2 3 Search index Search index Search index 2 2 N73 N N 1. 1.5 1. 1 0. 0.5 0. M1 M2 M 0 4t 4t 4t 1 2 1 2 3 1 2 Search index Search index Search index N82 N N 0. 0.4 0. 0. 0.2 0. M M1 M 0 3t 3t 1 2 3t 1 2 1 Search index Search index Search index 3 3 N95 N N 2 1 M M1 M 0 t t t 1 2 1 2 3 1 2 Рис. 5.7. Идентификация интереса пользователей Google к моделям телефонов компании Nokia. Серый цвет – исходные данные, черная толстая линия – идентифицированная модель, черная тонкая линия – прогноз Другим примером могут быть путеводители – класс товаров с дли тельным жизненным циклом. В последнее время, в связи с развитием ин тернет-технологий (интерактивных карт и др.), спрос на путеводители стал снижаться.

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

Для моделирования интереса пользователей системы Google к пу теводителям в целом (запрос «Travel guide»), так и к наиболее распро страненному изданию Lonely Planet применены выражения (5.12)-(5.14) (таблица 5.8, рис. 5.8).

Таблица 5. Идентификация интереса пользователей Google к путеводителям R 1 1 D1 D2 D3 С1 A1 B1 A2 B2 T2, % Travel Guide M4 1,67 -0,1343 0,488 0,092 -0,073 0,144 0,017 0,040 0,954 4, M5 1,547 -0,013 -0,0525 0,033 -0,612 0,284 0,483 0,099 -0,071 0,951 11, M6 1,694 -0,024 0,000130 -0,0283 0,480 0,243 -0,128 0,963 0,010 0,062 0,965 6, Lonely Planet M4 1,73 -0,1471 0,163 0,026 -0,042 0,490 0,041 -0,075 0,935 13, M5 1,513 -0,013 -0,0237 0,099 0,176 0,112 0,477 0,068 -0,051 0,936 17, M6 1,738 -0,025 0,000140 -0,0207 0,486 0,092 -0,129 0,963 0,021 0,074 0,953 6, Search index 2 Search index 2 Search index travel guide travel guide travel guide 1. 1. 1. 0. M4 M5 M 0. 0.5 6t 6t 6t 123 4 123 4 5 123 4 Search index Search index Search index 2 Lonely Planet Lonely Planet Lonely Planet 1.5 1. 1. 1 0.5 0. M4 M5 M 0. 0 6t 6t 6t 1 2 3 4 1 2 3 4 5 1 2 3 4 Рис. 5.8. Идентификация интереса пользователей Google к путеводителям.

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

Приведем еще один пример, когда модель (5.9) с успехом приме нена для описания ЖЦП при рассмотрении индекса запросов пользова телей поисковой системы Google на две модели телефонов Nokia (рис. 5.9).

Индекс поисковых запросов на Nokia n95 описывается моделью:

Yk = ( 0,269k 11,5) e0,034k + ( 0,045k + 11,8) e0,007 k + k с коэффициентом детерминации R 2 = 0,9182 (рис. 5.10 а).

Индекс поисковых запросов на Nokia 5300 описывается моделью:

Yk = ( 0,164k + 4,7 ) e 0,017 k + ( 0,080k 5,0 ) e 0,136 k + k с коэффициентом детерминации R 2 = 0,9209 (рис. 5.10 б).

а) б) Рис. 5.9. Индекс запросов на поиск Nokia n95 (a) и Nokia 5300 (б) в системе Google а) б) Рис. 5.10. Результат идентификации интереса пользователей системы Google к телефонам Nokia 5300 (a) и Nokia n95 (б) Сравнивая частоты построенных моделей интереса пользователей системы Google к телефонам Nokia, можно отметить, что Nokia вышла на рынок быстрее, чем Nokia n95, но медленнее уходила с рынка.

В данном параграфе моделирование тренда ЖЦП моделями (5.9) и (5.10) представляло как самостоятельный интерес, так и могло быть со ставной частью моделирования при параметрической итерационной де композиции моделями (5.11) и (5.12) с колебательной компонентой в исходных данных.

Рассмотрение идентификации моделей (5.9) и (5.10) при помощи ARMA-моделирования показало уже достаточно высокую точность, во всяком случае, при небольших значениях соотношения отношения мощностей помехи и полезного сигнала ( K n / s 10% ).

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

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

Случайным образом генерируется конечный набор пробных решений – «первое поколение». Оценка приспособленности решений текущего по коления осуществляется исходя из заданного критерия. В случае иден тификации параметров моделируемых функции (5.9) и (5.10) использо ван, как и для решения в базисах Гребнера, тот же критерий минимума суммы квадратов отклонений модельной функции от реальных данных.

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

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

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

Для идентификации модели трендов (5.9) и (5.10) был использован генетический алгоритм, реализованный функциями gatool системы МATLAB.

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

Таблица 5. Диапазон параметров модельных задач исследования алгоритмов идентификации моделей M1 (1) и M2 (2) 1 C1 D1 C2 D M -15 -0,3 3 -0, Min -3 -0,1 15 -0, Max M -15 -1,0 -0,3 3 0,2 -0, Min -3 -0,2 -0,1 15 1,0 -0, Max Генерировались программно выборки из 24 точек «полезных сиг налов» (5.9) и (5.10) при различных сочетаниях данных из таблицы 5.9.

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

На рисунке 5.11 представлены полученные для модели (5.9) гра фики критерия точности моделировании (коэффициента детерминации R 2 ) и критерия точности прогнозирования (второго коэффициента Тей ла K T 2 ) с горизонтом прогноза 20% от длины исходного ряда.

Результат с полным основанием можно считать репрезентатив ным: каждое значение графиков R 2 и K T 2 в функции коэффициента от ношения мощностей K n / s помехи k и полезного сигнала (соответст вующим моделям (5.9) и (5.10)) является средним из 500 случайных реализаций шума и полезного сигнала (во всем динамическом диапазо не сочетаний его параметров).

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

На рисунке 5.12 представлены соответствующие графики для мо дели (5.10).

T2, % R21 0. ARMA algorithm Gin algorithm Gin algorithm 0. 0. ARMA algorithm 0.92 KN/S, % KN/S, % 0 10 20 0 10 Рис. 5.11. Коэффициент детерминации R 2 и второй коэффициент Тейла K T в случае идентификации модели (5.9) при различных сочетаниях дисперсий шума и сигнала K n / s T2, % R 0.95 ARMA algorithm Gin algorithm 0. ARMA algorithm 0. Gin algorithm 0.8 K,% K,% 0 10 20 30 n/s 0 10 20 30 n/s Рис. 5.12. Коэффициент детерминации R 2 и второй коэффициент Тейла K T в случае идентификации модели (5.10) при различных сочетаниях дисперсий шума и сигнала K n / s Из рисунков 5.11 и 5.12 видим, что оба метода идентификации дают довольно высокую точность моделирования и прогнозирования.

При этом генетический алгоритм дает существенно меньшую погреш ность прогнозирования и несколько лучшую точность моделирования на модельных выборках при уровне «сигнал-шум» более 2%.

Интересно, что при малом шуме (до 2%) преимущество в модели ровании и в прогнозировании имеет ARMA-модель.

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

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

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

Объем продаж Тренд Первич Повторный ный цикл цикл Время Рис. 5.13. Кривая жизненного цикла товара с повторным циклом Более того, зачастую реальные временные ряды жизненных цик лов продуктов содержат, кроме тренда с повторным циклом, еще и ко лебательную компоненту.

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

( ) Yk = YkT 1 + YkS + k, (5.26) где YkT – параметрическая модель тренда, YkS – параметрическая модель колебательной компоненты, k – стохастическая компонента, удовле творяющая принятым условиям Гаусса-Маркова.

Для идентификации модели (5.26) целесообразно применить метод параметрической итерационной декомпозиции тренда и колебательной компоненты. Рассмотрим последовательно три модели гладких трендов с повторным циклом.

В качестве первой выберем новую модель:

YkT = C1e1k + C2e 2k + C3e3k + C4e 4k + k. (5.27) Вторую модель можно рассмотреть как сумму колоколообразных гауссовых кривых или моделей Хабберта [103, 116]:

2 1 ( k µ1 ) 2 ( k µ2 ) 2 + k.

Yk = C1e + C2e (5.28) Третью модель тренда жизненного цикла предложим в виде сум мы двух дробно-рациональных функций:

YkT = Tk1 + Tk2 + k, (5.29) P0 + Pt P0 + Pt где Tk1 =, Tk2 = 1, причем суммирование произво 1 + Q1t + Q2t 1 + Q1t + Q2t дится в случае неотрицательных значений каждого из циклов тренда. В случае если Tki 0 в модели считается, что Tki = 0.

Колебательную компоненту будем моделировать двумя моделями.

Первая – достаточно традиционная, состоит из трех гармоник ряда Фу рье:

Yk = A1 sin (1k + 1 ) + A2 sin ( 21k + 2 ) + A3 sin ( 31k + 3 ) + k. (5.30) Вторая модель позволяет рассмотреть более общий случай не кратных частот:

Yk = A1 sin (1k + 1 ) + A2 sin ( 2 k + 2 ) + A3 sin (3k + 3 ) + k. (5.31) Идентификацию модели (5.27) возможно осуществлять с помо щью предложенных моделей ARMA-моделей.

ARMA-модель тренда в этом случае имеет вид:

YkT = (1 + 2 + 3 + 4 ) YkT (12 + 13 + 14 + 23 + 24 + 34 ) YkT2 + (5.32) + (123 + 124 + 134 + 234 ) YkT3 1234YkT4 + k, j где j = e, k – линейная форма того же вида, что и для YkT в (5.32) в уровнях и от 1, 2, 3, 4.

k Применение МНК для идентификации параметров ARMA модели (5.32) приводит к необходимости решения системы четырех полиноми альных уравнений от переменных 1, 2, 3, 4. Использование базисов Гребнера, как следует из таблицы 3.1 при таком числе независимых пе ременных требует большой разрядной сетки компьютера, может при вести к неустойчивости решения. Можно использовать в этом случае метод гомотопического продолжения, который значительно улучшен в последние десятилетия коллективом во главе с J. Verschelde [117], или показанный выше генетический алгоритм.

Для идентификации модели тренда (5.28) был реализован генети ческий алгоритм оптимизации функцией gatool системы MATLAB.

Идентификации модели (5.29) проводилась при помощи взвешен ного МНК.

Покажем моделирование кривых жизненного цикла на примере мыла торговой марки «Absolut», выпускаюшегося одним из предпри ятий г.о. Самара.

В качестве рабочей части для моделирования временного ряда ис пользованы данные продаж за 2001-2009 годы. Для контрольной (про гнозной) части ряда, расчета MAPE-оценки прогноза и второго коэффи циента Тейла K T 2 взяты данные за 2010 год и с января по май 2011 года.

На рисунках 5.14, 5.15 и 5.16 представлены исходные данные про даж и графики моделей временного ряда, построенные на ежемесячных данных, а также на сгруппированных данных по кварталам и на сгруп пированных данных по полугодиям.

Сезонная компонента (5.30) использовалась для всех моделей тренда, а сезонная компонента (5.31) – только для модели (5.29) тренда, который оказался наиболее точным по значению коэффициента детер минации R 2 для выборки.

Y Y Trend Model Trend Model 400 Trend 200 200 Trend phase phase Data Data 0 а) б) 2002 2004 2006 2008 2010 t 2002 2004 2006 2008 2010 t Y Y 600 Trend Trend Model Model 400 Trend Trend 200 Trend 200 Trend phase 1 phase phase 2 phase Data Data 0 в) г) 2002 2004 2006 2008 2010 t 2002 2004 2006 2008 2010 t Рис. 5.14. Моделирование ежемесячных продаж: а)-в) тренды (5.27), (5.28), (5.29) и сезонная компонента (5.30);

г) тренд (5.29) и сезонная компонента (5.31) Y Y Trend Trend Model 1500 Model 1000 Trend Data phase 2 Trend 500 Data phase 0 а) б) 2002 2004 2006 2008 2010 t 2002 2004 2006 2008 2010 t Y Y Trend Trend Model Model 1500 1000 Trend Data Trend Data Trend Trend 500 phase 2 500 phase phase 1 phase 0 в) г) 2002 2004 2006 2008 2010 t 2002 2004 2006 2008 2010 t Рис. 5.15. Моделирование ежеквартальных продаж: а)-в) тренды (5.27), (5.28), (5.29) и сезонная компонента (5.30);

г) тренд (5.29) и сезонная компонента (5.31) Y Y Data Trend 3000 Model 2000 Trend Data Model Trend phase 2 Trend 1000 phase 0 t 0 t а) б) 2002 2004 2006 2008 2010 2002 2004 2006 2008 Y Y Trend Trend 3000 Model Model 2000 Data Data Trend Trend 1000 Trend Trend phase 1 phase phase 2 phase 0 t 0 t в) г) 2002 2004 2006 2008 2010 2002 2004 2006 2008 Рис. 5.16. Моделирование продаж по полугодиям: а)-в) тренды (5.27), (5.28), (5.29) и сезонная компонента (5.30);

г) тренд (5.29) и сезонная компонента (5.31) Оценки точности построенных моделей приведены в таблице 5. и на рисунке 5.17.

Из таблицы 5.10 видим, что модели тренда (5.28) и (5.29) приводят к наилучшему коэффициенту детерминации R 2, а модели тренда (5.27) и (5.29) дают лучший прогноз. Другим достоинством моделей тренда (5.28) и (5.29) является наличие явно выраженных тенденций для перво го и второго циклов продаж, что позволяет сравнивать модели этих циклов с продажами отдельных марок товарной линейки.

Таблица 5. Точность моделирования и прогнозирования моделей Модель R2 Mape T Модель колебат.



Pages:     | 1 |   ...   | 4 | 5 || 7 |
 





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

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