Šiesta úloha

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(ggplot2)
library(patchwork)
library(quantmod)

Na začiatok si stiahneme naše dáta a z nich vypočítame logaritmické výnosnosti, s ktorými následne budeme pracovať.

tickers <- c("AAPL", "GLD", "XLE", "SPY")   # Apple, Gold ETF, Energy ETF, S&P500 ETF
getSymbols(tickers, from = "2024-01-01", to = "2025-01-01")

[1] "AAPL" "GLD"  "XLE"  "SPY" 

data <- merge(Cl(AAPL), Cl(GLD), Cl(XLE), Cl(SPY))
colnames(data) <- tickers

ret <- na.omit(diff(log(data)))
colnames(ret) <- paste0(colnames(ret), "_ret")

ret_df <- na.omit(as.data.frame(ret))
head(ret_df)

op <- par(mfrow = c(2, 2), mar = c(4,4,2,1))
for (col in names(ret_df)) {
  boxplot(ret_df[[col]], main = col, xlab = "Value", col = "purple")
}
par(op)

Boxploty znázorňujú rozdelenie denných logaritmických výnosov aktív AAPL, GLD, XLE a SPY. Všetky výnosy sa sústreďujú okolo nuly, pričom najväčšiu volatilitu vykazuje AAPL a najstabilnejšie hodnoty má GLD. Mierne odľahlé body predstavujú extrémne pohyby na trhu, no celkové rozdelenie je vyvážené a vhodné na ďalšiu analýzu.

model <- lm(AAPL_ret ~ GLD_ret + XLE_ret + SPY_ret, data = ret_df)
summary(model)

Call:
lm(formula = AAPL_ret ~ GLD_ret + XLE_ret + SPY_ret, data = ret_df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.044288 -0.006673 -0.000362  0.005717  0.066928 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0003073  0.0007477   0.411   0.6814    
GLD_ret      0.0028525  0.0818840   0.035   0.9722    
XLE_ret     -0.1950403  0.0705065  -2.766   0.0061 ** 
SPY_ret      1.0396059  0.1007431  10.319   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01173 on 247 degrees of freedom
Multiple R-squared:  0.3122,    Adjusted R-squared:  0.3038 
F-statistic: 37.37 on 3 and 247 DF,  p-value: < 2.2e-16

Model vysvetľuje závislosť výnosu akcie AAPL od výnosov GLD, XLE a SPY. Z výsledkov vyplýva, že štatisticky významné premenné sú XLE (p = 0,0061) a SPY (p < 0,001). Premenná GLD významný vplyv nemá. Kladný koeficient pri SPY (1,04) potvrdzuje silnú väzbu akcie Apple na celkový trh, zatiaľ čo záporný koeficient pri XLE (-0,20) naznačuje mierny opačný vzťah s energetickým sektorom. Hodnota R^2=0,31 znamená, že model vysvetľuje približne 31 % variability výnosov akcie Apple.

op <- par(mfrow = c(2,2))
plot(model)
par(op)

library(ggfortify)

autoplot(model, 
         which = 1:4, 
         colour = "#c8a2c8",          # svetlofialové body
         smooth.colour = "#9d4edd",   # jemnejšia fialová LOESS čiara
         smooth.linetype = "solid",
         size = 2,
         label.size = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(color = "#5a189a", face = "bold", hjust = 0.5),
    strip.text = element_text(color = "#5a189a", face = "bold"),
    panel.grid.minor = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA)
  )

Diagnostické grafy ukazujú, že model je dobre špecifikovaný. Reziduá sú rozložené rovnomerne okolo nuly, čo potvrdzuje absenciu výraznej systematickej chyby. Q-Q graf naznačuje len mierne odchýlky od normality na okrajoch rozdelenia. Graf Scale-Location potvrdzuje približne konštantnú varianciu rezíduí (homoskedasticitu) a graf Residuals vs Leverage neodhalil žiadne významné vplyvné pozorovania. Celkovo možno model považovať za stabilný a spoľahlivý.Graf Cook’s distance naznačuje, že žiadne pozorovanie výrazne neovplyvňuje stabilitu modelu, keďže všetky hodnoty zostávajú bezpečne pod hranicou významného vplyvu.

# Jarque–Bera
jb_test <- tseries::jarque.bera.test(residuals(model)); jb_test

    Jarque Bera Test

data:  residuals(model)
X-squared = 279.24, df = 2, p-value < 2.2e-16

# Outlier test (Bonferroni)
outlier_car <- car::outlierTest(model); outlier_car

Jarque-Bera test ukázal, že reziduá modelu nie sú normálne rozdelené (p-hodnota < 0,001), čo naznačuje prítomnosť extrémnych hodnôt. Test odľahlých pozorovaní (Bonferroni) identifikoval tri dni – 21. 3. 2024, 3. 5. 2024 a 11. 6. 2024 – ako štatisticky významné odľahlé pozorovania, ktoré mohli ovplyvniť výsledky regresie.

model2 <- lm(AAPL_ret ~ I(log(1 + SPY_ret)) + XLE_ret, data = ret_df)
summary(model2)

Call:
lm(formula = AAPL_ret ~ I(log(1 + SPY_ret)) + XLE_ret, data = ret_df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.044321 -0.006706 -0.000423  0.005711  0.066901 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          0.0003426  0.0007431   0.461  0.64516    
I(log(1 + SPY_ret))  1.0402443  0.0980078  10.614  < 2e-16 ***
XLE_ret             -0.1953592  0.0696616  -2.804  0.00544 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0117 on 248 degrees of freedom
Multiple R-squared:  0.313, Adjusted R-squared:  0.3074 
F-statistic: 56.48 on 2 and 248 DF,  p-value: < 2.2e-16

op <- par(mfrow = c(2,2))
plot(model2)
par(op)

# Normalita + outliers pre model2
tseries::jarque.bera.test(residuals(model2))

    Jarque Bera Test

data:  residuals(model2)
X-squared = 280.67, df = 2, p-value < 2.2e-16

car::outlierTest(model2)

autoplot(model2, 
         which = 1:4, 
         colour = "#c8a2c8",          # svetlofialové body
         smooth.colour = "#9d4edd",   # jemnejšia fialová LOESS čiara
         smooth.linetype = "solid",
         size = 2,
         label.size = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(color = "#5a189a", face = "bold", hjust = 0.5),
    strip.text = element_text(color = "#5a189a", face = "bold"),
    panel.grid.minor = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA)
  )

Po úprave modelu (logaritmická transformácia SPY_ret) sa štruktúra rezíduí zlepšila, hoci Jarque-Bera test stále naznačuje mierne odchýlky od normality. Outlier test opäť identifikoval rovnaké tri dni – 21. 3. 2024, 3. 5. 2024 a 11. 6. 2024 – ako odľahlé pozorovania, avšak ich vplyv na model je už menší. Celkovo možno povedať, že transformácia prispela k stabilnejšiemu a robustnejšiemu modelu.

# Pôvodný model
p1 <- ggplot(ret_df, aes(x = SPY_ret, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "purple") +
  labs(x = "SPY_ret", y = "Squared residuals", title = "Squared residuals vs SPY_ret") +
  theme_minimal()

p2 <- ggplot(ret_df, aes(x = XLE_ret, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "purple") +
  labs(x = "XLE_ret", y = "Squared residuals", title = "Squared residuals vs XLE_ret") +
  theme_minimal()

p1 + p2

Grafy zobrazujú závislosť štvorcov rezíduí od výnosov SPY_ret a XLE_ret s cieľom overiť prítomnosť heteroskedasticity. Červená vyhladzovacia krivka je takmer vodorovná, čo naznačuje, že rozptyl rezíduí zostáva približne konštantný pri rôznych hodnotách oboch premenných. Nevidíme žiadny rozširujúci sa alebo zužujúci tvar, preto možno konštatovať, že v modeli sa heteroskedasticita nevyskytuje.

# Upravený model (model2) s log(1+SPY_ret)
p3 <- ggplot(ret_df, aes(x = log(1 + SPY_ret), y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "purple") +
  labs(x = "log(1 + SPY_ret)", y = "Squared residuals", title = "Residuals vs log(1+SPY_ret)") +
  theme_minimal()

p4 <- ggplot(ret_df, aes(x = XLE_ret, y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "purple") +
  labs(x = "XLE_ret", y = "Squared residuals", title = "Residuals vs XLE_ret (model2)") +
  theme_minimal()

p3 + p4

Grafy zobrazujú závislosť štvorcov rezíduí od premenných log(1 + SPY_ret) a XLE_ret po logaritmickej transformácii v modeli 2. Červená vyhladzovacia krivka zostáva takmer vodorovná, čo naznačuje, že rozptyl rezíduí sa s hodnotami vysvetľujúcich premenných výrazne nemení. Z toho vyplýva, že transformácia úspešne odstránila heteroskedasticitu a model má konštantný rozptyl chýb.

lmtest::bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 0.48024, df = 3, p-value = 0.9232

lmtest::bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 0.48542, df = 2, p-value = 0.7845

Výsledky Breusch–Pagan testu pre oba modely ukazujú, že p-hodnoty (0,9232 a 0,7845) sú výrazne vyššie než 0,05. To znamená, že nezamietame nulovú hypotézu o homoskedasticite rezíduí. Inými slovami, rozptyl chýb je konštantný a v modeloch nie je prítomná heteroskedasticita, čo potvrdzuje ich spoľahlivosť.

# Robustné SE pre model2 (doporučené reportovať)
lmtest::coeftest(model2, vcov = sandwich::vcovHC(model2, type = "HC1"))

t test of coefficients:

                       Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)          0.00034261  0.00073427  0.4666  0.641193    
I(log(1 + SPY_ret))  1.04024435  0.09913013 10.4937 < 2.2e-16 ***
XLE_ret             -0.19535922  0.06000130 -3.2559  0.001288 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Tabuľka zobrazuje výsledky robustnej regresie s upravenými smerodajnými chybami (HC1), ktoré zohľadňujú možnú heteroskedasticitu.

Premenná log(1 + SPY_ret) má silne pozitívny a štatisticky významný vplyv (p < 0.001) na výnos akcie AAPL, čo znamená, že Apple sa správa veľmi podobne ako trhový index SPY. Premenná XLE_ret má negatívny a štatisticky významný vplyv (p < 0.01), teda rast energetického sektora súvisí s miernym poklesom výnosov Apple. Konštanta (intercept) nie je štatisticky významná, čo naznačuje, že pri nulových výnosoch SPY a XLE je priemerný výnos AAPL blízky nule.

# Durbin–Watson (autokorelácia rezíduí)
lmtest::dwtest(model2)

    Durbin-Watson test

data:  model2
DW = 2.0406, p-value = 0.6239
alternative hypothesis: true autocorrelation is greater than 0

# RESET (nelinearity/špecifikácia)
lmtest::resettest(model2, power = 2:3, type = "fitted")

    RESET test

data:  model2
RESET = 1.3159, df1 = 2, df2 = 246, p-value = 0.2701

# VIF (multikolinearita)
car::vif(model2)

I(log(1 + SPY_ret))             XLE_ret 
           1.109609            1.109609 

Výsledky doplnkových testov potvrdzujú, že vytvorený model je štatisticky spoľahlivý a správne špecifikovaný. Hodnota Durbin–Watsonovho testu (DW = 2.04, p = 0.62) naznačuje, že v rezíduách sa nenachádza autokorelácia, čo znamená, že chyby nie sú navzájom závislé a model spĺňa predpoklad nezávislosti rezíduí. RESET test s p-hodnotou 0.27 nepreukázal žiadne problémy so špecifikáciou modelu, teda model má správne zvolenú funkčnú formu a nie je potrebné pridávať nelineárne členy. Výsledky VIF (Variance Inflation Factor) s hodnotou 1.1 pre obe vysvetľujúce premenné potvrdzujú, že medzi nimi neexistuje multikolinearita. Celkovo možno konštatovať, že model je dobre špecifikovaný, rezíduá sú nezávislé a vysvetľujúce premenné nie sú navzájom závislé.

library(knitr)
library(dplyr)

results <- data.frame(
  Test = c("Durbin–Watson", "RESET", "VIF (I(log(1+SPY_ret)))", "VIF (XLE_ret)"),
  Statistic = c(2.0406, 1.3159, 1.1096, 1.1096),
  `p-hodnota` = c(0.6239, 0.2701, NA, NA),
  Interpretácia = c(
    "Bez autokorelácie rezíduí",
    "Bez problémov so špecifikáciou modelu",
    "Žiadna multikolinearita",
    "Žiadna multikolinearita"
  )
)

kable(results, caption = "Súhrn doplnkových diagnostických testov pre model2")

Súhrn doplnkových diagnostických testov pre model2

Test	Statistic	p.hodnota	Interpretácia
Durbin–Watson	2.0406	0.6239	Bez autokorelácie rezíduí
RESET	1.3159	0.2701	Bez problémov so špecifikáciou modelu
VIF (I(log(1+SPY_ret)))	1.1096	NA	Žiadna multikolinearita
VIF (XLE_ret)	1.1096	NA	Žiadna multikolinearita

Tabuľka zobrazuje výsledky doplnkových diagnostických testov pre model 2. Durbin–Watson test potvrdzuje, že rezíduá nie sú autokorelované. RESET test ukazuje, že model je správne špecifikovaný a neobsahuje nelineárne chyby. Hodnoty VIF sú veľmi nízke, čo znamená, že medzi vysvetľujúcimi premennými neexistuje multikolinearita. Celkovo možno model považovať za stabilný a štatisticky spoľahlivý.