14.10.2022 15:48
ANazarov 2 1800 Источник

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

Критерий Дарбина-Уотсона (Durbin–Watson statistic) - один из самых распространенных критериев для проверки автокорреляции.

Данный критерий входит в стандартный инструментарий python:

  • присутствует в таблице выдачи результатов регрессионного анализа модуля линейной регрессии Linear Regression;

  • может быть рассчитан с помощью функции statsmodels.stats.stattools.durbin_watson.

К сожалению, стандартные инструменты python не позволяют получить табличные значения статистики критерия Дарбина-Уотсона, нам предлагается воспользоваться методом грубой оценки: считается, что при расчетном значении статистики критерия в интервале [1; 2] автокорреляция отсутствует (см. Durbin–Watson statistic). Однако, для качественного статистического анализа такой подход неприемлем.

Представляет интерес реализовать в полной мере критерий Дарбина-Уотсона средствами python, добавив этот важный критерий в инструментарий специалиста DataScience.

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

Применение пользовательских функций

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

Вот перечень данных функций:

  • graph_plot_sns_np - функция строит линейный график средствами seaborn;

  • graph_regression_plot_sns - функция строит график регрессионной модели и график остатков средствами seaborn;

  • regression_error_metrics - функция возвращает ошибки аппроксимации регрессионной модели;

  • graph_hist_boxplot_probplot_sns - функция позволяет визуализировать исходные данные для одной переменной путем одновременного построения гистограммы, коробчатой диаграммы и вероятностного графика средствами seaborn; имеется возможность выбирать, какие графики строить (h - hist, b - boxplot, p - probplot);

  • norm_distr_check - проверка нормальности распределения исходных данных с использованием набора из нескольких статистических тестов.

  • Goldfeld_Quandt_test, Breush_Pagan_test, White_test - проверка гетероскедастичности с использование тестов Голдфелда-Квандта, Бриша-Пэгана и Уайта соответственно;

  • graph_regression_pair_predict_plot_sns - прогнозирование: построение графика регрессионной модели (с доверительными интервалами) и вывод расчетной таблицы с данными для заданной области значений X.

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

Основы теории

Информацию о критерии Дарбина-Уотсона можно почерпнуть в [1, с.659], [2, с.117], [3, с.239], [4, с.188], а также:

Итак, предположим, мы рассматриваем регрессионную модель:

y_i = b_0 + \displaystyle\sum_{j=1}^{m} b_j x_{ij} + \varepsilon_i, \quad i=\overline{1,n}\widehat{y}_i = b_0 + \displaystyle\sum_{j=1}^{m} b_j x_{ij}, \quad i=\overline{1,n}

или в матричном виде:

\textbf{Y} = \textbf{X}b + \varepsilon\widehat{\textbf{Y}} = \textbf{X}b

Критерий Дарбина-Уотсона применяется в ситуации, когда регрессионные остатки связаны автокорреляционной зависимостью 1-го порядка [2, с.111]:

\varepsilon_i = \rho \varepsilon_{i-1} + \delta_i, \quad i=\overline{2,n}

где \rho - некоторое число (|\rho|<1), а случайные величины \delta_i удовлетворяют требованиям, предъявляемым к регрессионным остаткам классической модели (т.е. равенство нулю среднего значения, постоянство дисперсии и некоррелированность между собой):

E(\delta_i)=0E(\delta_i\delta_j)=\begin{cases}{\sigma_0}^2, & \quad i=j \\0, & \quad i \neq j\end{cases}

Проверяется нулевая гипотеза об отсутствии автокорреляции:

\begin{aligned} \qquad H_0: \rho = 0 \end{aligned}

Альтернативной гипотезой может быть:

  • существование отрицательной автокорреляции (левосторонняя критическая область):

\begin{aligned} \qquad H_1: \rho < 0 \end{aligned}

  • существование положительной автокорреляции (правосторонняя критическая область):

\begin{aligned} \qquad H_1: \rho > 0 \end{aligned}

  • существование автокорреляции вообще (двусторонняя критическая область):

\begin{aligned} \qquad H_1: \rho \neq 0 \end{aligned}

Расчетное значение статистики критерия Дарбина-Уотсона имеет вид:

DW_{calc} = \frac {\displaystyle\sum_{i=2}^{n} (\varepsilon_i - \varepsilon_{i-1})^2} {{\displaystyle\sum_{i=1}^{n} {\varepsilon_i}^2}}

где \varepsilon_i = y_i - \widehat{y}_i - остатки (невязки) регрессионной модели.

По таблицам (см. [1, с.659], [2, с.402], [3, с.291]) в зависимости от уровня значимости \alpha_{level} (5%, 2.5%, 1%), числа параметров регрессионной модели m (кроме свободного члена b_0) (от 1 до 5) и объема выборки n (от 15 до 100) определяются критические значения статистики Дарбина-Уотсона: нижний d_L и верхний d_U предел.

Правила принятия гипотез по критерию Дарбина-Уотсона выглядят довольно своеобразно - критические значения образуют пять областей различных статистических решений (причем критические границы принятия H_0 и непринятия H_1 не совпадают):

Значение DW_{calc}

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

Вывод

0 \leq DW_{calc} < d_L

отвергается H_0, принимается H_1: \rho > 0

есть положительная автокорреляция

d_L \leq DW_{calc} \leq d_U

неопределенность

d_U < DW_{calc} < 4-d_U

принимается H_0

автокорреляция отсутствует

4-d_U \leq DW_{calc} \leq 4-d_L

неопределенность

4-d_L < DW_{calc} \leq 4

отвергается H_0, принимается H_1: \rho < 0

есть отрицательная автокорреляция

Есть очень удачная мнемоническая схема, приведенная в [3, с.240]:

Особенности критерия Дарбина-Уотсона:

  1. Критические значения критерия табулированы для объема выборки от 15 до 100, аппроксимаций мне обнаружить не удалось. При меньших значениях критерий применять нельзя, при больших - очевидно, приходиться пользоваться грубым оценочным правилом: при расчетном значении статистики критерия в интервале [1; 2] автокорреляция отсутствует (см. https://en.wikipedia.org/wiki/Durbin–Watson_statistic).

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

  3. Критерий построен в предположении, что регрессоры \textbf{X} и ошибки \varepsilon не коррелированы, поэтому его нельзя применять, в частности, для моделей авторегрессии [4, с.191].

  4. Критерий не подходит для моделей без свободного члена b_0.

  5. Критерий имеет зону неопределенности, когда нет оснований ни принимать, ни отвергать нулевую гипотезу.

  6. Между статистикой критерия и коэффициентом автокорреляции существует приближенное соотношение:

DW_{calc} \approx 2 \cdot (1-\rho)

Существуют и другие критерии для проверки автокорреляции (тест Бройша-Годфри, Льюнга-Бокса и пр.).

Как было указано выше, большой проблемой является отсутствие табличных значений статистики критерия Дарбина-Уостона в стандартном инструментарии python. Для реализации возможностей данного критерия в полном объеме нам потребуется оцифровка весьма объемных таблиц критических значений.

Оцифровка табличных значений статистики критерия Дарбина - Уотсона

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

Однако, если мы хотим выполнять качественный статистический анализ, неизбежно придется работать с большим количеством статистических критериев и далеко не все из них реализованы в python. Критерий Дарбина-Уотсона - это только один из многих. Количество критериев, рассматриваемых в литературе по прикладной статистике в последние годы постоянно увеличивается. Специалисту придется реализовывать многие критерии самостоятельно и одна из проблем, с которой придется столкнуться - это таблицы критических значений. Далеко не все табличные значения имеют аппроксимации, а значит придется каким-то образом оцифровывать эти таблицы. Небольшие таблицы можно сохранить в файлах вручную, а вот такой подход с объемными таблицами (как в нашем случае) - это слишком непроизводительно и нерационально.

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

Замечу сразу, что я не являюсь глубоким специалистом в области анализа и обработки изображений и текстов на python - это не совсем мой профиль. Профессионалы в этой области, возможно, раскритикуют то, как решается поставленная задача и предложат более удачное решение. Если будет так - то заранее спасибо. Я же эту задачу старался решить наиболее простым и рациональным способом, доступным для широкого круга специалистов. На всякий случай могу процитировать Давоса Сиворта из "Игры престолов": "Простите за то, что увидите".

Алгоритм действий:

Для оцифровки я использовал таблицы, приведенные в [3, с.290-292].

  1. Сканируем таблицы, сохраняем в виде jpg-файлов (Durbin_Watson_test_1.jpg, Durbin_Watson_test_2.jpg, Durbin_Watson_test_3.jpg) в папке text_processing, расположенной внутри папки с рабочим .ipynb-файлом:

  1. Распознаем текст (я воспользовался онлайн-сервисом https://convertio.co/), полученные текстовые файлы Durbin-Watson-test-1.ocr.txt, Durbin-Watson-test-2.ocr.txt, Durbin-Watson-test-3.ocr.txt также помещаем в папке text_processing.

  2. Откроем файлы, запишем содержимое файлов в переменные, каждая из которых соответствует одной странице:

with open('text_processing\Durbin-Watson-test-1.ocr.txt') as f1:
    Durbin_Watson_test_1 = f1.readlines()
display(Durbin_Watson_test_1, type(Durbin_Watson_test_1), len(Durbin_Watson_test_1))

С остальными файлами - действуем аналогично:

with open('text_processing\Durbin-Watson-test-2.ocr.txt') as f2:
    Durbin_Watson_test_2 = f2.readlines()
display(Durbin_Watson_test_2, type(Durbin_Watson_test_2), len(Durbin_Watson_test_2))

with open('text_processing\Durbin-Watson-test-3.ocr.txt') as f3:
    Durbin_Watson_test_3 = f3.readlines()
display(Durbin_Watson_test_3, type(Durbin_Watson_test_3), len(Durbin_Watson_test_3))

Видим, что переменные представляют собой списки, элементами которых является строки.

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

Durbin_Watson_test = [Durbin_Watson_test_1, Durbin_Watson_test_2, Durbin_Watson_test_3]

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

  1. Исключаем все строки, которые начинаются не с цифр; при этом воспользуемся алгоритмом перезаписи списка:

# создаем новый список
Durbin_Watson_test_new = list()

# удаляем строки
for page in Durbin_Watson_test:
    page_temp = list()    # временная страница
    for line in page:
        if line[0].isdigit():
            page_temp.append(line)    # перезаписываем список
    Durbin_Watson_test_new.append(page_temp)

  1. Исключаем из текста управляющие символы (\t, \n) - с помощью регулярных выражений (regex) (модуль re):

# задаем шаблон для удаления символов
pattern = r'[\t+\n+]'    

# выполняем обработку
Durbin_Watson_test_new = [
    [re.sub(pattern, '  ', elem) for elem in page] 
    for page in Durbin_Watson_test_new]

  1. Удаляем все символы, кроме цифр, точек, запятых и пробелов:

# задаем шаблон для удаления символов
pattern = r'[^0-9,. ]'    

# выполняем обработку
Durbin_Watson_test_new = [
    [re.sub(pattern, '', elem) for elem in page] 
    for page in Durbin_Watson_test_new]

  1. Заменяем запятые на точки:

# задаем шаблон для удаления символов
pattern = r'[,]'

# выполняем обработку
Durbin_Watson_test_new = [
    [re.sub(pattern, '.', elem) for elem in page] 
    for page in Durbin_Watson_test_new]

  1. Разделяем строки:

# задаем шаблон
pattern = r'[ ]+'

# выполняем обработку
Durbin_Watson_test_new = [[re.split(pattern, elem) for elem in page] 
                          for page in Durbin_Watson_test_new]

  1. Сохраняем данные в DataFrame - для этого создадим список Durbin_Watson_list_df, элементами которого являются отдельные DataFrame, каждый из которых соответствует отдельной странице:

# создаем новый список
Durbin_Watson_list_df = list()
for page in Durbin_Watson_test_new:
    Durbin_Watson_list_df.append(pd.DataFrame(page))

  1. Исправляем вручную отдельные аномалии, возникшие при распознавании отсканированных данных - к сожалению, работы вручную совсем избежать не удается.

  • Корректируем DataFrame, соответствующий 1-й странице:

# создаем временный DataFrame
temp_df = Durbin_Watson_list_df[0]
# удаляем последние столбцы
temp_df = temp_df.drop(columns=[11, 12])
# корректируем вручную отдельные ошибки
temp_df.loc[[4],] = [19, 1.18, 1.40, 1.08, 1.53, 0.97, 1.68, 0.86, 1.85, 0.75, 2.02]
temp_df.loc[[8],[3]] = 1.17
temp_df.loc[[10],[3]] = 1.21
temp_df.loc[[17],[9]] = 1.11
temp_df.loc[[21],[4]] = 1.59
temp_df.loc[[25],[5]] = 1.34
temp_df.loc[[31],[10]] = 1.77
# записываем изменения
Durbin_Watson_list_df[0] = temp_df

  • Корректируем DataFrame, соответствующий 2-й странице:

# создаем временный DataFrame
temp_df = Durbin_Watson_list_df[1]
# удаляем последние столбцы
temp_df = temp_df.drop(columns=[11, 12])
# корректируем вручную отдельные ошибки
temp_df.loc[[2],[8]] = 1.77
temp_df.loc[[10],[9]] = 0.86
temp_df.loc[[10],[10]] = 1.77
temp_df.loc[[14],[9]] = 0.96
temp_df.loc[[17],[10]] = 1.71
temp_df.loc[[34],[10]] = 1.71
# записываем изменения
Durbin_Watson_list_df[1] = temp_df

  • Корректируем DataFrame, соответствующий 3-й странице:

# создаем временный DataFrame
temp_df = Durbin_Watson_list_df[2]
# удаляем последние столбцы
temp_df = temp_df.drop(columns=[11, 12])
# корректируем вручную отдельные ошибки
temp_df.loc[[2],[9]] = 0.48
temp_df.loc[[13],] = [28, 1.10, 1.24, 1.04, 1.32, 0.97, 1.41, 0.90, 1.51, 0.83, 1.62]
temp_df.loc[[20],[3]] = 1.14
temp_df.iloc[21:26, 7] = [1.04, 1.06, 1.07, 1.09, 1.10]
temp_df.loc[[26],[9]] = 1.11
temp_df.loc[[35],] = [90, 1.50, 1.54, 1.47, 1.56, 1.45, 1.59, 1.43, 1.61, 1.41, 1.64]
# записываем изменения
Durbin_Watson_list_df[2] = temp_df

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

11. Преобразуем значения из строкового в числовой тип:

for elem_df in Durbin_Watson_list_df:
    for col in elem_df.columns:
        elem_df[col] = pd.to_numeric(elem_df[col], errors='ignore')

  1. Корректируем структуру DataFrame:

  • меняем индекс - индексом теперь будет объем выборки n

  • каждый DataFrame снабжаем мультииндексом по столбцам (подробнее см. [7, с.169])

# меняем индекс
Durbin_Watson_list_df = [
    elem_df.set_index([0])
    for elem_df in Durbin_Watson_list_df]

# добавляем мультииндекс по столбцам
multi_index_list = ['p=0.95', 'p=0.975', 'p=0.99']    # список, содержащий значения для верхней строки мульииндекса
for i, elem_df in enumerate(Durbin_Watson_list_df):
    elem_df.index.name = 'n'
    elem_df.columns = pd.MultiIndex.from_product(
        [[multi_index_list[i]],
         ['m=1', 'm=2', 'm=3', 'm=4', 'm=5'],
         ['dL','dU']])

  1. Объединяем отдельные DataFrame в один:

Durbin_Watson_test_df = Durbin_Watson_list_df[0].copy()

for i, elem_df in enumerate(Durbin_Watson_list_df):
    if i > 0:
        Durbin_Watson_test_df = Durbin_Watson_test_df.join(elem_df)

display(Durbin_Watson_test_df)        
Durbin_Watson_test_df.info()

Итак, мы сформировали DataFrame с оцифрованными данными таблиц критических значений статистики Дарбина-Уотсона. Получить доступ к данным теперь очень просто - например, нам требуется вывести табличные значения статистики критерия при объеме выборки n=40, доверительной вероятности p=0.95 и числе параметров регрессионной модели m=2:

n = 40
p = 0.95
m=2

Durbin_Watson_test_df.loc[[n], (f'p={p}', f'm={m}')]

  1. Построим график табличных значений.

График получился весьма объемным - 3х5 элементов - однако он необходим: на графике можно увидеть те ошибки (пики и впадины), которые мы могли пропустить при ручной обработке ранее (некорректно отсканированные и распознанные цифры), тогда придется вернуться к этапу 10.

# меняем настройки Mathplotlib
plt.rcParams['axes.titlesize'] = 10      # шрифт заголовка
plt.rcParams['legend.fontsize'] = 9      # шрифт легенды
plt.rcParams['xtick.labelsize'] = 8      # шрифт подписей меток
plt.rcParams['ytick.labelsize'] = 8

fig = plt.figure(figsize=(297/INCH, 420/INCH))

ax_1_1 = plt.subplot(5,3,1)
ax_2_1 = plt.subplot(5,3,2)
ax_3_1 = plt.subplot(5,3,3)
ax_1_2 = plt.subplot(5,3,4)
ax_2_2 = plt.subplot(5,3,5)
ax_3_2 = plt.subplot(5,3,6)
ax_1_3 = plt.subplot(5,3,7)
ax_2_3 = plt.subplot(5,3,8)
ax_3_3 = plt.subplot(5,3,9)
ax_1_4 = plt.subplot(5,3,10)
ax_2_4 = plt.subplot(5,3,11)
ax_3_4 = plt.subplot(5,3,12)
ax_1_5 = plt.subplot(5,3,13)
ax_2_5 = plt.subplot(5,3,14)
ax_3_5 = plt.subplot(5,3,15)

fig.suptitle('Табличные значения статистики критерия Дарбина-Уотсона', fontsize = 16)

(Ymin, Ymax) = (0.3, 2.2)
x = Durbin_Watson_test_df.index

title_fontsize = 10

name_1_1 = ['p=0.95', 'm=1']
ax_1_1.set_title(name_1_1[0] + ' ' + name_1_1[1])
ax_1_1.plot(x, Durbin_Watson_test_df[tuple(name_1_1 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_1_1 + ['dU'])])

name_1_2 = ['p=0.95', 'm=2']
ax_1_2.set_title(name_1_2[0] + ' ' + name_1_2[1])
ax_1_2.plot(x, Durbin_Watson_test_df[tuple(name_1_2 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_1_2 + ['dU'])])

name_1_3 = ['p=0.95', 'm=3']
ax_1_3.set_title(name_1_3[0] + ' ' + name_1_3[1])
ax_1_3.plot(x, Durbin_Watson_test_df[tuple(name_1_3 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_1_3 + ['dU'])])

name_1_4 = ['p=0.95', 'm=4']
ax_1_4.set_title(name_1_4[0] + ' ' + name_1_4[1])
ax_1_4.plot(x, Durbin_Watson_test_df[tuple(name_1_4 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_1_4 + ['dU'])])

name_1_5 = ['p=0.95', 'm=5']
ax_1_5.set_title(name_1_5[0] + ' ' + name_1_5[1])
ax_1_5.plot(x, Durbin_Watson_test_df[tuple(name_1_5 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_1_5 + ['dU'])])

name_2_1 = ['p=0.975', 'm=1']
ax_2_1.set_title(name_2_1[0] + ' ' + name_2_1[1])
ax_2_1.plot(x, Durbin_Watson_test_df[tuple(name_2_1 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_2_1 + ['dU'])])

name_2_2 = ['p=0.975', 'm=2']
ax_2_2.set_title(name_2_2[0] + ' ' + name_2_2[1])
ax_2_2.plot(x, Durbin_Watson_test_df[tuple(name_2_2 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_2_2 + ['dU'])])

name_2_3 = ['p=0.975', 'm=3']
ax_2_3.set_title(name_2_3[0] + ' ' + name_2_3[1])
ax_2_3.plot(x, Durbin_Watson_test_df[tuple(name_2_3 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_2_3 + ['dU'])])

name_2_4 = ['p=0.975', 'm=4']
ax_2_4.set_title(name_2_4[0] + ' ' + name_2_4[1])
ax_2_4.plot(x, Durbin_Watson_test_df[tuple(name_2_4 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_2_4 + ['dU'])])

name_2_5 = ['p=0.975', 'm=5']
ax_2_5.set_title(name_2_5[0] + ' ' + name_2_5[1])
ax_2_5.plot(x, Durbin_Watson_test_df[tuple(name_2_5 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_2_5 + ['dU'])])

name_3_1 = ['p=0.99', 'm=1']
ax_3_1.set_title(name_3_1[0] + ' ' + name_3_1[1])
ax_3_1.plot(x, Durbin_Watson_test_df[tuple(name_3_1 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_3_1 + ['dU'])])

name_3_2 = ['p=0.99', 'm=2']
ax_3_2.set_title(name_3_2[0] + ' ' + name_3_2[1])
ax_3_2.plot(x, Durbin_Watson_test_df[tuple(name_3_2 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_3_2 + ['dU'])])

name_3_3 = ['p=0.99', 'm=3']
ax_3_3.set_title(name_3_3[0] + ' ' + name_3_3[1])
ax_3_3.plot(x, Durbin_Watson_test_df[tuple(name_3_3 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_3_3 + ['dU'])])

name_3_4 = ['p=0.99', 'm=4']
ax_3_4.set_title(name_3_4[0] + ' ' + name_3_4[1])
ax_3_4.plot(x, Durbin_Watson_test_df[tuple(name_3_4 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_3_4 + ['dU'])])

name_3_5 = ['p=0.99', 'm=5']
ax_3_5.set_title(name_3_5[0] + ' ' + name_3_5[1])
ax_3_5.plot(x, Durbin_Watson_test_df[tuple(name_3_5 + ['dL'])],
            x, Durbin_Watson_test_df[tuple(name_3_5 + ['dU'])])

ax_1_1.set_ylim(Ymin, Ymax)
ax_2_1.set_ylim(Ymin, Ymax)
ax_3_1.set_ylim(Ymin, Ymax)
ax_1_2.set_ylim(Ymin, Ymax)
ax_2_2.set_ylim(Ymin, Ymax)
ax_3_2.set_ylim(Ymin, Ymax)
ax_1_3.set_ylim(Ymin, Ymax)
ax_2_3.set_ylim(Ymin, Ymax)
ax_3_3.set_ylim(Ymin, Ymax)
ax_1_4.set_ylim(Ymin, Ymax)
ax_2_4.set_ylim(Ymin, Ymax)
ax_3_4.set_ylim(Ymin, Ymax)
ax_1_5.set_ylim(Ymin, Ymax)
ax_2_5.set_ylim(Ymin, Ymax)
ax_3_5.set_ylim(Ymin, Ymax)

legend = (r'$d_L$', r'$d_U$')
ax_1_1.legend(legend)
ax_2_1.legend(legend)
ax_3_1.legend(legend)
ax_1_2.legend(legend)
ax_2_2.legend(legend)
ax_3_2.legend(legend)
ax_1_3.legend(legend)
ax_2_3.legend(legend)
ax_3_3.legend(legend)
ax_1_4.legend(legend)
ax_2_4.legend(legend)
ax_3_4.legend(legend)
ax_1_5.legend(legend)
ax_2_5.legend(legend)
ax_3_5.legend(legend)

plt.show()

# возвращаем настройки Mathplotlib
plt.rcParams['axes.titlesize'] = f_size + 10      # шрифт заголовка
plt.rcParams['legend.fontsize'] = f_size + 6      # шрифт легенды
plt.rcParams['xtick.labelsize'] = f_size + 4      # шрифт подписей меток
plt.rcParams['ytick.labelsize'] = f_size + 4

  1. Сохраняем полученный DataFrame в csv-файл, помещаем его в папку table, расположенную внутри папки с рабочим .ipynb-файлом (в которой папку table у нас хранятся файлы с данными из статистических таблиц):

Durbin_Watson_test_df.to_csv(
    path_or_buf='table\Durbin_Watson_test_table.csv',
    mode='w+',
    sep=';',
    index_label='n')

Табличные значения статистики критерия Дарбина-Уотсона у нас теперь имеются, можем приступать к созданию пользовательской функции.

Создание пользовательской функции для реализации критерия Дарбина - Уотсона

Рассчитать статистику критерия Дарбина-Уотсона мы можем с помощью функции statsmodels.stats.stattools.durbin_watson.

Создадим пользовательскую функцию Durbin_Watson_test для проверки гипотезы об автокорреляции:

def Durbin_Watson_test(
    data,
    m = None,
    p_level: float=0.95):
    
    a_level = 1 - p_level
    data = np.array(data)
    n = len(data)
    
    # расчетное значение статистики критерия
    DW_calc = sms.stattools.durbin_watson(data)
    
    # табличное значение статистики критерия
    if (n >= 15) and (n <= 100):
        # восстанавливаем структуру DataFrame из csv-файла
        DW_table_df = pd.read_csv(
            filepath_or_buffer='table/Durbin_Watson_test_table.csv',
            sep=';',
            #index_col='n'
            )
                            
        DW_table_df = DW_table_df.rename(columns={'Unnamed: 0': 'n'})
        DW_table_df = DW_table_df.drop([0, 1, 2])
        
        for col in DW_table_df.columns:
            DW_table_df[col] = pd.to_numeric(DW_table_df[col], errors='ignore')
            
        DW_table_df = DW_table_df.set_index('n')

        DW_table_df.columns = pd.MultiIndex.from_product(
            [['p=0.95', 'p=0.975', 'p=0.99'],
            ['m=1', 'm=2', 'm=3', 'm=4', 'm=5'],
            ['dL','dU']])    
        
        # интерполяция табличных значений
        key = [f'p={p_level}', f'm={m}']
        f_lin_L = sci.interpolate.interp1d(DW_table_df.index, DW_table_df[tuple(key + ['dL'])])
        f_lin_U = sci.interpolate.interp1d(DW_table_df.index, DW_table_df[tuple(key + ['dU'])])
        DW_table_L = float(f_lin_L(n))
        DW_table_U = float(f_lin_U(n))
                   
        # проверка гипотезы
        Durbin_Watson_scale = {
            1: DW_table_L,
            2: DW_table_U,
            3: 4 - DW_table_U,
            4: 4 - DW_table_L,
            5: 4}
        
        Durbin_Watson_comparison = {
            1: ['0 ≤ DW_calc < DW_table_L',                   'H1: r > 0'],
            2: ['DW_table_L ≤ DW_calc ≤ DW_table_U',          'uncertainty'],
            3: ['DW_table_U < DW_calc < 4 - DW_table_U',      'H0: r = 0'],
            4: ['4 - DW_table_U ≤ DW_calc ≤ 4 - DW_table_L',  'uncertainty'],
            5: ['4 - DW_table_L < DW_calc ≤ 4',               'H1: r < 0']}
        
        r_scale = list(Durbin_Watson_scale.values())
        for i, elem in enumerate(r_scale):
            if DW_calc <= elem:
                key_scale = list(Durbin_Watson_scale.keys())[i]
                comparison = Durbin_Watson_comparison[key_scale][0]
                conclusion = Durbin_Watson_comparison[key_scale][1]
                break
           
    elif n < 15:        
        comparison = '-'
        conclusion = 'count less than 15'
    else:
        comparison = '-'
        conclusion = 'count more than 100'
    
    
    # формируем результат            
    result = pd.DataFrame({
        'n': (n),
        'm': (m),
        'p_level': (p_level),
        'a_level': (a_level),
        'DW_calc': (DW_calc),
        'ρ': (1 - DW_calc/2),
        'DW_table_L': (DW_table_L if (n >= 15) and (n <= 100) else '-'),
        'DW_table_U': (DW_table_U if (n >= 15) and (n <= 100) else '-'),
        'comparison of calculated and critical values': (comparison),
        'conclusion': (conclusion)
        },
        index=['Durbin-Watson_test'])
    
    
    return result

Протестируем созданную функцию - будем моделировать временные ряды с различными свойствами и выполнять проверку автокорреляции:

y_func = lambda x, b0, b1: b0 + b1*x

N = 30    # число наблюдений
(mu, sigma) = (0, 25)    # параметры моделируемой случайной компоненты (среднее и станд.отклонение)

  1. Смоделируем временной ряд с трендом, без автокорреляции остатков:

# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, 5)    # параметры моделируемого временного ряда
Y_model = np.array(y_func(T_model, b0, b1)) + np.random.normal(mu, sigma, N)

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))
# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, -5)    # параметры моделируемого временного ряда
Y_model = np.array(y_func(T_model, b0, b1)) + np.random.normal(mu, sigma, N)

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))

  1. Смоделируем временной ряд без тренда:

# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, 0)    # параметры моделируемого временного ряда
Y_model = np.array(y_func(T_model, b0, b1)) + np.random.normal(mu, sigma, N)

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))
# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (0, 0)    # параметры моделируемого временного ряда
Y_model = np.array(y_func(T_model, b0, b1)) + np.random.normal(mu, sigma, N)

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))

  1. Смоделируем временной ряд с трендом, с положительной автокорреляцией:

# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, 5)    # параметры моделируемого временного ряда

E = np.array([np.random.normal(mu, sigma, 1)])
r = 0.9
for i in range(1, N):
    elem = r*E[i-1] + np.random.normal(mu, sigma, 1)
    E = np.append(E, elem)
Y_model = np.array(y_func(T_model, b0, b1)) + E

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))

  1. Смоделируем временной ряд с трендом, с отрицательной автокорреляцией:

# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, 5)    # параметры моделируемого временного ряда

E = np.array([np.random.normal(mu, sigma, 1)])
r = -0.9
for i in range(1, N):
    elem = r*E[i-1] + np.random.normal(mu, sigma, 1)
    E = np.append(E, elem)
Y_model = np.array(y_func(T_model, b0, b1)) + E

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))

  1. Смоделируем временной ряд без тренда, с положительной автокорреляцией:

# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, 0)    # параметры моделируемого временного ряда

E = np.array([np.random.normal(mu, sigma, 1)])
r = 0.9
for i in range(1, N):
    elem = r*E[i-1] + np.random.normal(mu, sigma, 1)
    E = np.append(E, elem)
Y_model = np.array(y_func(T_model, b0, b1)) + E

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))
# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (0, 0)    # параметры моделируемого временного ряда

E = np.array([np.random.normal(mu, sigma, 1)])
r = 0.9
for i in range(1, N):
    elem = r*E[i-1] + np.random.normal(mu, sigma, 1)
    E = np.append(E, elem)
Y_model = np.array(y_func(T_model, b0, b1)) + E

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))

  1. Смоделируем временной ряд без тренда, с отрицательной автокорреляцией:

# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (100, 0)    # параметры моделируемого временного ряда

E = np.array([np.random.normal(mu, sigma, 1)])
r = -0.9
for i in range(1, N):
    elem = r*E[i-1] + np.random.normal(mu, sigma, 1)
    E = np.append(E, elem)
Y_model = np.array(y_func(T_model, b0, b1)) + E

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))
# моделирование
T_model = np.linspace(1, N, N)    # независимая переменная - номер наблюдения
(b0, b1) = (0, 0)    # параметры моделируемого временного ряда

E = np.array([np.random.normal(mu, sigma, 1)])
r = -0.9
for i in range(1, N):
    elem = r*E[i-1] + np.random.normal(mu, sigma, 1)
    E = np.append(E, elem)
Y_model = np.array(y_func(T_model, b0, b1)) + E

# визуализация
axes = sns.jointplot(
    x=T_model, y=Y_model,
    kind='reg')
plt.show()

# проверка автокорреляции
display(Durbin_Watson_test(Y_model, m=1, p_level=0.95))

Конечно, данный вычислительный эксперимент не может претендовать на всеобъемлемость, однако определенный любопытный предварительный вывод можно сделать: при наличии любого тренда (даже если этот тренда представляет собой равенство постоянной величине Y=b_0 + \varepsilon) критерий Дарбина-Уотсона выдает нам наличие положительной автокорреляции (даже если в модели автокорреляция не заложена нет или она отрицательная). Такой вывод нужно исследовать более глубоко, но это не входит в цель данного обзора. Специалист должен помнить об особенностях критерия Дарбина-Уотсона.

Теперь мы можем перейти к практическим примерам.

Пример 1: проверка автокорреляции модели временного ряда

Формирование исходных данных

В качестве исходных данных рассмотрим динамику показателей индексов пересчета сметной стоимости проектно-изыскательских работ в РФ. Эти показатели ежеквартально публикует Министерство строительства и ЖКХ РФ, а все проектные и изыскательские организации используют эти показатели при составлении смет на свои работы.

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

Исходные данные содержаться в файле Ежеквартальные индексы ПИР.xlsx, который помещен в папку data.

Прочитаем xlsx-файл:

data_df = pd.read_excel('data/Ежеквартальные индексы ПИР.xlsx', sheet_name='БД')  
#display(data_df)
display(data_df.head(), data_df.tail())
data_df.info()

Не будем подробно останавливаться на содержимом файла и его первичной обработке - это выходит за пределы данного обзора. Специалисты, причастные к сфере строительства и проектирования, поймут, а для остальных специалистов эти цифры можно воспринимать по аналогии с индексами инфляции Росстата и Минэкономразвития.

Прочитаем из этого файла интересующие нас данные - индексы изменения сметной стоимости проектных работ к уровню цен на 01.01.2001 г.:

Ind_design_2001 = np.array(data_df['Ипроект2001'])
print(Ind_design_2001, '\n', type(Ind_design_2001), len(Ind_design_2001))

Сохраним также вспомогательные (технические) переменные, необходимые при анализе временных рядов - дату (Date) и номер наблюдения (T):

# Дата показателя
Date = np.array(data_df['ДАТА'])

# Номер наблюдения
T = np.array(data_df['N'])

Для удобства дальнейшей работы сформируем сформируем отдельный DataFrame:

dataset_df = pd.DataFrame({
    'T': T,
    'Date': Date,
    'Ind_design_2001': Ind_design_2001})
display(dataset_df.head(), dataset_df.tail())

Визуализация

Настройка заголовков:

# Общий заголовок проекта
Task_Project = "Анализ динамики индексов изменения сметной стоимости проектно-изыскательских работ в РФ"
# Заголовок, фиксирующий момент времени
AsOfTheDate = "за 2008-2022 гг."
# Заголовок раздела проекта
Task_Theme = ""
# Общий заголовок проекта для графиков
Title_String = f"{Task_Project}\n{AsOfTheDate}"
# Наименования переменных
Variable_Name_T_month = "Ежемесячные данные"
Variable_Name_Ind_design_2001 = "Индекс изменения сметной стоимости проектных работ к уровню цен на 01.01.2001 г."

# Границы значений переменных (при построении графиков):
(X_min_graph, X_max_graph) = (0.0, max(T))
(Y_min_graph, Y_max_graph) = (2.0, 6.0)
graph_plot_sns_np(
    Date, Ind_design_2001,
    Ymin_in=Y_min_graph, Ymax_in=Y_max_graph,
    color='orange',
    title_figure=Title_String, title_figure_fontsize=12,
    title_axes=Variable_Name_Ind_design_2001, title_axes_fontsize=15,
    x_label=Variable_Name_T_month, label_fontsize=12)

Построение и анализ регрессионной модели

Построим линейную регрессионную модель и проведем ее экспресс-анализ:

model_linear_ols_1 = smf.ols(formula='Ind_design_2001 ~ T', data=dataset_df)
result_linear_ols_1 = model_linear_ols_1.fit()
print(result_linear_ols_1.summary2())

Формализация модели:

# Функция линейной регрессионной модели (SLRM - simple linear regression model)
SLRM_func = lambda x, b0, b1: b0 + b1*x

# параметры модели
b0 = result_linear_ols_1.params['Intercept']
b1 = result_linear_ols_1.params['T']

# уравнение модели
regr_model_linear_ols_1_func = lambda x: SLRM_func(x, b0, b1)

График модели:

R2 = round(result_linear_ols_1.rsquared, DecPlace)
legend_equation = f'линейная регрессия ' + r'$Y$' + f' = {b0:.4f} + {b1:.5f}{chr(183)}' + r'$X$' if b1 > 0 else \
                  f'линейная регрессия ' + r'$Y$' + f' = {b0:.4f} - {abs(b1):.5f}{chr(183)}' + r'$X$'

# Пользовательская функция
graph_regression_plot_sns(
    T, Ind_design_2001,
    regression_model=regr_model_linear_ols_1_func,
    #Xmin=X_min_graph, Xmax=X_max_graph,
    Ymin=Y_min_graph, Ymax=Y_max_graph,
    display_residuals=True,
    title_figure=Variable_Name_Ind_design_2001, title_figure_fontsize=16,
    title_axes = 'Линейная регрессионная модель',
    x_label=Variable_Name_T_month,
    #y_label=Variable_Name_Ind_design_2001,
    label_legend_regr_model = legend_equation + '\n' + r'$R^2$' + f' = {R2}',
    s=60)

Ошибки аппроксимации модели:

(model_error_metrics, result) = regression_error_metrics(model_linear_ols_1, model_name='linear_ols')
display(result)

Проверка нормальности распределения остатков:

res_Y_1 = np.array(result_linear_ols_1.resid)

# Пользовательская функция
graph_hist_boxplot_probplot_sns(
    data=res_Y_1,
    data_min=-0.25, data_max=0.25,
    graph_inclusion='hbp',
    data_label=r'$ΔY = Y - Y_{calc}$',
    #title_figure=Task_Project,
    title_axes='Остатки линейной регрессионной модели', title_axes_fontsize=16)    
norm_distr_check(res_Y_1)

Проверка гетероскедастичности:

Goldfeld_Quandt_test_df = Goldfeld_Quandt_test(result_linear_ols_1, p_level=0.95, model_name='linear_ols')
Breush_Pagan_test_df = Breush_Pagan_test(result_linear_ols_1, p_level=0.95, model_name='linear_ols')
White_test_df = White_test(result_linear_ols_1, p_level=0.95, model_name='linear_ols')

heteroscedasticity_tests_df = pd.concat([Breush_Pagan_test_df, White_test_df, Goldfeld_Quandt_test_df])
display(heteroscedasticity_tests_df)

Проверка автокорреляции:

sms.stattools.durbin_watson(res_Y_1)

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

display(Durbin_Watson_test(res_Y_1, m=1, p_level=0.95))

Выводы по результатам анализа модели:

Итак, мы провели статистический анализ регрессионной модели и установили:

  1. Регрессионная модель хорошо аппроксимирует фактические данные.

  2. Остатки модели имеют нормальное распределение (хотя результаты тестов противоречивы).

  3. Коэффициент детерминации значим; модель объясняет 97% вариации независимой переменной.

  4. Коэффициенты регрессии значимы.

  5. Обнаружена гетероскедастичность.

  6. Тест критерия Дарбина-Уотсона свидетельствует о наличии значимой положительной автокорреляции остатков.

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

  1. На графике модели хорошо заметна точка излома, которая говорит о смене тенденции (существуют специальные статистические тесты для проверки гипотез о смене тенденции, например, тест Чоу, но мы в данном обзоре рассматривать их не будем).

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

  3. Противоречивость тестов проверки нормальности распределения остатков.

  4. Наличие гетероскедастичности.

  5. Наличие автокорреляции. Явление автокорреляции может возникать в случае смены тенденции [5, с.118].

Тот факт, что распределение остатков признается нормальным по результатам таких тестов как Шапиро-Уилка, Эппса-Палли, Андерсона-Дарлинга может иметь разные причины, например, мы можем иметь дело со смесью двух распределений. Этот вопрос требует отдельного тщательного исследования.

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

graph_regression_pair_predict_plot_sns(
    model_fit=result_linear_ols_1,
    regression_model_in=regr_model_linear_ols_1_func,
    Xmin=X_min_graph, Xmax=X_max_graph+12, Nx=25,
    Ymin_graph=2.0, Ymax_graph=Y_max_graph,
    title_figure=Variable_Name_Ind_design_2001, title_figure_fontsize=16,
    title_axes='Линейная регрессионная модель', title_axes_fontsize=14,
    #x_label=Variable_Name_X,
    #y_label=Variable_Name_Y,
    label_legend_regr_model=f'линейная регрессия Y = {b0:.3f} + {b1:.4f}*X',
    s=50,
    result_output=False)

Нет, такой прогноз нам не нужен.

Пример 2: проверка автокорреляция регрессионной модели

Формирование исходных данных

Рассмотрим пример множественной линейной регрессионной модели, приведенный в источнике [6, с.192].

В качестве исходных данных рассматриваются ряд макроэкономических показателей США за 1960-1985 гг. (в сопоставимых ценах 1982 г., млрд.долл):

  • DPI - годовой совокупный располагаемый личный доход;

  • CONS - годовые совокупные потребительские расходы;

  • ASSETS - финансовые активы населения на начало календарного года.

Предполагается, что между переменной CONS и регрессорами DPI, ASSETS имеется линейная регрессионная связь.

Исходные данные содержаться в файле Macroeconomic_indicators_USA_1960_1985.csv, который помещен в папку data.

Прочитаем csv-файл:

data_df = pd.read_csv(filepath_or_buffer='data/Macroeconomic_indicators_USA_1960_1985.csv', sep=';')
display(data_df)
#display(data_df.head(), data_df.tail())
data_df.info()

Визуализация

fig, axes = plt.subplots(figsize=(297/INCH, 210/INCH))

title_figure = 'Анализ макроэкономических показателей США за 1960-1985 гг.'
fig.suptitle(title_figure, fontsize = 18)

sns.lineplot(
    x = data_df['YEAR'], y = data_df['DPI'],
    linewidth=3,
    legend=True,
    label='DPI',
    ax=axes)

sns.lineplot(
    x = data_df['YEAR'], y = data_df['CONS'],
    linewidth=3,
    legend=True,
    label='CONS',
    ax=axes)

sns.lineplot(
    x = data_df['YEAR'], y = data_df['ASSETS'],
    linewidth=3,
    legend=True,
    label='ASSETS',
    ax=axes)

axes.set_xlabel('Year')
axes.set_ylabel('US$ billion')

plt.show()

Построение и анализ регрессионной модели

Построим линейную регрессионную модель и проведем ее экспресс-анализ:

y = data_df['CONS']
X = data_df[['DPI', 'ASSETS']]
X = sm.add_constant(X)
model_linear_ols_2 = sm.OLS(y, X)
result_linear_ols_2 = model_linear_ols_2.fit()
print(result_linear_ols_2.summary2())

График модели:

fig = plt.figure(figsize=(297/INCH, 420/INCH/1.5))
ax1 = plt.subplot(2,1,1)
ax2 = plt.subplot(2,1,2)

title_figure = 'Анализ макроэкономических показателей США за 1960-1985 гг.'
fig.suptitle(title_figure, fontsize = 18)

fig = sm.graphics.plot_fit(
    result_linear_ols_2, 'DPI',
    vlines=True,    # это параметр отвечает за отображение доверительных интервалов для Y
    ax=ax1)
ax1.set_ylabel('CONS (US$ billion)', fontsize = 12)
ax1.set_xlabel('DPI (US$ billion)', fontsize = 12)
ax1.set_title('Fitted values vs. DPI', fontsize = 15)

fig = sm.graphics.plot_fit(
    result_linear_ols_2, 'ASSETS',
    vlines=True,    # это параметр отвечает за отображение доверительных интервалов для Y
    ax=ax2)
ax2.set_ylabel('CONS (US$ billion)', fontsize = 12)
ax2.set_xlabel('ASSETS (US$ billion)', fontsize = 12)
ax2.set_title('Fitted values vs. ASSETS', fontsize = 15)

plt.show()

Ошибки аппроксимации модели:

(model_error_metrics, result) = regrpy
ession_error_metrics(model_linear_ols_2, model_name='linear_ols')
display(result)

Проверка нормальности распределения остатков:

res_Y_2 = np.array(result_linear_ols_2.resid)

# Пользовательская функция
graph_hist_boxplot_probplot_sns(
    data=res_Y_2,
    data_min=-60, data_max=60,
    graph_inclusion='hbp',
    data_label=r'$ΔY = Y - Y_{calc}$',
    #title_figure=Task_Project,
    title_axes='Остатки линейной регрессионной модели', title_axes_fontsize=16)    

norm_distr_check(res_Y_2)

Прроверка гетероскедастичности:

Goldfeld_Quandt_test_df = Goldfeld_Quandt_test(result_linear_ols_2, p_level=0.95, model_name='linear_ols')
Breush_Pagan_test_df = Breush_Pagan_test(result_linear_ols_2, p_level=0.95, model_name='linear_ols')
White_test_df = White_test(result_linear_ols_2, p_level=0.95, model_name='linear_ols')

heteroscedasticity_tests_df = pd.concat([Breush_Pagan_test_df, White_test_df, Goldfeld_Quandt_test_df])
display(heteroscedasticity_tests_df)

Проверка автокорреляции:

display(Durbin_Watson_test(res_Y_2, m=1, p_level=0.95))

Выводы по результатам анализа модели:

Как видим, в целом результаты расчетов совпадают с результатами из первоисточника [6], в части выявления автокорреляции аналогично.

Информация к размышлению.

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

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

В дальнейшем автор при анализе остатков модели [6, с.198] выявляет структурный сдвиг (обусловленный мировым топливно-энергетическим кризисом в 1973 г.) и вводит в модель фиктивные переменные, учитывающие этот структурный сдвиг

Итоги

Итак, подведем итоги:

  • мы рассмотрели способы реализации в полной мере критерия Дарбина-Уотсона средствами python, создали пользовательскую функцию, уменьшающую размер кода;

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

Исходный код находится в моем репозитории на GitHub.

Надеюсь, данный обзор поможет специалистам DataScience в работе.

Литература

  1. Кобзарь А.И. Прикладная математическая статистика. Для инженеров и научных работников. - М.: ФИЗМАТЛИТ, 2006. - 816 с.

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

  3. Фёрстер Э., Рёнц Б. Методы корреляционного и регрессионного анализа / пер с нем. - М.: Финансы и статистика, 1983. - 302 с.

  4. Магнус Я.Р. и др. Эконометрика. Начальный курс - М.: Дело, 2004. - 576 с.

  5. Тихомиров Н.П., дорохина Е.Ю. Эконометрика. - М.: Экзамен, 2003. - 512 с.

  6. Носко В.П. Эконометрика. Кн.1. Ч.1, 2. - М.: Издательский дом "Дело" РАНХиГС, 2011. - 672 с.

  7. Вандер Плас Дж. Python для сложных задач: наука о данных и машинное обучение. - СПб: Питер, 2018. - 576 с.

Комментарии (2)


  1. eltardowut
    14.10.2022 16:03
    #24820674
    +2

    Отличная статья.


  1. slava_k
    15.10.2022 00:48
    #24823164

    Спасибо за статью!