CC BY
333
№ 2(32) 2011
А. Н. Порунов, канд. экон. наук, лаборатория стратегических исследований и операционного проектирования при Самарском государственном техническом университете
Методика приведения ненормально распределенного ряда к нормальному распределению и оценка методической ошибки
В статье рассматриваются методика реализации в среде Mathcad ненормально распределенного макроэкономического ряда к нормально распределенному на основе преобразования Бокса-Кокса и возникающие при этом ошибки в оценке «нормальности» распределения.
Введение
В большинстве случаев экономисту-аналитику приходится иметь дело со статистическими данными, которые по тем или иным причинам не проходят тест на нормальность. В этой ситуации есть два выхода: либо обратиться к непараметрическим методам, что весьма проблематично для экономиста, поскольку требует особой математической подготовки, либо воспользоваться специальными методами, позволяющими преобразовать исходную «ненормальную статистику» в «нормальную», что само по себе также непросто.
Широко распространено мнение, что если данных много (например, п >100) или исследуются переменные, значения которых определяются бесконечным числом независимых факторов, то не имеет смысла использовать непараметрические статистики, и в этой ситуации лучше обратиться к методам трансформации ненормально распределенных данных в нормально распределенные. Среди множества таких методов преобразований одним из лучших (при неизвестном типе распределения)считается Бокса-Кокса преобразование.
Авторы этого преобразования — известные статистики Джордж Эдвард Пелхэм Бокс (George Edward Pelham Box), профессор Висконсинского университета в г. Мэдисон (США), и Дэвид Роксби Кокс (Sir David
РохЬее Сох), профессор колледжа Бирбека Лондонского университета. Впервые суть предлагаемого метода была изложена ими в 1964 г. в журнале Королевского статистического общества (йВ) [1]. Практические аспекты Бокса-Кокса (БК) преобразования на сегодняшний день достаточно подробно изложены в специальной англоязычной литературе [2-6], чего нельзя сказать об отечественной. Рассмотрим, насколько «всемогуще» БК преобразование в борьбе с «ненормально» распределенным макроэкономическим рядом и какие иллюзии могут возникнуть у исследователя-экономиста в зависимости от степени его «статистической испорченности» при оценке согласия функций эмпирического и теоретического распределений.
Бокса-Кокса преобразование
Пусть некоторая совокупность X представлена вектором непрерывных данных х(, / е 1,... N. Бокса-Кокса (БК) преобразование определяется следующим образом:
х(Х) =
х1
X
, X ф о
(1)
ln(À), Х = 0
Выражение (1) — это универсальное параметрическое семейство преобразований, которое экономисты часто используют в ал-
№ 2(32) 2011
горитмах сезонной (циклической) корректировки для того, чтобы сезонная составляющая преобразованного динамического ряда стала (хотя бы в первом приближении) неэволюционирующей по амплитуде, что упрощает ее последующую идентификацию [3]. Тиражируемые в литературе по экономической статистике и по этой причине популярные среди экономистов, логарифмическое и степенное преобразования — лишь частный случай преобразования БК. Так, например, в зависимости от значений X получаем: при X = 0 — логарифмическое, при Х<>2 — степенное.
Один из способов выбрать оптимальное значение X — использование X, максимизирующей логарифм функции правдоподобия. Логарифм функции правдоподобия следующий:
f (х, X) = -1 • In
£ (X, (X) - X (X))2
i=1
N
(2)
+(*-1) -Xln<
1 n
где x(X) = — • ^x( (X) — среднеарифметиче-
N ,=1
екая БК преобразованных данных.
Поскольку изначально БК преобразование было ориентировано только на положительные величины, проблему учета отрицательных значений данных снимают, добавляя к исходным значениям некоторое смещение, переводящее все отрицательные величины в положительную область1:
х(Х) =
(х + смещение f -1
-, Хф О
х -■-■ О)
ln( х + смещение), Х = О
При этом должно выполняться условие (х + смещение) > 0 Vx( е X.
Доверительная оценка X (с использованием статистики отношения правдоподобия) может быть произведена следующим образом:
f(х, X) > f(х, ¡V) - 0.5%^
(4)
где X — оценка максимального правдоподобия для X;
%2а>1— верхняя 100(1 - а) процентиль хи-квадрат распределения с 1 степенью свободы.
Практическая реализация
Для иллюстрации процедуры БК преобразования в среде Mathcad2 использовался макроэкономический ряд ВВП РФ — ряд X (табл. 1).
С помощью мастера импорта данных (Data Import Wizard) файл исходных данных, созданный в Excel, импортируется в Mathcad:
N = rows{X); i = 1... N\ N = 127,
где rows (X) — встроенная mathcad-функция, возвращающая число строк N (элементов) в векторе X (рис. 1).
Для нахождения уравнения тренда (в случае экспоненциальной зависимости) воспользуемся стандартной, встроенной
1 Таким образом, получается двухпараметрическое семейство преобразований, которое называется преобразованием Бокса-Кокса.
2 В большинстве современных математических пакетов сдвиг на константу (смешение) не предусмотрен, т.е. используется алгоритм более простого однопара-метрического преобразования.
ПРИКЛАДНАЯ ИНФОРМАТИКА /-
1 № 2(32) 2011
Таблица 1
Динамика уровней ВВП РФ в 1885-2009 гг.
1 1 1 1
1885 76 1917 143 1949 301 1981 1440
1886 73 1918 116 1950 374 1982 1423
1887 80 1919 92 1951 440 1983 1477
1888 86 1920 77 1952 453 1984 1664
1889 79 1921 74 1953 476 1985 1661
1890 75 1922 69 1954 483 1986 1668
1891 65 1923 64 1955 536 1987 1666
1892 93 1924 82 1956 569 1988 1754
1893 92 1925 98 1957 610 1989 1763
1894 95 1926 121 1958 616 1990 1784
1895 106 1927 146 1959 692 1991 1745
1896 93 1928 162 1960 721 1992 1735
1897 105 1929 173 1961 691 1993 1716
1898 94 1930 152 1962 789 1994 1518
1899 89 1931 175 1963 830 1995 1282
1900 90 1932 166 1964 818 1996 1267
1901 87 1933 171 1965 849 1997 1141
1902 86 1934 208 1966 958 1998 1119
1903 99 1935 242 1967 970 1999 1156
1904 95 1936 293 1968 1020 2000 1068
1905 114 1937 289 1969 1062 2001 1111
1906 98 1938 295 1970 1047 2002 1173
1907 88 1939 333 1971 1086 2003 1332
1908 89 1940 359 1972 1203 2004 1585
1909 108 1941 382 1973 1273 2005 1746
1910 111 1942 344 1974 1218 2006 1860
1911 123 1943 225 1975 1253 2007 2001
1912 107 1944 202 1976 1349 2008 2271
1913 118 1945 217 1977 1420 2009 2074
1914 134 1946 194 1978 1469
1915 158 1947 225 1979 1466
1916 160 1948 280 1980 1424
* Таблица составлена автором с учетом современных границ РФ.
№ 2(32) 2011
тах(Х) 2000
X
1000
1900
1950 Г
2000
Рис. 1. Динамический ряд X и тренд
в МаИпсас!3 функцией ехрX,д). Она возвращает вектор, содержащий три коэффициента экспоненциальной кривой вида:
а • ехр(Ь • х) + с
и наилучшим образом аппроксимирует данные в векторах Г и у. Необязательный вектор д содержит начальное приближение для этих трех коэффициентов:
9 =
'0.0014 0.001 ч0.001у
с = ехр Ш (Г, X, д)\ '0.00000000075 с = 0.0143518673
-508.15140440551, X?гепс! = с1 • ехр(с2 • Г) + с3.
Для приведения ряда к стационарному виду вычитают из ряда X найденный тренд Хйепб и определяют ряд остатков ДЯ (рис. 2):
ДЯ = х - XIгепб.
Для проверки близости распределения ряда остатков к нормальному построим гистограмму распределения Н (рис. 3):
Н = ¡п&одгат^гипс^ -1), ДЯ),
где 11Шодгат(Кипс(4^ -1), ДЯ) —функция, возвращающая матрицу Н из двух столбцов, содержащую средние точки Гшлс>/-1 подынтервалов. Результирующая матрица содержит иипсуЩ -1 строк. Причем 1гипс4ы -1— функция, возвращающая целую часть выражениял/м -1, удаляя дробную часть.
тах( ДЯ)
ДЯ
тю( ДЯ)
-500
-1000
Рис. 2. Динамика ряда остатков ДЯ
3 Использовалась последняя модифицированная версия пакета МаШсасМД М-035.
№ 2(32) 2011
1 ЭЙ
Рис. 3. Гистограмма распределения ряда остатков
Как видно из гистограммы на рис. 3, характер распределения ряда остатков далек от нормального. Как показывает практика, может оказаться, «...что преобразование квадратного корня еще слабое (не поджимает справа «хвост» распределения), а логарифмическое — уже слишком сильное («хвостик» появляется слева). Раньше пришлось бы выбирать из этих двух, но преобразование Бокса-Кокса в данном случае {X между 0 и 0,5) найдет промежуточное решение. Поэтому, если истинное нормализующее преобразование неизвестно, преобразование Бокса-Кокса считается лучшим» [7].
Поскольку БК преобразование применяется только к положительным уровням ряда, выберем величину смещения так, чтобы (АИ+смещение) > 0 при любых значениях ряда остатков ЛЯ. Примем величину смещения несколько большей (для наглядности — на 20%) минимального значения в ряду остатков ЛЯ: смещение = 1,2 • 1Гнп(ЛЯ).
Тогда новый ряд остатков ЛЯд с учетом смещения будет равен:
ЛЯд = ЛЯ -1,2 • 1Гнп(ЛЯ),
где 1Гнп(ЛЯ) — функция, возвращающая наименьшее из значений ЛЯ.
Пусть показатель степени изменяется в пределах Х = -1,-1 + 0,1...15 с шагом 0,1, тогда логарифмическую функцию правдо-
подобия РР (ЛЯд, X) можно определить следующим образом:
РР(ЛЯдД) =
■ 1п
(ЛЯд, 1 _ А (АЯд, -1 X N 'Ь X
N
+(Х-1) -ХИЛЯд.
Для того чтобы найти оптимальное значение Хор„ итеративно подставляем значения X, при которых логарифмическая функция правдоподобия РР (ЛЯ, X) достигает максимума. Ориентируясь по графику логарифмической функции правдоподобия (рис. 4), выберем «вилку» из значений:
РР(ЛЯд, 1,48) = - 682,903;
РР(ЛЯд, 1,49) = - 682,902;
РР(ЛЯд, 1,50) = - 682,903.
Промежуточное значение РР (ЛЯ, 1,49) соответствует максимуму функции РР (ЛЯ, X), т.е. в данном случае Xopt = ^,49.
Тогда преобразованный ряд остатков ВС будет определяться по формуле:
ВС =
ЛЯд149 -1 1,49
№ 2(32) 2011
-500
РР(ДЯд, X)
-1000
-500 -600 -700 -800 -900 -1x103
1, i9
10
15
Рис. 4. График логарифмической функции правдоподобия
Определим еще один ряд ВЭ, получаемый в результате сортировки ряда остатков ВС:
ВБ = эоП{ВС),
где зой{ВС) — функция, возвращающая вектор со значениями из ВС, упорядоченными по возрастанию.
Это позволит автору отразить кривую плотности нормального распределения на гистограмме (рис. 5):
Н = ¡11з1одгат{1гипс{4ы) -1, ВС).
в принятых обозначениях будет иметь следующий вид:
Ып(ВЗ) =
1
•ехр
-1
2
'ВЭ - теап(ВЗ)л2
81с1еу(ВЗ)
где теап{ВЗ) — функция, возвращающая арифметическое среднее (среднее значение) элементов ряда бв; Stdev(BS) —функция, возвращающая сред-неквадратическое отклонение совокупности элементов Вв.
На гистограмме, изображенной на рис. 5, видно, что характер распределения ос-Классическая форма функции плотно- татков после преобразования по методу сти нормального распределения (гаусин) Бокса-Кокса близок к нормальному. «За-
тах(Н<2>)+10 30
н<2> 20 0,45-Л/-Л/п(Бб)
10000 20000 т/'п(Н<1>) Н<1>, Вв тах(Н<л>)
Рис. 5. Гистограмма ряда остатков после БК преобразования
8
№ 2(32) 2011
I
N • •
cnorm(BSn)
-3-2-10 1 2 3
-2,4 BSn, 2,1
Рис. 6. Графики эмпирической и теоретической функций распределения
быв» о критериях согласия4, оценим ряд остатков на нормальность распределения на основе показателей эксцесса и асимметрии. Коэффициент асимметрии:
skew (ВС) = -0,0334,
где skew(BC)—функция, возвращающая асимметрию элементов ВС.
Коэффициент эксцесса:
kurt{BC) = -0,01163,
где kurt(BC) — функция, возвращающая асимметрию элементов ВС.
Рассчитаем вспомогательные величины оА и аЕ:
а А =
6 • (N - 2) (N +1) • (N + 3)
= 0,2123;
оЕ =
24 • N ■ (N - 2) • (N - 3) (N +1)2 • (N + 3) • (N + 5)
= 0,4099.
Для ряда с распределением, близким к нормальному, должны выполняться следующие условия:
\skew{BC)\ = 0,0334 < 1,5 • о • /4 = 0,3185
kurt{BC) -
N+1
= 0,16 < 1,5-о-Е = 0,6149.
В данном случае условия выполняются.
Продолжим проверку. С этой целью проведем очень популярный сегодня у экономистов визуальный анализ нормальности. Стандартизируем сортированный ранее ряд остатков Вв, предполагая, что справедлива гипотеза о нормальности ряда:
BSn =
BS - mean(BS) Stdev(BS)
4 Такой и аналогичные примеры «забывчивости» в последние годы нередки в кандидатских (и не только) диссертациях, подробнее об этом см. на сайте: http:// www.biometrica.tomsk.ru/konf/index2.htm.
Построим эмпирическую функцию распределения — и сравним ее с теоретическим распределением (рис. 6), используя встроенную mathcad-функцию enorm (BSn). Эта функция возвращает кумулятивное распределение вероятностей со средним, равным 0, и дисперсией, равной 1.
На графиках рис. 6 видна близость кривых распределения ^ и cnorm(BSn). На основе mathcad-функции qnorm(F,ц,о), возвращающей обратное кумулятивное нормальное распределение ряда F с заданными средним р и среднеквадратическим отклонением а, построим еще один график зависимости BSn (Nr,), (рис. 7). Пред-
№ 2(32) 2011
варительно определим: / = 1...N -1; ^ = —; стандартизированного ряда остатков:
Щ = цпогт(Р1,0,1). м
На первый взгляд может показаться, что и рис. 7 не дает оснований для беспокойства — большая часть точек стандартизированного ряда остатков ВЭп располагается очень близко к прямой, и поэтому распределение ряда можно считать нормальным. Подобные заключения нередки в работах, посвященных исследованию макроэкономических рядов. Но огорчает то, что множатся случаи, когда этим и ограничивается процедура проверки гипотезы о нормальности распределения. Между тем использование старого «доброго» критерия согласия Пирсона (в данном случае при N =127 это оправдано), критерия Колмогорова или омега-квадрат говорит об обратном. Покажем, такли это. Тем более, что МаШсас! позволяет сделать это достаточно просто (для понимания) и наглядно.
Для начала рассчитаем критерий Пирсона. Определим размах вариации
Я = ВБпы - ВБп, = 4,5.
Проведем группировку ряда, число групп:
К = 1гипс{4ы -1) = 10, к = 1... К.
Величинаинтервала: Л = — = 0,45. Середины интервалов:
тк = ВБпк +1 + (к -1) • Л;
тк = (-2,2 -1,71 -1,18 - 0,51 -0,04 0,56 1,02 2 2,48 2,97).
Рассчитаем теоретические частоты f¡í\
I = Л ■ N ■-
1 )2
^ = (2 5 11 20 23 20 14 3 1 0).
№ 2(32) 2011
Эмпирические частоты (используем определенные ранее данные для построения гистограммы (рис. 5):
Н{г) = (5 2 10 16 22 25 23 6 8 10).
Тогда расчетный критерий Пирсона %2 будет равен:
Г'* -К*) I2
При уровне значимости а=0,05 и числе степеней свободы в = К - 3 = 7 табличное значение критической точки правосторонней критической области %гкр = 14,2. Таким образом, эмпирические и теоретические частоты отличаются значимо.
Далее определим значения статистики Колмогорова:
KBSn, =
— - cnorm{BSni)
где спогт(ВЗп1) — таИ"юас1-функция, возвращающая кумулятивное распределение вероятностей со средним, равным 0, и дисперсией, равной 1.
Статистика Колмогорова:
0 = та х(КВЗп) = 0,06.
Расчетное значение статистики К? = ^• О = 0,71 при выбранном уровне значимости а=0.05 превышает табличное
1 36
значение^= = 0,12, это означает, что нулевую гипотезу следует отвергнуть, т. е. характер распределения ряда остатков далек от нормального, несмотря на проведенное ранее его БК преобразование.
Заключение
Практика статистических исследований показывает, «...что распределения реальных данных никогда не входят в какое-либо параметрическое семейство» [8]. На сегодняшний день в статистической литерату-
ре есть немало примеров, показывающих, || что распределения ошибок измерений поч- ^ ти всегда отличаются от нормальных [9]. Та- ^ кие семейства — лишь возможные прибли- ^ жения, заведомо неточные. Представленный ранее анализ конкретных данных приводит к аналогичному заключению.
Поэтому нельзя не согласиться с мнением одного из авторитетных отечественных статистиков профессора А. И. Орлова, что необходимо переходить от методов параметрической статистики к непараметрическим и робастным методам [8]. В первую очередь, по мнению автора, это относится к исследованию макроэкономических рядов. И экономистам об этом нужно помнить.
Описок литературы
1. Box G. Е. P., Сох D. R. An analysis of transformations. (With discussion) J. Roy. Statist. Soc. Ser. В 26. 1964. P. 211-252. http://www.ams.org/mathscinet-getitem?mr=192611.
2. Box-Cox Transformations: An Overview. Pengfei Li. Department of Statistics, University of Connecticut. Apr 11, 2005. http://www.stat.uconn.edu/~student-journai/index_fiies/pengfi_s05.pdf.
3. Carroll R. J., and Ruppert D. On prediction and the power transformation family. Biometrika 68: 609-615.
4. Box-Cox Transformation. http://www-stat.stanford. du/~olshen/manuscripts/selenite/node6.html.
5. Davidson, Russell, andJames G. MacKinnon. Estimation and Inference in Econometrics. Oxford University Press. 1993.
6. Definition of Box-Cox Transformation. http://eco-nomics.about.com/cs/economicsglossary/g/box_ cox.htm.
7. Приведение данных к нормальному распределению: преобразование Бокса-Кокса. Тематический форум, http://molbiol.ru/forums/index.php? showtopic=201368.
8. ОрловА. И. О критериях согласия с параметрическим семейством, http://www.newtech.ru/~orlov/ kritsogl.htm.
9. Мирвапиев М., Никулин М. Критерии согласия типа хи-квадрат // Заводская лаборатория. 1992. Т. 58. №3. С. 52-58.
11
