Введение

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

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

Подгружаем ggplot2

Если у вас еще не стоит этот замечательный пакет, вы должны установить его с помощью

install.packages("ggplot2")

Когда он уже стоит, достаточно его подгрузить

library(ggplot2)

Строим гомоскедастические данные

(то есть данные без гетероскедастичности)

Сначала задаем параметры наших будущих данных.

N <- 150     # число наблюдений(точек)
x.max <- 100 # максимальное значение X
noise.level <- 20  # насколько сильно отклонение от линейности
a <- 2  
b <- 0.0  # коэффициенты уравнения прямой y = ax + b

Затем строим точки вдоль идеальной прямой

perfect.x <- runif(N, min=0, max=x.max)
perfect.y <- a*perfect.x + b

Добавляем немного шума (распределение шума будет нормальным):

noisy.y   <- perfect.y + rnorm(N, mean=0, sd=noise.level)

Создаем data.frame со всеми нашими данными и выводим его заголовок, чтобы убедиться, что получили таблицу такой формы, на которую рассчитывали:

noisy.data <- data.frame(x=perfect.x, y=noisy.y)

head(noisy.data) # выводим первые несколько строк
##           x        y
## 1 50.056318 49.28416
## 2 25.473245 55.87474
## 3 38.426436 96.33662
## 4  4.592379 22.77633
## 5 48.080299 78.40616
## 6 19.627918 39.17321

Строим график со всеми данными, а также линейным приближением:

ggplot (noisy.data, aes(x=x, y=y)) +
  geom_point(shape=1) +    
  geom_smooth(method=lm, se=FALSE)

Такие данные и называются гомоскедастическими - в них уровень отклонения реальных значений от идеальной зависимости не зависит от X или Y.

Строим гетероскедастические данные

Воспользуемся уже заготовленными perfect.x и perfect.y, но в этот раз добавим шум, пропорциональный X:

noise.factor <- 0.4
heterosc.y <- perfect.y + 
  rnorm(N, mean=0, sd=noise.factor) * perfect.x 
  # строчкой выше мы поэлементно умножаем случайный шум каждого i-го 
  # наблюдения на X[i]
heterosc.data <- data.frame(x=perfect.x, y=heterosc.y)

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

ggplot (heterosc.data, aes(x=x, y=y)) +
  geom_point(shape=1) +    
  geom_smooth(method=lm, se=FALSE)

Что только что произошло

Я думаю, все видят опреденную разницу в данных на двух графиках. Действительно, линия, которой эти данные фитятся в ходе линейной регрессией, в обоих случаях выглядит почти одинаково. Но во втором случае уровень отклонения Y от прямой зависит от величины X или соответственно величины Y (так как X и Y сами очень зависимые друг от друга величины).

Есть ли в этом физический смысл? Да, определенно. Представим, что X - это квадрат роста человека, а Y - его вес. На самом деле все будет не идеально линейно. Но во-первых, предположим, что линейно. А во-вторых, гетероскедастичность существует не только при использовании линейной модели, но и при любой(?) другой, например, полиномиальной.

Итак, положим, что у нас есть два отклонения. Один человек отклоняется от общей формулы на 50 кг и при росте в 2 м весит не 90, как должен, а 140. Сильное ли это отклонение? Ну допустим да. Другой человек отклоняется от общей формулы на те же 50кг и при росте 0.8 м весит не 12 кг, как должен, а 62 кг. Или -38. Согласитесь же, что это отклонение будет НАМНОГО сильнее.

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

И в чем проблема?

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

Что делать?

Для начала убедиться, что гетероскедастичность есть, что это не наши буйные фантазии. Один из вариантов - тест Бройша-Пагана (Breusch-Pagan test).

Подгружаем пакет lmtest(его тоже придется сначала поставить):

library(lmtest)

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

my.model <- lm(perfect.x ~ heterosc.y)
bp <- bptest(my.model)
bp
## 
##  studentized Breusch-Pagan test
## 
## data:  my.model
## BP = 5.0567, df = 1, p-value = 0.02453

p-value очень маленькое, поэтому гипотеза выполняется:)

Теперь поробуем то же на наших первых данных с равномерным шумом:

my.model <- lm(perfect.x ~ noisy.y)
bp <- bptest(my.model)
bp
## 
##  studentized Breusch-Pagan test
## 
## data:  my.model
## BP = 0.43404, df = 1, p-value = 0.51

p-value велико, и действительно, гипотеза о гетероскедастичности не подтвердилась.

И что делать, если подтвердилась?

В лекции говорилось, что в таких случаях мы используем взвешенный МНК, то есть чем больше X, тем меньше мы учитываем отклонение соответствующего Y. Функция lm имеет дополнительный аргумент weight для таких случаев. Можно задать его совсем агрессивно, на уровне

new.model <- lm(perfect.x ~ heterosc.y, weights=1/perfect.x)

Но может быть, это не лучшее решение.

(И может быть, продолжение следует)

Вопросы и замечания

Я не математик, не экономист и не R-программист, поэтому я мог понять что-то не совсем правильно. С удовольствием приму ваши замечания и уточнения.

Спасибо!