Модель прогнозирования объектов с цикличностью развития

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

Если считать, что временной ряд состоит из 3-х аддитивных составляющих

, (5.1)

где g_t – тренд;

s_t – циклические колебания;

e_t – случайная составляющая

и обозначить операцию вычисления скользящего среднего символом Т, то [23]

и

(при условии идеального выделения тренда ).

Это означает, что члены и могут исказить истинные значения остаточного ряда, в том числе индуцировать ложное колебательное движение (эффект Слуцкого-Юла).

Так, например, если , то операция вычисления простого скользящего среднего (с равными весами) k последовательных членов приведет к синосуидальному ряду с тем же периодом и фазой, но с амплитудой, равной [23]

.

Если выбрать кратным (длина отрезка усреднения равна или кратна периоду колебания), то и искажение отсутствует.

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

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

Однако полученные при этом данные относятся к интервалам с серединами между кварталами, а не к серединам интервалов (15 февраля, 15 мая, 15 августа, 15 ноября). Поэтому необходимо «центрировать» среднее. Наиболее просто это можно сделать, взяв средние последовательных пар средних, вычисленных по 4-м элементам. Такой процесс эквивалентен вычислению средних пяти элементов с весами [ ].

Если средний индекс сезонности за 12 месяцев или 4 квартала не равен единице, производится выравнивание индексов сезонности – деление всех индексов на их средний индекс.

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

, (t=1,…,n ; q=1,…,4), (5.2)

где s_q – индекс сезонности.

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

, (5.3)

где i – номер гармоники;

t – номер точки временного ряда;

– амплитуды гармоник.

Следует отметить, что гармоники обладают свойством ортогональности

( ),

( ), (5.4)

Это позволяет получить оценки параметров ( a₀, и соответственно) по методу наименьших квадратов

Если n – четное, то по этим формулам рассчитываются коэффициенты для гармоник. Для последней гармоники вычисляется

и в этом случае модель сезонной волны принимает вид

(5.5)

Если считать что периоды s_t задаются числами n/k_j, j=1,…,q, где (k₁,…,k_q) – подмножество последовательности целых чисел 1, 2, … ,(n-1)/2, если n – нечетное, или 1, 2, … n/2-1, если n – четное, то приведенные выше формулы примут вид

(5.6)

(для нечетных n)

(5.7)

(для четных n).

Иногда сезонность может быть выражена настолько явно, что нет необходимости доказывать ее существование. Однако в общем случае необходим статистический тест проверки гипотезы об отсутствии циклического слагаемого с заданным периодом n/k_j.

Выдвинем нулевую гипотезу

H₀:

или

При нулевой гипотезе статистика

(5.8)

имеет распределение Фишера с 2 и n-p степенями свободы.

Здесь p – число оцениваемых коэффициентов (2q+1 или 2q+2);

;

S² – оценка дисперсии случайной составляющей модели, которая определяется по формулам

(при нечетных n)

или

(при четных n).

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

Равномерно наиболее мощный критерий проверки гипотезы Н₀ (все ) состоит в принятии гипотезы Н₀, если

( j=1,…,q). (5.9)

При этом константа g выбирается таким образом, чтобы вероятность ошибки первого рода, которая определяется по формуле

(5.10)

где не превосходила заданную.

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

( . При этом

.

Имеет место приближенное соотношение

.

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

С учетом вышесказанного предложенная нами математическая модель, описывающая процессы такого характера, имеет вид [15]

k

t - t₀- S T_i

i=1

y_t = g(t) + a(t) Cos ( ------------------ 2p ) + e_t , (5.11)

T_k+1

где g(t) – основная тенденция (тренд);

a(t) – амплитуда циклических колебаний;

t₀– момент первого пика;

T_i– промежуток времени между i-ым и i+1 пиками;

k+1 – количество пиков, наблюдавшихся до момента времени t;

e_t – случайная составляющая с нулевым математическим

ожиданием.

При этом в качестве моделей g(t) и a(t) будем использовать класс полиномов

g(t) = g₀+ g₁t + g₂t²+ ... + g_ntⁿ, (5.12)

a(t) = a₀+ a₁t + a₂t²+ ... + a_mt^m.

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

считать, что в интервале между первым и последним наблюдавшимися пиками значения T_i( i = 1,…, k) известны так же, как и значение t₀. Поэтому для всех точек временного ряда, принадлежащих этому интервалу, величина

k

t - t₀- S T_i

i=1

cos ( ------------------ 2p )

T_k+1

может быть рассчитана, что позволяет получить одновременно оценки параметров модели ₀, ₁ , ₂ , ... , _n, ₀, ₁, ₂ , ... , _m и оценку дисперсии случайной составляющей по методу наименьших квадратов. Кроме того, МНК позволяет получить оценки их дисперсий и ковариаций: { _i} (i=0,…,n) ; { _i} (i=0,…,m) ; Cov { _i, _j} ( i,j = 0,…,n) ;

Cov { _i, _j} ( i,j = 0,…,m) ; Cov { _i, _j} ( i=0,…,n ; j = 0,..,m) .

Формула для краткосрочного прогноза показателя (на период до ближайшего нового пика) с использованием модели (5.11) будет иметь вид

k

t_пр- t₀- S T_i

i=1

(t_пр) = (t_пр) + (t_пр) cos ( -------------------------- 2p ) =

_k+1

= ₀+ ₁t_пр+ ₂t_пр²+ ... + _nt_прⁿ + (5.13)

k

t_пр- t₀- S T_i

i=1

+( ₀+ ₁t_пр+ ₂t_пр²+ ... + _mt_пр^m) Cos ( -------------------------------- 2p ) ,

_k+1

где t_пр– прогнозный момент времени;

_k+1 – оценка длительности периода между последним наблюдавшимся и вновь ожидаемым пиком.

В качестве оценки _k+1можно использовать среднее значение длительности промежутков между пиками по данным предыдущей статистики

k

_k+1= = S T_i / k . (5.14)

i=1

При этом оценка дисперсии случайной величины _k+1

1 k

{ _k+1} = ------- S (T_i- )² , (5.15)

k (k-1) i=1

а оценка дисперсии ошибки прогноза показателя может быть рассчитана по формуле [7]

A(t_пр)

{ (t_пр) - y(t_пр)} = + { (t_пр)} + cos²( ------------- ) { (t_пр)}+

_k+1

A(t_пр)

+ ( (t_пр) sin ( ------------ ) A(t_пр) ^-2_k+1 )². (5.16)

_k+1

A(t_пр)

{ _k+1 - T_k+1} + 2 cos ( -------------- ) Cov { (t_пр), (t_пр) } ,

_k+1

где

k

A(t_пр) = ( t_пр- t₀- S T_i) 2p ;

i=1

k+1 k

{ _k+1 - T_k+1 } = ------------------ S (T_i- )² ;

k (k-1) i=1

n n-1 n

{ (t_пр) } = S t_пр²ⁱ { _i} + 2 S S t_пр^i+jCov ( _i, _j) ;

i=0 i=0 j=i+1

m m-1 m

{ (t_пр) } = S t_пр²ⁱD { _i} + 2 S S t_пр^i+jCov ( _i, _j);

i=0 i=0 j=i+1

n m

Cov { (t_пр) , (t_пр) } = S S t_пр^i+jCov ( _i, _j) .

i=0 j=0

Пример. Данные о динамике валового национального продукта (ВНП) США за 1960–1985 гг. приведены в табл.10 и на рис.5.

Т а б л и ц а 10. ВНП США (трл. $)

Год	ВНП	Год	ВНП	Год	ВНП
	1,665		2,423		3,115
	1,708		2,416		3,192
	1,799		2,484		3,187
	1,873		2,608		3,248
	1,973		2,744		3,166
	2,087		2,729		3,279
	2,208		2,695		3,501
	2,271		2,826		3,318
	2,365		2,958

Из этих данных следует, что наблюдаемые пики в динамике ВНП США имели место в 1969, 1973, 1979 и 1984 гг., а периоды циклов составили соответственно 4, 6 и 5 лет (выброс 1981 г. по всей видимости, является чисто случайным).

Визуальный анализ рис. 5 говорит о линейности тренда, поэтому нами были апробированы линейные модели

,

.

Однако коэффициент ₁ оказался статистически незначимым (при a=0,05 и 17-5=12 cтепенях свободы). Поэтому нами была изменена модель, а именно исключена зависимость амплитуды колебаний от времени.

Проведенные на ЭВМ расчеты позволили получить оценки параметров ₀, ₁, ₀ по МНК, которые составили

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

Рис.5. Динамика ВНП США

Т а б л и ц а 11. Прогнозы и дисперсии ошибок прогнозов

Год прогноза	Прогноз ВНП	Дисперсия ошибки прогноза	Реальный ВНП	Отклонение
	3,45	0,45	3,72	0,27
	3,53	0,59	3,85	0,32
	3,68	0,82	4,02	0,34
	3,74	0,50	4,14	0,40

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