Analiza portfela kryptowalut w oparciu o modele GARCH

Streszczenie

Celem badania jest analiza porównawcza ryzyka rozumianego jako oszacowanie funkcji warunkowej wariancji w modelach klasy GARCH związanego z portfelem kryptowalut zbudowanym na bazie czterech aktywów: Algorand, Avalanche, PancakeSwap oraz PolkaDot w oparciu o metodologię market-cap-weighted-portfolio. Badanie obejmuje estymację warunkowej wariancji oraz wyznaczenie wartości narażonej na ryzyko (Value at Risk, VaR) dla wybranych modeli: GARCH(1,1), eGARCH(1,1) oraz GARCH-t(1,1). Analiza prowadzona jest na danych dziennych w przedziale od października 2020 do czerwca 2025, z wyodrębnieniem próby in-sample i out-of-sample. Dodatkowo porównano kwantyle teoretyczne i empiryczne, uwzględniając leptokurtyczność oraz niestacjonarność danych. Całość uzupełnia analiza wrażliwości VaR w zależności od długości okresu prognozy.

Wstęp

Kryptowaluty to coraz popularniejsze narzędzie inwestycyjne wykorzystywane zarówno przez uczestników rynków finansowych, jak i inwestorów prywatnych. Często uznawane za “pieniądz przyszłości”. Nazwa ta wzięła się z wynikających różnić nad tradycyjnymi walutami fiducjarnymi emitowanymi przez banki centralne. Wymienia się tutaj brak presji inflacyjnej (co jest charakterystyczne dla wspomnianych walut), wygodę użytkowania, szybkośc transakcji, a także przejrzystość i bezpieczeństwo, jakie zapewnia technologia blockchain (Marszałek P., 2019).

By lepiej zrozumieć istotę kryptowalut, warto przywołać pojęcie “waluty wirtualnej”. Europejski Bank Centralny w 2015 r. zdefiniował ją jako cyfrową reprezentację wartości, która nie jest emitowana przez bank centralny, instytucję kredytową ani instytucję pieniądza elektrycznego, a w określonych warunkach może być wykorzystywana jako substytut tradycyjnych środków płatniczych (Marszałek P., 2019).

Nawiązując jednakże do strony technologicznej kryptowaluty funkcjonują jako zdecentralizowane systemy księgowe. Informacje o stanie ich posiadania zapisywane są w postaci umownych jednostek wartości, opierając na kryptografii - technologii blockchain. Dostęp do środków posiada wyłącznie właściciel potrfela, który dysponuje odpowiednim kluczem prywatnym. Dzięki mechanizmowi sieci peer-to-peer i adekwatnym zabezpieczeniom technologicznym możliwe jest eliminowanie tzw. podwójnego wydatkowania tych samych jednostek. Transakcje realizowane są w oparciu o specjalne aplikacje, a cała sieć odpowiada za ich weryfikację i dystrybucję informacji (Marszałek P., 2019).

Z tych względów kryptowaluty postrzegane są jako narzędzie, które w przyszłości może przyczynić się do reformy systemów monetarnych. Niemniej jednak, zyskują one także opinie kontrowersyjnych aktywów. Spowodowane jest to niejasnym statusie prawnym, trudnego do zrozumienia dla przeciętnego użytkownika mechanizmu ich emisji i transferu, jak i występującej wysokiej zmienności ich wartości.

Wahania są istotnym czynnikiem ryzyka inwestycyjnego, a badaniu zmienności, która występuje w świecie kryptowalut oparto cel niniejszego raportu. Wykorzystano analizę porównawczą ryzyka inwestycyjnego mierzonego poprzez modelowanie funkcji warunkowej wariancji przy zastosowaniu modeli z rodziny GARCH. Badanie dotyczy portfela składającego się z czterech kryptowalut: Algorand, Avalanche, PancakeSwap oraz PolkaDot. Portfel został skonsturowany zgodnie z metodą market-cap-weighted portfolio, co oznacza, że udział każdej kryptowaluty w portfelu jest proporcjonalny do jej kapitalizacji rynkowej.

Celem badania jest analiza porównawcza ryzyka rozumianego jako oszacowanie funkcji warunkowej wariancji w modelach klasy GARCH. Dotyczy ona portfela składającego się z czterech kryptowalut: Algorand, Avalanche, PancakeSwap oraz PolkaDot. Zbudowano go w oparciu o metodę market-cap-weighted portfolio, co oznacza, że udziały procentowe każdej kryptowaluty są proporcjonalne do jej kapitalizacji rynkowej.

Dane

Dzienne ceny i kapitalizacje zostały pobrane ze strony CoinGeko. Każda kryptowaluta sięgała innego okresu, dlatego w celu ujednolicenia wybrano notowania zaczynające się od 21-10-2020, a kończące na 24-06-2025. Dało to łącznie 1707 obserwacji. Budując portfel wyznaczono udziały każdej kryptowaluty w dziennej kapitalizacji rynku, a następnie logarytmiczne stopy zwrotu. Zwrot z portfela ustalono, jako sumę iloczynów zrebalansowanych dziennych wag (przesuniętych o jeden dzień wstecz) oraz zlogarytmizowanej stopy zwortu z kryptowaluty.

# Obliczanie łącznej kapitalizacji i wagi
merged <- merged %>%
  mutate(
    total_cap = MCap_ALG + MCap_AVA + MCap_PAN + MCap_DOT,
    w_alg = MCap_ALG / total_cap,
    w_ava = MCap_AVA / total_cap,
    w_pan = MCap_PAN / total_cap,
    w_dot = MCap_DOT / total_cap
  )

# Obliczanie logarytmicznych dziennych zwrotów z cen
merged <- merged %>%
  mutate(
    r_alg = c(NA, diff(log(Close_ALG))),
    r_ava = c(NA, diff(log(Close_AVA))),
    r_pan = c(NA, diff(log(Close_PAN))),
    r_dot = c(NA, diff(log(Close_DOT)))
  )

# Przesunięcie wag o jeden dzień (rebalansowanie codzienne)
merged <- merged %>%
  mutate(
    w_alg_lag = lag(w_alg),
    w_ava_lag = lag(w_ava),
    w_pan_lag = lag(w_pan),
    w_dot_lag = lag(w_dot)
  )

# Obliczanie portfelowego zwrotu
merged <- merged %>%
  mutate(
    r_portfolio = w_alg_lag * r_alg +
      w_ava_lag * r_ava +
      w_pan_lag * r_pan +
      w_dot_lag * r_dot
  ) %>%
  drop_na()

portfolio_returns <- merged %>%
  select(Date, r_portfolio)

algorand <- algorand %>%
  filter(Date >= as.Date("2020-10-21")) # zerowa kapitalizacja w PancakeSwap

avalanche <- avalanche %>%
  filter(Date >= as.Date("2020-10-21")) 

pancake <- pancake %>%
  filter(Date >= as.Date("2020-10-21")) 

polkadot <- polkadot %>%
  filter(Date >= as.Date("2020-10-21")) 

alograndTS <- xts(algorand$Close_ALG, algorand$Date)
avalancheTS <- xts(avalanche$Close_AVA, avalanche$Date)
pancakeTS <- xts(pancake$Close_PAN, pancake$Date)
polkadotTS <- xts(polkadot$Close_DOT, polkadot$Date)

Powyższe wykresy prezentują ceny zamknięcia, można zauważyć, że najwyższe wartości kryptowaluty przyjmują w okresie 2021 - 2022, gdzie lata te charakteryzują się wysokimi wzrostami, a wraz z końcem 2022 następują gwałtowne spadki wycen. Późniejsze okresy ukazują względny spokój, a nawet stabilizację. Najwyżej wycenianą kryptowalutą jest avalanche, a najniżej algorand.

Powyższy wykres przedstawia zmiany udziałów kryptowalut w kapitalizacji całego portfela na przestrzeni lat 2021 - 2025. Wraz z początkiem okresu portfel był zdominowany przez DOT (PolkaDot), której udział przekraczał nawet 75%. Z czasem jego znaczenie systematycznie malało, a coraz bardziej zaczęł odznaczać się Avalanche, który obecnie stanowi najwyżej kapitalizowaną kryptowalutę w tym portfelu. Jego udział stanowi ok. 50% całości portfela.

Z kolei Algorand (ALG) oraz PancakeSwap (PAN) utrzymują stosunkowo niski i stabilny udział w strukturze portfela. ALG zazwyczaj oscyluje w granicach 5-15%, a PAN nie przekracza 5%.

Analizowanie zmienności

Typowe modele szeregów czaswoych, takie jak ARMA czy ARIMA, opierają się na założeniu homoskedastyczności. Oznacza to, że wariancja warunkowa jest stała w czasie. Założenie to chociaż upraszcza analizę, to w kontekście danych finansowych okazuje się trudne do utrzymania. Obserwując szeregi dot. instrumentów finansowych można zauważyć kilka powtarzających się własności, które w literaturze przedmiotu znane sa jako “stylizowane fakty”. Dane te bowiem cechują się zjawiskiem grupowania zmienności - okresy wysokiej wariancji mają tendencję do występowania w klastrach, podobnie jak okresy stabilności. W konsekwencji, momenty podwyższonej zmienności mogą sugerować podwyższoną warunkową wariancję również w kolejnych obserwacjach. Stylizowane fakty to m.in.

1. Leptokurtyczność rozkładów - rozkłady stóp zwrotów zwrotu z aktywów, w porównaniu z rozkładem normalnym mają “grube ogony” i wyższy szczyt funkcji gęstości.

2. Grupowanie warinacji - zarówno małe jak i duże zmiany kursów akcji mają tendencję do występowania seriami.

3. Efekt dźwigni - tendencja wariancji do większego wzrostu na skutek duże spadku cen, lecz jednocześnie niższego wzrostu w przypadku wzrostu cen o tą samą wielkość.

Z pomocą przychodzą modele GARCH (Generalized AutoRegressive Conditional Heteroskedasticity), które zostały zaprojektowane, aby uchwycić tego typu zjawiska. Jego podstawą są modele ARCH (AutoRegressive Conditional Heteroskedasticity), w którym bieżąca warunkowa wariancja zależy od przeszłych wartości reszt modelu. W GARCH koncept ten jest rozwijany, umożliwiając jednoczesne uwzględnienie zarówno opóźnionych składników wariancji, jak i wcześniejszych kwadratów reszt. Pozwala to na bardziej elastyczne odwzorowanie struktur zmienności obserwowanych w rzeczywistych szeregach finansowych. Tym samym GARCH umożliwia analizę ryzyka, prognozowanie zmienności oraz szacowanie wartości zagrożonej (VAR), która w późniejszej części pracy zostanie przedstawiona.

Normalność rozkładu

Wykres przedstawia porónwanie empirycznego rozkładu dziennych logarytmicznych dla poszczególnych składników portfela. Najbardziej powszechnym zjawiskiem jest tutaj leptokurtyczność, czyli występowanie wyższego szczytu i grubszych ogonów niż przy rozkładzie normalnym, który prezentowany jest zieloną linią. Szczególnie wyraźnie odchylenie od normalności można zaobserwować w przypadku PAN, podobne lecz mniej nasilone cechy widoczne są w przypadku AVA i ALG. DOT jest natomiast najbardziej zbliżony do rozkładu normalnego, aczkolwiek nadal występują dosyć grube ogony. Taka struktura wskazuje na wyższe prawdopodobieństwo wystąpienie ekstremalnych zmian cen.

Na wykresie przedstawiono empiryczny rozkład dziennych logarytmicznych stóp zwrotu dla całego portfela oraz jego porównanie z teoretycznym rozkładem normalnym. Widoczne jest, że empiryczna krzywa gęstości znacznie odbiega od rozkładu normalnego, podobnie jak w przypadku zwrotów z poszczególnych składników badanego portfela. Tym samym również wykazywana jest leptokurtyczność.

basicStats(portfolio_returns$r_portfolio)

##             X..portfolio_returns.r_portfolio
## nobs                             1707.000000
## NAs                                 0.000000
## Minimum                            -0.449900
## Maximum                             0.317012
## 1. Quartile                        -0.026505
## 3. Quartile                         0.026355
## Mean                               -0.000006
## Median                              0.001020
## Sum                                -0.009536
## SE Mean                             0.001244
## LCL Mean                           -0.002445
## UCL Mean                            0.002434
## Variance                            0.002641
## Stdev                               0.051391
## Skewness                           -0.306192
## Kurtosis                            7.114675

Przyglądając się podstawowym statystykom stóp zwrotu z portfela, można zauważyć, że średnia stopa zwrotu portfela jest zbliżona do zera. Mediana (0.001) natomiast przewyższa średnią, co może sugerować niewielką lewostronną asymetrię.

Rozstęp wartości wskazuje na bardzo wysoką zmienność. Największy dzienny spadek wyniósł aż -44,99%, a największy wzrost 31.70%. Są to dosyć ekstremalne wartości, charakterystyczne dla aktywów cechujących się dużą niestabilnością.

Dodatkowo, kurtoza wynosi 7.11, co wksazuje na wyraźnie grubsze ogony niż rozkład normalny (kurtoza = 3). Oznacza to wyższe prawdopodobieństwo występowania dużych strat, jak i zysków. Skośność -0.306, podobnie jak w przypadku wyższej mediany niż średniej wskazuje, że rozkład jest lewostronnie skośny.

qqnorm(portfolio_returns$r_portfolio)
qqline(portfolio_returns$r_portfolio, col = 2)

options(scipen = 999)
jarque.bera.test(na.omit(portfolio_returns$r_portfolio))

## 
##  Jarque Bera Test
## 
## data:  na.omit(portfolio_returns$r_portfolio)
## X-squared = 3639, df = 2, p-value < 0.00000000000000022

Na końcu przeprowadzono test Jarque-Bera, który silnie odrzuca hipotezę zerową mówiąca o rozkładzie normalnym. Potwierdza to wcześniejsze obserwacje z histogramów i wykresów gęstości.

Grupowanie wariancji

ggplot(portfolio_returns, aes(x = Date, y = r_portfolio)) +
  geom_line(linewidth = 0.3, color = "#411A1BFF") +
  labs(
    x = "Data",
    y = "Stopa zwrotu (%)",
    title = "Stopy zwrotu z całego portfela inwestycyjnego"
  ) +
  scale_x_date(date_minor_breaks = "1 year") +
  geom_vline(
    xintercept = as.Date("2024-01-01"),
    color = "#ff0000", linetype = "longdash",
    linewidth = 0.9
  )

Wykres przedstawia dzienne stopy zwrotu z portfela, zauważalne są wyraźne epizody zwiększonej zmienności, szczególnie w pierwszej części próby (2020 - 2022). Później natomiast obserwowana jest stopniowa stabilizacja. Widoczne są również okresy grupowania wariancji, kiedy to obserwowana jest wyższa zmienność. Zauważalna pionowa przerywana linia, to wprowadzony podział odnośnie do późniejszych badań, gdzie obserwacje na lewo stanowią próbę in-sample, a obserwacje na prawo próbę out-of-sample.

ACF i PACF

# wykres wartości ACF dla zwrotów
acf(portfolio_returns$r_portfolio, lag.max = 36, na.action = na.pass,
    ylim = c(-0.05, 0.2), # zmieniamy zakres wartości dla osi pionowej
    col = "darkblue", lwd = 7,
    main = "ACF zwrotów portfela")

# wykres wartości ACF dla kwadratów zwrotów
acf(portfolio_returns$r_portfolio^2, lag.max = 36, na.action = na.pass,
    ylim = c(-0.05, 0.7), # zmieniamy zakres wartości dla osi pionowej
    col = "darkblue", lwd = 7,
    main = "ACF kwadratów zwrotów portfela")

pacf(portfolio_returns$r_portfolio, lag.max = 36, na.action = na.pass,
    ylim = c(-0.05, 0.2), # zmieniamy zakres wartości dla osi pionowej
    col = "darkblue", lwd = 7,
    main = "PACF zwrotów portfela")

Wykres ACF dziennych stóp zwortu wskazuje, że zależności pomiędzy kolejnymi obserwacjami są relatywnie słabe. Jedynie 5, 8 i 9 opóźnienie przekracza granice istotności statystycznej. Może to sugerować pewną obecność efektów autoregresyjnych w danych, jednak nie są one dominujące.

Mimo to, w kontekście modelowania warunkowej wariancji bardziej interesujące są obecności zależności w drugich momentach rozkładów. Na drugim wykresie znacznie bardziej wyraźna jest struktura korelacyjna kwadratów stóp zwrotów. Wysoka i długo utrzymująca się autokorelacja świadczy o obecności efektów ARCH. Sugeruje to, że wariancja nie jest stała w czasie, lecz wykazuje zależności od własnej przeszłości, co uzasadnia stosowanie GARCH.

ArchTest(portfolio_returns$r_portfolio, lags = 5)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  portfolio_returns$r_portfolio
## Chi-squared = 152.3, df = 5, p-value < 0.00000000000000022

Test ARCH-LM jednoznacznie potwierdza, że w danych z portfela występują niestacjonarności w wariancji (heteroskedastyczność), tym samym obecne są efekty ARCH. Uzasadnia to dalsze modelowanie zmienności za pomocą modeli GARCH, eGARCH lub GARCH-t.

Modelowanie zmienności

Modele ARCH

# ARCH(1)
k.arch1 <- garchFit(formula = ~ garch(1, 0),
                    data = na.omit(portfolio_returns$r_portfolio),
                    cond.dist = "norm", # rozkład warunkowy reszt
                    trace = FALSE) 

# podsumowanie wyników i kilka testów diagnostycznych
summary(k.arch1)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 0), data = na.omit(portfolio_returns$r_portfolio), 
##     cond.dist = "norm", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 0)
## <environment: 0x7fce9b52e6a0>
##  [data = na.omit(portfolio_returns$r_portfolio)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##            mu          omega         alpha1  
## -0.0000056055   0.0021184161   0.2157934130  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##            Estimate   Std. Error  t value             Pr(>|t|)    
## mu     -0.000005606  0.001164516   -0.005                0.996    
## omega   0.002118416  0.000097508   21.726 < 0.0000000000000002 ***
## alpha1  0.215793413  0.042998419    5.019           0.00000052 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2682.607    normalized:  1.571533 
## 
## Description:
##  Sat Jun 28 14:11:55 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                    Statistic              p-Value
##  Jarque-Bera Test   R    Chi^2  3036.0648972 0.000000000000000000
##  Shapiro-Wilk Test  R    W         0.9433305 0.000000000000000000
##  Ljung-Box Test     R    Q(10)    31.8493028 0.000424258369119412
##  Ljung-Box Test     R    Q(15)    39.5362946 0.000532610689546131
##  Ljung-Box Test     R    Q(20)    43.9588263 0.001524154062077177
##  Ljung-Box Test     R^2  Q(10)   108.1086336 0.000000000000000000
##  Ljung-Box Test     R^2  Q(15)   135.9750660 0.000000000000000000
##  Ljung-Box Test     R^2  Q(20)   167.4520014 0.000000000000000000
##  LM Arch Test       R    TR^2     83.2058958 0.000000000001005973
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -3.139551 -3.129986 -3.139557 -3.136011

# stała w równaniu średniej jest nieistotna, więc możemy ją wykluczyć z modelu
k.arch1 <- garchFit(formula = ~ garch(1, 0),
                    data = na.omit(portfolio_returns$r_portfolio),
                    include.mean = F,
                    cond.dist = "norm",
                    trace = FALSE)

# podsumowanie wyników i kilka testów diagnostycznych
summary(k.arch1)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 0), data = na.omit(portfolio_returns$r_portfolio), 
##     cond.dist = "norm", include.mean = F, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 0)
## <environment: 0x7fce8eef7078>
##  [data = na.omit(portfolio_returns$r_portfolio)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##     omega     alpha1  
## 0.0021184  0.2157906  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value             Pr(>|t|)    
## omega  0.0021184   0.0000975   21.728 < 0.0000000000000002 ***
## alpha1 0.2157905   0.0429942    5.019          0.000000519 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2682.606    normalized:  1.571533 
## 
## Description:
##  Sat Jun 28 14:11:55 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                    Statistic              p-Value
##  Jarque-Bera Test   R    Chi^2  3036.3243729 0.000000000000000000
##  Shapiro-Wilk Test  R    W         0.9433296 0.000000000000000000
##  Ljung-Box Test     R    Q(10)    31.8492620 0.000424265011933378
##  Ljung-Box Test     R    Q(15)    39.5360483 0.000532656060662884
##  Ljung-Box Test     R    Q(20)    43.9584706 0.001524321155485731
##  Ljung-Box Test     R^2  Q(10)   108.1019873 0.000000000000000000
##  Ljung-Box Test     R^2  Q(15)   135.9641445 0.000000000000000000
##  Ljung-Box Test     R^2  Q(20)   167.4320655 0.000000000000000000
##  LM Arch Test       R    TR^2     83.2000707 0.000000000001008527
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -3.140722 -3.134345 -3.140725 -3.138362

k.arch3 <- garchFit(formula = ~ garch(3, 0),
                    data = na.omit(portfolio_returns$r_portfolio),
                    include.mean = F,
                    cond.dist = "norm",
                    trace = FALSE)
summary(k.arch3)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(3, 0), data = na.omit(portfolio_returns$r_portfolio), 
##     cond.dist = "norm", include.mean = F, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(3, 0)
## <environment: 0x7fce8e67ee10>
##  [data = na.omit(portfolio_returns$r_portfolio)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##     omega     alpha1     alpha2     alpha3  
## 0.0015646  0.1234058  0.2756343  0.0108288  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value             Pr(>|t|)    
## omega  0.00156462  0.00009147   17.105 < 0.0000000000000002 ***
## alpha1 0.12340578  0.02838100    4.348     0.00001372700050 ***
## alpha2 0.27563431  0.04042796    6.818     0.00000000000924 ***
## alpha3 0.01082884  0.02668125    0.406                0.685    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2759.233    normalized:  1.616422 
## 
## Description:
##  Sat Jun 28 14:11:55 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic         p-Value
##  Jarque-Bera Test   R    Chi^2  992.1917120 0.0000000000000
##  Shapiro-Wilk Test  R    W        0.9644916 0.0000000000000
##  Ljung-Box Test     R    Q(10)   22.0862762 0.0146710602440
##  Ljung-Box Test     R    Q(15)   27.7695973 0.0230565851566
##  Ljung-Box Test     R    Q(20)   30.6013253 0.0606784579718
##  Ljung-Box Test     R^2  Q(10)   37.8685729 0.0000400036910
##  Ljung-Box Test     R^2  Q(15)   55.8924819 0.0000012634665
##  Ljung-Box Test     R^2  Q(20)   65.6816412 0.0000009084538
##  LM Arch Test       R    TR^2    38.3436318 0.0001349044076
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -3.228158 -3.215405 -3.228169 -3.223438

Początkowo wypróbowano model ARCH(1) oraz ARCH(3). Testy Ljung-Boxa wskazują na istnienie autokorelacji w resztach oraz kwadratach reszt, tym samym mimo, że modele są statystycznie istotne, to nie uchwycają w pełni struktury zmienności portfela.

plot(k.arch3, which = 10)

plot(k.arch3, which = 11)

Przypatrując się wykresowi ACF standaryzowanych reszt oraz ich kwadratowi, potwierdza się, że autokorelacja w resztach ciągle występuje, dlatego modele te nie są adekwatne do tej analizy.

Dostosowanie modelu GARCH

k.garch11 <- garchFit(formula = ~ garch(1, 1),
                      data = na.omit(portfolio_returns$r_portfolio),
                      include.mean = F,
                      cond.dist = "norm",
                      trace = FALSE)
k.garch12 <- garchFit(formula = ~ garch(1, 2),
                      data = na.omit(portfolio_returns$r_portfolio),
                      include.mean = F,
                      cond.dist = "norm",
                      trace = FALSE)

k.garch21 <- garchFit(formula = ~ garch(2, 1),
                      data = na.omit(portfolio_returns$r_portfolio),
                      include.mean = F,
                      cond.dist = "norm",
                      trace = F)

k.garch22 <- garchFit(formula = ~ garch(2, 2),
                      data = na.omit(portfolio_returns$r_portfolio),
                      include.mean = F,
                      cond.dist = "norm",
                      trace = F)

Następnie wyestymowano modele GARCH(1,1), GARCH(1,2), GARCH(2,1) oraz GARCH(2,2) i za pośrednictwem kryterium informacyjnych wybrano najbardziej porządany zestaw parametrów.

compare.ICs(c("k.arch1", "k.arch3", "k.arch3", "k.garch11", "k.garch12", "k.garch21", "k.garch22"))

##         AIC       BIC       SIC      HQIC     model
## 1 -3.140722 -3.134345 -3.140725 -3.138362   k.arch1
## 2 -3.228158 -3.215405 -3.228169 -3.223438   k.arch3
## 3 -3.228158 -3.215405 -3.228169 -3.223438   k.arch3
## 4 -3.296415 -3.286850 -3.296421 -3.292875 k.garch11
## 5 -3.295262 -3.282508 -3.295273 -3.290541 k.garch12
## 6 -3.295264 -3.282511 -3.295275 -3.290544 k.garch21
## 7 -3.294466 -3.278524 -3.294483 -3.288565 k.garch22

Spośród powyższych modeli, najniższym Akaike odznacza się model GARCH (1,1) oraz GARCH (2,1).

summary(k.garch11)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = na.omit(portfolio_returns$r_portfolio), 
##     cond.dist = "norm", include.mean = F, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7fce8e8b1078>
##  [data = na.omit(portfolio_returns$r_portfolio)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##       omega       alpha1        beta1  
## 0.000045167  0.094012746  0.893457825  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value             Pr(>|t|)    
## omega  0.00004517  0.00001493    3.026              0.00248 ** 
## alpha1 0.09401275  0.01437755    6.539       0.000000000062 ***
## beta1  0.89345783  0.01581368   56.499 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2816.49    normalized:  1.649965 
## 
## Description:
##  Sat Jun 28 14:11:56 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic    p-Value
##  Jarque-Bera Test   R    Chi^2  438.3381721 0.00000000
##  Shapiro-Wilk Test  R    W        0.9739353 0.00000000
##  Ljung-Box Test     R    Q(10)   19.7662685 0.03154253
##  Ljung-Box Test     R    Q(15)   24.9872294 0.05011520
##  Ljung-Box Test     R    Q(20)   27.9415375 0.11079054
##  Ljung-Box Test     R^2  Q(10)   12.9917679 0.22413244
##  Ljung-Box Test     R^2  Q(15)   17.6067054 0.28390614
##  Ljung-Box Test     R^2  Q(20)   21.4263338 0.37243193
##  LM Arch Test       R    TR^2    15.3858947 0.22100819
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -3.296415 -3.286850 -3.296421 -3.292875

W modelu GARCH(1,1) wszystkie współczynniki są statystycznie istotne, a suma parametrów jest bliska 1, co oznacza wysoką trwałość szoków zmienności. Zmienność ma tendencję do klastrowania się w czasie. Testy diagnostyczne wskazują, że brak jest autokorelacji kwadratach reszt oraz nie występują efekty ARCH. Można zatem stwierdzić, że model został dobrany poprawnie.

summary(k.garch21)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(2, 1), data = na.omit(portfolio_returns$r_portfolio), 
##     cond.dist = "norm", include.mean = F, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(2, 1)
## <environment: 0x7fce9e21f5c0>
##  [data = na.omit(portfolio_returns$r_portfolio)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##       omega       alpha1       alpha2        beta1  
## 0.000045648  0.092633224  0.002034932  0.892683336  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value             Pr(>|t|)    
## omega  0.00004565  0.00001680    2.717             0.006585 ** 
## alpha1 0.09263322  0.02645538    3.501             0.000463 ***
## alpha2 0.00203493  0.03157723    0.064             0.948617    
## beta1  0.89268334  0.01972591   45.254 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2816.508    normalized:  1.649975 
## 
## Description:
##  Sat Jun 28 14:11:56 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic    p-Value
##  Jarque-Bera Test   R    Chi^2  436.6627337 0.00000000
##  Shapiro-Wilk Test  R    W        0.9739555 0.00000000
##  Ljung-Box Test     R    Q(10)   19.7149998 0.03206651
##  Ljung-Box Test     R    Q(15)   24.9242264 0.05097036
##  Ljung-Box Test     R    Q(20)   27.8733179 0.11243228
##  Ljung-Box Test     R^2  Q(10)   12.9801911 0.22478147
##  Ljung-Box Test     R^2  Q(15)   17.6180390 0.28327831
##  Ljung-Box Test     R^2  Q(20)   21.4375828 0.37179113
##  LM Arch Test       R    TR^2    15.4083571 0.21986043
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -3.295264 -3.282511 -3.295275 -3.290544

Drugim modelem pod względem wysokości kryterium Akaike jest GARCH (2,2). Ponownie parametry są bardzo bliskie 1, jednakże alpha2 jest nieistotna statystycznie, co oznacza, że dodatnie tego składnia nie poprawia modelu. Autokorelacja kwadratów reszt, jak i efekty ARCH nie występują, jednakże GARCH(2,1) nie przynosi istotnych korzyści względem GARCH(1,1), dlatego preferowany jest model mniejszy.

plot(k.garch11, which = 10)

plot(k.garch11, which = 11)

Wykresy wystandaryzowanych reszt potwierdzają wyniki uzyskane w testach. Nie ma zatem podstaw, aby twierdzić, że model GARCH(1,1) jest nieodpowiedni.

GARCH (1,1) oraz ARIMA(2,2)

Chcąc możliwie najlepiej zastosować modele rodziny GARCH, przeprowadzono estymację GARCH (1,1) z zastosowaniem odpowiedniego modelu średniej ARIMA. Do wyboru parametrów p i q posłużono się algorytmem, który dopasowuje odpowiednie modele do danych, gdzie za pomocą kryterium informacyjnego Akaike można wybrać najlepszy model.

cl <- makePSOCKcluster(10)
AC <- autoarfima(as.numeric(portfolio_returns$r_portfolio),
                 ar.max = 2, ma.max = 2,
                 criterion = 'AIC',
                 method = 'partial',
                 arfima = FALSE,
                 include.mean = NULL,
                 distribution.model = 'norm',
                 solver = 'solnp',
                 cluster = cl)

head(AC$rank.matrix)

##   AR MA Mean ARFIMA       AIC converged
## 1  2  2    0      0 -3.104624         1
## 2  2  2    1      0 -3.103452         1
## 3  1  0    0      0 -3.098112         1
## 4  0  1    0      0 -3.098055         1
## 5  0  2    0      0 -3.097633         1
## 6  2  0    0      0 -3.097560         1

Na podstawie kryterium informacyjnego AIC jako najlepszy model został wyłoniony ARIMA(2,2).

spec <- ugarchspec(
  # równanie warunkowej wariancji
  variance.model = list(model = "sGARCH",
                        garchOrder = c(1, 1)),
  # równanie warunkowej wartości oczekiwanej
  mean.model = list(armaOrder = c(2, 2),
                    include.mean = F),
  # zakladany rozklad warunkowy reszt
  distribution.model = "norm") # rozklad normalny

k.garch11_2 <- ugarchfit(spec=spec, data=na.omit(portfolio_returns$r_portfolio))
k.garch11_2

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## ar1     1.412483    0.005698   247.8834 0.000000
## ar2    -0.977260    0.004922  -198.5416 0.000000
## ma1    -1.410286    0.001851  -761.8747 0.000000
## ma2     0.985686    0.000087 11269.9526 0.000000
## omega   0.000045    0.000015     2.9599 0.003077
## alpha1  0.094441    0.014583     6.4760 0.000000
## beta1   0.893028    0.016156    55.2767 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## ar1     1.412483    0.005246   269.2245 0.000000
## ar2    -0.977260    0.005175  -188.8538 0.000000
## ma1    -1.410286    0.000828 -1702.4472 0.000000
## ma2     0.985686    0.000178  5548.0989 0.000000
## omega   0.000045    0.000026     1.7298 0.083659
## alpha1  0.094441    0.026515     3.5618 0.000368
## beta1   0.893028    0.027992    31.9027 0.000000
## 
## LogLikelihood : 2821.024 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.2970
## Bayes        -3.2747
## Shibata      -3.2971
## Hannan-Quinn -3.2888
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.4052  0.5244
## Lag[2*(p+q)+(p+q)-1][11]    6.6507  0.1411
## Lag[4*(p+q)+(p+q)-1][19]   12.9769  0.1112
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.04257  0.8365
## Lag[2*(p+q)+(p+q)-1][5]   1.98844  0.6212
## Lag[4*(p+q)+(p+q)-1][9]   5.15597  0.4069
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2506 0.500 2.000  0.6167
## ARCH Lag[5]    2.9608 1.440 1.667  0.2955
## ARCH Lag[7]    4.5708 2.315 1.543  0.2716
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.0551
## Individual Statistics:              
## ar1    0.10680
## ar2    0.40043
## ma1    0.12472
## ma2    0.05707
## omega  0.32291
## alpha1 0.35009
## beta1  0.30337
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.3870 0.6988    
## Negative Sign Bias  0.3384 0.7351    
## Positive Sign Bias  0.3857 0.6997    
## Joint Effect        0.6029 0.8958    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic  p-value(g-1)
## 1    20     73.30 0.00000002574
## 2    30     87.36 0.00000009232
## 3    40    102.00 0.00000015388
## 4    50    110.43 0.00000123121
## 
## 
## Elapsed time : 0.554914

Model GARCH z częścia ARMA(2,2) wykazuje trochę lepsze dopasowanie niż GARCH(1,1). Również kryterium AIC wskazuję, że dodatnie ARMA(2,2) polepsza model. Sugerując się jednak testem Pearsona, można dojść do wniosku, że rozkład reszt nie jest w pełni zgodny z założeniem normalności, więc alternatywnie warto rozważyć inne rozkład, jak t-Studenta.

Model EGARCH(1,1)

Istotnym zjawiskiem empirycznym w analize szeregów czasowych jest tzw. efekt dźwigni, polegający na asymetrycznym wpływie szoków pozytywnych i negatywnych na zmienność – przy czym spadki cen aktywów częściej generują silniejszy wzrost wariancji warunkowej niż analogiczne zmiany wzrostowe. Patrząc na rozkład stóp zwrotów oraz cen zamknięcia składowych portfela, można dojść do wniosku, że również w przypadku analziowanych danych takie zjawisko występuje. Jednym z rozszerzeń modelu GARCH, które w sposób relatywnie trafny uwzględnia tę właściwość stylizowaną, jest model E-GARCH.

Model ten charakteryzuje się zastosowaniem logarytmu wariancji warunkowej, co umożliwia eliminację konieczności narzucania warunku nieujemności parametrów, przy jednoczesnym zachowaniu interpretowalności i stabilności procesu. Kluczowym komponentem konstrukcyjnym modelu E-GARCH jest asymetryczna funkcja reakcji zmienności na nowe informacje – tzw. krzywa wpływu – której lewe ramię (odpowiadające szokom negatywnym) cechuje się większą stromością niż ramię prawe, co odzwierciedla silniejsze reakcje rynku na niekorzystne wiadomości.

W ramach niniejszego podejścia zdecydowano się na estymację modelu E-GARCH w wariancie (1,1), odpowiadającym standardowej specyfikacji używanej również w klasycznym modelu GARCH. Dodatkowo, w celu uwzględnienia struktury zależności w średniej, model wariancji uzupełniono o komponent ARMA(2,2).

# EGARCH

spec = ugarchspec(
  # rownanie warunkowej wariancji
  variance.model = list(model ="eGARCH",
                        garchOrder = c(1, 1)),
  # rownanie warunkowej sredniej
  mean.model = list(armaOrder = c(2, 2), include.mean = F),
  # zakladany rozklad warunkowy reszt
  distribution.model = "norm")

k.egarch11 <- ugarchfit( spec = spec,
                         data = na.omit(portfolio_returns$r_portfolio))
k.egarch11

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : eGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## ar1     1.415760    0.007812  181.2326 0.000000
## ar2    -0.971744    0.004792 -202.7706 0.000000
## ma1    -1.411412    0.005234 -269.6431 0.000000
## ma2     0.980997    0.000282 3478.6065 0.000000
## omega  -0.142913    0.026430   -5.4072 0.000000
## alpha1  0.023146    0.014069    1.6452 0.099935
## beta1   0.974913    0.004277  227.9505 0.000000
## gamma1  0.210996    0.023917    8.8222 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## ar1     1.415760    0.007498  188.8177 0.000000
## ar2    -0.971744    0.004506 -215.6518 0.000000
## ma1    -1.411412    0.006383 -221.1039 0.000000
## ma2     0.980997    0.000300 3274.0661 0.000000
## omega  -0.142913    0.035785   -3.9936 0.000065
## alpha1  0.023146    0.018750    1.2345 0.217029
## beta1   0.974913    0.005521  176.5914 0.000000
## gamma1  0.210996    0.044709    4.7193 0.000002
## 
## LogLikelihood : 2823.797 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.2991
## Bayes        -3.2736
## Shibata      -3.2992
## Hannan-Quinn -3.2897
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.1993 0.65527
## Lag[2*(p+q)+(p+q)-1][11]    6.8185 0.09124
## Lag[4*(p+q)+(p+q)-1][19]   13.3544 0.08815
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      0.652  0.4194
## Lag[2*(p+q)+(p+q)-1][5]     4.094  0.2425
## Lag[4*(p+q)+(p+q)-1][9]     7.762  0.1433
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.08386 0.500 2.000  0.7721
## ARCH Lag[5]   3.02244 1.440 1.667  0.2865
## ARCH Lag[7]   4.84816 2.315 1.543  0.2403
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.0931
## Individual Statistics:              
## ar1    0.33922
## ar2    0.78167
## ma1    0.06895
## ma2    0.20844
## omega  0.38690
## alpha1 0.13410
## beta1  0.36738
## gamma1 0.13925
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.44520 0.6562    
## Negative Sign Bias 0.08396 0.9331    
## Positive Sign Bias 0.29163 0.7706    
## Joint Effect       0.89015 0.8278    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic   p-value(g-1)
## 1    20     79.68 0.000000002113
## 2    30     97.62 0.000000002341
## 3    40    113.39 0.000000003486
## 4    50    132.28 0.000000001386
## 
## 
## Elapsed time : 0.3666849

W analizowanym przypadku, podobnie jak w uprzednio estymowanych modelach, suma parametrów modelu GARCH pozostaje bliska jedności. Wartość p dla statystyki Weighted Ljung-Box obliczonej dla kwadratów wystandaryzowanych reszt nie wykazuje istotności statystycznej, co oznacza brak autokorelacji w tych resztach. Dodatkowo, wyniki testu Weighted ARCH LM wskazują na brak efektów ARCH w resztach modelu. Jednakże nieistotny jest parametr alpha1, co nie potwierdza występowania efektu dźwigni. Dlatego lepszym rozwiązaniem jest pozostanie przy modelu GARCH(1,1)

plot(k.egarch11, which = 12)

plot(k.egarch11, which = 11)

plot(k.egarch11, which = 9)

Zwracając uwagę na wykres News Impact Curve, jej kształt jest dosyć nietypowy - pozytywne informacje wywołują większe reakcje zmienności niż negatywne. Odnośnie do QQ plot, reszty dalej nie wykazują rozkładu normalnego.

Model GARCH-t(1,1)

Ostatnim badanym modelm będzie GARCH-t(1,1), który zakłada, że reszty podlegają rozkładowi t-studenta, który lepiej modeluje leptokurtyczność. Biorąc pod uwagę poprzednie wyniki testów, wykresów QQ plot orz histogramów rozkładu stóp zwrotu, faktyczne odwołanie się do rozkładu t-studenta powinno wywołać pozytywne działanie na specyfikę modelu zmienności.

spec <- ugarchspec(
  # równanie warunkowej wariancji
  variance.model = list(model = "sGARCH",
                        garchOrder = c(1, 1)),
  # równanie warunkowej wartości oczekiwanej
  mean.model = list(armaOrder = c(2, 2),
                    include.mean = F),
  # zakladany rozklad warunkowy reszt
  distribution.model = "std") # rozklad t-Studenta

k.garcht11 <- ugarchfit(spec=spec, data=na.omit(portfolio_returns$r_portfolio))
k.garcht11

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## ar1    -0.394248    0.004629   -85.1766 0.000000
## ar2    -0.993218    0.001260  -788.0966 0.000000
## ma1     0.387032    0.001159   333.9806 0.000000
## ma2     0.997831    0.000066 15062.0991 0.000000
## omega   0.000035    0.000017     1.9866 0.046965
## alpha1  0.082087    0.019417     4.2275 0.000024
## beta1   0.911021    0.020387    44.6869 0.000000
## shape   4.644576    0.536003     8.6652 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## ar1    -0.394248    0.007053  -55.8986 0.000000
## ar2    -0.993218    0.002708 -366.7341 0.000000
## ma1     0.387032    0.001861  207.9803 0.000000
## ma2     0.997831    0.000209 4776.0815 0.000000
## omega   0.000035    0.000025    1.3962 0.162662
## alpha1  0.082087    0.026388    3.1107 0.001866
## beta1   0.911021    0.027713   32.8734 0.000000
## shape   4.644576    0.499925    9.2905 0.000000
## 
## LogLikelihood : 2890.379 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.3771
## Bayes        -3.3516
## Shibata      -3.3772
## Hannan-Quinn -3.3677
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.120 0.28982
## Lag[2*(p+q)+(p+q)-1][11]     7.402 0.01397
## Lag[4*(p+q)+(p+q)-1][19]    13.950 0.05992
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      0.637  0.4248
## Lag[2*(p+q)+(p+q)-1][5]     3.353  0.3463
## Lag[4*(p+q)+(p+q)-1][9]     6.933  0.2048
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.04863 0.500 2.000  0.8255
## ARCH Lag[5]   3.03472 1.440 1.667  0.2848
## ARCH Lag[7]   4.92103 2.315 1.543  0.2326
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.1381
## Individual Statistics:              
## ar1    0.05422
## ar2    0.24878
## ma1    0.08014
## ma2    0.18471
## omega  0.37036
## alpha1 0.21416
## beta1  0.25648
## shape  0.20811
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.15748 0.8749    
## Negative Sign Bias 0.03059 0.9756    
## Positive Sign Bias 0.34176 0.7326    
## Joint Effect       0.35072 0.9502    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     31.31      0.03729
## 2    30     40.96      0.06940
## 3    40     46.88      0.18058
## 4    50     44.46      0.65732
## 
## 
## Elapsed time : 0.333324

Wszystkie wspólczynniki AR, MA, alpha1, beta1 oraz shape są istotne statystycznie. Rozkład t-Studenta z shape = 4.64 dobrze odpowiada na leptokurtyczność danych. Wartość p dla statystyki Weighted Ljung-Box obliczonej dla kwadratów wystandaryzowanych reszt nie wykazuje istotności statystycznej, co oznacza brak autokorelacji w tych resztach. Dodatkowo, wyniki testu Weighted ARCH LM wskazują na brak efektów ARCH w resztach modelu.

plot(k.garcht11, which = 11)

plot(k.garcht11, which = 9)

plot(k.garcht11, which = 12)

Podsumowanie - wybór modelu

# Wyciągnięcie kryteriów informacyjnych z poszczególnych modeli
kryteria_info <- data.frame(
  infocriteria(k.garcht11),
  infocriteria(k.garch11_2),
  infocriteria(k.egarch11)
)

# Nadanie nazw kolumn
colnames(kryteria_info) <- c("GARCH-t(1,1)", "GARCH(1,1) ARMA(2,2)", "eGARCH(1,1)")

# Zaokrąglenie i estetyczna tabela
kable(round(kryteria_info, 4))

	GARCH-t(1,1)	GARCH(1,1) ARMA(2,2)	eGARCH(1,1)
Akaike	-3.3771	-3.2970	-3.2991
Bayes	-3.3516	-3.2747	-3.2736
Shibata	-3.3772	-3.2971	-3.2992
Hannan-Quinn	-3.3677	-3.2888	-3.2897

Podsumowując, model GARCH-t(1,1) wykazuje bardzo dobre dopasowanie do danych. Wszystkie kluczowe parametry są statystycznie istotne. Model dobrze odwzorowuje zmienność warunkową, a zastosowanie rozkładu t-Studenta skutecznie ujęło leptokurtyczność stóp zwrotu. Wykorzystanie rozkład t-Studenta wskazało, że model ten jest najelpszy spośród dotychczasowych pod względem wszystkich kryterium informacyjnych.

Prognozowanie

Wariancja

W kontekście modeli GARCH istotnym elementem opisu zmienności jest wariancja bezwarunkowa, oznaczający długookresowy poziom wariancji w czasie. W przypadku podstawowego modelu GARCH(1,1), wyraża się ją wzorem:

\[ E\varepsilon_t^2 = \frac{\alpha_0}{1 - \alpha_1 - \beta_1} \]

params <- coef(k.garch11_2)
var_uncond_2 <- params["omega"] / (1 - params["alpha1"] - params["beta1"])
names(var_uncond_2) <- "unconditional variance GARCH(1,1)"
var_uncond_2

## unconditional variance GARCH(1,1) 
##                       0.003588407

params <- coef(k.egarch11)
var_uncond_3 <- params["omega"] / (1 - params["alpha1"] - params["beta1"])
names(var_uncond_3) <- "unconditional variance EGARCH(1,1)"
var_uncond_3

## unconditional variance EGARCH(1,1) 
##                          -73.63111

params <- coef(k.garcht11)
var_uncond_1 <- params["omega"] / (1 - params["alpha1"] - params["beta1"])
names(var_uncond_1) <- "unconditional variance GARCH-t(1,1)"
var_uncond_1

## unconditional variance GARCH-t(1,1) 
##                         0.005007468

Następnie na podstawie dopasowanych wariantów modeli GARCH, możliwe jest prognozowanie wariancji warunkowej na kolejne okresy. Wykorzystuje się w tym celu funkcję ugarchforecast z paketu rugarch, która pozwala wygenerować prognozę zmienności na zadaną liczbę kroków w przód. Wynik następnie podnoszony jest do kwadratu w celu uzyskania prognozy.

fore500_2 <- ugarchforecast(k.garch11_2, n.ahead = 500)
sigma_forecast500_2 <- sigma(fore500_2)^2
cat("Prognozowana wariancja warunkowa – GARCH(1,1):\n")

## Prognozowana wariancja warunkowa – GARCH(1,1):

print(head(as.numeric(sigma_forecast500_2)))

## [1] 0.002025593 0.002045176 0.002064514 0.002083610 0.002102466 0.002121086

fore500_3 <- ugarchforecast(k.egarch11, n.ahead = 500)
sigma_forecast500_3 <- sigma(fore500_3)^2
cat("\nPrognozowana wariancja warunkowa – eGARCH(1,1):\n")

## 
## Prognozowana wariancja warunkowa – eGARCH(1,1):

print(head(as.numeric(sigma_forecast500_3)))

## [1] 0.002150132 0.002174300 0.002198123 0.002221599 0.002244728 0.002267508

fore500 <- ugarchforecast(k.garcht11, n.ahead = 500)
sigma_forecast500 <- sigma(fore500)^2
cat("\nPrognozowana wariancja warunkowa – GARCH-t(1,1):\n")

## 
## Prognozowana wariancja warunkowa – GARCH-t(1,1):

print(head(as.numeric(sigma_forecast500)))

## [1] 0.001932069 0.001953263 0.001974312 0.001995216 0.002015976 0.002036592

Uzyskane wektory reprezentują prognozowaną ścieżkę wariancji warunkowej na 500 przyszłych okresów.

# Wspólny wykres prognozowanej wariancji z 3 modeli
plot(sigma_forecast500_2, type = "l", col = "#1b9e77", lwd = 1.5,
     main = "Prognoza warunkowej wariancji z różnych modeli GARCH",
     ylab = "Warunkowa wariancja", xlab = "Dni naprzód",
     ylim = c(0.002, 0.0055))  # << rozszerzenie osi Y

lines(sigma_forecast500, col = "#d95f02", lwd = 1.5)        # GARCH-t(1,1)
lines(sigma_forecast500_3, col = "#7570b3", lwd = 1.5)      # eGARCH(1,1)

legend("topright",
       legend = c("GARCH(1,1)", "GARCH-t(1,1)", "eGARCH(1,1)"),
       col = c("#1b9e77", "#d95f02", "#7570b3"),
       lty = 1, lwd = 1.5)

Wykres przedstawia prognozę warunkowej wariancji dla trzech modeli GARCH. Wszystkie trzy modele wykazują wzrost wariancji w krótkim okresie, który następnie się stabilizuje. Najszybsza stabilizacja (około 30 dnia) zauważalna jest w przypadku modelu eGARCH, który uwzględnia asymetrię reakcji na dodatnie i ujemne szoki. Najpóźniej stabilizuje się model bazujący na rozkładzie t-Studenta, GARCH-t przewiduje wyraźnie wyższy poziom wariancji w długim okresie niż pozostałe modele.

In-sample

ggplot(portfolio_returns, aes(x = Date, y = r_portfolio)) +
  geom_line(linewidth = 0.3, color = "#411A1BFF") +
  labs(
    x = "Data",
    y = "Stopa zwrotu (%)",
    title = "Stopy zwrotu z całego portfela inwestycyjnego"
  ) +
  scale_x_date(date_minor_breaks = "1 year") +
  geom_vline(
    xintercept = as.Date("2024-01-01"),
    color = "#ff0000", linetype = "longdash",
    linewidth = 0.9
  )

Wyznaczenie wartości narażonej na ryzyko (Value at Risk, VaR) sprowadza się do oszacowania maksymalnej wielkości potencjalnych strat, jakie mogą wystąpić w portfelu inwestycyjnym przy zadanym poziomie istotności oraz w określonym horyzoncie czasowym. Analizę należy rozpocząć od wydzielenia przedziałów in-sample oraz out-of-sample. Zakres notowań portfela zaczyna się od 22-10-2020 r., a kończy na 24-06-2025 r. Na wykresie stóp zwrotu zauważalne są wyraźne epizody zwiększonej zmienności, szczególnie w pierwszej części próby (2020 - 2022). Później natomiast obserwowana jest stopniowa stabilizacja. W przypadku tej analizy, postanowiono, że próba out-of-sample zacznie się wraz 01-01-2024 r. (przerywana czerwona linia na wykresie). Przed tym okresem widoczne są dosyć duże wahania zmienności, których uwzględnienie w próbie in-sample pozwoli lepiej oszacować parametry modelu. Próba out-of-sample zawiera względną stabilizacje, jednakże nadal z występującymi wahaniami zmienności. Jest wystarczająco długa, dlatego umożliwia odpowiednią ocenę zdolności predykcyjnej modelu.

# wydzielenie probki out of sample
portfolio2 <- portfolio_returns[portfolio_returns$Date >= as.Date("2024-01-01"), ]
# probka in sample
portfolio1 <- portfolio_returns[portfolio_returns$Date < as.Date("2024-01-01"), ]

# standaryzacja zwrotow

portfolio1$rstd <- (portfolio1$r_portfolio - mean(portfolio1$r_portfolio, na.rm = TRUE)) /
  sd(portfolio1$r_portfolio, na.rm = TRUE)

# Obliczenie 1% kwantyla empirycznego
q01_empiryczny <- quantile(portfolio1$rstd, 0.01, na.rm = TRUE)
cat("Kwantyl empiryczny\n")

## Kwantyl empiryczny

as.numeric(q01_empiryczny)

## [1] -2.570963

# Obliczenie 1% kwantyla standardowego rozkładu normalnego
q01_normalny <- qnorm(0.01, mean = 0, sd = 1)
cat("Kwantyl teoretyczny\n")

## Kwantyl teoretyczny

as.numeric(q01_normalny)

## [1] -2.326348

Dla próby in-sample dokonano standaryzacji logarytmicznych stóp zwrotu portfela inwestycyjnego, tj. przekształcenia ich do postaci o średniej równej zero i odchyleniu standardowym równym jeden. Umożliwiło to porównanie empirycznego rozkładu zwrotów z rozkładem normalnym – standardowym założeniem w klasycznych modelach VaR.

Nastepnie wyznaczono empiryczny 1-procentowy kwantyl rozkładu standaryzowanych stóp zwrotu portfela, który wyniósł -2.571. Oznacza to, że w okresie estymacyjnym 1% najgorszych dziennych zwrotów było mniejszych niż -2.571 odchyleń standardowych od średniej.

Dla porównania, przy założeniu normalności rozkładu, teoretyczny 1-procentowy kwantyl wynosi -2.326. Otrzymana wartość jest znacznie bardziej ujemna niż teoretyczna, co wskazuje na leptokurtyczność rozkładu zwrotów portfela. Oznacza to, że ekstremalne spadki (czyli rzadkie, bardzo niekorzystne zdarzenia) występowały częściej, niż sugerowałby klasyczny model oparty na rozkładzie normalnym.

# GARCH(1,1) z rozkładem normalnym
spec_garchn <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
  distribution.model = "norm"
)

in.garch11 <- ugarchfit(spec = spec_garchn, data = na.omit(portfolio1$r_portfolio))

# eGARCH(1,1)
spec_egarch <- ugarchspec(
  variance.model = list(model = "eGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
  distribution.model = "norm"
)

in.egarch11 <- ugarchfit(spec = spec_egarch, data = na.omit(portfolio1$r_portfolio))


# GARCH-t(1,1)
spec_garcht <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
  distribution.model = "std"
)

in.garcht11 <- ugarchfit(spec = spec_garcht, data = na.omit(portfolio1$r_portfolio))

Oszacowano modele GARCH(1,1), eGARCH(1,1) oraz GARCH-t(1,1) na danych w okresie in-sample, a następnie wyznaczono wartość narażoną na ryzyko (VAR) jako iloczyn prognozowanej zmienności portfela oraz 1-procentowego kwantyla empirycznego rozkładu standaryzowanych stóp zwrotu.

portfolio1$VaR_garch11  <- q01_empiryczny * in.garch11@fit$sigma
portfolio1$VaR_egarch11 <- q01_empiryczny * in.egarch11@fit$sigma
portfolio1$VaR_garcht11 <- q01_empiryczny * in.garcht11@fit$sigma

plot(portfolio1$Date, portfolio1$r_portfolio, type = "l", col = "red", lwd = 1,
     ylim = c(-0.5, 0.5), xlab = "Data", ylab = "Zwrot / VaR",
     main = "Zwroty portfela vs. VaR z trzech modeli")

abline(h = 0, lty = 2)

# Dodajemy linie VaR
lines(portfolio1$Date, portfolio1$VaR_garch11,   col = "#1b9e77", lwd = 1.2)
lines(portfolio1$Date, portfolio1$VaR_egarch11,  col = "#d95f02", lwd = 1.2)
lines(portfolio1$Date, portfolio1$VaR_garcht11,  col = "#7570b3", lwd = 1.2)

# Legenda
legend("bottomleft",
       legend = c("Zwrot portfela", "VaR GARCH(1,1)", "VaR eGARCH(1,1)", "VaR GARCH-t(1,1)"),
       col = c("red", "#1b9e77", "#d95f02", "#7570b3"),
       lwd = 2,
       cex = 0.8,
       bty = "n")

Powyższy wykres przedstawia oszacowane wartości narażone na ryzyko wraz z szeregiem stóp zwrotu dla każdego z wybranych wcześniej modeli. Można zaobserwować, że wszystkie trzy modele dosyć dobrze oddają strukturę ryzyka w czasie. Linie VAR przebiegają poniżej zwrotów, wyznaczając granicę, której przekroczenie oznacza nieoczekiwane, wysokie straty.

n_g11 <- sum(portfolio1$r_portfolio < portfolio1$VaR_garch11, na.rm = TRUE)
p_g11 <- n_g11 / length(na.omit(portfolio1$VaR_garch11))

n_eg11 <- sum(portfolio1$r_portfolio < portfolio1$VaR_egarch11, na.rm = TRUE)
p_eg11 <- n_eg11 / length(na.omit(portfolio1$VaR_egarch11))

n_gt11 <- sum(portfolio1$r_portfolio < portfolio1$VaR_garcht11, na.rm = TRUE)
p_gt11 <- n_gt11 / length(na.omit(portfolio1$VaR_garcht11))

# Wyniki
cat("GARCH:", n_g11, "(", round(100 * p_g11, 2), "%)\n")

## GARCH: 19 ( 1.63 %)

cat("eGARCH:", n_eg11, "(", round(100 * p_eg11, 2), "%)\n")

## eGARCH: 17 ( 1.46 %)

cat("GARCH-t:", n_gt11, "(", round(100 * p_gt11, 2), "%)\n")

## GARCH-t: 16 ( 1.37 %)

Z punktu widzenia liczby przekroczeń - sytuacji, w których rzeczywista strata okazała się większa niż VAR, najwięcej wystąpiło w modelu GARCH(1,1), było ich aż 19, co stanowi 1,63$ przypadków. Najmniej przekroczeń uzyskano w modelu GARCH-t - 16, co stanowiło 1,37% przypadków.

Out-of-sample

#GARCH(1,1)
# Zakres próby out-of-sample
start  <- which(portfolio_returns$Date == as.Date("2024-01-01"))
finish <- nrow(portfolio_returns)
portfolio2 <- portfolio_returns[start:finish, ]

VaR_garch <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]

  # Standaryzacja
  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)

  # Specyfikacja modelu
  spec <- ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "norm"
  )

  tmp.garch11 <- ugarchfit(
    spec = spec,
    data = na.omit(tmp.data$r_portfolio),
    solver = "hybrid",
    solver.control = list(trace = 0)
  )

  sigma.forecast <- ugarchforecast(tmp.garch11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]

  VaR_garch[k - start + 1] <- q01 * sigma.forecast2
}

# Przypisanie wyników
portfolio2$VaR_garch <- VaR_garch

#EGARCH(1,1)

VaR_egarch <- rep(NA, times = finish - start + 1)


for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]

  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)

  spec <- ugarchspec(
    variance.model = list(model = "eGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "norm"
  )

  tmp.egarch11 <- ugarchfit(
    spec = spec,
    data = na.omit(tmp.data$r_portfolio),
    solver = "hybrid",
    solver.control = list(trace = 0)
  )

  sigma.forecast <- ugarchforecast(tmp.egarch11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR_egarch[k - start + 1] <- q01 * sigma.forecast2
}

portfolio2$VaR_egarch <- VaR_egarch

#GARCH-t
VaR_garcht <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]

  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)

  spec <- ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "std"
  )

  tmp.garcht11 <- ugarchfit(
    spec = spec,
    data = na.omit(tmp.data$r_portfolio),
    solver = "hybrid",
    solver.control = list(trace = 0)
  )

  sigma.forecast <- ugarchforecast(tmp.garcht11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR_garcht[k - start + 1] <- q01 * sigma.forecast2
}

portfolio2$VaR_garcht <- VaR_garcht

Oszacowano modele GARCH(1,1), eGARCH(1,1) oraz GARCH-t(1,1) na danych w okresie out-of-sample, a następnie przypisano wartość narażoną na ryzyko (VAR) do odpowiedniego zbioru.

n_g11 <- sum(portfolio2$r_portfolio < portfolio2$VaR_garch, na.rm = TRUE)
p_g11 <- n_g11 / sum(!is.na(portfolio2$VaR_garch))
cat("GARCH(1,1):", n_g11, "przypadków (", round(100 * p_g11, 2), "%)\n")

## GARCH(1,1): 5 przypadków ( 0.92 %)

n_eg11 <- sum(portfolio2$r_portfolio < portfolio2$VaR_egarch, na.rm = TRUE)
p_eg11 <- n_eg11 / sum(!is.na(portfolio2$VaR_egarch))
cat("eGARCH(1,1):", n_eg11, "przypadków (", round(100 * p_eg11, 2), "%)\n")

## eGARCH(1,1): 5 przypadków ( 0.92 %)

n_gt11 <- sum(portfolio2$r_portfolio < portfolio2$VaR_garcht, na.rm = TRUE)
p_gt11 <- n_gt11 / sum(!is.na(portfolio2$VaR_garcht))
cat("GARCH-t(1,1):", n_gt11, "przypadków (", round(100 * p_gt11, 2), "%)\n")

## GARCH-t(1,1): 5 przypadków ( 0.92 %)

Następnie wyznaczono liczbę przekroczeń - sytuacji, w których rzeczywista strata okazała się większa niż VAR.

# Wykres prognoz VaR i rzeczywistych zwrotów
plot(portfolio2$Date, portfolio2$r_portfolio, col = "black", lwd = 1, type = 'l',
     ylim = c(-0.2, 0.2), main = "Zwroty i prognozy VaR (1%) – GARCH(1,1), eGARCH(1,1), GARCH-t(1,1)",
     ylab = "Zwrot / VaR", xlab = "Data")

# Linie dla prognoz VaR z różnych modeli
lines(portfolio2$Date, portfolio2$VaR_garch, col = "darkgreen", lwd = 1.2)
lines(portfolio2$Date, portfolio2$VaR_egarch, col = "red", lwd = 1.2)
lines(portfolio2$Date, portfolio2$VaR_garcht, col = "purple", lwd = 1.2)

# Legenda
legend("bottomleft",
       legend = c("Zwrot portfela", "VaR – GARCH(1,1)", "VaR – eGARCH(1,1)", "VaR – GARCH-t(1,1)"),
       col = c("black", "darkgreen", "red", "purple"),
       lty = 1, lwd = 1.5, cex = 0.8, bty = "n")

Wszystkie trzy modele prognozują bardzo zbliżone poziomy VAR, linie często się pokrywają. W okresach zwiększonej zmienności obserwowane jest wyższe pogłębienie w przypadku GARCH(1,1). Wskaźnik liczbt przekroczeń wyniósł dokładnie 5 przypadków dla każdego z modeli, co stanowi 0,92% wszystkich obserwacji. Wynik jest zbliżony do przyjętego poziomu istotności (1%).

Analiza wrażliwości

Na etapie końcowym zdecydowano się przeprowadzić analizę wrażliwości, której celem było ukazanie wpływu przyjętych założeń początkowych. W szczególności długości okresu prognozy na estymację parametrów modelu oraz końcowe wyniki analizy ryzyka. W tym kontekście, wykorzystano wyłoniony uprzednio na podstawie wartości kryteriów informacyjnych – model GARCH-t(1,1). Został on ponownie oszacowany dla czterech wariantów długości próby out-of-sample: 9 miesięcy, 6 miesięcy, 2 miesiące oraz 2 tygodnie. 

# 9 miesięcy
start_date <- as.Date("2024-09-24")
end_date   <- as.Date("2025-06-24")

start  <- which(portfolio_returns$Date == start_date)
finish <- which(portfolio_returns$Date == end_date)

portfolio_9m <- portfolio_returns[start:finish, ]
VaR_garcht_9m <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]
  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)
  spec <- ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "std"
  )
  tmp.garcht11 <- ugarchfit(spec = spec, data = na.omit(tmp.data$r_portfolio), solver = "hybrid", solver.control = list(trace = 0))
  sigma.forecast <- ugarchforecast(tmp.garcht11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR_garcht_9m[k - start + 1] <- q01 * sigma.forecast2
}
portfolio_9m$VaR_garcht <- VaR_garcht_9m

# 6 miesięcy
start_date <- as.Date("2024-12-24")
end_date   <- as.Date("2025-06-24")

start  <- which(portfolio_returns$Date == start_date)
finish <- which(portfolio_returns$Date == end_date)

portfolio_6m <- portfolio_returns[start:finish, ]
VaR_garcht_6m <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]
  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)
  spec <- ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "std"
  )
  tmp.garcht11 <- ugarchfit(spec = spec, data = na.omit(tmp.data$r_portfolio), solver = "hybrid", solver.control = list(trace = 0))
  sigma.forecast <- ugarchforecast(tmp.garcht11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR_garcht_6m[k - start + 1] <- q01 * sigma.forecast2
}
portfolio_6m$VaR_garcht <- VaR_garcht_6m

# 2 miesiące
start_date <- as.Date("2025-04-24")
end_date   <- as.Date("2025-06-24")

start  <- which(portfolio_returns$Date == start_date)
finish <- which(portfolio_returns$Date == end_date)

portfolio_2m <- portfolio_returns[start:finish, ]
VaR_garcht_2m <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]
  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)
  spec <- ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "std"
  )
  tmp.garcht11 <- ugarchfit(spec = spec, data = na.omit(tmp.data$r_portfolio), solver = "hybrid", solver.control = list(trace = 0))
  sigma.forecast <- ugarchforecast(tmp.garcht11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR_garcht_2m[k - start + 1] <- q01 * sigma.forecast2
}
portfolio_2m$VaR_garcht <- VaR_garcht_2m

# 2 tygodnie
start_date <- as.Date("2025-06-10")
end_date   <- as.Date("2025-06-24")

start  <- which(portfolio_returns$Date == start_date)
finish <- which(portfolio_returns$Date == end_date)

portfolio_2w <- portfolio_returns[start:finish, ]
VaR_garcht_2w <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- portfolio_returns[1:(k - 1), ]
  tmp.data$rstd <- (tmp.data$r_portfolio - mean(tmp.data$r_portfolio, na.rm = TRUE)) /
    sd(tmp.data$r_portfolio, na.rm = TRUE)
  q01 <- quantile(tmp.data$rstd, 0.01, na.rm = TRUE)
  spec <- ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
    mean.model = list(armaOrder = c(2, 2), include.mean = FALSE),
    distribution.model = "std"
  )
  tmp.garcht11 <- ugarchfit(spec = spec, data = na.omit(tmp.data$r_portfolio), solver = "hybrid", solver.control = list(trace = 0))
  sigma.forecast <- ugarchforecast(tmp.garcht11, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR_garcht_2w[k - start + 1] <- q01 * sigma.forecast2
}
portfolio_2w$VaR_garcht <- VaR_garcht_2w

# Wykres
plot(portfolio_9m$Date, portfolio_9m$r_portfolio, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.2, 0.2), main = "VaR – GARCH-t(1,1), 9 miesięcy",
     ylab = "Zwrot", xlab = "Data")
abline(h = 0, lty = 2)
lines(portfolio_9m$Date, portfolio_9m$VaR_garcht, type = 'l', col = "purple", lwd = 1.2)

Pierwszy wykres przedstawia okno długości 9 msc

# Wykres
plot(portfolio_6m$Date, portfolio_6m$r_portfolio, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.2, 0.2), main = "VaR – GARCH-t(1,1), 6 miesięcy",
     ylab = "Zwrot", xlab = "Data")
abline(h = 0, lty = 2)
lines(portfolio_6m$Date, portfolio_6m$VaR_garcht, type = 'l', col = "purple", lwd = 1.2)

Następnie próbę skrócono do 6 msc.

plot(portfolio_2m$Date, portfolio_2m$r_portfolio, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.2, 0.2), main = "VaR – GARCH-t(1,1), 2 miesiące",
     ylab = "Zwrot", xlab = "Data")
abline(h = 0, lty = 2)
lines(portfolio_2m$Date, portfolio_2m$VaR_garcht, type = 'l', col = "purple", lwd = 1.2)

Kolejny wykres prezentuje okno czasowe dot. 2 msc.

plot(portfolio_2w$Date, portfolio_2w$r_portfolio, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.2, 0.2), main = "VaR – GARCH-t(1,1), 2 tygodnie",
     ylab = "Zwrot", xlab = "Data")
abline(h = 0, lty = 2)
lines(portfolio_2w$Date, portfolio_2w$VaR_garcht, type = 'l', col = "purple", lwd = 1.2)

Ostatni wykres przedstawia okno czasowe 2 tyg.

# 9 miesięcy
n_exceed_9m <- sum(portfolio_9m$r_portfolio < portfolio_9m$VaR_garcht, na.rm = TRUE)
p_exceed_9m <- n_exceed_9m / sum(!is.na(portfolio_9m$VaR_garcht))
cat("Przekroczenia VaR GARCH-t(1,1) – 9 miesięcy:", n_exceed_9m, "(", round(100 * p_exceed_9m, 2), "%)\n")

## Przekroczenia VaR GARCH-t(1,1) – 9 miesięcy: 4 ( 1.46 %)

# 6 miesięcy
n_exceed_6m <- sum(portfolio_6m$r_portfolio < portfolio_6m$VaR_garcht, na.rm = TRUE)
p_exceed_6m <- n_exceed_6m / sum(!is.na(portfolio_6m$VaR_garcht))
cat("Przekroczenia VaR GARCH-t(1,1) – 6 miesięcy:", n_exceed_6m, "(", round(100 * p_exceed_6m, 2), "%)\n")

## Przekroczenia VaR GARCH-t(1,1) – 6 miesięcy: 3 ( 1.64 %)

# 2 miesiące
n_exceed_2m <- sum(portfolio_2m$r_portfolio < portfolio_2m$VaR_garcht, na.rm = TRUE)
p_exceed_2m <- n_exceed_2m / sum(!is.na(portfolio_2m$VaR_garcht))
cat("Przekroczenia VaR GARCH-t(1,1) – 2 miesiące:", n_exceed_2m, "(", round(100 * p_exceed_2m, 2), "%)\n")

## Przekroczenia VaR GARCH-t(1,1) – 2 miesiące: 0 ( 0 %)

# 2 tygodnie
n_exceed_2w <- sum(portfolio_2w$r_portfolio < portfolio_2w$VaR_garcht, na.rm = TRUE)
p_exceed_2w <- n_exceed_2w / sum(!is.na(portfolio_2w$VaR_garcht))
cat("Przekroczenia VaR GARCH-t(1,1) – 2 tygodnie:", n_exceed_2w, "(", round(100 * p_exceed_2w, 2), "%)\n")

## Przekroczenia VaR GARCH-t(1,1) – 2 tygodnie: 0 ( 0 %)

Wydłużając okres prognozowania uzyskiwana jest większa liczba przekroczeń VAR. Odsetki wraz z rozwojem długości próby rosną. Jednakże w porównaniu do całego poprzedniego okresu out-of-sample, który rozpoczął się w 01-01-2024 r. (czyli 18 msc.) wyniki przekroczeń są gorsze, wtedy odsetek ich wynosił (0.92%). Pokazuje to, że odpowiednie manipulowanie próbą może znacząco zmieniać wyniki, w przypadku bardzo krótkich okresów jak 2 msc, czy 2 tyg, odsetek stanowi 0%. Okres 9 msc, patrząc na wykres poniżej obejmował największe wahania zmienności, podobnie jak z okresem 6 msc. Skuteczność prognoz VAR zależy zatem w dużej mierze od zmienności i struktury danych w danym okresie.

plot(portfolio_9m$Date, portfolio_9m$r_portfolio, col = "black", lwd = 1, type = 'l',
     ylim = c(-0.2, 0.2), main = "Zwroty i prognozy VaR – GARCH-t(1,1)",
     ylab = "Zwrot / VaR", xlab = "Data")
abline(h = 0, lty = 2)

# Linie prognoz z różnych horyzontów
lines(portfolio_9m$Date, portfolio_9m$VaR_garcht, col = "blue", lwd = 1.5)
lines(portfolio_6m$Date, portfolio_6m$VaR_garcht, col = "darkgreen", lwd = 1.5)
lines(portfolio_2m$Date, portfolio_2m$VaR_garcht, col = "orange", lwd = 1.5)
lines(portfolio_2w$Date, portfolio_2w$VaR_garcht, col = "purple", lwd = 1.5)

# Dodanie linii z pełnego okresu
lines(portfolio2$Date, portfolio2$VaR_garcht, col = "red", lwd = 1.5, lty = 2)

legend("bottomleft",
       legend = c("Zwrot portfela", 
                  "VaR – 9 miesięcy", 
                  "VaR – 6 miesięcy", 
                  "VaR – 2 miesiące", 
                  "VaR – 2 tygodnie",
                  "VaR – cały okres"),
       col = c("black", "blue", "darkgreen", "orange", "purple", "red"),
       lty = c(1, 1, 1, 1, 1, 2), lwd = 1.5, cex = 0.8, bty = "n")

Podsumowanie

W ramach przeprowadzonej analizy dokonano oszacowania ryzyka portfela inwestycyjnego składającego się z 4 kryptowalut (Algorand, Avalanche, PancakeSwap oraz PolkaDot) w okresie od października 2020r. do czerwaca 2025 r. Szczególnie uwzględniono modelowanie zmienności oraz wartości narażonej na ryzyko (VAR).

Przeprowadzono weryfikację istotnych efektów zmienności, w tym efektów ARCH i struktury autokorelacji. Następnie analiza objęła modelowanie przy pomocy klasycznego modelu GARCH oraz jego rozszerzeń eGARCH i GARCH-t. Następnie wykonano prognozowanie warunkowej wariancji na okresie 18 msc, po czym sprawdzono wrażliwość próbki odwołując się do okien czasowych 9 msc, 6 msc, 2 msc oraz 2 tyg.

Model GARCH-t(1,1) okazał się najlepiej dopasowanym do empirycznej struktury danych, dobrze uwzględniając leptokurtyczność analizowanych stóp zwrotu. Jednakże w samym prognozowaniu, modele GARCH, eGARCH oraz GARCH-t prezentowały podobną jakość, liczba przekroczeń VAR oraz jego odsetek był na identycznym poziomie.