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

Язык труда и переводы:
УДК:
519.257
Дата публикации:
14 декабря 2022, 15:14
Категория:
Математическое моделирование физических процессов и технических систем
Авторы
Масягин Михаил Михайлович
МГТУ им. Н.Э. Баумана
Горяинов Владимир Борисович
МГТУ им. Н.Э. Баумана
Аннотация:
Цель работы — сравнение методов наименьших квадратов и наименьших модулей с точки зрения их эффективности при оценивании параметров экспоненциальной авторегрессионной модели. Она имеет огромное значение как в теоретических, так и прикладных математических задачах, поскольку, с одной стороны, позволяет с высокой точностью описывать даже сложные нелинейные процессы, а с другой, обладает удовлетворительной вычислительной сложностью для использования ее в компьютерных расчетах. Последнее особенно относится к модели первого порядка. В рамках исследования проведено масштабное испытание указанных методов оценивания параметров с использованием вычислительного эксперимента. В качестве тестовых распределений вероятности обновляющего процесса будут рассмотрены распределения Гаусса, Тьюки, Лапласа, Стьюдента, Коши и логистическое распределение.
Ключевые слова:
экспоненциальная авторегрессионная модель, метод наименьших квадратов, метод наименьших модулей, алгоритм Левенберга — Марквардта
Основной текст труда

На сегодняшний день экспоненциальная авторегрессионная модель [1] является одной из наиболее удачных и активно используемых на практике моделей для анализа временных рядов. С ее помощью может быть представлено множество физических, биологических и экономических нелинейных процессов. Например, качка корабля и последовательности импульсов в коре головного мозга [2], колебания курсов валют и риски изменения банками процентных ставок [3].

Обобщенная экспоненциальная авторегрессионная модель описывается следующим рекуррентным соотношением:

X_{t}=(a+\sum _{i=0}^{N-1}b_{i}X_{t-1}^{i}e^{-cX_{t-1}^{2}})X_{t-1}+\varepsilon _{t},                                                 (1)

где a\in \mathbb {R} , b_{i}\in \mathbb {R} ,\ i\in [0,\ N-1] , c\in \mathbb {R} — параметры модели;

\varepsilon _{t},\ t\in \mathbb {N} — последовательность независимых случайных величин, удовлетворяющая следующим условиям:

  • {\mathit {E\varepsilon }}_{t}=0 , т. е. она имеет нулевое математическое ожидание;
  • {\mathit {D\varepsilon }}_{t}=E\varepsilon _{t}^{2}=\sigma ^{2}<+{\infty } — конечная дисперсия.

Данную последовательность принято называть «обновляющим процессом».

С ростом порядка модели, увеличивается и вычислительная сложность компьютерного эксперимента, в котором она используется. В связи с этим на практике наиболее часто применяется модель первого порядка. При N=1 соотношение (1) принимает следующий вид:

X_{t}=(a+be^{-cX_{t-1}^{2}})X_{t-1}+\varepsilon _{t}.                                               (2)

Предположим, что имеется процесс X_{t} , удовлетворяющий уравнению (2) и имеются n его наблюдений: X_{1,}X_{2,}...,X_{n-1},X_{n}.  В таком случае параметры уравнения можно оценить методами наименьшних квадратов (МНК) и наименьшних модулей (МНМ) ({\dot {a}},{\dot {b}},{\dot {c}}) , представляющих из себя решения задач минимизации функций (3) и (4) соответственно:

g_{\mathit {l.s.m}}(a,b,c)=\sum _{t=2}^{n}(X_{t}-(a+be^{-cX_{t-1}^{2}})X_{t-1})^{2};                                                             (3)

g_{\mathit {l.a.m.}}(a,b,c)=\sum _{t=2}^{n}\left|X_{t}-(a+be^{-cX_{t-1}^{2}})X_{t-1}\right|.                                                        (4)

Точность оценки {\widehat {\theta }} скалярного параметра \theta распределения вероятности произвольной случайной величины логично измерять как абсолютную величину её отклонения от оцениваемого параметра: \left|{\widehat {\theta }}-\theta \right| . С учетом того, что величина \left|{\widehat {\theta }}-\theta \right| случайна вследствие случайности самой оценки {\widehat {\theta }} , то имеет смысл заменить её на E\left|{\widehat {\theta }}-\theta \right| .

На практике удобнее работать со всюду дифференциируемыми функциями, поэтому E\left|{\widehat {\theta }}-\theta \right| в свою очередь может быть заменено на E({\widehat {\theta }}-\theta )^{2} .

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

{\underset {{\widehat {\theta _{i}}},\ i\in \mathbb {N} }{\mathit {argmin}}}E({\widehat {\theta _{i}}}-\theta )^{2},

где {\widehat {\theta _{i}}},\ i\in \mathbb {N} — оценки исследуемого скалярного параметра.

Таким образом, вводится возможность качественного сравнения оценок. Количественное же сравнение оценок может быть выражено соотношением (5). Будем называть его относительной эффективностью {\widehat {\theta _{i}}} по отношению к {\widehat {\theta _{j}}} :

e({\widehat {\theta _{i}}},{\widehat {\theta _{j}}})={\frac {E({\widehat {\theta _{j}}}-\theta )^{2}}{E({\widehat {\theta _{i}}}-\theta )^{2}}},\ i\neq j.                                                           (5)

В качестве тестовых распределений \varepsilon _{t} используем распределения Гаусса (6), Тьюки (7), Лапласа (8), Стьюдента (9), Коши (10) и логистическое распределение (11). Конкретные их значения получим с помощью датчика случайных чисел библиотеки NumPy языка программирования Python.

f(x)={\frac {1}{\sqrt {2\pi }}}e^{\frac {-x^{2}}{2}};                                                                (6)

f(x)=(1-\gamma ){\frac {1}{\sqrt {2\pi }}}e^{\frac {-x^{2}}{2}}+\gamma {\frac {1}{\sqrt {2\pi }}}e^{\frac {-x^{2}}{2\tau ^{2}}},\ \tau >1,\ 0<\gamma <1;                                         (7)

f(x)={\frac {1}{{\sqrt {2}}\sigma }}e^{\frac {-{\mathit {sqrt2}}\left|x\right|}{\sigma }};                                                 (8)

f(x)={\frac {\Gamma ({\frac {m+1}{2}})}{{\sqrt {m\pi }}\Gamma ({\frac {m}{2}})(1+{\frac {x^{2}}{m}})^{\frac {m+1}{2}}}};                                                         (9)

f(x)={\frac {1}{\pi (1+x^{2})}};                                                         (10)

f(x)={\frac {e^{-x}}{(1+e^{-x})^{2}}}.                                                                   (11)

Вектор наблюдений X_{1,}X_{2,}...,X_{n-1},X_{n} сгенерируем 1000 раз со следующими начальным условием и параметрами:

X_{0}=0,\ a=-0,3,\ b=-0,8,\ c=1,\ n=50,\ 100,\ 200,\ 500.

Важно отметить, что подобный выбор параметров a , b и c , и распределений вероятности обновляющего процесса не совсем случаен. Он позволяет согласно [4] считать, что модель (2) является стационарной, так как выполняются условия:

  • \left|a\right|<1 ;
  • c>0 ;
  • {\mathit {D\varepsilon }}_{t}=\sigma ^{2}<+{\infty } .

Минимизацию функций оценок параметров будем производить с помощью оптимизационного алгоритма Левенберга-Марквардта [5], представляющего из себя удачную комбинацию двух других методов: Гаусса-Ньютона и наивного градиентного спуска. Направление его поиска определяется следующим матричным уравнением:

где J({\vec {x}}_{k}) — якобиан минимизируемой функции {\vec {f}}({\vec {x}}_{k}) ,

\lambda _{k}\geq 0 — параметр регуляризации, индивидуальный для каждой итерации алгоритма,

I — единичная матрица,

{\vec {p}}_{k} — вектор изменения {\vec {x}}_{k} : {\vec {x}}_{k+1}={\vec {x}}_{k}+{\vec {p}}_{k} .

Для проведения вычислительного эксперимента была использована реализация метода Левенберга-Марквардта из библиотеки SciPy языка программирования Python.

В силу того что оценки и ({\dot {a}},{\dot {b}},{\dot {c}}) векторные, их сравнение происходило независимо покоординатно. Так, например, относительная эффективность параметра b рассчитывалась следующим образом:

e_{N}({\widetilde {b}},{\dot {b}})={\frac {E({\dot {b}}-b)^{2}}{E({\widetilde {b}}-b)^{2}}}={\frac {\sum _{i=1}^{N}({\dot {b_{i}}}-b)^{2}}{\sum _{i=1}^{N}({\widetilde {b_{i}}}-b)^{2}}}.

Для краткости приведем таблицу относительной эффективности только для параметра b .

 

Оценка e_{N}({\widetilde {b}},{\dot {b}}) относительной эффективности МНК по отношению к МНМ при разлчиных распределениях \varepsilon _{t} и числе шагов n

Распределение ε(t)

Число шагов n

50

100

200

500

Гаусса

1,667

1,553

1,473

1,531

Тьюки:

γ = 0,01, τ = 3

γ = 0,01, τ = 10

γ = 0,1, τ = 3

γ = 0,1, τ = 10

 

1,581

0,772

1,017

0,169

 

1,561

0,810

1,011

0,163

 

1,452

0,779

0,952

0,159

 

1,435

0,879

1,052

0,173

Лапласа

0,676

0,634

0,609

0,561

Стьюдента:

m = 15

m = 10

m = 5

m = 4

 

1,314

1,329

1,023

0,900

 

1,473

1,371

1,080

0,908

 

1,379

1,401

1,103

0,879

 

1,366

1,313

0,989

0,911

Коши

3,551 · 10–5

1,831 · 10–4

5,333 · 10–5

1,191 · 10–5

Логическое

1,215

1,323

1,185

1,301

 

В таблице демонстрируется, что МНК однозначно превосходит МНМ в случаях нормального и логистического распределений. Также он оказывается лучше МНМ для распределения Стьюдента при , , и для распределения Тьюки при . Из этого может быть сделан вывод, что МНК эффективнее МНМ для распределения Стьюдента при больших значениях параметра () и для распределения Тьюки при малых значениях параметра . В остальных случаях МНМ превосходит МНК. При этом погрешность методики тестирования уменьшается пропорционально {\frac {1}{\sqrt {N}}},  что в совокупности с погрешностями алгоритма Левенберга — Марквардта обуславливает высокую степень разброса значений в таблице.

В следующей работе планируется дополнить сравнение методом наименьших функций, который может быть выражен следующей формулой:

g_{\mathit {l.f.}}(a,b,c)=\sum _{t=2}^{n}\rho (X_{t}-(a+be^{-cX_{t-1}^{2}})X_{t-1}),

где \rho , например, является функцией потерь Хьюбера для устойчивой регрессии:

\rho _{\delta }(x)=\left\{{\begin{matrix}{\frac {1}{2}}x^{2,}\ \ \ \ \ \ \ \ {\mathit {if}}\left|x\right|\leq \delta ;\\\delta (\left|x\right|-{\frac {1}{2}}\delta ),\ {\mathit {if}}\left|x\right|>\delta .\end{matrix}}\right.

Литература
  1. Ozaki T. Non-linear time series models and dynamical systems. Elsevier Science Publishers B.V., 1985, pp. 25–83.
  2. Ozaki T. Time series modeling of neuroscience data. CRC Press, 2012. DOI: https://doi.org/10.1201/b11527
  3. Olugbode M., El-Masry A., Pointon J. Exchange rate and interest rate exposures of UK industries using first-order autoregressive exponential GARCH-in-mean (EGARCH-M) approach. The Manchester School, 2014, vol. 82, iss. 4, pp. 409–464. DOI: https://doi.org/10.1111/manc.12029
  4. Allal J., El Melhaoui S. Optimal detection of exponential component in autoregressive models. J Time Ser Anal, 2011, vol 27, iss. 6, pp. 793–810.
  5. Rhinehart R.R. Nonlinear regression modeling for engineering applications: modeling, model validation, and enabling design of experiments. Wiley, 2016.
Ваш браузер устарел и не обеспечивает полноценную и безопасную работу с сайтом.
Установите актуальную версию вашего браузера или одну из современных альтернатив.