Инструментарий моделирования колебательной компоненты в колоколообразных кривых жизненного цикла продукта Текст научной статьи по специальности «Математика»

В. К. Семёнычев, Е. И. Куркин, Е. В. Семёнычев, А. А. Данилова

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

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

ключевые слова: жизненный цикл продукта; нефть; газ; уголь; колоколообразный тренд; колебательные компоненты; генетический алгоритм; ARMA-модель. JEL classification: C22; L71; D91.

1. введение

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

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

Альтернативой является построение феноменологических моделей, под которыми понимают эмпирически устанавливаемые закономерности. Такие модели применяют, опираясь на статистику данных в случаях, когда детальная картина явления сложна, а задачи моделирования и, главное, прогнозирования добычи продукта актуальны. Известно, что большое число феноменологических моделей для трендов объемов добычи продукта для отдельных месторождений, регионов и стран имеют колоколообразную форму (Хасанов и др., 2001; Мирзаджанзаде и др., 2004).

Одной из первых и широко используемых до настоящего времени феноменологических моделей является симметричная относительно t0 формула Хабберта (Hubbert, 1956; Bardi, Yaxley, 2005; Бажанов, Выскребенцев, 2007):

T 2

Y (t)= ™ (-^ + e(t), (1)

1 + ch (&(t -10))

где ch (•) — гиперболический косинус, t0 — абсцисса ЖЦП, соответствующая максимуму тренда Tmax траектории ЖЦП, e(t) — стохастическая компонента наблюдений ЖЦП, а параметр о определяет наклоны кривой роста и падения добычи.

Симметричная формула С. П. Капицы (Капица, 1999) была предложена впервые для моделирования динамики численности народонаселения Земли, а более поздние исследования (Хасанова и др., 2001) показали, что она успешно моделирует динамику ЖЦП месторождений добычи нефти, газа, угля и других невозобновляемых ресурсов:

T ■ 02

Y(t) = -^b-r + e(t) . (2)

(t -to) +о2

В качестве третьей, также симметричной кривой ЖЦП примем модель Гаусса (Bartlett, 2000; Brandt, 2007):

Y (t )= Tmax )2/ 02 +e(t) .

Известна и несимметричная модель Хаммонда и Маккея (Hammond, Mackay, 1993):

где

Norm = t0 o(i> to e"o(i >t0.

Y (t) = ta(t} to e—°(t' +e(t), Norm

Применение известных моделей к динамике тренда продукта отдельных месторождений, регионов и стран показало, что в большинстве случаев колоколообразные траектории трендов несимметричны: длительность этапа спада ЖЦП дольше длительности этапа роста (Bartlett, 2000; Brandt, 2007; Hammond, Mackay, 1993; Мирзаджанзаде и др., 2004). При этом известным приемом трансформации симметричных моделей в несимметричные является задание параметра о в моделях Хабберта и Гаусса в виде логистической функции Ферхюль-

ста о(t) = о1 + (о2 — о1 ))(l + e—it—t)/0т ), изменяющейся во времени от значения о1 на этапе роста (t < t0) до значения о2 на этапе падения ЖЦП (t > t0) (Brandt, 2009).

Предложенное расширение законов изменения о^) на функцию Гомперт-ца о = о1 + (о2 —о1 )exp{— 0.7e—(t—)/от}, Ричардса о = о1 +(о2 —о1 )(1 + e~(t—t°)/от )) (в нее дополнительно заложена возможность регулировки скорости изменения наклона логистических моделей при помощи параметра оT ), а также Рам-сея о = о1 + (о2 — о 1 )[1 + (1 + (t —10) / оT)e—(t—10)/от ] для t > t* и о = о1 при t < tl (где to = t0 — 1.678оT), приведенное в (Семёнычев и др., 2012а, 2012б), показало, что большое число моделирований трендов можно с высокой точностью реализовать на таких моделях при соответствующем выборе логистической функции, обеспечивающей асимметрию.

CQ

Вместе с тем, анализ известных данных по траекториям ЖЦП показал, что во многих | случаях при довольно высокой точности моделирования трендов не обеспечивается удовле- § творительная точность прогнозирования при наличии в наблюдениях выборки колебательной компоненты, к которой отнесены сезонное и циклическое колебания траектории ЖЦП вокруг уровней тренда, меньшие тренда по амплитуде, но более динамичные.

Системного подхода к решению данной проблемы, отражающего многообразие возмож- |

•О

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

рубежной литературе авторам обнаружить не удалось. ®

ад ui

2. Модели взаимодействия колебательной компоненты с трендом Ц

i

Моделирование взаимодействия колебательной компоненты с трендом будем осуществ- ^ лять на основе параметрического подхода, который в большей мере ориентирован на задачу прогнозирования и обычно требует для своей реализации меньшего объема выборок (Айва- I зян, 2001), чем непараметрический, подробно описанный в работах (Berndstrup, Hylleberg, & 2002; Семёнычев, Семёнычев, 2011). $

Наиболее просто считать колебательную компоненту не взаимодействующей с трендом (независимой) и входящей аддитивно в структуру ряда динамики имеющихся на практике наблюдений

Y = Tk + SA +ek, (3)

где Yk — уровни определяемого параметра ЖЦП, Tk — уровни тренда, Sf — уровни аддитивной колебательной компоненты, ek — уровни стохастической компоненты, k = 1,..., n — номера наблюдений ряда динамики при шаге опроса A, n — объем выборки.

На практике зачастую имеет место и более сложное взаимодействие колебательной компоненты с трендом ЖЦП.

Например, встречаются случаи пропорционально-мультипликативного взаимодействия, когда колебательная компонента ряда пропорциональна уровням тренда (Семёнычев, Семёнычев, 2011):

Yk = Tk (1 + SM ) + % = Tk +TX +Ek, (4)

где TkSkf — мультипликативная колебательная компонента (при условии |SM|<1).

Более общим случаем является присутствие в структуре ряда одновременно аддитивной и пропорционально-мультипликативной колебательных компонент:

Yk = Tk (1 + SM ) + SA + ek. (5)

Для модели (5) известен алгоритм идентификации параметров (Семёнычев и др., 2010а) лишь в случае представления тренда в виде полинома (степени NT) с линейными параметра-

NT

ми D-: Tk = ^D++1 (kA)Nt ' , а колебательных компонент, как это обычно делается (Pollock,

1=0

1993) — рядами гармоник с частотами w и Qr:

SM

NM na

= 2[ Aqsin (qkA) + Bq cos (wqkA)], SA = 2[E sin (Qr£A) + Fr cos (QrkA)].

q=1 r=1

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

Y = T + sr + e, где sr = ^A, ■ 11-g + g ^ •sin(A + ).

T

(6)

max /

Веса g в модели (6) соответствуют частотам w¿ колебательной компоненты, а параметр Tmax = max{T, k = 1,...,n} обеспечивает нормировку для возможности сравнения амплитуд аддитивной и мультипликативной колебательных компонент в одном диапазоне значений. Модель (6) обобщает структуры (3) (при g = 0) и (4) (при g =1). Каждая гармоника (с час-

T

тотой w¿) в этой модели будет суммой мультипликативной S^ = A —~sin (kА + <) и ад-

Tmax

дитивной SA = A sin (w,kА + < t) составляющих, взятых с весами g и 1 — g.

Наглядная графическая иллюстрация наложения аддитивной модели (3), пропорционально-мультипликативной модели (4) и взвешенной по амплитуде аддитивно-мультипликативной модели (6) колебательных компонент на колоколообразный тренд ЖЦП с произвольно выбранными параметрами приведена на рис. 1.

Y

1

0.8 0.6 0.4 0.2 0

S = S'

Y

10 20 30 40

50

1

0.8 0.6 0.4 0.2 t 0

SM

S =

Y

10 20 30 40

50

1

0.8 0.6 0.4 0.2 t0

-S = S

AM _

= (1-Y)SA+ySm

Y = 0.7

10 20 30 40

50

Рис. 1. Демонстрация взвешенной по амплитуде аддитивно-мультипликативной

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

В модели (4) при значениях уровня тренда Тк , близких к нулю, пропорционально-мультипликативная колебательная компонента Тк (1 + ) также будет близка к нулю, становясь соизмеримой с текущими значениями стохастической компоненты ек, что может существенно снизить точность идентификации параметров Бк . Модели (5) и (6) избавлены от этого недостатка, поэтому их можно интерпретировать как метод повышения точности (вычислительной устойчивости) идентификации компоненты Бк4 для прогнозирования на этапе спада ЖЦП. Значения у; определяют «вес» мультипликативного вхождения каждой частоты юi независимо от других частот, что приводит к более наглядному разделению амплитуд, частот и фаз.

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

=l

чем при малых значениях тренда. Именно такой характер изменения частоты колебательной компоненты наблюдался в некоторых реальных выборках ЖЦП.

Для моделирования таких случаев введем параметр в в формулу для мгновенной частоты:

/ T \~в T

\ max /

(7)

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

Тогда фаза колебаний переменной частоты колебательной компоненты определится как интеграл от мгновенной частоты (Сергиенко, 2002):

ik ik F k = fw(s )ds + Po = wf

T (s )

T

to \ max /

' + Po

(8)

где Т(я), «(я) — значения тренда и частоты в момент времени 5, ^ — момент времени, соответствующий фазе Фк .

О

0

1

I

!ф §

ф

0 со iu

1

I

S: iu

I

!ф §

ф

о Щ

Учитывая выражения (7), (8) и (6), можно в данном случае говорить о взвешенной по частоте колебательной компоненте:

Q

5Ql = 2 A

i+g,

T

Jjl.-i

T

\ max /

81П

•"■Q

5Q = 2 A

T

T

81П

тах /

w

k fk

f

'o

T (s )

IM ^

T

:+<p,

T

+p

(9) (10)

Частота и амплитуда в (9) и (10) могут быть не зависимыми от тренда (при у = 0, в = 0), а могут и зависеть от него (рис. 2).

15 10 5 0 -5

fj /\f J / \ \ д

-..-г.. ..чЧ/ S N .........vr 4 \ * N 1 / * V t Is 4r -

5 10 15 20 25 30 0 5 10 15 20 25

а) б)

Рис. 2. Вид колебательной компоненты: а) с постоянной частотой в = 0; б) с переменной частотой в = - 0.3 (пунктир — тренд и колебательная компонента, сплошная линия — общий вид модели)

30

I

\

1=1

1=1

Y

Y

t

t

0

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

3. Идентификация параметров моделей

Для идентификации параметров рассматриваемых сложных многопараметрических и нелинейных моделей используется генетический алгоритм, в котором осуществляется поиск решения путем подбора, комбинирования и вариации искомых параметров методом, напоминающим биологическую эволюцию. Такой алгоритм впервые был предложен в (Holland, 1975) и получил в последние годы существенное распространение (Хасанов и др., 2001), в том числе и на задачи идентификации параметров моделей временных рядов (Ursu, Turkman, 2012). Генетический алгоритм для данных моделей состоит из следующих основных стадий реализации:

1) случайным образом генерируется конечный набор пробных решений — первое поколение параметров модели;

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

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

4) в случае продолжения алгоритма генерируется новое поколение параметров посредством операторов скрещивания и мутаций, затем осуществляется переход к п. 2 генетического алгоритма.

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

Для идентификации моделей, содержащих тренд и колебательные компоненты моделей (3), (4), (6), (9) и (10), оправданно воспользоваться методом параметрической итеративной декомпозиции рядов с колебательной компонентой (Семёнычев и др., 2010б). Согласно этому методу на первой итерации генетическим алгоритмом находят параметрическую модель тренда (Хасанов и др., 2011). Затем после детрендирования ряда идентифицируют колебательную компоненту, затем осуществляют десезонализацию, вновь идентифицируют тренд и т. д. На реальных выборках обычно достаточно трех-четырех итераций.

SAM - SA =2 A sin (w k A + p¡).

QQ

applied econometrics_прикладная эконометрика | № 33 (1) 2014

Для идентификации параметров сезонной компоненты в моделях вида (3) и (4) целесообразно использовать метод обобщенных параметрических ARMA-моделей ряда динамики | (Семёнычев, Семёнычев, 2011). При этом сначала находят оценки частот w как параметров * соответствующих ARMA-моделей. Затем определяют фазы p¿ путем представления гармони- ^ ки в виде A sin (w¡t + p¿) = Bi sin (ww f) + Ci cos (wwt), где Bi = A cos p t, Ct = A sin p i. После «j;

n (O

i ^

этого идентифицируют B¡ и C¡, рассчитывают параметр p ¡ = arctg — , определяют ампли- í

B¡ =;g

туду Af и параметр g¡. Одновременное нахождение амплитуд, частот, фаз всех гармоник ко- -лебательных компонент, а также весовых параметров в моделях (6), (7), (9), (10) требует

решения задачи идентификации параметров в три (для моделей (3), (4)), четыре (для мо- ^

дели (6)) и пять (для моделей (9), (10)) раз большей размерности, чем число определяе- | мых гармоник. Последовательное определение параметров колебательной компоненты

позволяет в два-три раза сократить размерность задачи идентификации колебательной s;

компоненты, что особенно актуально при идентификации моделей колебаний со многи- "i,

QQ

ми частотами. Использование ARMA-моделей для идентификации частот колебательных | компонент оправданно и при идентификации параметров модели (6), если предположить,

что частоты взвешенной колебательной компоненты можно определить, считая ее ® аддитивной:

ф о

CQ

Достоинством примененных АЯМА-моделей, в сравнении с известными методами максимального правдоподобия или нелинейным МНК (Айвазян, 2001), является возможность применения этих моделей на относительно коротких выборках при малом числе (до трех-четырех) итераций. Использованные методы идентификации рассматриваемых колоколо-образных моделей ЖЦП с колебательными компонентами реализованы авторами в «Программе моделирования и прогнозирования уровней добычи нефти и газа "011_ИеП;"» (Семёнычев и др., 2012в).

¡=1

4. Методика оценки точности моделирования и прогнозирования

траектории жцп

В работе рассматривается решение задачи прогнозирования уровней ЖЦП после того, как точка пика ЖЦП была пройдена. Разделим исходную выборку на «рабочую часть» наблюдений, по которой строится модель ЖЦП, и на известную «контрольную часть» наблюдений, по которой будем оценивать точность полученного прогноза (рис. 3). Границу окончания рабочей части наблюдений выберем на расстоянии 11Аеп± от времени начала выборки и исследуем точность при различной глубине прогноза ^

Выбор td¡etií может осуществляться достаточно произвольно. Однако имеет смысл изучать выборки после прохождения пика добычи, например, с наблюдения tpe¡¡k + 3 < < tm¡í¡. — 1, т. е. выбором tiAetií из следующего диапазона: от трех лет после прохождения пика добычи нефти и до года перед окончанием известной выборки. Точку разделения данных на «рабочую» и «контрольную» части будем в процессе анализа смещать по оси времени для оценки достигаемой точности моделирования и прогнозирования.

ident t

14 12 10 8 6 4

5 10 15 20 25 30 35 40

Рис. 3. Демонстрация методики оценки точности прогноза

Затем проведем усреднение оценок параметров по всем выборкам, начинавшихся с различных tident (усредняя значения критериев по всем tident с одинаковым tpr). В дальнейшем для каждого временного ряда будем сравнивать усредненные критерии точности прогноза с разным горизонтом t

В качестве критерия точности моделирования траектории принят коэффициент детерми-

^ident ^ /Rident ^

нации R2 = 1 — ^ (Yk — YM°d ) ^ (Yk —Yk ) , а точность прогнозирования оценена с помо-

k=1 / k=1

щью второго коэффициента Тейла

T = 12

I N

2 ( -YMod) 2 Yk2 +2 Kd)2

=Nident+l / \ k=Nident +1 k=Nident +1

■100%,

который близок по своим значениям к известной МАРЕ-оценке прогнозирования (Айвазян, 2001), но обладает большей устойчивостью при значениях уровня ЖЦП, стремящихся к нулю (на стадии падения уровней добычи). В описании критериев и N — номера последних значений уровней ЖЦП, используемых для построения модели и оценки точности прогноза соответственно.

Y

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

В рамках настоящего исследования были проанализированы десятки реальных выборок ЖЦП добычи нефти, газа и угля. Для каждой из них выбиралась модель тренда с большим значением Я2. На них показана возможность повышения точности прогнозирования с присутствием колебательной компоненты в траектории ЖЦП.

Иногда учет в модели ЖЦП колебательной компоненты не изменял заметно значения критериев точности. Вместе с тем зафиксировано много случаев заметного улучшения точности моделирования и, главное, прогнозирования, при наличии в модели колебательной компоненты.

На рисунке 4 представлен результат, в соответствии с которым наиболее точной моделью тренда добычи нефти в Норвегии (по критерию коэффициента детерминации Я) оказалась модель в виде тренда Хабберта с асимметрией Рамсея, а большую точность дает учет двух гармоник аддитивных колебаний модели (3). Однако главный итог исследований заключается, пожалуй, в существенном снижении погрешности прогноза по критерию Тей-ла Т2 (рис. 4в). К примеру, при горизонте прогноза пять лет дополнение модели тренда колебательной компонентой уменьшило ошибку прогноза добычи нефти в Норвегии с 7.5% до 2% (рис. 4в).

На рисунках 4-10 представлены: серой пунктирной линией — исходные данные, черной линией — построенная модель (сплошная для точек, по которым проводится идентификация параметров, пунктир — для точек оценки прогноза).

4 000 3 000 2 000 1 000 0

Y=Trend

\\ \

t

4 000 3 000 2 000 1 000 0

Y=Trend+SA

T, %

\ -

/

t

/ / /

Trenc /

у Trend + >А

о

0

1

«i «i

t

!ф §

Ф

0

CQ ai

1 I

s:

ai

t

!ф §

Ф

0

01

1970 1980 1990 2000 2010

1970 1980 1990 2000 2010

а)

б)

в)

Рис. 4. Добыча нефти в Норвегии (тыс. барр. в день): а) модель Хабберта с асимметрией Рамсея;

б) модель тренда, дополненная моделью колебаний (3) с двумя гармониками;

в) оценка точности прогноза по критерию Тейла, R = 0.9955, R:^+S = 0.9960

На рисунке 5 представлены результаты, относящиеся к добыче газа в Евросоюзе. Они получены на основе модели Хабберта с асимметрией Рамсея и демонстрируют учет аддитивного взаимодействия двух моделей трендов Trendj и Trend2, а также аддитивной колебательной компоненты (3) с тремя гармониками. В данном примере решена задача моделирования тренда ЖЦП, состоящего из суммы двух несимметричных колоколообразных моделей Хабберта с асимметрией Рамсея (ранее было проведено моделирование лишь для суммы симметричных моделей Хабберта (Patzek, Croft, 2010)), а дополнительный учет трех колебательных компонент привел к существенному повышению точности прогнозирования.

На рисунке 6 показано, что более точную модель добычи нефти для стран Организации экономического сотрудничества и развития (OECD) R дает тренд из суммы трех моделей Капицы с асимметрией Гомперца (Trendj, Trend2 и Trend3) и дополнительный учет колебательной компоненты пропорционально-мультипликативной модели (4) с двумя гармониками.

При этом добавление в модель ЖЦП двух гармоник пропорционально-мультипликативных колебаний модели (4) не позволило существенно увеличить точность моделирования коэффициента RT+S, но улучшило точность прогнозирования, снизив максимальные значения критерия T2 при горизонте прогноза до шести лет с 3 до 1.3%.

8

6

4

t

2

1

2

3

4

5

250 200 150 100 50 0

Y=Trend

¿i7 ч.

Л

Tre ndj

Trend2

250 200 150 100 50 0

Y=Trend+S

A

T2' %

\ •

Trend J

Trend2

Trend

> / /

— — —

Trend + SA

1970 1980 1990 2000 2010 1970 1980 1990 2000 2010

а)

б)

в)

Рис. 5. Добыча газа в Евросоюзе (млрд куб. м в год): а) модель Хабберта с асимметрией Рамсея; б) модель тренда, дополненная аддитивной моделью колебаний (3) с тремя гармониками; в) оценка точности прогноза по критерию Тейла, Я2 = 0.9661, = 0.9814

Y=Trend

Y=Trend+SM

T2' %

20 000 15 000 10 000 5 000 0

1970 1980 1990 2000 2010 а)

20 000 15 000 10 000 5 000 0

3 t

1970 1980 1990 2000 2010 б)

0

Tre nd \ / / /

--- — —

Trend + SM

1 2 3 4 5 6

в)

Рис. 6. Добыча нефти в странах OECD (тыс. барр. в день): а) модель Капицы с асимметрией Гомперца;

б) модель тренда, дополненная моделью колебаний (4) с двумя гармониками;

в) оценка точности прогноза по критерию Тейла, R = 0.9679, R^+s = 0.9685

Описание динамики добычи тренда нефти в США моделью Хабберта с асимметрией Гомперца при колебательной компоненте с зависимыми от тренда частотой и амплитудой по модели (10) с одной гармоникой незначительно повысило точность моделирования, но существенно повысило точность прогнозирования (рис. 7).

Рисунок 8 иллюстрирует точностные характеристики моделирования и прогнозирования добычи газа в Великобритании трендом из суммы двух моделей Хабберта с асимметрией Гомпертца (Trendj и Trend2) и с взвешенной по амплитуде колебательной модели (6) с тремя гармониками.

Существенный выигрыш в точности получен и для добычи твердого угля в Германии (рис. 9). Здесь использована уже модель Капицы с асимметрией Ричардса и одна гармоника колебательной компоненты взвешенной по частоте модели (9).

На рисунке 10 приведен пример моделирования и прогнозирования для отдельного месторождения ОАО «НК «Роснефть». Дополнение тренда моделью колебаний (10) приводит к улучшению точности прогноза в 1.5-2 раза и повышает коэффициент детерминации, оценивающий точность моделирования, с 0.81 до 0.97.

6

5

4

3

t

t

t

2

1

2

3

4

3

2

t

t

10 000 8 000 6 000 4 000 2 000 0

Y=Trend

1900 1925 1950 1975 2000

а)

10 000 8 000 6 000 4 000 2 000 0

Y=Trend+S'

m

1900 1925 1950 1975 2000 б)

3.5 3 2.5 2 1.5 1

T2' %

T ren d —

— --' ->

T ren d + 2

123456789 10

в)

Рис. 7. Добыча нефти в США (тыс. барр. в день): а) модель Хабберта с асимметрией Гомперца;

б) модель тренда, дополненная моделью колебаний (10) с одной гармоникой;

в) оценка точности прогноза по критерию Тейла, Я2 = 0.9892, = 0.9949

120 100 80 60 40 20 0

Y=Trend

Y=Trend+SW

Л

Trend2

v

Trend1

T2' %

Trend

/ ! / !

/ / / f f Trend + SW1

1970 1980 1990 2000 2010 а)

1970 1980 1990 2000 2010 б)

3

в)

Рис. 8. Добыча газа в Великобритании (млрд куб. м в год): а) модель Хабберта с асимметрией Гомперца;

б) модель тренда, дополненная моделью колебаний (6) с тремя гармониками;

в) оценка точности прогноза по критерию Тейла, Я = 0.9954, = 0.9960

200 150 100 50

Y = Trend

0 1945

Г 4

1 .4

4; .

1965

1985 а)

2005

0 1945

1965

1985

t

2005

T2' %

—i

Trend

Trend+Sn1

б)

2468

в)

t 10

Рис. 9. Добыча твердого угля в Германии (млн т в год): а) модель Капицы с асимметрией Ричардса; б) модель тренда, дополненная моделью колебательной компонентой (9) одной гармоники; в) оценка точности прогноза по критерию Тейла, Я2 = 0.9679, = 0.9883

О

0

1

t

!ф §

ф

0 со ai

1

si

ui

t

:ф §

ф о

со

t

t

t

t

t

1

2

4

5

t

Y = Trend

15 12.5 10 7.5 5

2.5 2001

А"»

[j N v V

s\

___._

2005

t

2009 2012

15

12.5

10

7.5

5

2.5 2001

40 30 20 10 0

T2' %

Trend

Trend + Sn2

2005

2009 2012

а) б) в)

Рис. 10. Добыча нефти на месторождении ОАО «НК «Роснефть» (тыс. барр. в день): а) модель Хаммонда-Маккея без дополнительной асимметрии;

б) модель тренда, дополненная моделью колебаний (10) с одной гармоникой;

в) оценка точности прогноза по критерию Тейла, Я = 0.8107, = 0.9707

Y = Trend+S "2 1ю

t

t

2

3

4

5

6. Заключение

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

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

Айвазян С. А. (2001). Прикладная статистика. Основы эконометрики. М.: ЮНИТИ-ДАНА.

Бажанов А. В., Выскребенцев А. С. (2007). Адекватность кривых Хабберта для прогнозирования темпов добычи нефти. MPRA Paper 15117.

Капица С. П. (1999). Сколько людей жило, живет и будет жить на Земле. Очерк теории роста человечества. М.: Международная программа образования.

Мирзаджанзаде А. Х., Хасанов М. М., Бахтизин Р. Н. (2004). Моделирование процессов нефтегазодобычи. Нелинейность, неравновесность, неопределенность. Москва -Ижевск: Институт компьютерных исследований.

Семёнычев В. К., Семёнычев Е. В., Куркин Е. И. (2010а). Моделирование рядов экономической динамики полиномиальным трендом и одновременным вхождением аддитивной и мультипликативной колебательных компонент. Вестник Самарского муниципального института управления, 3 (14), 7-20.

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

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

applied econometrics_прикладная эконометрика | № 33 (1) 2014

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

Со

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

s

Семёнычев Е. В. (2012). Эконометрическое моделирование жизненного цикла продукта. Самара, САГМУ ^ Семёнычев Е. В., Куркин Е. И., Молостова П. А. (2012а). Выбор модели колоколообразной фор-

со1

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

комплексом, 8, 28-34. i

' ' ¡ф

Семёнычев В. К., Куркин Е. И., Семёнычев Е. В. (20126). Идентификация параметров импульс- ®

ад ui

Семёнычев В. К., Куркин Е. И., Семёнычев Е. В., Рязанцев С. В., Данилова А. А. (2012в). Про- § грамма моделирования и прогнозирования уровней добычи нефти и газа «Oil_Ident». Свидетельство

S: ai

о государственной регистрации программ для ЭВМ № 2012619424 от 18.10.2012. Сергиенко А. Б. (2002). Цифровая обработка сигналов. Спб.: Питер.

Хасанов М., Карачурин Н., Тяжев Е. (2001). Оценка извлекаемых запасов нефти на основе феноменологических моделей. Вестник инжинирингового центра ЮКОС, 2, 3-7.

Bardi U., Yaxley L. (2005). How general is the Hubbert curve? The case of fisheries. In: Proceedings of the 4th International ASPO Conference. Lisbon, Portugal. http://europe.theoildrum.com/files/bardiyax-leyaspo2005.pdf.

Bartlett A. A. (2000). An analysis of US and world oil production patterns using Hubbert-style curves.

Mathematical Geology, 32 (1), 1-17.

Brandt A. R. (2007). Testing Hubbert. Energy Policy, 35 (5), 3074-3088.

Brandt A. (2009). Methods of forecasting future oil supply. UKERC Review of Evidence for Global Oil Depletion, Technical Report 6.

Berndstrup B., Hylleberg S. (2002). Seasonality in economic models. http://mit.econ.au.dk/vip_htm/ monielsen/papers/seasonality.pdf.

Hammond G. P., Mackay R. M. (1993). Projections of UK oil and gas supply and demand to 2010. Applied Energy, 44 (2), 93-112.

Holland J. H. (1975). Adaptation in natural and artificial systems. University of Michigan Press, Ann Arbor.

Hubbert M. K. (1956). Nuclear energy and the fossil fuels. Amer. Petrol. Inst. Drilling & Production Practice. Proc. Spring Meeting, San Antonio, Texas, 1956, 7-25.

Patzek T. W., Croft G. D. (2010). A global coal production forecast with multi-Hubbert cycle analysis.

Energy, 35, 3109-3122.

Pollock D. S. G. (1993). Lectures in time-series analysis and forecasting. Queen Mary and Westfield College, The University of London. http://www.le.ac.uk/users/dsgp1/courses/tseries/2cycles.pdf.

Ursu E., Turkman K. F. (2012). Periodic autoregressive model identification using genetic algorithms.

Journal of Time Series Analysis, 33 (3), 398-405.

t

:ф §

ф о

CQ