Econometrics in R - cvičenie 9

1. Úvod a údaje

V tomto cvičení pracujeme s databázou filmov data_ekonometria.csv, ktorá obsahuje 250 najlepšie hodnotených filmov podľa IMDb (umiestnenie, názov, rok uvedenia, žáner, IMDb hodnotenie, dĺžku, réžiu, filmové štúdio a rozpočet).

Údaje o filmoch sú uložené v súbore csv, stĺpce sú oddelené bodkočiarkou ";" a používajú desatinnú čiarku. Súbor data_ekonometria.csv je uložený v tom istom priečinku ako tento RMarkdown súbor.

Databáza obsahuje najlepšie hodnotené filmy podľa IMDb s informáciou o roku uvedenia, žánri, dĺžke trvania, režisérovi, filmovom štúdiu a rozpočte. Nie všetky údaje budú použité, preto sa zameriame len na filmy jednej filmovej spoločnosti, aby sme mali časové usporiadanie – konkrétne na filmy štúdia Warner Bros. Pictures.

1.1 Úvod do problému, stanovenie hypotéz

Rozhodli sme sa modelovať IMDb hodnotenie filmu (premenná IMDb) v závislosti od troch vysvetľujúcich premenných:

dĺžka filmu v minútach (Dlzka_min),
rozpočet v dolároch (Rozpocet),
rok uvedenia filmu (Rok), ktorý zachytáva aj možné časové trendy v hodnoteniach.

Pracovná hypotéza je, že:

vyšší rozpočet by mohol viesť k lepšiemu filmu a vyššiemu IMDb hodnoteniu (očakávame kladné znamienko pri koeficiente rozpočtu),
primerane dlhší film (väčší priestor na rozvinutie deja) môže byť vnímaný lepšie, takže aj pri dĺžke očakávame skôr kladný vplyv,
pri roku uvedenia nepredpokladáme jednoznačný smer – ide skôr o kontrolnú premennú, ktorá zachytáva možné zmeny preferencií divákov v čase.

1.2 Príprava databázy, úprava údajov

Najprv načítame celý súbor, upravíme formát čísel a následne vyberieme len filmy štúdia Warner Bros. Pictures:

# kontrola, či súbor existuje v aktuálnom pracovnom priečinku
if (!file.exists("data_ekonometria.csv")) {
  stop("Súbor 'data_ekonometria.csv' sa nenašiel v aktuálnom priečinku. 
Skontroluj, prosím, že Rmd a CSV sú v tom istom priečinku.")
}

# načítanie údajov – bodkočiarka, desatinná čiarka
udaje <- read.csv2("data_ekonometria.csv",
                   header = TRUE,
                   stringsAsFactors = FALSE)

# premenovanie stĺpcov do "R-friendly" tvaru
colnames(udaje) <- c("Umiestnenie", "Nazov", "Rok", "Zaner",
                     "IMDb_raw", "Dlzka_min", "Rezia",
                     "Spolocnost", "Rozpocet_raw")

# úprava IMDb hodnotenia (z formátu "9,3" na numerické 9.3)
udaje$IMDb <- as.numeric(gsub(",", ".", udaje$IMDb_raw))

# úprava rozpočtu (odstránenie medzier a prevod na číslo)
udaje$Rozpocet <- as.numeric(gsub(" ", "", udaje$Rozpocet_raw))

# výber filmov Warner Bros. Pictures a relevantných premenných
udajeWB <- udaje[udaje$Spolocnost == "Warner Bros. Pictures",
                 c("Rok", "IMDb", "Dlzka_min", "Rozpocet")]

# zoradenie podľa roku
udajeWB <- udajeWB[order(udajeWB$Rok), ]

udajeWB

##      Rok IMDb Dlzka_min  Rozpocet
## 43  1942  8.5       162   1000000
## 146 1948  8.2       126   2500000
## 162 1954  8.2       105   1400000
## 238 1967  8.1       127   3200000
## 103 1971  8.3       136   1300000
## 226 1973  8.1       122  12000000
## 191 1975  8.1       185  12000000
## 61  1980  8.4       146  19000000
## 176 1982  8.1       117  30000000
## 80  1984  8.4       229  30000000
## 105 1987  8.3       116  30000000
## 17  1990  8.7       145  25000000
## 140 1992  8.2       130  14400000
## 111 1995  8.3       170  60000000
## 117 1997  8.2       138  35000000
## 16  1999  8.7       136  63000000
## 245 1999  8.0        86  50000000
## 127 2005  8.2       140 150000000
## 158 2005  8.2       132  53000000
## 39  2006  8.5       151  90000000
## 41  2006  8.5       130  40000000
## 137 2006  8.2       118  19000000
## 3   2008  9.0       152 185000000
## 171 2008  8.1       116  30000000
## 14  2010  8.8       148 160000000
## 181 2011  8.1       130 250000000
## 69  2012  8.4       164 230000000
## 25  2014  8.6       169 165000000
## 196 2014  8.1       122   4000000
## 198 2015  8.1       120 185100000
## 75  2019  8.4       122  69000000

2. Lineárna regresia v základnom tvare

Budeme odhadovať rovnicu:

\[ \text{IMDb}_t = \beta_0 + \beta_1 \text{Dlzka}_t + \beta_2 \text{Rozpocet}_t + \beta_3 \text{Rok}_t + e_t, \]

kde:

\(\text{IMDb}_t\) – IMDb hodnotenie vybraného filmu Warner Bros. v roku \(t\),
\(\text{Dlzka}_t\) – dĺžka filmu v minútach,
\(\text{Rozpocet}_t\) – rozpočet filmu v dolároch,
\(\text{Rok}_t\) – rok uvedenia filmu,
\(e_t\) – náhodná zložka.

# základný regresný model
model <- lm(IMDb ~ Dlzka_min + Rozpocet + Rok, data = udajeWB)
summary(model)

## 
## Call:
## lm(formula = IMDb ~ Dlzka_min + Rozpocet + Rok, data = udajeWB)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33337 -0.13538 -0.05900  0.09991  0.53793 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.269e+00  5.168e+00   1.213   0.2356  
## Dlzka_min   3.413e-03  1.583e-03   2.156   0.0402 *
## Rozpocet    6.779e-10  7.355e-10   0.922   0.3648  
## Rok         7.713e-04  2.586e-03   0.298   0.7678  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2257 on 27 degrees of freedom
## Multiple R-squared:  0.2249, Adjusted R-squared:  0.1388 
## F-statistic: 2.612 on 3 and 27 DF,  p-value: 0.07183

3. Autokorelácia rezíduí

V tejto časti sa pozrieme na ďalší dôležitý predpoklad klasického lineárneho regresného modelu – nezávislosť rezíduí. Pri údajoch usporiadaných v čase (napr. filmy jedného štúdia podľa roku uvedenia) sa môže stať, že chyba v čase \(t\) je systematicky spätá s chybou v čase \(t-1\). Vtedy hovoríme o autokorelácii rezíduí.

3.1 Čo je autokorelácia rezíduí?

Autokorelácia rezíduí je situácia, keď platí:

\[ \operatorname{Corr}(e_t, e_{t-k}) \neq 0 \quad \text{pre niektoré } k \neq 0. \]

Najčastejšie sa skúma autokorelácia prvého rádu:

\[ e_t = \rho e_{t-1} + \nu_t,\quad |\rho| < 1. \]

3.2 Dôsledky autokorelácie rezíduí

odhady koeficientov \(\hat{\beta}\) sú síce nestranné,
ale neefektívne,
štandardné chyby môžu byť podhodnotené → p-hodnoty sa javia menšie, než by mali byť,
môže dochádzať k falošnému zdaniu štatistickej významnosti,
teda t-testy a F-testy sú skreslené.

3.3 Detekcia autokorelácie

Grafická informácia

Najprv sa pozrieme na to, ako model vysvetlí vývoj IMDb hodnotenia filmov v čase. Pridáme k údajom vyrovnané (fitted) hodnoty z regresie a zobrazíme ich spolu s empirickými bodmi.

# pridáme vyrovnané hodnoty do dát
udajeWB$fitted <- fitted(model)

# scatterplot + vyrovnaná línia
ggplot(udajeWB, aes(x = Rok, y = IMDb)) +
  geom_point(size = 2) +
  geom_line(aes(y = fitted), size = 1) +
  labs(
    title = "IMDb hodnotenie filmov Warner Bros. Pictures: empirické vs. vyrovnané hodnoty",
    x = "Rok",
    y = "IMDb hodnotenie"
  ) +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Z grafu vidíme, že hodnotenia filmov Warner Bros. sú rozmiestnené okolo vyrovnanej línie modelu relatívne pravidelne a nevidíme dlhé súvislé úseky, kde by boli rezíduá iba kladne alebo iba záporne. Už vizuálne teda nemáme silné podozrenie na výraznú autokoreláciu.

# uložíme si reziduá z pôvodného modelu
res <- residuals(model)

ACF graf (Autocorrelation Function)

Táto funkcia priraďuje odhad korelácie medzi rezíduami v aktuálnom období a rezíduami posunutými o \(k\) období späť.

acf(res, lag.max = 4, main = "Autokorelačná funkcia rezíduí")

Na grafe je prerušovanou čiarou vyjadrený aj 95 % interval spoľahlivosti pre hodnotu autokorelačného koeficientu pri danom posune. V našich údajoch žiadny z odhadnutých autokorelačných koeficientov výrazne nevybočuje z tohto intervalu, takže ACF graf nenaznačuje výraznú sériovú koreláciu rezíduí.

Durbin–Watsonov test

Durbin–Watsonov koeficient (DW) je vypočítaný z rezíduí podľa vzorca

\[ DW = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^{2}}{\sum_{t=1}^{n} e_t^{2}}, \]

kde \(n\) je počet pozorovaní. Medzi koeficientom autokorelácie dvoch susedných rezíduí a DW platí približný vzťah

\[ \hat{\rho} \approx 1 - \frac{DW}{2}. \]

Hodnoty DW sú v intervale od 0 po 4.

hodnoty blízke 0 → silná pozitívna autokorelácia,
hodnoty blízke 4 → silná negatívna autokorelácia,
hodnoty okolo 2 → bez problémov s autokoreláciou prvého rádu.

V praxi sa často používa intuitívne pravidlo: ak sa DW nachádza v intervale približne 1,8 až 2,2, problém autokorelácie prvého rádu zväčša nemusíme riešiť.

dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.3437, p-value = 0.7856
## alternative hypothesis: true autocorrelation is greater than 0

V našom prípade sa Durbin–Watsonov koeficient pohybuje blízko hodnoty 2 a p-hodnota testu nenaznačuje štatisticky významnú autokoreláciu prvého rádu rezíduí.

DW test má svoje obmedzenia – regresory nesmú byť časovo posunuté a model nesmie obsahovať oneskorené hodnoty vysvetľovanej premennej ako regresory. Tieto obmedzenia nemá Breusch–Godfreyov test.

Breusch–Godfreyov test (BG test)

BG test je formálnym testom autokorelácie (sériovej korelácie) rezíduí regresného modelu až do zvoleného rádu \(p\).

Pre test autokorelácie až do rádu \(p\) platí:

\(H_0\): \(\rho_1 = \rho_2 = \cdots = \rho_p = 0\) – žiadna sériová korelácia,
\(H_1\): aspoň jedna z \(\rho_j \neq 0\) – sériová korelácia prítomná.

bgtest(model, order = 1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 1.381, df = 1, p-value = 0.2399

V našom prípade BG test rovnako ako Durbin–Watsonov test nepotvrdzuje prítomnosť autokorelácie prvého rádu – p-hodnota je relatívne vysoká, takže nulovú hypotézu o „žiadnej autokorelácii“ nezamietame. Napriek tomu si demonštratívne ukážeme, ako by bolo možné model dynamizovať pomocou Koyckovho prístupu.

4. Dynamizácia modelu – Koyckov model

Teoreticky vychádzame z modelu s distribuovaným oneskorením, kde efekty vysvetľujúcich premenných v čase klesajú geometricky. Koyckova transformácia vedie k modelu:

\[ \text{IMDb}_t = \alpha(1-\lambda) + \beta_0 \text{Rozpocet}_t + \gamma_0 \text{Dlzka}_t + \lambda \text{IMDb}_{t-1} + v_t, \]

kde \(\lambda\) reprezentuje rýchlosť prispôsobenia a \(\text{IMDb}_{t-1}\) je oneskorená hodnota vysvetľovanej premennej.

V našom prípade:

\(y_t\) = IMDb hodnotenie filmu Warner Bros. v roku \(t\),
vysvetľujúce premenné: rozpočet a dĺžka filmu.

4.1 Odhad Koyckovho modelu v R

udajeWB <- udajeWB %>%
  arrange(Rok) %>%
  mutate(
    IMDb_lag1 = dplyr::lag(IMDb)
  )

model_koyck <- lm(IMDb ~ Rozpocet + Dlzka_min + IMDb_lag1,
                  data = udajeWB)

summary(model_koyck)

## 
## Call:
## lm(formula = IMDb ~ Rozpocet + Dlzka_min + IMDb_lag1, data = udajeWB)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32467 -0.13251 -0.00605  0.07865  0.52386 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.671e+00  1.693e+00   5.713 5.18e-06 ***
## Rozpocet     8.675e-10  5.726e-10   1.515    0.142    
## Dlzka_min    2.258e-03  1.751e-03   1.289    0.209    
## IMDb_lag1   -2.071e-01  1.887e-01  -1.097    0.283    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2231 on 26 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2568, Adjusted R-squared:  0.171 
## F-statistic: 2.994 on 3 and 26 DF,  p-value: 0.04905

dwtest(model_koyck)

## 
##  Durbin-Watson test
## 
## data:  model_koyck
## DW = 1.986, p-value = 0.4658
## alternative hypothesis: true autocorrelation is greater than 0

V dynamizovanom modeli je koeficient pri oneskorenej IMDb (\(\text{IMDb}_{t-1}\)) v našich údajoch relatívne malý a štatisticky nevýznamný, čo je v súlade s predchádzajúcimi testami – výrazný zotrvačný efekt v hodnoteniach filmov sa tu nepotvrdil. Porovnanie ukazovateľa Adjusted R-squared medzi pôvodným a dynamickým modelom tiež ukazuje, že pridaním oneskorenej hodnoty IMDb sa kvalita modelu výrazne nezlepšila.

4.2 Newey–West robustné (alebo aspoň robustné) štandardné chyby

Aj keď sme v našich údajoch nenašli silný dôkaz o autokorelácii, v praxi je často vhodné pracovať s robustnými štandardnými chybami, ktoré sú robustné voči heteroskedasticite a určitej forme autokorelácie (Newey–West). Pri krátkych časových radoch však môže výpočet Newey–Westovho odhadu zlyhať. Preto použijeme bezpečnú verziu:

pokúsime sa vypočítať Newey–Westov odhad,
ak to zlyhá, použijeme aspoň heteroskedasticky robustný odhad typu HC0.

vcov_nw <- tryCatch(
  {
    # konzervatívny malý lag, aby to nepadalo na krátkych radoch
    NeweyWest(model, prewhite = FALSE, adjust = TRUE, lag = 1)
  },
  error = function(e) {
    message("Newey-West zlyhal, používam vcovHC (HC0). Chyba bola: ", e$message)
    sandwich::vcovHC(model, type = "HC0")
  }
)

coeftest(model, vcov = vcov_nw)

## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.2692e+00 3.4929e+00  1.7948  0.08388 .
## Dlzka_min   3.4126e-03 1.4549e-03  2.3455  0.02660 *
## Rozpocet    6.7791e-10 9.1655e-10  0.7396  0.46591  
## Rok         7.7129e-04 1.7633e-03  0.4374  0.66529  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5. Záver

Na dátach filmov štúdia Warner Bros. Pictures sme odhadli lineárny regresný model, v ktorom IMDb hodnotenie závisí od dĺžky filmu, rozpočtu a roku uvedenia. Následne sme sa zamerali na diagnostiku autokorelácie rezíduí pomocou grafického posúdenia, ACF grafu, Durbin–Watsonovho testu a Breusch–Godfreyho testu.

V našej databáze sa výrazná autokorelácia rezíduí nepotvrdila, čo je prirodzené aj vzhľadom na to, že ide o filmy zoradené v čase, nie o klasický ekonomický časový rad s ročnou frekvenciou makroekonomických ukazovateľov. Napriek tomu sme demonštratívne ukázali, ako možno model dynamizovať pomocou Koyckovho modelu a ako použiť Newey–Westove (alebo aspoň heteroskedasticky robustné) štandardné chyby.

Metódy odstraňovania dôsledkov autokorelácie (Koyckova transformácia, Almonove polynómy, Cochran–Orcutt a najmä ARIMA modely) majú v ekonometrii veľký význam, najmä pri „skutočných“ časových radoch (napr. HDP, inflácia, nezamestnanosť). Tu sme ich ilustrovali na jednoduchej filmovej databáze, ktorá je prispôsobená požiadavkám zadania.