Kvalita regresného modelu ceny ropy Brent
V zadaní 4 sme na údajoch z databázy Global Oil Prices & Geopolitical Events (2010–2026) (Kaggle, autor Kavya Dhyani, licencia CC BY 4.0) odhadli dva modely:
model1: cena Brent vysvetlená iba indexom geopolitického rizika GPR,model: cena Brent vysvetlená indexmi GPR, VIX a DXY.V zadaní 4 sme koeficienty iba interpretovali a zámerne sme sa nevenovali štatistickému testovaniu. V tomto zadaní na tento model nadviažeme a posúdime jeho kvalitu ako celku:
R.squared, Adj.R.squared.Sigma, F Stat., p.value, df, df.residual.library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(dplyr)
# Pomocná funkcia zo zadania 4: číslo so slovenskou desatinnou čiarkou.
sk_cislo <- function(x, digits = 3) {
formatC(x, format = "f", digits = digits, decimal.mark = ",", big.mark = " ")
}Údaje pripravíme rovnakým postupom ako v zadaní 4, aby sme pracovali s tou istou vzorkou a s tým istým modelom. Súbor oil_geopolitics_dataset_2010_2026.csv môže byť v podpriečinku udaje alebo priamo vedľa dokumentu.
cesta <- "udaje/oil_geopolitics_dataset_2010_2026.csv"
if (!file.exists(cesta)) cesta <- "oil_geopolitics_dataset_2010_2026.csv"
if (!file.exists(cesta)) {
stop("Vložte CSV do podpriečinka udaje alebo vedľa QMD dokumentu.")
}
ropa <- read.csv(cesta, sep = ",", dec = ".",
na.strings = c("", "NA"),
check.names = FALSE, stringsAsFactors = FALSE)
ropa$date <- as.Date(ropa$date)
# Rovnaké premenné ako v zadaní 4.
udaje <- ropa[, c("date", "brent_price", "gpr_index", "vix", "dxy_index")]
# Rovnaké ošetrenie chýbajúcich hodnôt ako v zadaní 4:
# riadky bez ceny Brent vynecháme, vysvetľujúce premenné doplníme mediánom.
udaje <- udaje[!is.na(udaje$brent_price), ]
for (premenna in c("gpr_index", "vix", "dxy_index")) {
chyba <- is.na(udaje[[premenna]])
udaje[[premenna]][chyba] <- median(udaje[[premenna]], na.rm = TRUE)
}
cat("Počet dní vo vzorke:", nrow(udaje), "\n")Počet dní vo vzorke: 4047 cat("Obdobie:", format(min(udaje$date)), "až", format(max(udaje$date)), "\n")Obdobie: 2010-02-17 až 2026-03-12 Odhadujeme ten istý viacnásobný model ako v zadaní 4:
\[ BrentPrice_t = \beta_0 + \beta_1 GPR_t + \beta_2 VIX_t + \beta_3 DXY_t + u_t \]
kde \(t\) označuje obchodný deň. Pre porovnanie odhadneme znova aj jednoduchý model model1 so samotným indexom GPR.
model1 <- lm(brent_price ~ gpr_index, data = udaje)
model <- lm(brent_price ~ gpr_index + vix + dxy_index, data = udaje)V zadaní 4 sme zobrazili iba odhady koeficientov. Teraz pridáme aj štandardné chyby, t-štatistiky, p-hodnoty a hviezdičky významnosti:
tidy(model) |>
mutate(
Significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.10 ~ ".",
TRUE ~ ""
)
) |>
kable(
digits = 4,
col.names = c("Premenná", "Odhad", "Std. Error", "t value",
"p-value", "Významnosť"),
caption = "Regresné koeficienty modelu ceny Brent"
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Odhad | Std. Error | t value | p-value | Významnosť |
|---|---|---|---|---|---|
| (Intercept) | 216.8578 | 3.0072 | 72.1120 | 0e+00 | *** |
| gpr_index | 0.2705 | 0.0102 | 26.6284 | 0e+00 | *** |
| vix | -0.1698 | 0.0420 | -4.0470 | 1e-04 | *** |
| dxy_index | -1.7637 | 0.0345 | -51.0791 | 0e+00 | *** |
Všetky tri vysvetľujúce premenné sú pri tejto veľkej vzorke štatisticky významné. Koeficient pri GPR je kladný: o jeden bod vyšší index geopolitického rizika sa pri nezmenených VIX a DXY spája s cenou Brent vyššou približne o 0,271 USD/barel. Koeficient pri DXY je záporný: silnejší dolár sa spája s nižšou cenou ropy vyjadrenou v USD, čo zodpovedá ekonomickej intuícii.
Súhrnné charakteristiky modelu ako celku získame funkciou glance(). Tabuľku zobrazíme celú (bez deviance), ale v ďalších častiach sa sústredíme iba na prvé dve skupiny charakteristík.
glance(model) |>
select(-deviance) |>
kable(
digits = 4,
col.names = c("R.squared", "Adj.R.squared", "Sigma", "F Stat.",
"p.value", "df", "logLik", "AIC", "BIC",
"df.residual", "N"),
caption = "Súhrnné charakteristiky modelu (glance)"
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| R.squared | Adj.R.squared | Sigma | F Stat. | p.value | df | logLik | AIC | BIC | df.residual | N |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.3962 | 0.3958 | 18.2894 | 884.4787 | 0 | 3 | -17502.33 | 35014.65 | 35046.18 | 4043 | 4047 |
# Hodnoty si uložíme pre ďalšie výpočty a text.
g <- glance(model)
g1 <- glance(model1)
n <- nobs(model)
k <- length(coef(model)) # počet koeficientov vrátane konštantyKoeficient determinácie vyjadruje, akú časť variability ceny ropy Brent vysvetľuje regresný model:
\[ R^2 = 1 - \frac{RSS}{TSS}, \qquad RSS = \sum_{t=1}^{n} e_t^2, \qquad TSS = \sum_{t=1}^{n} (y_t - \bar y)^2 . \]
Hodnotu si overíme aj ručným výpočtom:
y <- udaje$brent_price
RSS <- sum(residuals(model)^2)
TSS <- sum((y - mean(y))^2)
R2_rucne <- 1 - RSS / TSS
R2_adj_rucne <- 1 - (1 - R2_rucne) * (n - 1) / (n - k)
data.frame(
Ukazovatel = c("R²", "Adjusted R²"),
Rucny_vypocet = c(R2_rucne, R2_adj_rucne),
glance = c(g$r.squared, g$adj.r.squared)
) |>
kable(digits = 6, col.names = c("Ukazovateľ", "Ručný výpočet", "glance()"),
caption = "Overenie R² a upraveného R²") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Ukazovateľ | Ručný výpočet | glance() |
|---|---|---|
| R² | 0.396246 | 0.396246 |
| Adjusted R² | 0.395798 | 0.395798 |
Interpretácia: Model vysvetľuje približne 39,6 % variability dennej ceny ropy Brent v rokoch 2010–2026. Zvyšných približne 60,4 % variability zostáva nevysvetlených – ide o vplyv faktorov, ktoré model neobsahuje (napr. ponuka OPEC+, dopyt Číny, zásoby ropy, pandémia COVID-19).
\[ \bar R^2 = 1 - (1 - R^2)\,\frac{n-1}{n-k} \]
kde \(k\) je počet odhadovaných koeficientov vrátane konštanty (u nás \(k = 4\)).
Interpretácia: Upravený koeficient determinácie má hodnotu 0,3958 a je takmer zhodný s \(R^2\) = 0,3962. Je to preto, že máme veľmi veľa pozorovaní (\(n = 4047\)) a iba tri vysvetľujúce premenné, takže penalizácia za počet premenných je zanedbateľná.
Upravený koeficient determinácie je vhodný práve na porovnanie modelov s rôznym počtom premenných, ktoré majú rovnakú vysvetľovanú premennú a sú odhadnuté na rovnakých pozorovaniach – presne ako naše dva modely zo zadania 4.
porovnanie_modelov <- data.frame(
Model = c("model1: Brent ~ GPR", "model: Brent ~ GPR + VIX + DXY"),
R2 = c(g1$r.squared, g$r.squared),
R2_adj = c(g1$adj.r.squared, g$adj.r.squared),
Sigma = c(g1$sigma, g$sigma),
F_stat = c(g1$statistic, g$statistic),
df = c(g1$df, g$df),
df_res = c(g1$df.residual, g$df.residual)
)
porovnanie_modelov |>
kable(digits = 4,
col.names = c("Model", "R²", "Adj. R²", "Sigma", "F Stat.",
"df", "df.residual"),
caption = "Porovnanie jednoduchého a viacnásobného modelu") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | R² | Adj. R² | Sigma | F Stat. | df | df.residual |
|---|---|---|---|---|---|---|
| model1: Brent ~ GPR | 0.0031 | 0.0028 | 23.4959 | 12.5097 | 1 | 4045 |
| model: Brent ~ GPR + VIX + DXY | 0.3962 | 0.3958 | 18.2894 | 884.4787 | 3 | 4043 |
Samotný index geopolitického rizika vysvetľuje iba približne 0,31 % variability ceny Brent. Po pridaní VIX a DXY vzrástol upravený koeficient determinácie na 39,6 %. Rast \(\bar R^2\) potvrdzuje, že pridané premenné model skutočne zlepšili a nejde iba o mechanický nárast \(R^2\). Cena ropy teda súvisí oveľa silnejšie s makroekonomickým prostredím (najmä silou dolára) ako so samotným geopolitickým rizikom.
Vysoká hodnota \(R^2\) sama osebe nezaručuje správnu špecifikáciu modelu ani kauzálny vzťah medzi premennými.
\[ \hat\sigma = \sqrt{\frac{RSS}{n-k}} \]
Interpretácia: \(\hat\sigma\) = 18,29 USD/barel. Typická odchýlka skutočnej ceny Brent od ceny predikovanej modelom je teda približne 18,3 USD za barel. Pri priemernej cene 77,62 USD/barel to predstavuje približne 23,6 % priemernej ceny – model je teda na presnú predikciu dennej ceny ropy pomerne nepresný. V porovnaní s jednoduchým modelom (\(\hat\sigma\) = 23,50 USD/barel) sa však chyba znížila.
F-test overuje spoločnú významnosť všetkých vysvetľujúcich premenných. Testujeme hypotézy:
\[ H_0: \beta_1 = \beta_2 = \beta_3 = 0 \]
\[ H_1: \text{aspoň jeden z koeficientov } \beta_1, \beta_2, \beta_3 \neq 0 \]
Testovacia štatistika:
\[ F = \frac{R^2/(k-1)}{(1-R^2)/(n-k)} \sim F(k-1,\; n-k) \]
Všetky tieto hodnoty si overíme ručným výpočtom a doplníme aj kritickú hodnotu F-rozdelenia na hladine 5 %:
df1 <- k - 1
df2 <- n - k
sigma_rucne <- sqrt(RSS / df2)
F_rucne <- (R2_rucne / df1) / ((1 - R2_rucne) / df2)
p_rucne <- pf(F_rucne, df1, df2, lower.tail = FALSE)
F_krit <- qf(0.95, df1, df2)
data.frame(
Ukazovatel = c("Sigma", "F Stat.", "p.value", "df", "df.residual",
"Kritická hodnota F (α = 0,05)"),
Rucny_vypocet = c(sigma_rucne, F_rucne, p_rucne, df1, df2, F_krit),
glance = c(g$sigma, g$statistic, g$p.value, g$df, g$df.residual, NA)
) |>
kable(digits = 4, col.names = c("Ukazovateľ", "Ručný výpočet", "glance()"),
caption = "Overenie charakteristík presnosti a významnosti modelu") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Ukazovateľ | Ručný výpočet | glance() |
|---|---|---|
| Sigma | 18.2894 | 18.2894 |
| F Stat. | 884.4787 | 884.4787 |
| p.value | 0.0000 | 0.0000 |
| df | 3.0000 | 3.0000 |
| df.residual | 4043.0000 | 4043.0000 |
| Kritická hodnota F (α = 0,05) | 2.6071 | NA |
Interpretácia: F-štatistika má hodnotu 884,48, čo je mnohonásobne viac ako kritická hodnota 2,61. P-hodnota je prakticky nulová (R ju zobrazuje ako 0, lebo je menšia ako najmenšie zobraziteľné číslo). Nulovú hypotézu zamietame – aspoň jedna z premenných GPR, VIX, DXY má štatisticky významný vplyv na cenu Brent, model ako celok je štatisticky významný.
Aj jednoduchý model1 je štatisticky významný (F = 12,51, p = 0,000409), hoci vysvetľuje menej ako 1 % variability. Ide o dôležité poučenie: štatistická významnosť neznamená vysokú vysvetľovaciu schopnosť. Pri tisícoch pozorovaní vychádzajú významné aj veľmi slabé vzťahy.
P-hodnoty t-testov aj F-testu platia za klasických predpokladov (nezávislosť a normalita náhodných chýb, konštantný rozptyl). Pri denných časových radoch sú tieto predpoklady často porušené, čo ukáže grafická analýza nižšie. Významnosť preto treba brať s rezervou.
Ak hovoríme o vlastnostiach modelu ako celku, pozrime sa najskôr na štandardné diagnostické grafy funkcie plot(). Červená LOWESS krivka vzniká lokálnym vyhladzovaním pomocou váženej metódy najmenších štvorcov.
# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))
# Vykresliť všetky 4 diagnostické grafy modelu
plot(model)
# Resetovať layout
par(mfrow = c(1, 1))Keďže naše údaje sú denný časový rad, štandardné grafy doplníme o dva grafy v čase, ktoré problém modelu ukážu najzreteľnejšie.
graf_udaje <- data.frame(
date = udaje$date,
Skutocnost = udaje$brent_price,
Predikcia = fitted(model),
Reziduum = residuals(model)
)
ggplot(graf_udaje, aes(x = date)) +
geom_line(aes(y = Skutocnost, colour = "Skutočná cena"), linewidth = 0.5) +
geom_line(aes(y = Predikcia, colour = "Predikcia modelu"), linewidth = 0.5) +
scale_colour_manual(values = c("Skutočná cena" = "darkblue",
"Predikcia modelu" = "firebrick")) +
labs(x = "Rok", y = "Cena Brent (USD/barel)", colour = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top")
Model zachytáva hlavné vlny ceny ropy (napr. prepad koncom roka 2014 a vrchol na jar 2022), ale pri dlhších obdobiach sa od skutočnosti výrazne vzďaľuje: v rokoch 2015–2017 a počas pandémie v roku 2020 predikuje vyššiu cenu, než aká bola, a v rokoch 2011–2014 a 2022–2024 naopak nižšiu.
ggplot(graf_udaje, aes(x = date, y = Reziduum)) +
geom_hline(yintercept = 0, colour = "grey40") +
geom_hline(yintercept = c(-2, 2) * sigma(model),
colour = "firebrick", linetype = "dashed") +
geom_line(colour = "steelblue", linewidth = 0.4) +
labs(x = "Rok", y = "Rezíduum (USD/barel)") +
theme_minimal(base_size = 12)
# Korelácia rezíduí dvoch po sebe idúcich dní
autokor <- cor(residuals(model)[-1], residuals(model)[-n])
Prerušované čiary vyznačujú hranicu \(\pm 2\hat\sigma\). Rezíduá nie sú náhodne rozhádzané okolo nuly, ale tvoria dlhé obdobia, keď model systematicky podhodnocuje (napr. 2012–2013, 2022–2023) alebo nadhodnocuje (2015–2017, 2020) cenu. Korelácia rezíduí dvoch po sebe idúcich dní je 0,994, teda takmer 1.
Ide o silnú autokoreláciu rezíduí, ktorá porušuje predpoklad nezávislosti náhodných chýb. Dôsledkom je, že štandardné chyby koeficientov sú pravdepodobne podhodnotené a p-hodnoty z tabuliek vyššie vychádzajú „optimistickejšie“, než v skutočnosti sú.
Model preto považujeme za užitočný opisný model súvislosti ceny ropy s makroekonomickým prostredím, nie za spoľahlivý prognostický model. Riešením by mohlo byť modelovanie denných výnosov (brent_return) namiesto cien, prípadne zahrnutie oneskorených hodnôt (brent_lag_1), ktoré databáza obsahuje.
Obsah pamäte si uložíme, aby sme ho mohli v ďalšom zadaní načítať príkazom load():
save.image("pamat_ropa.Rdata")