ИДЕНТИФИКАЦИЯ МОДЕЛЕЙ ЖИЗНЕННОГО ЦИКЛА ПРОДУКЦИИ НА ОСНОВЕ МОДЕЛЕЙ АВТОРЕГРЕССИИ - СКОЛЬЗЯЩЕГО СРЕДНЕГО И БАЗИСОВ ГРёБНЕРА Текст научной статьи по специальности «Математика»

№ 1 (25) 2012

В. К. Семёнычев, Е. И. Куркин, Е. В. Семёнычев

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

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

Ключевые слова: модель жизненного цикла продукции; авторегрессия; МНК; базисы Грёбнера; автомобиль; сотовый телефон; путеводитель.

JEL classification: C22; D91.

1. Введение

Жизненный цикл продукции (товара или услуги) характеризуется динамикой объемов продаж и получаемой прибыли от момента начала разработки новой продукции до ее ухода с рынка. Исследование жизненного цикла осуществляется для принятия решений по минимизации затрат, увеличению прибыли, продлению наиболее прибыльных его этапов путем изменения товарной политики, объемов производства и сбыта, уровня цен, методов продвижения товаров и др. Типичный жизненный цикл продукции состоит из этапов вывода на рынок, роста, зрелости и упадка (Котлер, 2007). В общем случае можно подбирать модель сразу для всех этапов жизненного цикла, но интерес могут представлять и модели отдельных этапов эволюции: например, роста и/или спада, нахождение моментов смены видов или параметров моделей для оценки принятых управленческих или маркетинговых решений, для оценки воздействия внешней среды. Реализация моделирования жизненного цикла продукта является актуальной и далеко не простой задачей в силу сложности и нелинейности требуемых обычно аналитических выражений, а также из-за необходимости использования зачастую малых (менее тридцати наблюдений) выборок для «молодых» (инновационных) процессов и мониторинга эволюции моделей.

Известно много моделей жизненного цикла продукции (Оразбаев и др., 2007; Мазуркин, 2002; Соколов, 2006; Brody, 1945; Aman, Ohkochi, 2010; Hook et al., 2011; Grubler, 1990), как правило, специализированных по отдельным видам продукции. Однако исследования по их возможной области применения (по динамическому диапазону параметров моделей и по мощности стохастической компоненты в рядах динамики), а также по учету колебательной компоненты, часто фиксируемой в соответствующих временных рядах, не обнаружены.

122 J--=

Теория и методология •

№ 1(25) 2012

2. Модели жизненного цикла продукции

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

Yk = Qe aik D+ C2 ea*hD+ek,

(1)

где, здесь и далее, величины ек удовлетворяют условиям Гаусса-Маркова (Айвазян, 2001), чтобы при идентификации с помощью метода наименьших квадратов (МНК) получать оптимальные оценки параметров модели (1), к = 1,2,..., п — номер наблюдения, п — объем выборки, А — период дискретизации.

Уравнение (1) исходит из предпосылки, что одна из экспонент модели уровней ряда Ук , например Yl = С1еакА, формирует в большей мере тренд кривой жизненного цикла на этапе роста, а другая (Yк2 = С2еакА) — этап спада (рис. 1).

1500 1000 500 0

-500 -1000 -1500

Sales

/

0 ф г

1 'ф

§

ф

0

со ui

1 I

s: ui

(о Ф г

I 'ф

§

ф О

со

Рис. 1. Иллюстрация моделирования жизненного цикла продаж автомобилей Мустанг (тыс. шт.) соотношением (1)

Источник: (Котлер, 2007).

Известны примеры применения модели тренда с использованием экспонент (В1^у, 1945), в виде произведения степенной функции с дробным показателем и экспоненты (Ораз-баев и др., 2007; Соколов, 2006). Известен также принцип формирования феноменологических моделей в виде сумм и/или произведений некоторых других функций (Мазуркин, 2002), описывающих этапы роста и спада.

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

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

Yk = (( + D2kА + ... + Dn (к А)" ) + 4 , (2)

Yk = ( + ЦкА)е а1кА + (С2 + D2 к А)е а кА + ек, (3)

№ 1 (25) 2012

э-

р

Yk = C1ea<kA + C2eakD + CfakA ( sin (o^kA) + B1 cos (o^kA)) + + C2eakA (( sin (o2kA) + B2 cos(o2kA)) + ek, Yk = C1e a'k A + C1e a'kA ( sin (o1k A) + B1 cos (o1k A)) +

+ A2 sin (o2 k A) + B2 cos (o2 k A) + ek, Yk = D1 + D2kA + eakA ( sin (o1kA) + B1 cos (o1kA)) + + A2 sin (o

2 k A) + B2 cos (o 2 k A) + ek, Yk = D1 + D2kA + D3 (kA)2 + eakA (A1 sin (o1kA) + B1 cos(o1kA)) + + A2 sin (o 2k A) + B2 cos (o 2 k A) + ek.

(4)

(5)

(6) (7)

<u

Ц Отметим, что уравнения (4) - (7) могут, в зависимости от конкретных приложений, мо-

¡5 делировать сезонные компоненты временного ряда жизненного цикла гармониками с не-0

g кратными частотами, а также применяться к задаче моделирования экономических циклов,

^ рассмотренных в работах (Beveridge, Nelson, 1981; Hodrick, Prescott, 1997; Enders, 2004).

| Модели (1) и (3) можно рассматривать как развитие известных моделей (Brody, 1945; Aman,

a Ohkochi, 2010), а (2) — как обобщение моделей работы (Оразбаев и др., 2007).

га

« 3. Использование АЯМА-моделей для жизненного цикла продукции

§

¡5 Идентификация нелинейных по параметрам моделей (1) - (7) является довольно сложной

си

§ задачей, особенно на коротких выборках. Ее решение можно получить на основе обобщенна ных параметрических моделей авторегрессии - скользящего среднего (АЯМА-моделей) ря-^ да наблюдений (Семёнычев, Семёнычев, 2011). Идентификация при этом осуществляется § в два этапа. На первом строится (с помощью 2-преобразования) АЯМА-модель из наблю-

дений Ук , в которой при помощи МНК идентифицируются ее коэффициенты. Затем через Ц них рассчитываются некоторые нелинейные параметры моделей: в (1) - (3) это показатели экспонент, а в (4) - (7) — показатели экспонент и частоты гармоник. На втором этапе, ис-

| пользуя результаты первого этапа и исходные данные, с помощью МНК идентифицируют

¡э оставшиеся параметры моделей — амплитуды гармоник и коэффициенты полиномов.

■с Например, модели (2) соответствует (при к > п) АЯМА-модель п-го порядка

I п+1

8 ¥к=2(-1)с+1 К,--у- )+Хк, (8)

®

| где = еа'Л, СП+1 — биномиальные коэффициенты, а Хк будет АЯМА-моделью того же ви-

§ да, что и (8), но записанной в величинах стохастической компоненты ек. §

а-

| Соотношениям (1) и (3) соответствуют АЯМА-модели: (9) — второго порядка и (10) —

■fr четвертого порядка соответственно:

I Y = (h +h2 )Yk-1 -h1 h2Yk-2 +Xk , при k > 2, (9)

124 у =

Теория и методология •

№ 1(25) 2012

Yk = 2 ( +h2 K-! -(hf + 4hih2 + h2 )-2 + 2hi h2 (hi +h2 )Yk-3 -hf h2Yk-4 + Xk , IS

, . (10) 3

при к > 4 ,

is

где, здесь и далее, hi = ea'A, h2 = ea'A. ^

ai

Очевидно, что Хк, описываемые моделями (9) и (10), но в величинах стохастической g компоненты ек, также удовлетворяют условиям Гаусса-Маркова, за исключением наличия автокорреляции, которая может быть скомпенсирована с помощью прореживания исход- ^ ной выборки (Семёнычев, Семёнычев, 2011). Прореживание заключается в использовании uj выборок через одно, два и т. д. наблюдения и выбора тех идентифицируемых оценок коэф- ® фициентов ARMA-моделей, для которых окажется лучше принятый критерий точности мо- £ делирования. Исследования показали, что автокорреляция существенно уменьшается при § шаге прореживания, равном порядку ARMA-моделей. ^

Для моделей (4) - (7) целесообразно применить метод параметрической итерационной щ декомпозиции тренда и колебательной компоненты (Семёнычев и др., 2010). Он, в отличие от известного непараметрического метода сезонной декомпозиции (Айвазян, 2001), на каждой итерации использует оценки компонент ряда при помощи ARMA-моделей для гармоник и экспонент. При этом сначала идентифицируется тренд, а затем, после его вычитания из уровней временного ряда — колебательная компонента. Далее находят уточненную оценку тренда, потом — уточненную оценку колебательной компоненты и т. д. Прекращение итерирования происходит при стабилизации меры точности моделирования. Как показали исследования, при этом удается повысить вычислительную точность за счет уменьшения порядков конструируемых ARMA-моделей на каждой итерации и размерности соответствующих нормальных уравнений, получая высокую точность параметрического моделирования уже при трех-четырех итерациях.

Так, для модели (4) тренд описывается соотношением (1), а колебательная компонента — выражением

YkS = СхеакА ( sin (wikА) + B cos (wikА)) + С2ea'kA (( sin (w2kA) + B2 cos (w2кД)) + Xk. (11)

Коэффициенты hi, h2 при параметрической итерационной декомпозиции будут известны по результатам идентификации тренда ARMA-моделью (9) (через величины hi, h2 можно рассчитать ai, а2). Линейные параметры Ci и С2 модели (9) определит МНК. ARMA-модель колебательной компоненты в (11) после параметрической итерационной декомпозиции будет иметь при к > 4 вид

Yk =2((hi h2)Yk-i -(hi2+h2+4mmhih2)+2(hfh2m+hih2m)-hi2h2Yk-4 +xk. (12)

Здесь и далее mi = cos (wi A), m2 = cos (w2 A) . В (12) коэффициенты hi, h2 находят итерациями после идентификации тренда, а неизвестными являются коэффициенты mi, m2 (или, с учетом принятых обозначений, параметры модели wi, w2).

Модели (5) - (7) при применении параметрической итерационной декомпозиции состоят из трендов

YT = СхеакА+гк, (13)

YT = Di + D2kA + 8k , (14)

№ 1 (25) 2012

YT = D1 + D2 k A + D3 (k A)2 + ek (15)

и колебательной компоненты

YkS = CieakA (A sin (wjkA) + Bi cos(WjkA)) + A2 sin (w2kA) + B2 cos(w2kA) + Xk. (16) ARMA-модель тренда (13) имеет при k >1 вид

7k = hiYk-i +Xk, (17)

а тренды (14) и (15), линейные относительно коэффициентов полиномов, легко находятся с помощью МНК.

ч АЯМА-модель колебательной компоненты (16) после параметрической итерационной ^ декомпозиции при k > 4 имеет вид:

Yk = 2(/¿2+mihi) -(i+hi2 + 4mmhi)+2( + 2,^2h)-hiYk-+X. (18)

15 В случае модели (5) параметр ^ может быть найден из идентификации тренда, а в случае

I

моделей (6) и (7) неизвестны все три параметра в модели (18): ^, т, т2. Заменой с = т^, с2 = степень многочлена модели (18) может быть понижена:

ф

Ук = _ (у_2 + Тк_4) с + 2¥к_3т2С + 2 (-1 + у- ) - 4Гк_2С1 т + 2У- т - У- + Хк. (19)

р

со МНК-идентификация коэффициентов модели (8) приводит к нормальному алгебраическому уравнению одной переменной, а (17) — к нормальному линейному (относитель-| но ) уравнению. В остальных случаях будут более сложные нормальные нелинейные § системы полиномиальных уравнений. Распространенные итерационные численные мето-§ ды решения таких систем (к примеру, метод Ньютона) имеют существенные недостатки: о их сходимость и набор найденных корней зависят от выбора начального приближения или ^ диапазона допустимых значений корней, существует также и проблема разделения близких | корней. Поэтому в данной работе предлагается использовать свободное от этих недостатков £ решение в базисах Грёбнера.

| 4. Базисы Грёбнера

а-

е

§ Метод базисов Грёбнера сводит систему нормальных уравнений к треугольному виду, а задачу ее решения — к решению последовательности алгебраических уравнений, что легко | сделать программным путем. Базисы Грёбнера позволяют найти все корни системы и не требуют для них начального приближения (Ершов, 2007). Метод базисов Грёбнера имеет, одна-| ко, свои недостатки. Например, при его реализации не должны происходить округления при § операциях с коэффициентами. Поэтому вместо вещественных коэффициентов применяются | целые, которые быстро растут в процессе вычислений из-за наличия в алгоритме функций | вычисления наименьшего общего кратного. Скорость и результат вычисления базиса Грёб-нера могут меняться в зависимости от различных факторов: упорядочивания переменных, ■ё использования различных методов (в частности, критериев Бухбергера тривиальности кри-§ тических пар полиномов) и т. д. (Ершов, 2000).

126

(25) 2012

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

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

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

0 ф г

1

¡ф §

ф

0 со ui

1 I

s: ui

(о Ф г

I 'ф

§

ф О

со

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

Величина коэффициентов

1- -10 50-100 1000 -10000 10000 -100000

2 3 4 5 6 7 2 3 4 5 6 2 3 4 5 6 2 3 4 5

2 1 1 1 1 2 3 1 111 3 1 1 113 11 13

3 1 1 1 3 1 13 1 1 3 12 3

4 1 1 3 1 23 1 3 13

5 1 2 3 1 3 1 3 23

6 1 3 1 3 1 3 23

7 1 3 1 3 1 3 33

8 2 3 2 3 3 3 33

9 2 3 2 3 3 3 33

10 2 3 3 3 3 3 33

Примечание. N (по столбцам) — количество переменных, М (по строкам) — их максимальная степень. Результаты вычислений: 1 — базис Грёбнера применен успешно; 2 — базис Грёбнера вычисляется, но его корни не совпадают с корнями системы из-за накопления больших вычислительных погрешностей; 3 — базис Грёб-нера не вычисляется.

Из таблицы 1 следует, что базисы Грёбнера могут быть использованы для решения систем полиномиальных уравнений не более чем с 5-ю переменными, входящими не более чем в 7-ой степени (выделено в таблице серым фоном).

Для обеспечения точности идентификации моделей экономической динамики необходимо как минимум 3- 4 значащих цифры, поэтому лучшие результаты применение базиса Грёбнера дает для систем с 2-3 переменными, входящими в степени не выше 4.

№ 1 (25) 2012

5. Идентификация моделей с помощью базисов Грёбнера

■ф

е-

со о о

S

и

5

О ф

t

6

2 ф ■J

£

и

л §

о

S

5

о

о ф

6 ф

а о

>s

ф

§ §

ф eg

0

1 о

0

W

1 5 S

1

w

Si

S

a

2

0

1 I

ф

I

§

>s

ф

§ § §

as

¡5 i ф

(20)

В качестве примера идентификации моделей с помощью базисов Грёбнера покажем, что для модели (11) МНК приводит к системе нормальных нелинейных систем полиномиальных уравнений вида

11,1^14 й2 + ¿1,2 йЗ +11,3^14 й + Л,4 ^ +11,5^13 й2 + 11,6^13 + + 11,7^2 й2 +11,8^3 й + Л,9 й й2 + ¿1,10 +^1,11^2 + ¿1,12 ^ й2 + + ^1,13^й + 11,14й +\15йЗ + ^1,16й2 + ¿1,17 й й +11,18й + ¿1,19й = 0, ^2,1^2 йЗ +¿2,2 й2 йЗ + ¿2,3 й2 й + ¿2,4 й2 + ¿2,5 йЗ й2 + ¿2,6 йЗ й2 + +12,7 ^3й2 + ¿2,8йЗй + ¿2,9й й2 + ¿2,10^3 + ^2,11^2 + ¿2,12^+ +¿2,13 йЗ й + ¿2,14 й й2 + ¿2,15 йЗ + ¿2,16 й2 + ^2,17 й й + ¿2,18 й + ^2,19 й = 0.

Например, при идентификации модели (3), построенной для первых восьми лет выпуска автомобилей Мустанг, (20) примет вид:

-222 +1169^2 _ 2256^2 + 1915^2 _ 587^ + 871^ _ 4512^й2 +

■ 6666hjh2 +1895hj

1149^?+5746^2h2

+ 8474hjh2

+7251ft2h2 —1906h?h? + 518h3 - 2349h3h2 + 3790h3h2

- 9999h?h2 + ■2541h3 h? + 662h3 h4 = 0,

(21)

—222 + 1169h1 — 2256h? + 1915h3 — 587 ft4 + 871h2 — 4512h1 h2 + + 8474h? h? — 6666h? h3 + 1895h? h4 — 1149h? + 5746h? h1 — 9999h? h? + +7251h? h3 — 1906h? hi + 518h? — 2349h? h1 + 3790h? h? — 2541h? h3 + 662h? hi = 0.

В (20) переменные h1, h2 входят в уравнения в степени не выше четвертой, что удовлетворяет требованиям табл. 1, поэтому может быть применен метод базиса Грёбнера.

Для вычисления базисов Грёбнера использовалась функция Groebner [Basis] системы Maple. Такой метод сводит систему (20) к уравнению 25-ой степени относительно одной переменной h1 и линейному уравнению относительно h2 , в котором h1 выступает как параметр.

В силу симметрии модели (3) уравнения, определяющие h1 и h2, совпадают и имеют вид:

juh? + jl 2h?4 + ... + jl 25h1 + jl 26 = 0.

(22)

В рассматриваемом примере идентификации модели (3) первое уравнение базиса Грёбнера примет следующий вид (разрядность коэффициентов в (23) высокая, например, коэффициент при й25 — сотни квадриллионов (см. (23)).

Действительные корни уравнения (23) равны - 0.2745883443, 0.3280396491 и 0.6259083348. Переходя к показателям экспонент а1 = 1п (й1) и исключая первый корень (как лишенный смысла), получим, что а1 =-1.11 или а1 =-0.46 . Поскольку модель (3) симметрична, те же вычисления справедливы и для а2. Рассматривая возможные варианты сочетаний решений а, а2, получим, что лучшее моделирование (определяемое значением показателем точности моделирования) достигается при а1 = -1.11 и а2 =- 0.46.

128

№ 1(25) 2012

33061088041489760^25—449904585018573406^24+288071460651089971Ц23— |

— 11473089010404936366^22+31845164619533836137^21—65778533047042063560^+

+ 106345160863694313906^9—142026093076653314179^+168076326537052466587^ -

—190105771871483745187^6+212085302717314069442^5—223859608455635825467(23) 3 +207827724106180827459^3—159400586765441931405^+95711270696912673117—

—41716792177369901470^1°+10770145559766427305^9+236628362932400559^— ^

— 514407175897667360^+601949179387679673^—77137194321734878^— ^

— Ш77549243240366^4+8105011388113752^3-629506244676756^2— | —145645410637800^ +23484803022888=°. |

Применение МНК к АЯМА-моделям (9) - (11) и (19) приводит к системам полиноми- щ альных уравнений, похожих на (20), причем количество неизвестных в моделях, а также их максимальные степени удовлетворяют условиям табл. 1, поэтому алгоритм построения базиса Грёбнера реализуется аналогично.

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

6. Методика оценки точности и области применения идентификации

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

Предложим более общий подход: будем рассматривать точность идентификации моделями (1) - (7) в широких диапазонах значений их параметров и при соотношении мощностей «шум/полезный сигнал» от 0 до 30%. Исследуем тестовые выборки в диапазоне параметров, указанных в табл. 2 (для моделей (1) - (4), обозначаемых далее как М1 - М4), и в табл. 3 (для моделей (5) - (7), обозначаемых далее как М5 - М7 ).

Соотношения дисперсий «шум/сигнал» задавались программно от 0% до 30% для временного ряда, состоящего из 24 точек. На рисунках 2 и 3 представлены графики рассчитанных характеристик точности: для моделирования — коэффициента детерминации Я а для прогнозирования (с горизонтом прогноза в 20% от длины исходного ряда) — второго коэффициента Тейла Т2, вычисляемых по формулам:

№ 1 (25) 2012

-у)2

R2 = k=

■ф

£ о о о S п

\о

5

о ф

Í

6

S ф

■J £

л §

о

S

5

о

о ф

6 ф

а о

>s

ф

S §

ф

Щ

о Í о

0

W í 5 S

1

w

Sí

s

a-

2

0 t t ф t n

1

>s

Ф

s §

§

as ■fr

¡5

I

ф

=1__t = k=l_

N ' T 2 l l

s (y -Y )2 s Yk

k=1 | k=1 k=1

-YM0d)2

100%,

где Y^od — модельные значения ряда динамики.

Выбор принятых простых критериев точности моделирования и прогнозирования объясняется необходимостью частого использования малых выборок при моделировании жизненного цикла. При наличии больших выборок можно обращаться к более точным (и, естественно, более сложным) критериям, например, (Diebold, Mariano, 1995).

Таблица 2. Диапазон параметров модельных задач для исследования устойчивости алгоритмов идентификации моделей Mj - М4 к шуму

Q А «1 А B1 Wj С2 А «2 л2 B2 w2

M,

min -15 -0.3 3 -0.10

max -3 -0.1 15 -0.03

M3

min -15 -1.0 -0.3 3 0.2 -0.10

max -3 -0.2 -0.1 15 1.0 -0.03

M4

min -15 -0.3 0.01 0.01 0.3 3 -0.10 0.01 0.01 0.3

max -3 -0.1 0.10 0.10 1.0 15 -0.03 0.10 0.10 1.0

М2

D, А А А А1 «1

min 10 -10 10 1 50 -0.6

max 50 -1 50 5 150 -0.4

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

алгоритмов идентификации моделей М5 - М7 к шуму

А А А с, «1 Л B1 w2 л2 B2

М5

min 20 -0.10 1.2 0.1 0.1 0.6 10 10

max 50 -0.01 1.5 0.7 0.7 0.9 30 30

М6

min 20 5 -0.10 1.2 10 10 0.6 10 10

max 50 10 -0.01 1.5 30 30 0.9 30 30

М7

min 20 5 0.1 -0.10 1.2 10 10 0.6 10 10

max 50 10 1.0 -0.01 1.5 30 30 0.9 30 30

130

№ 1(25) 2012

21

И

0.95 0.9 0.85 0.8

V

\ \ _

м3>ч м4

ч

7"2,%

10

20

30

30252015105

0?

¿..•Л»"*'

• м

* М1

// м2

//

//

10

20

30

Рис. 2. Коэффициент детерминации Я2 и второй коэффициент Тейла Т2 в случае идентификации моделей М1 - М4 при различных сочетаниях дисперсий шума и сигнала Км

КН!8'% 0

КЛ//3'%

Рис. 3. Коэффициент детерминации Я2 и второй коэффициент Тейла Т2 в случае идентификации моделей М5 - М7 при различных сочетаниях дисперсий шума и сигнала К

■т

131

о ф г

и 'ф

§ ф

0

Щ

ни

1 I

5;

иц

(о Ф г

и 'ф

§ ф

о

03

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

Всего для построения рис. 2 и 3 было рассчитано 17500 реализаций модельных задач, параметры которых выбирались из диапазонов, указанных в табл. 2, 3. Модель М2 исследована для случая многочлена 3-ей степени. Представленные на рис. 2, 3 значения коэффициентов Я , Т2 говорят о высокой точности моделирования и прогнозирования моделей М1, М2, М4, М6, М7 при уровнях шума < 10% (коэффициент детерминации > 0.95, а коэффициент Тейла <22%) и о довольно высокой точности моделирования и прогнозирования моделей М1 -М3, М5, М6 при уровнях шума < 30% (коэффициент детерминации > 0.85, а коэффициент Тейла < 30%). Точность идентификации моделей М3, М5 и прогнозирования с помощью М5 оказалась несколько ниже, чем точность других моделей. Таким образом, область возможного применения предложенных моделей и методов их идентификации оказывается достаточно широкой.

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

№ 1 (25) 2012

■ф

е-

со о о S

п

5 о

L. Ф

t

6

2 ф ■J s

л §

о

5

5

о о

ф

6

ф а о

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

Проведенное исследование на основе более 2000 генераций модельных задач показало, что наиболее точно определяются оценки значений временных рядов: коэффициент детерминации Я и оценка прогноза Т2. Вариация Я в большинстве случаев составляет 1-2% и не превышает 6%, а вариация коэффициента Тейла — менее 10%.

Вариации отдельных параметров моделей оказались несколько выше. Так, показатели экспонент а1, а2 и частоты <, <2 имеют вариацию 10 - 15%, а амплитуды колебательных компонент С1, С2, Л1, Л2, В1, В2 — до 24%. Невозможность использования длинных реальных выборок при этом «заменяется» усреднением по множеству модельных выборок.

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

>s

ф §

§

ф о о t о

0

W t S S

1

«I §

S

а

S

0

1

I

ф

I

I

>s

ф

§ § §

а*

s ■&■

Й

I

ф

Коэффициенты моделей М1, М3, М4 идентификации жизненного цикла автомобилей Мустанг производства «Форд» (Котлер, 2007), а также коэффициенты Я и Т2 представлены в табл. 4, а результаты построения моделей — на рис. 4. Оказалось, что наилучшее приближение жизненного цикла продаж автомобилей Мустанг и одновременно наименьшую ошибку прогноза на 9-ый год дает модель М4.

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

Q А «1 A1 B1 w1 С2 D2 «2 A2 B2 w2 R2 A %

M, -1341 -0.8341 1604 -0.3734 0.979 22.10

Мз -688 -97 -1.1133 953 280 -0.4619 0.979 22.09

М4 -1130 -0.8972 0.085 -0.126 1.501 1394 -0.3482 -0.020 -0.100 1.979 1.000 10.39

600 400 200

Sales

/ \ M usta ng

/ \ 4 N

M, Л

600.

Sales

400

200

) s M usts ng

/ 1 4 N v

M, N

600

400

200

Sales

/ -- \ M usta rng

/ > 4

M N

123456789

123456789

123456789

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

черная тонкая линия — прогноз)

132

0

0

(25) 2012

Для жизненного цикла добычи нефти из пласта Турнейского яруса использование АЯМА-моделей для М2 с полиномом третьей степени (Семёнычев, Куркин, 2010) позволяет получить существенно более точные результаты по сравнению с (Оразбаев и др., 2007) — см. рис. 5 и табл. 5.

300

200

100

a)

300

200

100

б)

10

20

30

Рис. 5. Моделирование уровня годовой добычи нефти моделями М2 с использованием многочленов: а) 1-ой и 2-ой степени; б) 3-ей (пунктир), 4-ой и 5-ой степени

Таблица 5. Идентификация жизненного цикла месторождения моделями М2

Степень модели Dj d2 D3 D4 D5 D6 aj Aj R2

1 -250.4 135.2 -0.2224 134.5 0.705

2 -113.8 -20.6 28.4 -0.2521 101.0 0.909

3 -123.1 -1.0 -53.6 20.1 -0.4313 125.6 0.947

4 -125.6 -94.7 42.9 20.4 5.9 -0.5369 129.4 0.943

5 -134.6 487.3 -1144.8 682.6 -165.4 15.6 -0.7812 138.1 0.922

0

ф

г

1 'ф

§ ф

0

CQ Ui

1 I

s: ui

(о Ф г

I 'ф

§

ф О

оа

0

Параметры модели М2 и тренда (13), с учетом использования ARMA-моделей (8) и (17), а также параметры трендов (14), (15) можно найти в рамках системы линейных уравнений и решения алгебраического уравнения одной переменной.

Рассмотрим теперь ситуацию, когда жизненный цикл продукта отражен не только данными реальных продаж, но и косвенно — динамикой интереса пользователей системы Google к сотовым телефонам марок Nokia N70, N73, N82 и N95. Для этого используем данные о запросах пользователей, которые еженедельно агрегируются в статистике Google Trends (2010).

Смоделируем с помощью Mj, М3 и М4 ежемесячную динамику индекса запросов (табл. 6 и рис. 6). Усредняя коэффициенты детерминации и Тейла по всем маркам телефонов, получим, что Mj описывает индекс запросов с R2 = 0.88, T2 = 12.3%, М3 — с R2 = 0.95, T2 = 14.6%, М4 — с R2 = 0.97, T2 = 5.3%. Наилучшее приближение, как видим, дает модель М4.

Особенностью рынка путеводителей является сезонность, поскольку летом люди путешествуют чаще. Для моделирования интереса пользователей системы Google к путеводителям как в целом (запрос «Travel guide»), так и к наиболее распространенному изданию Lonely Planet применены модели М5 - М7 (табл. 7 и рис. 7). На рисунках серый цвет соответствует исходным данным, черная толстая линия — идентифицированной модели, черная тонкая линия — прогнозу. Худшую точность дала модель М6, модели М5 и М7 дали близкую точность, причем у М7 она несколько лучше.

№ 1 (25) 2012

Таблица 6. Идентификация интереса пользователей Google к моделям телефонов Nokia

■ф

е-

со о о s

п

\о s

о ф

t

&

s ф

s

л

4 §

о s

5

о

о ф

е-

ф а о

>s

ф

s

s

ф eg О Ï о

0

W

1 5 S

I s

W

Si

s

aS

0 ï ï ф ï

1

>s

ф

s

s §

a*

s ■fr

¡5 i ф

Ci Di ai А, в, w, C2 D2 a2 A2 B2 w2 R2 T2,%

N70

М, -12.68 -0.0643 13.11 -0.0517 0.920 8.98

М, 1.25 0.59 -0.3620 -0.91 0.25 -0.0638 0.961 12.69

М4 -3.59 -0.0813 -0.089 -0.078 0.2736 4.06 -0.0402 0.012 -0.080 0.1166 0.975 7.36

N73

М, -6.81 -0.0896 6.97 -0.0463 0.945 2.12

М, -3.11 0.04 0.0121 3.39 0.35 -0.0398 0.966 16.28

М4 -5.76 -0.0957 -0.006 -0.101 0.1531 5.92 -0.0439 -0.006 -0.082 0.1342 0.965 3.98

N82

М, -2.93 -0.1178 2.82 -0.0735 0.877 7.30

М, -0.80 -0.35 -0.3229 0.82 0.09 -0.0963 0.950 12.90

М4 -4.25 -0.1130 0.006 -0.040 0.9326 4.15 -0.0817-0.005 -0.015 0.2379 0.985 4.35

N95

М, -6.55 -0.0909 6.71 -0.0345 0.787 30.94

М, -12.69 -1.09 -0.1365 13.18 -0.20 -0.0317 0.922 16.42

М4 -2.73 -0.1374-0.033 -0.295 0.7671 3.07 -0.0195 0.270 -0.194 0.7671 0.950 5.66

N 70

M 3

1 2 3 4 i

1 2 3 4 i Search index

2 1.5 1

0.5 0

0.5!

0 L S

1.5 -01

f r, V N 70

r Л S s

4 s 'S m

1 2 3 4 i Search index

/ J N7

r V <

f V «

1 2 3 4t ' 1 2 3 4t ' 1 2 3 4t

Search index_ Search index_ Search index

N82

II Л

N82

M3

N82

I Л4

1 2 Search index

3 t 1 2

Search index

-С N 95

/ s ■s.

s

3 t 1 2

Search index

* h N 95

A Л S »V

r M3 >

h * 14 N 95

ï Л Ч »V

Л 4 Ь.

12 3t 123t 12 3t

Рис. 6. Идентификация интереса пользователей Google к моделям телефонов Nokia (серый цвет — исходные данные, черная толстая линия — идентифицированная модель,

черная тонкая линия — прогноз)

134

1.5

1.5

1.5

0.5

0.5

0

0

0.4

0.4

0.2

0.2

0

3 t

3

3

№ 1(25) 2012

Таблица 7. Идентификация интереса пользователей Google к путеводителям

А А А Q «1 w1 A1 Bi w1 a2 B2 R2 T2, %

Travel Guide

М5 1.67 -0.1343 0.488 0.092 -0.073 0.144 0.017 0.040 0.954 4.28

М6 1.547 -0.013 -0.0525 0.033 -0.612 0.284 0.483 0.099 -0.071 0.951 11.22

М7 1.694 -0.024 0.00013 -0.0283 0.480 0.243 -0.128 0.963 0.010 0.062 0.965 6.80

Lonely Planet

М5 1.73 -0.1471 0.163 0.026 -0.042 0.490 0.041 -0.075 0.935 13.98

М6 1.513 -0.013 -0.0237 0.099 0.176 0.112 0.477 0.068 -0.051 0.936 17.50

М7 1.738 -0.025 0.00014 -0.0207 0.486 0.092 -0.129 0.963 0.021 0.074 0.953 6.50

Search index

Search index

1.5

1

0.5

2 1.5 1

0.5 0

1 2 3 Search index

4 5 6 t

Lonely Planet

M 15

1 2 3 4 5 6 t

2 1.5 1

0.5 0

2 1.5 1

0.5 0

1 2 3 4 5 6 Search index

2 1.5

1

0.5

Search index

Щ

M7

travel guide

V

1 2 3 4 5 6 t Search index

1 2 3 4 5 6 t

1.5 1

0.5

1 2 3 4 5 6 t

135

О

ф

г

Ü 'ф

S ф

0

CQ ai

1 I

s;

ai

to ф

г

Ü 'ф

S ф О

ü<

оа

Рис. 7. Идентификация интереса пользователей Google к путеводителям

Точность оценки параметров моделей, представленных в табл. 4 -7, может быть оценена величиной их вариации. Исходные данные, лежащие в основе построения моделей, измерены один раз, и величина содержащейся в них (в результате влияния случайных процессов) стохастической компоненты не определяется из условий наблюдений. Дисперсия стохастической компоненты исходных данных оценивается путем сравнения исходных временных рядов с построенными моделями. Соотношение дисперсий стохастической компоненты и сигнала KN/S в рассмотренных случаях лежит в диапазоне 1-14%, за исключением описания интереса пользователей системы Google к телефону Nokia N95 моделью Mj, где уровень KN/S = 42%. Отметим, что при описании этого же случая моделью M3 уровень

KN/S = 7%.

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

Наиболее точно определяются в рассмотренных случаях оценки модельных значений временных рядов: коэффициент детерминации R2 и оценка прогноза Тейла T2. Вариация R2

2

2

№ 1 (25) 2012

в большинстве случаев составляет 1-2% и не превышает 6%, а вариация коэффициента Тей-ла не превышает 10%. Близкие значения временного ряда Ук могут быть описаны существенно различающимся набором параметров моделей. Поэтому вариации параметров перечисленных выше моделей — показатели экспонент а1, а2 и частоты , «2 — имеют вариацию а 10 -15%, а параметры амплитуд и фаз колебательных с-, с 2, А, , в-, в2

<а риации до 24%.

8. Выводы

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

На тестовых и реальных выборках показана целесообразность использования для Ц их идентификации подхода на основе обобщенных параметрических АЯМА-моделей, методов параметрической итеративной идентификации, а для сложных моделей — базисов

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

1« Грёбнера. Для многих тестовых и восьми реальных выборок (из трех отраслей экономи-

I

0 возможность применения предложенных моделей для широких диапазонов значений па-¡1 раметров и сочетаний дисперсий шума и сигнала. Для каждого из реальных примеров уда-а лось подобрать модель с коэффициентами детерминации Я2 > 0.947 и оценкой прогноза

1 Тейла Г2 < 10.4%.

ГС 2

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

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

5

Айвазян С. А. (2001). Прикладная статистика. Основы эконометрики. М.: ЮНИТИ-ДАНА. | Бессонов В. А. (2005). Проблемы анализа российской макроэкономической динамики переходно-о го периода. М.: ИЭПП.

§ Ершов А. Г. (2000). Использование алгоритма построения базисов Грёбнера в рамках концепции программирования в ограничениях. В кн.: Тезисы докладов Четвертого сибирского конгресса

| по прикладной и индустриальной математике (ИНПРИМ-2000). Новосибирск: Институт математики СО РАН. IV, 105-106.

с;

¡^ Ершов А. Г. (2007). Алгоритмы и программные средства для геометрических задач параметри-

О

§ ческого проектирования. Диссертация на соискание ученой степени кандидата физ.-мат. наук. Но-§

¡5 Мазуркин П. М. (2002). Статистическая теория реальных циклов. В кн.: Материалы Четвертой

® Международной конференции «Циклы». Ч. 1. Ставрополь: СевКавГТУ 140 -155. §

восибирск.

Котлер Ф. (2007). Основы маркетинга. Краткий курс. М.: Издательский дом «Вильяме».

№ 1(25) 2012

Оразбаев Б. Б., Курмангазиева Л. Т., Кабылхамит Ж. М. (2007). Задачи прогноза и идентификации щ

нефтедобычи, математические методы и алгоритмы их решения. Электронная библиотека Атырау- |

ского института нефти и газа. http://www.aing.isd.kz/e-lib/okk.pdf. <ф

Семенов С. П., Славский В. В., Татаринцев П. Б. (2004). Системы компьютерной математики. о

Учебное пособие для студентов математического факультета АГУ. Барнаул: Изд-во Алт. ун-та. ®

Семёнычев В. К., Семёнычев Е. В. (2011). Параметрическая идентификация рядов динамики: t

§

структуры, модели, эволюция. Самара: Изд-во СамНЦ РАН. а

Семёнычев В. К., Семёнычев Е. В., Коробецкая А. А. (2010). Исследование точности метода мо- ^ делирования и прогнозирования экспоненциальной тенденции на основе обобщенных параметриче- uj ских ARMA-моделей. Вестник Самарского муниципального института управления, 1, 36-45. §

Семёнычев В. К., Куркин Е. И. (2010). ARMA-моделирование уровня годовой добычи нефти | из пласта и оценка геологического риска инвестиций в нефтегазодобывающей промышленности. § Вестник Самарского муниципального института управления, 2, 7-14. ^

Соколов В. А. (2006). Эволюционные уравнения как феноменологическая модель разработ- ® ки нефтегазовых месторождений. Электронный научный журнал «Нефтегазовое дело», 2/2006. http://www.ogbus.ru/authors/SokolovVA/SokolovVA_1.pdf.

Четыркин Е. М. (1977). Статистические методы прогнозирования. М.: Статистика.

Aman H., Ohkochi T. (2010). An application of growth curve model for predicting core churn in open source development. In: Proceedings of the Joint Conference on Knowledge-Based Software Engineering (JCKBSE'10). A. Caplinskas. H. Pranevicius, T. Nakatani (Eds.). Kaunas: Technologija. 32-38. http://www.hpc.cs.ehime-u.ac.jp/~aman/pdf/jckbse2010.pdf.

Beveridge S., Nelson C. (1981). A new approach to decomposition of economic time series into permanent and transitory components with particular attention to measurement of the business cycle. Journal of Monetary Econometrics, 7, 151-174.

Brody S. (1945). Bioenergetics and growth. New York: Reinhold publishing.

Diebold F., Mariano R. (1995). Comparing predictive accuracy. Journal of Business and Economic Statistics, 13, 253 - 263.

Enders W. (2004). Applied econometric time series. 2nd edition. John Wiley & Sons.

Google Trends. (2010). Статистика запросов пользователей системы Google. http://trends.google.com.

Grubler A. (1990). The rise and fall of infrastructures. Dynamics of evolution and technological change in transport. Physica-Ver^.

Hook M., Li J., Oba N., Snowden S. (2011). Discriptive and predictive growth curves in energy system analysis. Natural Resources Research, 20 (2), 103 - 116.

Hodrick R., Prescott E. (1997). Postwar business cycles. Journal of Money, Credit and Banking, 19, 1-16.

137