# Spustite raz v konzole, iba ak tieto balíky chýbajú.
install.packages(c("knitr", "rmarkdown", "broom", "kableExtra",
"ggplot2", "dplyr"))Zadanie 4 – Lineárna regresia
Zoznámenie s lineárnou regresiou na údajoch o cenách ropy
1 Predmet výskumu
V tomto zadaní použijeme údaje z databázy Global Oil Prices & Geopolitical Events (2010–2026) z portálu Kaggle (autor Kavya Dhyani, licencia CC BY 4.0). Databáza obsahuje denné ceny ropy Brent a WTI, makroindikátory (index dolára DXY, index strachu VIX), index geopolitického rizika GPR aj anotácie významných geopolitických udalostí.
Budeme postupovať od jednoduchej otázky k regresnému modelu:
- Súvisí vyššia hodnota indexu geopolitického rizika (GPR) s vyššou cenou ropy Brent?
- Grafické znázornenie tejto závislosti.
- Odhad regresnej priamky a jej interpretácia.
- Pridanie ďalších vysvetľujúcich premenných (VIX a DXY).
- Predikcia.
V tomto zadaní sa nevenujem diagnostickým testom.
Na rozdiel od cvičenia s dĺžkou života, kde jeden riadok predstavoval jednu krajinu, tu jeden riadok predstavuje jeden obchodný deň. Regresiu na ňom najprv odhadneme jednoducho, akoby išlo o nezávislé pozorovania; predpoklady, ktoré to vyžaduje, si overím v neskorších cvičeniach.
2 Príprava údajov
2.1 Umiestnenie súboru
Súbor oil_geopolitics_dataset_2010_2026.csv je uložení v podpriečinku udaje. Kód umožňuje aj uloženie CSV priamo k dokumentu. Ak existujú obe kópie, použije sa kópia v udaje.
Modely odhadneme funkciou lm(). Výsledky prevedieme na tabuľky pomocou balíka broom a zobrazíme pomocou známych funkcií kable() a kable_styling(). Na renderovanie dokumentu je potrebné prostredie Quarto s balíkmi knitr a rmarkdown.
library(broom)
library(dplyr)
library(ggplot2)
library(knitr)
library(kableExtra)
# Pomocná funkcia: číslo na slovenský formát (desatinná čiarka).
sk_cislo <- function(x, digits = 3) {
formatC(x, format = "f", digits = digits, decimal.mark = ",")
}
# Pomocná funkcia: člen rovnice so znamienkom (napr. "+ 0,600" alebo "- 0,600").
sk_clen <- function(x, digits = 3) {
znamienko <- ifelse(x < 0, "- ", "+ ")
paste0(znamienko,
formatC(abs(x), format = "f", digits = digits, decimal.mark = ","))
}Pri formátovaní použijeme zápis |>, ktorý odovzdá výsledok predchádzajúceho príkazu ako prvý argument nasledujúcej funkcie. Napríklad tabulka |> kable() zobrazí objekt tabulka pomocou kable().
2.2 Načítanie CSV
V našom súbore sú stĺpce oddelené čiarkou a desatinné čísla používajú bodku. Zo všetkých stĺpcov databázy vyberieme tie, ktoré budeme v modeli potrebovať.
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)
# Dátum prevedieme z textu na triedu Date.
ropa$date <- as.Date(ropa$date)
# Vyberieme stĺpce potrebné pre regresiu.
udaje <- ropa[, c("date", "brent_price", "gpr_index", "vix", "dxy_index")]
prehled <- data.frame(
Pocet_riadkov = nrow(udaje),
Pocet_stlpcov = ncol(udaje),
Prvy_den = min(udaje$date, na.rm = TRUE),
Posledny_den = max(udaje$date, na.rm = TRUE)
)
prehled |>
kable(format = "html", row.names = FALSE,
col.names = c("Počet riadkov", "Počet stĺpcov",
"Prvý deň", "Posledný deň"),
caption = "Základné informácie o vybraných údajoch") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Počet riadkov | Počet stĺpcov | Prvý deň | Posledný deň |
|---|---|---|---|
| 4047 | 5 | 2010-02-17 | 2026-03-12 |
head(udaje) |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Prvých šesť dní vo vybraných údajoch") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| date | brent_price | gpr_index | vix | dxy_index |
|---|---|---|---|---|
| 2010-02-17 | 76.27 | 80.73 | 21.72 | 80.38 |
| 2010-02-18 | 77.78 | 80.73 | 20.63 | 80.40 |
| 2010-02-19 | 78.19 | 80.73 | 20.02 | 80.64 |
| 2010-02-22 | 78.61 | 80.73 | 19.94 | 80.51 |
| 2010-02-23 | 77.25 | 80.73 | 21.37 | 80.85 |
| 2010-02-24 | 78.09 | 80.73 | 20.27 | 80.85 |
Jeden riadok predstavuje jeden obchodný deň. Nemodelujeme prierezovú vzorku, ale vývoj cien a rizika v čase takzvané “časové rady”.
| Premenná | Význam |
|---|---|
| date | Dátum obchodného dňa. |
| brent_price | Cena ropy Brent v USD za barel; vysvetľovaná premenná. |
| gpr_index | Index geopolitického rizika GPR; vyjadruje napätie v svete na základe spravodajského pokrytia. |
| vix | Index VIX; trhové označenie "index strachu", meria očakávanú volatilitu akciového trhu. |
| dxy_index | Index amerického dolára DXY; sila dolára voči košu hlavných mien. |
2.3 Chýbajúce hodnoty
Chýbajúcu hodnotu R označuje ako NA. Najprv zistíme ich počty:
chybajuce <- data.frame(
Premenna = names(udaje),
Pocet_NA = colSums(is.na(udaje))
)
chybajuce |>
kable(format = "html", col.names = c("Premenná", "Počet chýbajúcich hodnôt"),
row.names = FALSE, caption = "Chýbajúce hodnoty pred doplnením") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Počet chýbajúcich hodnôt |
|---|---|
| date | 0 |
| brent_price | 0 |
| gpr_index | 0 |
| vix | 0 |
| dxy_index | 0 |
Aby sme sa mohli sústrediť na regresiu, ponecháme jednoduchý postup z predlohy: chýbajúce hodnoty vysvetľujúcich premenných doplníme mediánom. Je to didaktické zjednodušenie, ktoré môže ovplyvniť výsledky. Chýbajúcu vysvetľovanú premennú takto nenahrádzame; deň s neznámou cenou Brent z analýzy vynecháme.
# Zachováme pôvodné údaje pre prípadnú neskoršiu prácu.
udaje.povodne <- udaje
# Ak chýba vysvetľovaná premenná, taký riadok nepoužijeme.
udaje <- udaje[!is.na(udaje$brent_price), ]
for (premenna in c("gpr_index", "vix", "dxy_index")) {
median_premennej <- median(udaje[[premenna]], na.rm = TRUE)
if (!is.finite(median_premennej)) stop("Premennú nemožno doplniť mediánom.")
chyba <- is.na(udaje[[premenna]])
udaje[[premenna]][chyba] <- median_premennej
}
premenne <- c("brent_price", "gpr_index", "vix", "dxy_index")
opis <- data.frame(
Premenna = premenne,
Minimum = sapply(udaje[premenne], min),
Median = sapply(udaje[premenne], median),
Priemer = sapply(udaje[premenne], mean),
Maximum = sapply(udaje[premenne], max)
)
opis |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Premenná", "Minimum", "Medián", "Priemer", "Maximum"),
caption = "Opisná štatistika po doplnení chýbajúcich hodnôt") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Minimum | Medián | Priemer | Maximum |
|---|---|---|---|---|
| brent_price | 19.33 | 75.17 | 77.62 | 127.98 |
| gpr_index | 58.42 | 93.52 | 102.60 | 318.95 |
| vix | 9.14 | 16.65 | 18.39 | 82.69 |
| dxy_index | 72.93 | 95.31 | 92.91 | 114.11 |
3 Graf
Očakávame, že dni s vyšším geopolitickým rizikom majú spravidla aj vyššiu cenu ropy Brent, pretože napätie vo svete zvyčajne zvyšuje obavy z prerušenia dodávok. Ide zatiaľ o pracovnú hypotézu.
graf_riziko <- ggplot(udaje, aes(x = gpr_index, y = brent_price)) +
geom_point(color = "steelblue", size = 2, alpha = 0.4) +
labs(x = "GPR index (geopolitické riziko)",
y = "Cena ropy Brent (USD/barel)") +
theme_minimal(base_size = 13)
graf_riziko
Body v grafe naznačujú rastúcu závislosť: dni s vyšším indexom GPR majú spravidla vyššiu cenu ropy Brent. Zároveň vidíme, že dni s podobným GPR môžu mať dosť odlišnú cenu; presný tvar závislosti preto posúdime odhadom regresnej priamky v ďalšej kapitole.
4 Jednoduchá lineárna regresia
4.1 Zápis modelu
Najprv použijeme iba jednu vysvetľujúcu premennú:
[ BrentPrice_t = _0 + _1 GPR_t + u_t. ]
- (t) označuje obchodný deň.
- (_0) je konštanta, teda priesečník regresnej priamky s osou Y.
- (_1) je smernica priamky: vyjadruje rozdiel predikovanej ceny Brent v USD/barel pri jednotkovom rozdiele indexu GPR.
- (u_t) zachytáva tú časť ceny ropy, ktorú tento jednoduchý model nevysvetľuje.
4.2 Odhad pomocou lm()
model1 <- lm(brent_price ~ gpr_index, data = udaje)
vysledky1 <- tidy(model1)
vysledky1[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: cena Brent a index GPR") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 73.377 |
| gpr_index | 0.041 |
Funkcia tidy(model1) z balíka broom prevedie výsledky modelu na tabuľku: jeden riadok zodpovedá jednému koeficientu. Stĺpec term obsahuje názov člena modelu a estimate jeho odhad. Označenie (Intercept) znamená konštantu.
tidy() pripravuje údaje o výsledkoch; kable() vytvára tabuľku a kable_styling() upravuje jej vzhľad. Odhad modelu sa formátovaním nemení.
Znak ~ čítame ako „vysvetľujeme pomocou“. Vľavo je vysvetľovaná premenná a vpravo vysvetľujúca premenná. Konštantu R zahrnie automaticky.
Pri našich údajoch získame odhadnutú rovnicu
[ _t = 73,377 + 0,041 , GPR_t. ]
Interpretácia smernice: Rozdiel jedného bodu indexu GPR sa v tomto jednoduchom modeli spája s rozdielom približne 0,04 USD/barel v predikovanej cene ropy Brent.
Interpretácia konštanty: Pri GPR = 0 by priamka predikovala cenu približne 73,38 USD/barel. Nula však môže ležať mimo pozorovaného rozsahu indexu GPR v tejto vzorke, preto konštantu chápeme predovšetkým ako parameter polohy priamky.
Poznámka
Ani zistená závislosť sama o sebe nedokazuje kauzálny účinok. Dni sa líšia aj v mnohých ďalších charakteristikách (ponuka, dopyt, sezónnosť, makroekonomická situácia), ktoré model nezahŕňa.
4.3 Regresná priamka v grafe
K bodovému grafu pridáme priamku pomocou geom_abline(). Jej konštantu (intercept) a smernicu (slope) prevezmeme priamo z odhadnutého modelu model1.
graf_riziko +
geom_abline(intercept = coef(model1)[1],
slope = coef(model1)[2],
color = "firebrick", linewidth = 1)
5 Viacnásobná lineárna regresia
Cena ropy môže súvisieť s viacerými charakteristikami súčasne. Rozšírime model o index VIX (očakávaná trhová volatilita a strach) a index dolára DXY (ropu na svetových trhoch oceňujeme v USD, takže sila dolára môže cenu ovplyvňovať):
[ BrentPrice_t = _0 + _1 GPR_t + _2 VIX_t + _3 DXY_t + u_t. ]
model <- lm(brent_price ~ gpr_index + vix + dxy_index,
data = udaje)
vysledky <- tidy(model)
vysledky[c("term", "estimate")] |>
kable(format = "html", digits = 4, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: GPR, VIX a DXY") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 216.8578 |
| gpr_index | 0.2705 |
| vix | -0.1698 |
| dxy_index | -1.7637 |
Pri našich údajoch dostaneme približne:
[ _t = 216,858 + 0,271 , GPR_t - 0,170 , VIX_t - 1,764 , DXY_t. ]
Koeficient GPR sa teraz interpretuje pri rovnakých hodnotách VIX a DXY. Rozdiel jedného bodu GPR sa pri nezmenených ostatných premenných spája s rozdielom približne 0,27 USD/barel v predikovanej cene Brent.
Všimnite si: Koeficient GPR sa po pridaní ďalších premenných zmenil v porovnaní s jednoduchým modelom. Jednoduchý a viacnásobný model odpovedajú na odlišné otázky: druhý zohľadňuje aj ďalšie charakteristiky obchodných dní. Ak sa znamienko niektorého koeficientu líši od očakávania (napríklad pri DXY by silnejší dolár mohol znamenať lacnejšiu ropu v USD), premennú nevyradíme iba preto, že znamienko nezodpovedá intuícií; otázke vyvrátenia alebo potvrdenia vzťahu sa budeme venovať pri štatistickom testovaní.
6 Stručný súhrn modelu pomocou glance()
Zatiaľ čo tidy() poskytuje riadok pre každý koeficient, glance() poskytuje jeden riadok so súhrnnými charakteristikami celého modelu. Zvolíme počet pozorovaní, (R^2) a odhad štandardnej odchýlky chyby.
suhrn <- glance(model)
suhrn[c("nobs", "r.squared", "sigma")] |>
kable(format = "html", digits = c(0, 3, 3), row.names = FALSE,
col.names = c("Počet dní", "R²", "Odhad SD chyby (USD/barel)"),
caption = "Základné charakteristiky viacnásobného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Počet dní | R² | Odhad SD chyby (USD/barel) |
|---|---|---|
| 4047 | 0.396 | 18.289 |
nobsje počet obchodných dní použitých pri odhade.r.squaredje podiel výberovej variability ceny Brent zachytený modelom; v našom modeli je to približne 39,6 %.sigmaje odhad štandardnej odchýlky modelovej chyby, tu vyjadrený v USD/barel.
Neznamená to, že daný podiel predikcií je správnych, ani to nedokazuje kauzalitu. Úplný textový súhrn je stále dostupný príkazom summary(model). tidy(model) okrem názvov a odhadov obsahuje aj std.error, statistic a p.value. Ich vysvetlenie ponechám na hodinu venovanú štatistickému testovaniu. Preto ich v úvodných tabuľkách nezobrazujeme.
7 Skutočná hodnota, predikcia a rezíduum
Pre každý deň môžeme porovnať skutočnú cenu Brent s hodnotou na regresnej ploche:
[ e_t = y_t - y_t. ]
Tento rozdiel nazývame rezíduum. Kladné rezíduum znamená, že skutočná cena je vyššia než predikcia; záporné znamená opak. Rezíduum vieme vypočítať, zatiaľ čo skutočnú modelovú chybu (u_t) nepoznáme.
porovnanie <- data.frame(
Datum = format(udaje$date, "%Y-%m-%d"),
Skutocnost = udaje$brent_price,
Predikcia = fitted(model),
Reziduum = residuals(model)
)
head(porovnanie, 8) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Dátum", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Skutočné a predikované ceny Brent (USD/barel)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Dátum | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|
| 2010-02-17 | 76.27 | 93.24 | -16.97 |
| 2010-02-18 | 77.78 | 93.39 | -15.61 |
| 2010-02-19 | 78.19 | 93.07 | -14.88 |
| 2010-02-22 | 78.61 | 93.31 | -14.70 |
| 2010-02-23 | 77.25 | 92.47 | -15.22 |
| 2010-02-24 | 78.09 | 92.66 | -14.57 |
| 2010-02-25 | 76.29 | 92.79 | -16.50 |
| 2010-02-26 | 77.59 | 93.65 | -16.06 |
Vypracovaný príklad pre prvý riadok tabuľky: Skutočnosť 76,27 − Predikcia 93,24 = -16,97, čo sa zhoduje s rezíduom -16,97 v prvom riadku. Kladné rezíduum by znamenalo, že skutočná cena je vyššia než predikcia; záporné znamená opak.
Funkcia coef(model) vracia koeficienty, fitted(model) vyrovnané hodnoty a residuals(model) rezíduá. Ak budeme neskôr potrebovať maticu vysvetľujúcich premenných, získame ju príkazom model.matrix(model).
8 Predikcia pre zadanú hodnotu
Vráťme sa k jednoduchému modelu s indexom GPR. Ako zadané hodnoty nezvolíme ľubovoľné čísla, ale kvartile nášho súboru: dolný kvartil, medián a horný kvartil GPR. Ide o modelové scenáre („čo by priamka predikovala za takéhoto rizika“), nie o predpoveď budúcej ceny ropy.
nove_udaje <- data.frame(
Scenar = c("Dolný kvartil GPR", "Medián GPR", "Horný kvartil GPR"),
gpr_index = as.numeric(quantile(udaje$gpr_index,
c(0.25, 0.5, 0.75), na.rm = TRUE))
)
nove_udaje$Predikcia <- predict(model1, newdata = nove_udaje)
nove_udaje |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Scenár", "GPR index", "Predikcia ceny Brent (USD/barel)"),
caption = "Predikcie jednoduchého regresného modelu") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Scenár | GPR index | Predikcia ceny Brent (USD/barel) |
|---|---|---|
| Dolný kvartil GPR | 81.39 | 76.75 |
| Medián GPR | 93.52 | 77.25 |
| Horný kvartil GPR | 116.72 | 78.21 |
Ručná kontrola pre medián: dosadíme hodnotu GPR = 93,52 do odhadnutej rovnice a výsledok porovnáme s funkciou predict(). Pri ručnom dosadení pracujeme s nezaokrúhlenými koeficientami, takže oba postupy sa musia presne zhodovať.
gpr_test <- median(udaje$gpr_index)
b <- coef(model1)
pred_rucne <- unname(b[1] + b["gpr_index"] * gpr_test)
pred_funkcia <- unname(predict(model1,
newdata = data.frame(gpr_index = gpr_test)))
# Kontrola, že oba postupy dávajú rovnaký výsledok.
stopifnot(isTRUE(all.equal(pred_rucne, pred_funkcia)))
kontrola <- data.frame(
Postup = c("Ručné dosadenie do rovnice", "Funkcia predict()"),
GPR = gpr_test,
Predikcia = c(pred_rucne, pred_funkcia)
)
kontrola |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Postup", "GPR index", "Predikcia ceny Brent (USD/barel)"),
caption = "Ručná kontrola predikcie pre medián GPR") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Postup | GPR index | Predikcia ceny Brent (USD/barel) |
|---|---|---|
| Ručné dosadenie do rovnice | 93.52 | 77.25 |
| Funkcia predict() | 93.52 | 77.25 |
Všimnite si, že rovnica v texte zobrazuje koeficienty zaokrúhlené na tri desatinné miesta; dosadenie takto zaokrúhlených hodnôt by preto dalo iba približný výsledok. Presná zhoda vyžaduje nezaokrúhlené koeficienty, ako v tabuľke vyššie.
Ide o bodové predikcie, nie o isté výsledky. Skutočný deň s podobným GPR môže mať odlišnú cenu ropy.
9 Zhrnutie
- Odhadli sme jednoduchú regresiu ceny Brent na index geopolitického rizika GPR a rozšírili sme ju o indexy VIX a DXY.
- Smernicu pri GPR interpretujeme ako zmenu predikovanej ceny Brent v USD/barel pri rozdiele jedného bodu indexu; vo viacnásobnom modeli pri rovnakých hodnotách ostatných premenných.
- Viacnásobný model zachytáva približne 39,6 % výberovej variability ceny Brent; zvyšok pripadá na faktory, ktoré model nezahŕňa.
- Rezíduá ukazujú, o koľko sa skutočná cena líši od predikcie v jednotlivých dňoch.
Varovanie
Ide o opisnú regresiu denných časových radov. Stacionaritu údajov ani autokoreláciu rezíduí sme neoverovali, preto výsledky nepovažujeme za overený prognostický model a nevyvodzujeme z nich štatistickú významnosť ani kauzalitu.