Анализ затрат предприятия с помощью вейвлет-преобразований Текст научной статьи по специальности «Математика»

ЭКОНОМИКО-МАТЕМАТИЧЕСКОЕ

МОДЕЛИРОВАНИЕ

УДК 519.6:517.44+658.15

АНАЛИЗ ЗАТРАТ ПРЕДПРИЯТИЯ С ПОМОЩЬЮ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЙ*

А. А. МИЦЕЛЬ,

доктор технических наук,

профессор кафедры автоматизированных

систем управления, Томский государственный университет систем управления

и радиоэлектроники; профессор кафедры высшей математики и математической

физики, Национальный исследовательский

Томский политехнический университет

E-mail: maa@asu.tusur.ru

А.Н. ШЕМЯКИНА,

E-mail: alesya.as23@yandex.ru

Введение

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

* Статья подготовлена по материалам журнала «Экономический анализ: теория и практика». 2013. № 46 (349).

статье авторы используют доступные в Интернете экономические показатели с наименьшим шагом дискретизации - один месяц.

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

Проблема анализа экономических показателей в виде сложных нестационарных рядов встречается в источниках [2, 11, 16, 18, 19], где основная цель работ состоит в том, чтобы одномерный временной ряд представить в виде набора компонент, а именно, в выделении циклических составляющих. В этом случае формирование модели ряда состоит в его декомпозиции на соответствующие компоненты.

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

1-

0,5-

-0,5

лиз. В данной работе используются вейвлет-преобразования.

Вейвлет-анализ широко используется в различных областях науки и техники, таких как сейсмология, анализ речи, обработка изображений, изучение мульти-фрактальных объектов, разных сигналов в медицине (ЭКГ). Это можно встретить в работах [1, 8, 17].

В экономике применение вей-влет-преобразования представлено в работах [4, 15].

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

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

Постановка задачи

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

Пусть у = {уп, п = 0,1,..., L _ 1} - временной ряд с расстоянием между отсчетами ^ т. е. частота дискретизации составляет / = 1 / ^ . Будем считать, что L = 2], J е N.

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

Обозначим ^ = (Ь

0,0, °0,0, °1,С °1,1, °2,0" ■

2,3' '

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

Рис. 1. Вейвлет Добеши 4

временной ряд. Если положить Ь00 равным нулю, все остальные вейвлет-коэффициенты оставить без изменений и сделать обратное дискретное вейвлет-преобразование (ОДВП), то результат будет такой же, как если из ряда у вычесть полиномиальный тренд степени N — 1. Напомним также, что для широкого класса стохастических процессов с хорошей точностью можно считать, что вейвлет-коэффици-енты а независимы и распределены как М(0,2_^-о2).

Для удобства будем это обозначать записью ^ = Ж (у), а обратное ОДВП будем обозначать записью

у = ж-^у

Энергия ряда на уровне разрешения 3 определяется выражением

23 -1

3

(1)

Определим дискретный вейвлет-спектр (ДВС) как вектор размерности ] + 1, равный

Р(Ь20,р0, Р1, . . .Ру-1).

Здесь у-я компонента ДВС показывает, сколько энергии временного ряда содержится в указанном отрезке частот.

Предположим, необходимо разделить временной ряд длины Ь = 23, ] е N на составляющие (циклическую и трендовую). Для этого на основе вейвлет-коэффициентов ^ = W (у) подсчитывается ДВС. Пусть, например, ДВС имеет два пика на уровнях разрешения у1 и у2 > ] + 2. Разделим вектор коэффициентов V! на два вектора коэффициентов wl и одинаковой с вектором V длины по формулам:

3 = 0,1,... 3 +1 1; 3 = 3 + 2,..., J -1/'

■аНа.3,1!

(2)

а(2 =

0,

а

3 = 0,1,...,32 -1 1

31 = 32 - 1, 32,...,-1 - 1

2

а

Ъ(1) = Ь . ¿(2) = 0

"0,0 "0,0 > 0,0

После этого подсчитываем ОДВП векторов w1 и w2, т. е. выделяем две искомые составляющие у О) = Ж1^ (1)) и у (2) = Ж-1^ (2)). Трудности возникают, если, например, /2 = ] + 2, т. е. когда между пиками ДВС нет промежутка. Необходимо решить, какому именно из векторов ^ (1) или w (2)) приписать коэффициенты, у которых значение уровня разрешения равно ] = ] + 1. Предлагается использовать следующий подход.

Пусть . такой уровень разрешения, что энергия этого уровня разрешения Р. заключена между двумя пиками ДВС. Вейвлет-коэффициенты припишем векторам w (1) и w (2) по следующему правилу:

а ш = а -. ;

а + Ъ

П у п Ъп

а ® = —•— а -. ,

а. + Ъ ] ,п'

где коэффициенты апи Ъп, п = 0,1,..., 2] - 1 вычисляются по формулам:

ап = а,Г п . 2

Ъп = (а]+1,[2п] + а ]+1,[2п+1] ,(1)

)А

При таком подходе а(Щ + а(2) = а.п, т. е. векторы w (1) и w (2) остаются аддитивными: w (1) + w (2) = w. Это очень важно, так как выделяемые компоненты тоже будут аддитивными у (1)+у (2) = Ж ^ (1)) + + (2)) = у.

Стоит отметить, что масштабные коэффициенты Ъ.п несут информацию о трендовом, глобальном поведении сигнала, в то время как вейвлет-ко-эффициенты а.п несут информацию о локальном отклонении сигнала от тренда [4, 5].

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

Проведем декомпозицию временного ряда общих затрат предприятия. Исходными данными будут служить помесячные данные общих эксплуатационных затрат с 2008 по 2013 г. одного из филиалов ООО «Газпром трансгаз Томск». Выборка состоит из 72 значений затрат. Однако для дискретного вейвлет-преобразования

обязательным условием является представление ряда в виде 2' элементов, где J - целое число. Были выбраны 26 значений затрат предприятия. Этот временной ряд представлен на рис. 2.

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

Расчет вектора вейвлет и масштабных коэффициентов w с помощью прямого вейвлет-преобразо-вания проводился в пакете Mathcad 15. Был выделен дискретный вейвлет-спектр сигнала согласно формуле (1), который представлен на рис. 3.

Анализ рис. 3 показывает, что ДВС имеет пики при ] = 0 и ] = 3. Проведем разделение вектора на три вектора w (1), w (2), w (3) при . = 0, 2, 4 по формулам (2). Причем данные векторы имеют ту же длину, что и вектор вейвлет-коэф-фициентов w. Проведя обратное ДВП векторов (2) и w (3), получим составляющие сигнала

w

(1)

w

и 64

1IIIIIIIIII1 ООО 111 им

о сч л о сч « о сч л о сч л о сч «

га и 3 о 12 га и 3

X и

I I I I I I I I I I I I I I I I I I I I I I

и

и

О

и 64

и

и

О

и 64

Рис. 2. Затраты предприятия с января 2008 г. по апрель 2013г., млн руб. Р

Рис. 3. Дискретный вейвлет-спектр затрат предприятия

0

2

5

1 1 1111 ММ 1 II 1 11111 111 НИМИ 1111111 11 м МП 111 м 111 м 1 м м м 1

ОС ос ос © © © со со

© © © © © © © © © © © © © о © © ©

Л « Л Л Л Л Л Л Л Л 'И Л Л Л

£ ю та ш £ ю та ш £ Ю « та ш ю я та ш 3 ю я та ш и о.

£ и £ и и с <

и и и ¿3 ¿3

~"Исходный сигнал Трендовая составляющая

Рис. 4. Исходный ряд и трендовая составляющая затрат, млн руб.

Исходный сигнал

Циклическая составляющая

Рис. 5. Исходный ряд и циклическая составляющая затрат, млн руб.

ОС ОС © ОС © © © ©

© © © © ©

Л Л Л Л

£ и ю ^ ю

X Ё и и Ё и и

■ ■ 1 © 1 ■ ■ ■ 1111 |||| ||||| М11 1111 1 ■ ■ 1 со

© © © © © © © ©

Л Л Л Л « Л Л

л л <3 о. л та о. л

ю « та м £ ю я та м £ ю я та м

и Ё и ё и

и и и

Исходный сигнал

Шумовая составляющая

Рис. 6. Исходный ряд и шумовая составляющая затрат, млн руб.

Т = Ж-1^ (1)), 1 = 1, 2, 3. Графики этих векторов представлены на рис. 4-6.

Видно, что первая составляющая указывает на бизнес-цикл или трендо-вую составляющую (см. рис. 4). Это можно объяснить тем, что именно в V (1) приписан коэффициент Ь00, который, как отмечалось ранее, отвечает за тренд. В остальных случаях этот коэффициент был обнулен.

Составляющая Т (2) (см. рис. 5) представляет собой более высокочастотные (сезонные) колебания.

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

Это хорошо иллюстрирует график суммы компонент Т(1) и Т (2) и исходный временной ряд (рис. 7). Получена более сглаженная кривая, не включающая шумовые колебания.

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

Прогноз затрат предприятия

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

Прогноз трендовой составляющей. Прогнозная модель строится на основе коэффициентов трендовой составляющей вейвлет-преобразования.

В качестве регрессионной модели выбрано уравнение ^ = 29 890 + + 441,82? - 4,113?2.

Поиск коэффициентов уравнения регрессии осуществлялся с помощью метода наименьших квадратов. Процедура поиска была реализована в Mathcad.

Коэффициент детерминации построенной модели Я2 = 0,92, Р-статистика Фишера равна 332 при теоретическом значении, равном 3,15 (уровень значимости а = 0,95); значения ¿-статистики для коэффициентов модели равны соответственно: ¿1 = 15, ¿2 = 9 (теоретическое значение ¿-статистики Стьюдента при уровне значимости 0,95 ¿Т = 1,67).

Значимость коэффициента детерминации подтверждается по статистике Фишера - Снедекора - РБнабл = 350,75 при теоретическом значении 3,15 при уровне значимости 0,95.

После проверки адекватности модели необходимо проанализировать остатки

е = у " y,

где у - эмпирические данные.

Проверим гипотезу о нормальном распределении остатков с помощью критерия Пирсона. Остатки распределены нормально с вероятностью 95 %. Проверка реализована в программе Statistica.

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

На рис. 8 приведен прогноз трендовой составляющей затрат на 12 мес. (номера месяцев с 65 по 76). Здесь же приведен коридор ошибок. Видно, что с увеличением периода прогноза коридор ошибок растет.

Прогноз сезонной составляющей. Прогнозная модель строится на основе коэффициентов сезонной составляющей вейвлет-преобразования.

45 п....................................................................................................................................................................................................................

гп ТТЛ гттп ТТЛ II 1 ГТ тл- ГГГП гпт ггт 11111 ГТП 1 II И 1111111 гптт 11111

со со со сл сл ел о о о ,—! СО СО

о о о о о о •—1 •—1 •—1 •—1

<4 <ч

Л (X М « Л Л Л (X м 5К Л 1а А (X м 5К Л 1а Л (X м 5К Л л ю »1 Л (X м 5К Л л ю »1 Л (X м Л

1 1 1 1 1 ш л

Я и О Я и о Я и О Ё и О и о <

' Исходный сигнал Сумма составляющих Т*1-1 и Т®

Рис. 7. Сравнение суммы составляющих Т (1) и Т (2) с исходным рядом затрат, млн руб.

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

~к(Г -12,46)"

К = 0,735 - 5 079^

+10 65Ып

+4 7^т

-4 803sin

- 23,09) 7,975 к(Г - 29,15)'

4,536 к(Г - 37,37) 8,106

+4 435sin

-1617cos

- 2 077cos

- 2 665cos

8,146

к(Г - 25,25)

6,665 к(Г - 47,66) 11,66 . к(Г - 64,23)' 6,239

оО «Ооо о_____

40-

11 111 111 111 I 11 III 111 1111 11111 111111IIII11I11111 1111 11 1111 111 111 11111 111 111 11

1 4 7 10 16 22 28 34 40 46 52 58

— Линия регрессии ----Коридор ошибок

16 22 28 34

О Эмпирические данные •• Прогнозные данные

Номер месяца

Рис. 8. Прогноз трендовой составляющей затрат на 2013-2014 гг., млн руб.

- 51,36) 4,423

где t - время, отсчитываемое от января 2008 г., мес.

Модель прогноза сезонной составляющей на 12 мес. (с мая 2012 г. по апрель 2013 г.) приведена на рис. 9.

Коэффициент детерминации построенной модели Я2 = 0,89.

Оценка дисперсии ошибок модели: с = 2 870 при уровне значимости 0,9; Р-статистика Фишера для полученной модели равна 40,3 при теоретическом значении,

О Эмпирические данные Линия регрессии ...... Прогнозные данные ----Коридор ошибок

Рис. 9. Модель прогноза сезонной составляющейзатрат на 2013-2014 гг., млнруб.

равном 1,95 (уровень значимости а = 0,95); теоретическое значение ^-статистики Стьюдента при уровне значимости равно 0,95, ?Т=2,0003, а полученные значения ^-статистики для коэффициентов модели равны соответственно: ? = 3,5; ? „ = 4,3; ?, = 5,6; t . = 4;

п1 ' ' п2 ' ' п3 ' ' п4 '

п = 9,8; п = 5,б п = 4,8; ?л = 3; п=8,8; и=5,3

Значимость коэффициента детерминации подтверждается. Наблюдаемое значение ^^ бл = 42,882 при теоретическом значении 1,9 и уровне значимости 0,95. Таким образом, можно сделать вывод о значимости всех полученных коэффициентов и модели в целом. То есть использование модели в расчетах статистически обоснованно.

После проверки адекватности модели необходимо проанализировать остатки Б = У " У,

где у - эмпирические данные.

Гипотеза о нормальном распределении остатков также подтверждается с вероятностью 95 %.

Анализ результатов

Построив регрессионные модели составляющих затрат и спрогнозировав их, была построена прогнозная модель затрат предприятия на период с мая 2012 г. до мая 2014 г. Эта модель (рис. 10) имеет вид

# = 29 890 + 441,8? - 4,113?2 -

-5 079^

+1617^

%(? -12,46) 8,146 _ '%(? - 25,25) 6,665

-10 65Ыи

-4 7^т

%(? - 23,09)

7,975 к(? - 29,15)" 4,536

-2 077^ -2 665cos

%(? - 47,66) 11,66 %(? - 64,23) 6,239

- 4 803sin 4 435sin

%(? - 37,37) 8,106 %(? - 51,36) 4,423

Анализ рис. 10 показывает, что фактические значения операционных затрат в апреле и мае 2013 г. попадают в доверительный интервал прогноза, что подтверждает качество построенной модели. С вероятностью Р = 81 % прогноз верен.

Заключение

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

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

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

55 т

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

1. Астафьева Н. М. Вей-влет-анализ: основы теории и примеры применения // Успехи физических наук. 1996. Т. 166. Вып. 11. С. 1145-1170.

2. Бессонов В. А. Введение в анализ российской макроэкономической динамики переходного периода. М.: Институт экономики переходного периода, 2003. 152 с.

3. Бессонов В. А., Петроневич А. В. Сезонная корректировка как источник ложных сигналов. М.: ВШЭ, 2013.

4. Бурнаев Е. В. Применение вейвлет-преоб-разования для анализа экономических временных рядов / Математическое моделирование развивающихся экономических систем // В сб. научн. трудов летней школы по экономико-математическому моделированию ЭК0М0Д-2006. Киров: ВятГУ, 2006. С.95-170.

5. БурнаевЕ. В., ОленевН. Н. Мера близости для временных рядов на основе вейвлет-коэффициентов // Тр. XLVIII научн. конф. МФТИ. Долгопрудный: ФУПМ, 2005. С. 108-110.

6. Витязев В. В. Вейвлет-анализ временных рядов. СПб: СПбГУ, 2001. 58 с.

7. ВоскобойниковЮ. Е., Гочаков А. В., КолкерА. Б. Фильтрации сигналов и изображений: фурье и вейв-лет-алгоритмы (с примерами в Mathcad): монография. Новосибирск: НГАСУ (Сибстрин), 2010. 188 с.

8. Гончаров А. А. Использование вейвлетов при построении моделей объектов управления с нерегулярным измерением выхода // Моделирование систем. 2013. № 1.

9. Губанов В. А. Выделение нестационарной циклической составляющей из временных рядов. URL: http://www.ecfor.ru/pdf.php?id=books/ sa2006/07.

О Реальные данные — Прогнозные данные

VIII 2013

----Коридор ошибок

Рис. 10. Модель для прогноза затрат предприятия на 12 мес. (май 2012 г. - апрель 2013 г.), млн руб.

10. Губанов В. А. Выделение тренда из временных рядов макроэкономических показателей // Научые труды ИНПРАН. М.: МАКС Пресс, 2005.

11. Губанов В. А. Оценка и прогноз конъюнктурных циклов в трендах экономических временных рядов // Науч. тр. ИНПРАН. 2006. Т. 4. С. 154-175.

12. Дьяконов В. П. Вейвлеты. От теории к практике. М.: СОЛОН-Р, 2002. 448 с.

13. Дьяконов В. П. MATLAB. Обработка сигналов и изображений: спец. справочник. СПб: Питер, 2002. 608 с.

14. Короленко П. В., Рыжикова Ю. В. Моделирование и обработка случайных сигналов и структур: учеб. пособие. М.: МГУ, 2012. 67 с.

15. Огородов А. П. Применение теории вейвлет-преобразования в исследовании финансовых рядов. URL: http://fetmag. mrsu. ru/2009-3/pdf/Financial_ transient_series. pdf.

16. Садовникова Н. А, Шмойлова Р. А. Анализ временных рядов и прогнозирование: учеб. пособие. М.: МЭСИ, 2001. 67 с.

17. Солодчук А. А. Исследование суточного хода геоакустической эмиссии на озере Микижа в период 2006-2007 гг. URL: http://www. kscnet. ru/ivs/ publication/young_conf/2008/2/art15.pdf.

18. Butter F. A. G. den, Fase M. M. G. Seasonal Adjustment as a Practical Problem. N.-Y.: Elsevier Science Publisher, 1991.

19. Gomez V., Maravall V. Estimation, Prediction, and Interpolation for Nonstationary Series With the Kalman Filter // Journal of the American Statistical Association. June 1994. Vol. 87. No. 426.