ЭКОНОМИКА И МАТЕМАТИЧЕСКИЕ МЕТОДЫ, 2012, том 48, № 2, с. 108-115
СТАТИСТИЧЕСКИЕ МЕТОДЫ И ТЕОРИЯ ВЕРОЯТНОСТЕЙ
МЕТОДЫ ИДЕНТИФИКАЦИИ ЛОГИСТИЧЕСКОЙ ДИНАМИКИ И ЖИЗНЕННОГО ЦИКЛА ПРОДУКТА МОДЕЛЬЮ ВЕРХУЛСТА
© 2012 г. В.К. Семенычев, В.Н. Кожухова, Е.В. Семенычев
(Самара)
На тестовых выборках различного объема, в широких диапазонах сочетаний параметров модели и мощности помехи сравнивается точность известных аналитических методов идентификации модели логистической динамики Верхулста (Перла-Рида). Показано преимущество применения метода Левенберга-Марквардта для реализации метода наименьших квадратов, возможность его использования при моделировании жизненного цикла продукта типа "фетиш".
Ключевые слова: модель Верхулста, моделирование, прогнозирование, методы идентификации, точность, метод Левенберга-Марквардта, модель жизненного цикла продукта типа "фетиш".
Логистические тренды находят широкое применение в экономике при моделировании цен на товары, уровней спроса, емкости рынков, демографических процессов и т.д. Динамика подобных показателей является эволюционирующей: вначале следует стадия медленного роста, который затем ускоряется и снова замедляется в период зрелости, достигая затем насыщения.
При моделировании жизненного цикла продуктов обычно выделяют четыре стадии: внедрение на рынок, рост, насыщение и упадок. Одной из наиболее трудных для идентификации моделей жизненного цикла является модель "фетиш", в которой продукт (например, новогодние товары, компьютерные игры и др.) очень быстро достигает пика популярности (продаж) и так же быстро переходит к стадии упадка. Кумулятивную (накопленную) динамику уровней таких циклов (например, продаж или спроса) отражает логистическая модель и ее обобщения.
В настоящее время известно более десяти моделей логистической динамики (Семёнычев В., Семёнычев Е., 2006, гл. 5) уровней показателя Уь из которых наиболее распространенной можно считать нелинейную по параметрам модель Верхулста (Перла-Рида) для тренда и с аддитивной (ее часто называют основной) структурой стохастической компоненты (помехи) £к:
А 0
7к =-^ + £к, (!)
1+ А1 е
где А0, А1, а - параметры логистического тренда, к = 1, ..., N - номер наблюдения, N - объем выборки, А - период опроса (дискретизации).
В экономической практике выполнение условий Гаусса-Маркова (центрированность, го-москедастичность, отсутствие автокорреляции, нормальный закон распределения) для £к в (1) считается справедливым и достаточно оправданным требованием, которое позволяет получать оптимальные (несмещенные, эффективные и состоятельные) оценки параметров тренда (Боро-дич, 2001, гл. 5) с помощью метода наименьших квадратов (МНК).
Все известные методы идентификации модели Верхулста основываются на использовании тех или иных преобразований модели тренда, приводящих к моделям, линейным по параметрам. При этом для стохастической компоненты принимают "удобные" (для применения МНК) предположения о ее свойствах и/или о месте ее вхождения в структуру модели.
Статья посвящена актуальной, не решенной до настоящего времени задаче сравнения точности моделирования и прогнозирования известных методов идентификации модели (1) в широком диапазоне сочетаний параметров, при различной мощности помехи, а также предложению по повышению точности идентификации.
Наиболее распространенный в настоящее время метод идентификации модели (1) заключается в переходе к обратным значениям уровней ряда. При этом предполагаются "удобные" для реализации метода свойства модели: стохастическая компонента £к находится в знаменателе модели тренда, она мультипликативна по отношению к экспоненциальной функции и имеет логарифмически нормальное распределение. Кроме того, уровень насыщения А0 считается априорно известным (обычно он представляет самостоятельный интерес для определения). Такую преобразованную модель ряда можно логарифмировать, а затем применить МНК к полученной парной линейной регрессии с параметрами 1п А1, а и нормальным законом распределения помехи 1п £к (Мхитарян, 2008, гл. 7):
Yk =
A о
1 + A1 e
-аМ,
^-1= A i e Yk 1
-аакЬ,
ln(A 0/Yk -1) = In A1 - akD +ln sk.
Добавим, что МНК минимизирует в этом случае другую среднеквадратическую невязку - на логарифмах значений 1п(А0 /Тк - 1) и 1пА1.
Принятие указанных свойств можно считать недостатками данного метода идентификации, ограничивающими область его применения. Поэтому перейдем к обзору и сравнению других девяти методов идентификации модели (1).
В известном методе трех сумм исходный ряд динамики разбивается на три равных отрезка, затем вычисляются суммы значений ряда внутри каждого отрезка и определяются разности этих сумм (Четыркин, 1977, гл. 4). Параметры модели оцениваются по следующим формулам:
1
1
1
а0 = — | ln
m \ U Yk R Yk/
-lnj
1
1
R2 Yk R3 Yk ,
A0 = m
A? =
1 l(1- m e
где т - длина отрезка;
Методы Фишера и Готеллинга предполагают исследование дифференциального уравнения логистической функции dY/dk = Та - У2а/Л0. В обоих методах используются приближения для расчета производной.
В методе Фишера (Четыркин, 1977, гл. 4) приближенно вычисляются темпы прироста:
^ -1. 0,51п Тк+1, а0, В0 = а^штЦ(0,51п(Тк+1/Тк-1)-а +ВТк_х)\ А0 = а°/В°. dJ Тк Тк-1 - ^
а, B
k = 2
В методе Готеллинга (Суслов, Ибрагимов, Талышева, 2005, гл. III, разд. 11) производная рассчитывается как разность текущего и предшествующего значений показателя Yk: dY/dk-- DY/Dk = Yk - Yk-1.
Данные методы, как и приведенные ниже, предполагают аддитивное "включение" в модель стохастической компоненты, отвечающей условиям Гаусса-Маркова, уже после проведения линеаризующих преобразований.
k
В частности, для метода Готеллинга будем иметь AYk = Yk-1 a - Y2k-1 а/А0 + hk, где hk - "новая" (не равная £к в (1)) стохастическая компонента.
Оценки параметров a, B находятся с помощью МНК:
N
ao, Bo = argmin/ (Yk -Yk_, - Yk-X a + BY^)2, A0 = ao/Bo.
ao, Bo к = 2
Метод Юла сводит задачу идентификации параметров модели (1) к идентификации параметров регрессии темпов прироста (Yk+1 - Yk) / Yk на Yk+1:
(Yk+1 - Yk)/Yk = (e-aD -1) - Yk+1(e-aD -1) /Ao + ]k ,
где ]k - "новая" стохастическая компонента.
Нелинейно входящие в уравнение параметры заменяют линейными e-aD - 1 = Q, (e-aD - 1)/A0 = P и применяют МНК:
N-1
Qo, Po = argmin/ ((Yk+1 -Yk)/Yk - Q + PYk+1)2.
Q, p
k = 1
Метод Родса (Четыркин, 1977, гл. 4) основан на взятии разности между соседними обратными значениями ряда, при этом идентификация осуществляется относительно параметров
(1 - е-аА)/А0 и е-аА:
1 -(1-е )1г-аА ± + Р
е + рЪ
Ук+1 А 0 Ук где рк - "новая" стохастическая компонента.
Метод Нейра использует регрессию разности соседних обратных значений ряда на их сумму. Решение осуществляется относительно параметров (е-аА - 1)/(е-аА + 1) и [(е-аА - 1)/(е-аА + 1)] х
X
(2/Ао):
1 1 e-аД - 1 2 e-tA-1 J 1 ML, — —--■# — - —-# ITT" + — I+}k,
Yk+1 Yk e + 1 A о e-aA + 1 \Yk+1 Yk,
где }k - "новая" стохастическая компонента.
Все перечисленные методы предполагают последовательное вычисление оценок параметров модели: на первом этапе производится расчет оценок параметров a°, Ao0, на втором - оценки параметра A1
В работах (Четыркин, 1977, гл. 4; Суслов, Ибрагимов, Талышева, 2005, гл. III, разд. 11) предлагается оценивать значение параметра A1 методом моментов: А0/Yk - 1 = A1e-akA, ln A1 =
1 N 1 N
= akA + ln(A0 /Yk - 1), а затем находить средние значения — / lnA 1=— / akA +
1 N N N k=1 N k=1
+--/ ln(A0/Yk - 1). Учитывая, что / k = 0,5N(N +1), получаем оценку
N k = 1 k = 1
A1 = exp j-N f0,5a°AN(N +1) + / ln(A0/Yk -1)jj.
Однако при наличии в выборке значений, превышающих найденную оценку уровня насыщения логистической кривой A0, метод моментов становится неработоспособным, поскольку возникает отрицательное число под знаком логарифма. Казалось бы, оценка параметра A° может быть найдена с помощью перехода к обратным значениям уровней ряда и МНК:
N I 1 1 e-a°M\2
A1 = arg min / ( — --- A1 -I .
1 1 fri\ Yk A° 1 A0 I
Однако нужно учесть, что в этом случае получим гетероскедастическую стохастическую компоненту:
„ __А о , _ А о + в, (1+ А1 е )
¿-к —
■■ 1+ A 1 e-akD * 1+ A1 e^ 1+ A1 e -а*д 1 1+ A1 e -а*д £k (1+ A1 e -akD )
Y* Ao + £k(1+ A1 e-акд) Y* A0 Фк' {k' A0 Y*
Для компенсации гетероскедастичности - уменьшения неэффективности оценок параметров модели (Бородич, 2001, гл. 5) - можно рассмотреть возможность развития в этом направлении предложенного в работе (Семёнычев В., Семёнычев Е., 2006, гл. 5) метода обобщенных параметрических моделей авторегрессии - скользящего среднего (ARMA-моделей) для модели (1).
Для реализации этого метода осуществим замену переменных модели A0 = 1/A00, A1 = A 10 /A00, а затем с помощью прямого и обратного Z-преобразования сконструируем разностную схему (для к > 3):
D* = À(Dk-1 - Dk-2) + Dk-1, (2)
где À = е-аД, Dk - детерминированная часть модели, Dk = A00 + A10e-akD.
Первый способ реализации этого метода (ARMA I) учитывает соотношение Dk = 1/Yk - фк и приводит к следующей модели авторегрессии:
1/Yt = G* = À(Gk-1 - Gk-2) + Gk-1 + 6*, 6* - Фк-1 + À({*-2 - Фк-1), где 6* - гетероскедастическая стохастическая компонента.
Оценка параметра À (и с учетом обозначений в (2) параметра а) находится с помощью взвешенного МНК для компенсации гетероскедастичности. В качестве оценок весов wk можно использовать:
1) уровни ряда Gk, если помехи малы по сравнению с уровнями ряда;
2) МНК в два этапа - на первом этапе модель идентифицируется с помощью обычного МНК, на втором этапе полученные оценки D* применяются в качестве весов;
3) один из методов непараметрического сглаживания: вначале оцениваются значения 6* = Y* - Yk, где Yk* - полученные сглаженные значения ряда, а затем по соседним m точкам строится ряд оценок дисперсии 6k, при этом теряются m - 1 значения:
m -1 ! ^ m -1 \ 2
Sk = -Г" / Uk - i / {k-i I .
m -1 i = 0 \ mi = 0 I
Полученные оценки дисперсии используются в качестве весов wk:
N 1
À0 = argmin/ — (Gk - À(GM - G^) - Gk-1)2.
— 2
À k=3 Wk
Второй способ (ARMA II) учитывает соотношение Dk = 1/(Yk - £k):
(Yk-1 - Yk)Yk-2 - À(Yk-2 - Yk-1)Yk =
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.