Krótka analiza i prognoza premii za ryzyko S&P 500 (język R)

Powiedzmy, że chcemy oszacować premię za ryzyko amerykańskiej giełdy w roku 2025 i 2026. Określenie "premia za ryzyko" używam jako skrót na oznaczenie nadwyżki stopy zwrotu z indeksu (tu S&P 500) nad stopę zwrotu z 10-letnich obligacji skarbowych. Precyzyjnie przecież mówiąc o premii mamy na myśli oczekiwane lub wymagane wynagrodzenie za przyszłe ryzyko, czyli takie, którego okres jeszcze nie zaczął się. I może się zrealizować (wtedy albo stracimy, albo mniej zarobimy), albo nie zrealizować (zarobimy co najmniej tyle, ile wymagamy). Dodatnia premia za ryzyko wynika z tego, że inwestorzy są niechętni wobec ryzyka. To znaczy, że wyceniają niżej aktywa o wyższym ryzyku. Mechanizm ten pozwala łatwo zrozumieć, dlaczego oczekiwana (a więc i średnia) stopa zwrotu jest dodatnia, a nie zerowa, jak w hazardzie. Otóż cała tajemnica leży w tym, że to ludzie wyceniają aktywa, a nie natura. Jeżeli wielu inwestorów będzie wyceniać niżej akcje z powodu ich większego ryzyka, to przecież spółka nie przestaje zarabiać tyle co zarabiała, jej wartość księgowa się nie zmienia. Spadek ceny oznacza więc, że inwestorzy mogą więcej zarobić, sprzedając w przyszłości akcje. Ryzyko jest mimo wszystko miarą psychologiczną, a więc dość subiektywną - i to właśnie sprawia, że inwestowanie daje dużo większe możliwości niż hazard.

Historyczne roczne premie za ryzyko S&P 500 dla lat 1945-2024 biorę ze strony Damodarana. Średnia wyniosła 7,7%. To coś dużo biorąc pod uwagę, że średnia stopa wolna od ryzyka to ok. 5%. Porównałem więc te dane z S&P U.S. Equity Risk Premium Index. Obejmuje on zajęcie długiej pozycji w indeksie S&P 500 Futures Excess Return Index i krótkiej pozycji w indeksie S&P U.S. Treasury Bond Futures Excess Return Index [zob. S&P Factor Indices Methodology]. Od 2015 roku dostałem taki wykres:

A więc dane są poprawne, bo idealnie korelują. Oczywiście w tym okresie średnia była jeszcze wyższa i dla S&P U.S. Equity Risk Premium Index wyniosła 13,2%, a Damodarana 13,5%. 

Stabilność średniej

Pytanie brzmi czy średnia się zmieniała w ciągu tych 80 lat. Przez pierwsze 40 lat wyniosła ok. 8,6%, a kolejne niecałe 6,8%. Sprawdźmy czy możemy uznać to za coś więcej niż przypadek. 

Jeżeli dzielimy okres na dwie części, to najlepiej użyć testu Chowa. W tym artykule pokazałem, jak go wykonać w R.

if (require("strucchange")==FALSE) {

  install.packages("strucchange")

  library("strucchange")

}

sctest(premia ~ 1, type="Chow", point=40)

	Chow test

data:  premia ~ 1
F = 0.18, p-value = 0.7

Jak widać nie ma podstaw do odrzucenia hipotezy, że średnia jest stała. Wykonałem dodatkowo najróżniejsze testy (ADF, KPSS, Bai-Perrona, CUSUM, MOSUM) i nie wykryłem niestacjonarności czy niestabilności. 

Modelowanie

Jeżeli średnia się nie zmienia, to trzeba sprawdzić autokorelację. Najpierw ACF:

acf(premia, 15)

Do dziś nie mogę z tego, że R pokazuje bez sensu ten zerowy rząd. Głupota. W cząstkowej jest już normalnie od pierwszego rzędu: 

pacf(premia, 15)

W Dodatku do tego artykułu zrobiłem analogiczną analizę dla nadwyżki nad bonami. Tam PACF 4-tego rzędu jest na granicy istotności. Tutaj natomiast nie można mówić o jakiejkolwiek istotności. Z tego z kolei wynika, że obligacje skarbowe mają swój udział w cyklu 4-letnim (bo je tam dodaliśmy). 

Mimo wszystko sprawdzę model AR(4) i jego stabilność. W R kod będzie następujący:

  # Pierwszy model - arima

  premia_model_arma1 <- arima(

    premia,

    order           = c(4, 0, 0),          # ARMA(4, 0, 0)

    include.mean    = TRUE,                # estymuj μ

    fixed           = c(0, 0, 0, NA, NA),  # φ1=0, φ2=0, φ3=0, φ4=estymuj, μ=estymuj

    transform.pars  = FALSE                # aby fixed=0 było literalne

  )

  summary(premia_model_arma1)

Coefficients:
      ar1  ar2  ar3    ar4  intercept
        0    0    0  0.205      0.079
s.e.    0    0    0  0.111      0.026

sigma^2 estimated as 0.0344:  log likelihood = 21.16,  aic = -36.33

Gdyby zastosować taki model, dostalibyśmy rzeczywiście średnią 7,9% (+/- 2,6%), ale dodatkowo współczynnik 4-tego rzędu 0,205 (+/- 0,11).

Aby sprawdzić stabilność tego modelu, trzeba niestety zamienić na formę typu:

premia4 <- tail(zmienna, -4)

premia_4 <- head(zmienna, -4)

premia.model <- premia4 ~ premia_4

i przykładowo, żeby użyć testu CUSUM-Score, wpisać:

par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))

cus = efp(premia.model, type="Score-CUSUM")

plot(cus, functional=NULL)

Za jednym zamachem mamy tu sprawdzenie stabilności wszystkich współczynników, w tym wariancji. Brak jakichkolwiek wychyleń (byłyby zaznaczenia) świadczy o stabilności. 

To jednak tylko punkt wyjścia, bo znacznie lepszym sposobem jest użycie autoarfima() w rugarch lub auto.arima() w forecast. Wybrałem ten pierwszy, bo dostarcza wszystkie potrzebne statystyki. Ten argument  omówiłem w tym artykule. Sprawdziłem więc w autoarfima max ARMA(7, 7) w metodzie full. Wybrałem rozkład "std", bo dawał najlepsze rezultaty, jeśli chodzi o statystyki dopasowania. Wadą rugarch jest czasochłonność obliczeń, dlatego koniecznie trzeba używać argumentu cluster tam gdzie się da. Ponadto wybrałem solver "nlminb", ponieważ wydaje się dość efektywny. Użyłem funkcji dwukrotnie - wg kryterium BIC i AIC (HQIC nie używam, bo danych jest za mało, zob. tu). Okazało się, że oba dały ten sam model. Zwracam też uwagę, że nie stosuję tutaj modelu GARCH ze względu na małą liczbę danych. Efekt:

# Drugi model - autoarfima (rugarch)

    if (require("rugarch")==FALSE) {

      install.packages("rugarch")

    }

    library("rugarch")  

# Kryterium BIC

  klaster <- makeCluster(detectCores() - 1)

  premia_model_arma_bic <- autoarfima(data = premia, ar.max = 7, ma.max = 7, criterion = "BIC",

             method = "full", arfima = FALSE, include.mean = NULL,

             distribution.model = "std", cluster = klaster, external.regressors = NULL,

             solver = "nlminb")

  stopCluster(klaster)

  premia_model_arma_bic

  # Kryterium AIC

  klaster <- makeCluster(detectCores() - 1)

  premia_model_arma_aic <- autoarfima(data = premia, ar.max = 7, ma.max = 7, criterion = "AIC",

                                      method = "full", arfima = FALSE, include.mean = NULL,

                                      distribution.model = "std", cluster = klaster, external.regressors = NULL,

                                      solver = "nlminb")

  stopCluster(klaster)

  premia_model_arma_aic

  premia_model_arma2 <- premia_model_arma_aic

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(7,0,4)
Distribution	: std 

Optimal Parameters
------------------------------------
       Estimate  Std. Error    t value Pr(>|t|)
ar1     0.44901    0.000926    484.974        0
ar2     0.00000          NA         NA       NA
ar3     0.00000          NA         NA       NA
ar4     0.00000          NA         NA       NA
ar5    -0.22529    0.000465   -484.974        0
ar6     0.00000          NA         NA       NA
ar7     0.34287    0.000027  12785.413        0
ma1    -1.28980    0.000061 -21191.451        0
ma2     0.00000          NA         NA       NA
ma3     0.61606    0.000029  21191.451        0
ma4     0.44741    0.000021  21191.451        0
sigma   0.13607    0.005403     25.184        0
shape 100.00000    0.004585  21809.383        0

Robust Standard Errors:
       Estimate  Std. Error    t value Pr(>|t|)
ar1     0.44901    2.324983    0.19312  0.84686
ar2     0.00000          NA         NA       NA
ar3     0.00000          NA         NA       NA
ar4     0.00000          NA         NA       NA
ar5    -0.22529    1.166585   -0.19312  0.84686
ar6     0.00000          NA         NA       NA
ar7     0.34287    0.000385  891.10681  0.00000
ma1    -1.28980    0.005203 -247.88844  0.00000
ma2     0.00000          NA         NA       NA
ma3     0.61606    0.000766  804.54128  0.00000
ma4     0.44741    0.001619  276.38488  0.00000
sigma   0.13607    0.619235    0.21974  0.82608
shape 100.00000    1.325372   75.45050  0.00000

LogLikelihood : 46.12 

Information Criteria
------------------------------------
                     
Akaike       -0.95299
Bayes        -0.71479
Shibata      -0.97067
Hannan-Quinn -0.85749

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                    0.000611  0.9803
Lag[2*(p+q)+(p+q)-1][32] 11.645908  1.0000
Lag[4*(p+q)+(p+q)-1][54] 17.816677  0.9976

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4785  0.4891
Lag[2*(p+q)+(p+q)-1][2]    0.7811  0.5750
Lag[4*(p+q)+(p+q)-1][5]    2.3214  0.5447


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      1.519   2  0.4678
ARCH Lag[5]      2.647   5  0.7542
ARCH Lag[10]     9.429  10  0.4920

Nyblom stability test
------------------------------------
Joint Statistic:  231429
Individual Statistics:             
ar1   0.09969
ar5   0.05997
ar7   0.01355
ma1   0.03352
ma3   0.11951
ma4   0.09983
sigma 0.13154
shape 0.42456

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.89 2.11 2.59
Individual Statistic:	 0.35 0.47 0.75

BIC i AIC wskazały ten sam model: AR o rzędach 1, 5 i 7; MA o rzędach 1, 3 i 4. W nieodpornych statystykach (tj. przy założeniu danego rozkładu, tu t-Studenta) wszystkie są istotne. Jednak statystyki odporne (na odchylenia od założonego rozkładu reszt) wskazują jedynie istotność rzędu 7 dla AR, a dla MA wszystkie są istotne. Pozostałe statystyki są okej, wszystkie współczynniki są stabilne, chociaż znowu cały model mocno niestabilny (pewnie niezbyt idealny rozkład).

Test MCS

Aby się upewnić, że na pewno drugi model jest lepszy od pierwszego, użyjemy testu MCS, który omówiłem w poprzednim artykule. Tak jak wtedy za stratę uznam sytuację, gdy znak prognozy in-sample nie zgadza się ze znakiem rzeczywistej wartości.

  #MCS test

  Funkcja_strat <- function(rzeczywiste, prognoza) {

    roznica <- rzeczywiste - prognoza

    strata <- ifelse(sign(rzeczywiste) != sign(prognoza), abs(roznica), 0)

    return(strata)

  }

  # Obliczenie funkcji strat dla każdego modelu (in-sample)

  straty1 <- Funkcja_strat(premia, as.numeric(fitted(premia_model_arma1)))

  straty2 <- Funkcja_strat(premia, premia_model_arma2$fit@fit$fitted.values) 

  straty <- cbind(straty1, straty2)

  mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 500, nblock = 10, boot = "stationary")

  print(mcs_stationary)

$includedR
[1] 1 2

$pvalsR
     [,1]
[1,] 0.23
[2,] 1.00

$excludedR
NULL

$includedSQ
[1] 1 2

$pvalsSQ
     [,1]
[1,] 0.23
[2,] 1.00

$excludedSQ
NULL

Test MCS nie pozwolił wykluczyć żadnego z dwóch modeli, ale wskazuje dużo bardziej na zachowanie drugiego, tj. z autoarfima, więc go wybierzemy.

Wykres modelu

Wobec tego sprawdźmy jak te dopasowania dla najlepszego modelu wyglądają:

  # Wartości dopasowane

  dopas <- premia_model_arma2$fit@fit$fitted.values

  okres <- c(1945:2024)

  # Wykres najlepszego modelu z prawdziwymi obserwacjami

  plot(x=okres, y=premia, type="l", ylab = "")

  lines(x=okres, y=dopas, col="red")

  abline(h=0, col="grey")

  abline(h=mean(premia), col="blue", lwd=2)

  legend("topright", legend=c("premia za ryzyko S&P 500", "prognoza in-sample", "średnia premia za ryzyko"), 

         col=c("black","red", "blue"), lwd=2, cex=0.6)

Prognoza na 2025-2026

Niestety nie mam wiedzy, czy można od razu w rugarch robić prognozy mając już najlepszy model z autoarfima. jestem niemal pewny, że jest sposób, ale szybciej jest mi tak zrobić. Dlatego robię to "na piechotę". Obliczam arfimaspec(), z tego arfimafit i dopiero z tego arfimaforecast(). W arfimaspec() wpisujemy normalnie maksymalne rzędy ARMA, ale za to dodajemy argument fixed.pars = list(ar2=0, ar3=0), który właśnie kontroluje ustawienie konkretnych współczynników. 

  rzadAR <- premia_model_arma2$fit@model$modelinc[2]

  rzadMA <- premia_model_arma2$fit@model$modelinc[3]

premia_spec <- arfimaspec(mean.model = list(armaOrder = c(rzadAR, rzadMA), include.mean = TRUE,

                                            arfima = FALSE, external.regressors = NULL), distribution.model = "std",

                          start.pars = list() 

                          , fixed.pars = list(ar2=0, 

                                                                 ar3=0,

                                                                 ar4=0,

                                                                 ar6=0,

                                                                 ma2=0)

                          )

premia_fit <- arfimafit(spec = premia_spec, data = premia, solver = "nlminb")

prognoza <- arfimaforecast(premia_fit, n.ahead = 2)

prognoza

*----------------------------------*
*        ARFIMA Model Forecast     *
*----------------------------------*

Horizon: 2
Roll Steps: 0
Out of Sample: 0

0-roll forecast: 
    T+1     T+2 
0.05436 0.07451 

Na rok 2025 dostajemy premię niecałe 5,4% i na 2026 7,5%. Dodajemy nowe prognozy do wykresu:

  nowa_prognoza <- as.numeric(prognoza@forecast$seriesFor)

  

  # Wykres najlepszego modelu z prawdziwymi obserwacjami

  plot(x=c(okres, 2025, 2026), y=c(premia, NA, NA), type="l", ylab = "")

  lines(x=c(okres, 2025, 2026), y=c(dopas, NA, NA), col="red")

  lines(x=c(okres, 2025, 2026), y=c(rep(NA, length(premia)-1), tail(dopas, 1), nowa_prognoza), col="green", lwd=2, type="o", cex=0.4)

  abline(h=0, col="grey")

  abline(h=mean(premia), col="blue", lwd=2)

  legend("topright", legend=c("premia za ryzyko S&P 500", "prognoza in-sample", "prognoza poza próbą", "średnia premia za ryzyko"), 

         col=c("black","red", "green", "blue"), lwd=2, cex = 0.6)

Prognoza krocząca

Czerwone linie są tak naprawdę jedynie dopasowaniem do danych, a nie prognozą. Aby uzyskać bardziej realny obraz dobroci modelu, wykonam prognozy kroczące za pomocą arfimaroll(), czyli prognozy out-of-sample. Jest to odpowiednik ugarchroll bez procesu GARCH (zob. ten art.). Zobaczmy 10 ostatnich kroczących-nawijanych prognoz, tj. od roku 2015:

  # Prognoza krocząca-nawijana

  klaster <- makeCluster(detectCores() - 1)

  premia_prognoza_roll <- arfimaroll(spec = premia_spec, data = premia, forecast.length = 10,

                                     refit.every = 1, refit.window = "recursive",

                                     window.size = NULL, solver = "nlminb", fit.control = list(ar2=0, ar3=0),

                                     cluster=klaster)

  premia_prognoza_roll <- resume(premia_prognoza_roll)

  stopCluster(klaster)

  prognoza_roll <- premia_prognoza_roll@forecast$density$Mu

  okres_prognozy_roll <- c(2015:2024)

  names(prognoza_roll) <- okres_prognozy_roll

  prognoza_roll 

  

  # Wykres z prognozą kroczącą

  plot(x=c(okres, 2025, 2026), y=c(premia, NA, NA), type="l", ylab = "", xlim=c(2015, 2026))

  lines(x=c(okres, 2025, 2026), y=c(dopas, NA, NA), col="red")

  lines(x=c(okres, 2025, 2026), y=c(rep(NA, length(premia)-1), tail(prognoza_roll, 1), nowa_prognoza), col="green", lwd=2, type="o", cex=0.4)

  lines(x=c(NA, okres_prognozy_roll), y=c(tail(dopas, 1), prognoza_roll), col="green", lwd=2, type="o", cex=0.4)

  abline(h=0, col="grey")

  abline(h=mean(premia), col="blue", lwd=2)

  legend("topright", legend=c("premia za ryzyko S&P 500", "prognoza in-sample", "prognoza poza próbą", "średnia premia za ryzyko"), 

         col=c("black","red", "green", "blue"), lwd=2, cex = 0.6)

  

Prognoza sugeruje spadek premii w 2025 i może wzrost 2026. Zauważmy tylko, że spadek nie poniżej zera, tylko poniżej średniej. No i oczywiście nie mówi to nic na temat zmian w obligacjach. 

Rozkład częstości i gęstości

Skupiając się na prognozowaniu łatwo jednak zgubić szeroki obraz, tj. ekonomiczne znaczenie premii za ryzyko. Dlatego musimy mieć w pamięci, że prawdopodobnie tak będzie, ale faktycznie może się zdarzyć niespodzianka. Ostatnie lata przyzwyczajają nie tylko zwykłych ludzi, ale i ekspertów do tego, że co roku z indeksów USA jest dodatnia stopa zwrotu. Od 2009 r do 2024, czyli przez 16 lat zdarzyły się zaledwie 2 ujemne roczne (całkowite) stopy zwrotu! Tworzy to złudzenie, że ryzyka nie ma, bo wystarczy przeczekać parę lat. Jednak ekonomia robi swoje i musi to jakoś wyrównać. Dlatego muszą nie tylko przydarzać się szokowe, bardzo silne i szybkie spadki (prędzej czy później). Muszą występować dłuższe okresy, w których obligacje skarbowe przyniosą większe zwroty od akcji. Gdyby nie ich nie było, to rzeczywiście nie byłoby ryzyka.

Nawet historyczna ujemna skośność (średnia < mediana) i kurtoza w rozkładzie premii o tym świadczy, zobaczmy:

  

  # Zapisujemy stare ustawienia graficzne

  old_par <- par(no.readonly = TRUE)

  

  # Dajemy więcej miejsca po prawej (dół, lewo, góra, prawo)

  par(mar = c(3.5, 5, 4, 5) - 0.1, mgp = c(2, 1, 0))

  

  # 1. Obliczamy histogram, ale nie rysujemy

  h <- hist(premia,

            breaks = 10,

            plot   = FALSE)

  

  # 2. Estymujemy KDE (Kernel Density Estimation)

  d <- density(premia,

               kernel = "gaussian",

               bw     = "nrd0")

  

  # 3. Rysujemy histogram na lewą oś Y (counts)

  plot(h,

       freq   = TRUE,

       col    = "lightgray",

       border = "white",

       xlab   = "Premia",

       ylab   = "Częstość",

       xlim  = range(h$breaks, d$x),

       ylim   = c(0, max(h$counts) * 1.1),

       main   = "Histogram i gęstość rozkładu – lewa oś: częstość, prawa oś: gęstość")

  

  # 4. Nałóż KDE na ten sam wykres

  par(new = TRUE)

  plot(d$x, d$y,

       type  = "l",

       lwd   = 2,

       axes  = FALSE,

       xlab  = "",

       ylab  = "",

       ylim  = c(0, max(d$y) * 1.1))

  

  # 5. Dodajemy prawą oś Y (density)

  axis(side = 4,

       at   = pretty(c(0, max(d$y) * 1.1)))

  mtext("Gęstość (KDE)", side = 4, line = 3)

  

  # 7. Przywracamy poprzednie ustawienia

  par(old_par)

  

Jeżeli rynek jest efektywny w długim okresie, to powinniśmy się spodziewać, że podobny rozkład wystąpi dla okresów wieloletnich. Zobaczmy dla 4-letnich. Zagregujmy stopy zwrotu:

  grupy <- ceiling(1:length(premia)/4)

  premia_4lata <- exp(tapply(X = log(1+premia), INDEX = grupy, FUN = sum)) - 1

I powtarzamy kod dla wykresu rozkładu częstości i gęstości, zastępując zmienną premia przez premia_4lata. Dostajemy:

Rozkład gęstości 4-letnich premii zawiera niewiele danych, ale też występuje na nim skośność. 

Dodatek dla premii za ryzyko nad bonami skarbowymi 3-miesięcznymi

Powtórzmy całą procedurę dla nadwyżki ponad bony skarbowe. To znaczy premia zawiera również ryzyko długoterminowych obligacji skarbowych w postaci zmienności stopy procentowej.

 

Stabilność średniej

W całym okresie średnia wyniosła aż 9%. Odwrotnie niż dla obligacji pierwsze 40 lat dały mniej, bo 8%, a kolejne prawie 10%. Czy więc tym razem to istotna zmiana średniej? Test Chowa wskazuje na stabilność:

	Chow test

data:  premia ~ 1
F = 0.22, p-value = 0.6

ACF i PACF

acf(premia, 15)

pacf(premia, 15)

Okazuje się, że jedynie ta 4-tego rzędu jest na granicy istotności (przekracza 0,2). To potwierdzałoby istnienie 4-letnich cykli, chociaż wykorzystanie tego fenomenu na pewno nie jest proste. 

Modelowanie

Sprawdzamy zatem AR(4). Efekt:

Coefficients: ar1 ar2 ar3 ar4 intercept 0 0 0 0.236 0.090 s.e. 0 0 0 0.112 0.024 sigma^2 estimated as 0.0282: log likelihood = 29.05, aic = -52.11 Training set error measures: ME RMSE MAE MPE MAPE MASE ACF1 Training set 0.000006594 0.168 0.1362 -266.2 433.2 0.6704 -0.06222

Gdyby zastosować taki model, dostalibyśmy średnią 9% (+/- 2,4%) -  faktycznie, tyle co średnia arytmetyczna. Dodatkowo współczynnik 4-tego rzędu 0,24 (+/- 0,11).

To teraz autoarfima z tymi samymi ustawieniami.

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(4,0,7)
Distribution	: std 

Optimal Parameters
------------------------------------
       Estimate  Std. Error     t value Pr(>|t|)
mu     0.098811    0.000012   8073.8245  0.00000
ar1    0.000000          NA          NA       NA
ar2    0.000000          NA          NA       NA
ar3    0.945906    0.000100   9503.1496  0.00000
ar4    0.054088    0.000008   6695.8385  0.00000
ma1    0.279539    0.000033   8524.8266  0.00000
ma2    0.000000          NA          NA       NA
ma3   -1.685603    0.000062 -27299.6027  0.00000
ma4   -0.135260    0.000014  -9448.6750  0.00000
ma5   -0.118209    0.000014  -8582.6486  0.00000
ma6    0.000000          NA          NA       NA
ma7   -0.036416    0.000006  -5762.9626  0.00000
sigma  0.111497    0.009553     11.6718  0.00000
shape 16.523609   12.901423      1.2808  0.20028

Robust Standard Errors:
       Estimate  Std. Error      t value Pr(>|t|)
mu     0.098811    0.000566   174.510530  0.00000
ar1    0.000000          NA           NA       NA
ar2    0.000000          NA           NA       NA
ar3    0.945906    0.010638    88.917365  0.00000
ar4    0.054088    0.000072   750.823674  0.00000
ma1    0.279539    0.002076   134.669756  0.00000
ma2    0.000000          NA           NA       NA
ma3   -1.685603    0.002733  -616.737764  0.00000
ma4   -0.135260    0.000331  -408.478534  0.00000
ma5   -0.118209    0.000610  -193.739688  0.00000
ma6    0.000000          NA           NA       NA
ma7   -0.036416    0.000020 -1799.605712  0.00000
sigma  0.111497    0.093611     1.191060  0.23363
shape 16.523609 1287.342190     0.012835  0.98976

LogLikelihood : 62.2351 

Information Criteria
------------------------------------
                    
Akaike       -1.3059
Bayes        -1.0081
Shibata      -1.3327
Hannan-Quinn -1.1865

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                       3.101 0.07826
Lag[2*(p+q)+(p+q)-1][32]    17.518 0.04329
Lag[4*(p+q)+(p+q)-1][54]    27.653 0.46556

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.3235  0.5695
Lag[2*(p+q)+(p+q)-1][2]    0.6261  0.6373
Lag[4*(p+q)+(p+q)-1][5]    1.6151  0.7116


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]     0.6869   2  0.7093
ARCH Lag[5]     2.5554   5  0.7681
ARCH Lag[10]    5.7786  10  0.8335

Nyblom stability test
------------------------------------
Joint Statistic:  411.572
Individual Statistics:             
mu    0.03948
ar3   0.04534
ar4   0.02661
ma1   0.04062
ma3   0.15471
ma4   0.02460
ma5   0.03678
ma7   0.02549
sigma 0.13303
shape 0.08077

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 2.29 2.54 3.05
Individual Statistic:	 0.35 0.47 0.75

Dostałem model ARMA(4, 7), ale przy usuniętych ar1, ar2, ma2 i ma6.

Jednak całkowity model okazuje się niestabilny, bo Joint Statistic = 412. Za to poszczególne współczynniki są stabilne shape, czyli kształt rozkładu reszt, jest niestabilny. Wynika to po prostu z tego, że rozkład t-Studenta nie jest idealnym rozkładem tutaj.

Test MCS

$includedR
[1] 2

$pvalsR
     [,1]
[1,]    0
[2,]    1

$excludedR
[1] 1

$includedSQ
[1] 2

$pvalsSQ
     [,1]
[1,]    0
[2,]    1

$excludedSQ
[1] 1

Test MCS potwierdza jednoznacznie, że drugi model (uzyskany z autoarfima) jest lepszy, albowiem każe odrzucić pierwszy model.

Prognozy kroczące i dwie nowe

*----------------------------------*
*        ARFIMA Model Forecast     *
*----------------------------------*

Horizon: 2
Roll Steps: 0
Out of Sample: 0

0-roll forecast: 
   T+1    T+2 
0.0309 0.1279 

Podobnie jak poprzednio - rok 2025 spadkowa premia.

Rozkład częstości i gęstości

Roczne premie:

4-letnie premie:

Porównując ze sobą obydwa rodzaje premii (nad bonami i nad obligacjami), to co się rzuca w oczy, to wpływ obligacji skarbowych na 4-letnią cykliczność stopy zwrotu. Po odjęciu bonów premia wykazuje pewną autokorelację 4-tego rzędu. Po odjęciu obligacji skarbowych premia ta autokorelacja mieści się w granicach nieistotności, a więc przypadkowości (jednak zwracam uwagę, że to PACF to taka klasyczna metoda szukania cykli, raczej wstępna). Wydaje się, że taka dysproporcja wynika z tego, że długoterminowe obligacje odzwierciedlają oczekiwania koniunktury, a ta przecież porusza się właśnie w takich kilkuletnich cyklach. Dostajemy więc pewnego rodzaju "dowód", że indeks S&P 500 podąża w jakimś zakresie w kierunku zmian gospodarczych.

Modele w zbiorze ufności w rugarch (język R)

Jeśli powiem, że każda prognoza gospodarcza czy giełdowa jest okraszona niepewnością, to będzie to banał. Ale jeśli powiem, że celem prognozy powinno być nie tyle odgadnięcie przyszłej wartości, co minimalizacja jej niepewności, to już banałem nie jest. Dla uproszczenia założymy, że niepewność to przedział ufności dla prognozy naszego modelu. Im więcej informacji zawiera model, tym przedział ufności staje się mniejszy, zbliżając się tym samym do prognozy punktowej. Widzimy więc, że idea nie jest wcale banalna: należy zebrać jak najwięcej istotnych informacji, aby zawężać przedział prognozy, zmniejszając tym samym ryzyko.

Analogicznie do przedziału ufności, Hansen et al.* zaproponowali zbiór ufności dla wielu modeli: z całej puli testowanych modeli, poszukiwane są te najlepsze, które wpadają do zbioru ufności. Im więcej posiadamy informacji, tym mniejszy będzie ten zbiór, zbliżając się tym samym do jednego najlepszego modelu.

Standardową metodą wyboru najlepszego modelu jest użycie kryterium informacyjnego, np. AIC, BIC czy HQIC. Niestety nie uwzględniają one niepewności czy przypadkowości parametrów. Z tym się ciągle zmagamy. Przykładowo, mam zawsze wątpliwości co do używania modeli parametrycznych, ponieważ wydają się one sztuczne, realizując uśrednianie na siłę. Dlatego tak ważne są testy stabilności parametrów. Ale stabilność to dopiero początek, bo: 

(1) jeśli spośród wielu modeli różnica w kryterium informacyjnym jest niewielka, to może ona wynikać po prostu z przypadku,

(2) dodanie tylko jednej obserwacji bardzo często zmienia wybór modelu.

W rugarch jest funkcja mcsTest (Model Confidence Set - MCS), która testuje, czy dany model należy zachować w zbiorze ufności czy go odrzucić. Ma ona następujące argumenty:

- losses: macierz lub ramka danych z wartościami strat dla każdego modelu (każda kolumna odpowiada jednemu modelowi). To w jaki sposób definiujemy straty, zależy od nas. W przypadku rozkładów symetrycznych, najczęściej stosuje się kwadraty odchyleń od prognozy lub wartości absolutne tych odchyleń. Dla rozkładów niesymetrycznych stosować można semi-odchylenia albo miary oparte na wartości zagrożonej.  

- alpha: poziom istotności (domyślnie 0.05).

- nboot: liczba replikacji bootstrapowych (domyślnie 100). Bootstrap to procedura wielokrotnego losowania próbki ze zwracaniem. nboot to właśnie liczba tych losowań.

- nblock: długość bloku dla bootstrappu blokowego (domyślnie 1). Jeżeli obserwacje mogą być zależne od siebie, to bootstrap musi losować blok występujących po sobie x obserwacji. Jeżeli uważamy, że każda obserwacja jest niezależna od innej, to x = 1. Im rząd autokorelacji większy, tym blok powinien być naturalnie też większy.

- boot: rodzaj bootstrapu ("stationary" lub "block"). Chociaż w obydwu przypadkach mamy do czynienia z blokami, to metoda blokowa, "block", oznacza, że całą próbę dzielimy na stałą liczbę przedziałów, tak że każdy przedział składa się bloku obserwacji o długości wyznaczonej dokładnie przez nblock. Potem dopiero losujemy blok, powtarzamy to losowanie dokładnie nboot razy. Natomiast  metoda stacjonarna, "stationary", tworzy bloki o losowej długości ze średnią równą nblock. Jeśli występuje autokorelacja w szeregu czasowym, to blok nie może być za krótki. Najczęściej "stationary" jest lepszy od "block".

Co właściwie robi ten test? W skrócie, kiedy mamy macierz strat dla każdego modelu, liczymy różnice między każdym modelem (prognozą każdego modelu) i są one odpowiednio uśredniane. Tworzymy bloki bootstrapowe i dla każdego bloku obliczamy średnie różnice między średnią stratą w bloku, średnią stratą dla całego modelu i średnią stratą ze wszystkich modeli. Dla każdego modelu będą różne wartości tworzące rozkład różnic. Wartości zbyt odchylone na plus (tzn. prawostronne ogony rozkładu) zostają odrzucone i tak uzyskujemy zbiór ufności. 

Oczywiście najpierw sami budujemy funkcję strat. Ja użyłem definicji najbardziej dla mnie intuicyjnej: jeżeli znak prognozy różni się od znaku wartości rzeczywistej, to bierzemy moduł z różnicy między nimi. Jeżeli znak jest ten sam, to strata wynosi zero. 

Przetestujmy tę funkcję. Będziemy używać mcsTest na zbiorze uczącym. Reguły używania poniższych funkcji w rugarch omówiłem szczegółowo tutaj. 

Stwierdziłem wcześniej, że próba między 300 a 400 tygodniowych danych giełdowych dla GARCH wydaje się już w miarę odpowiednia. Zrobię symulację EGARCH o próbie równej 330. Powiedzmy, że prawdziwy proces to EGARCH(1, 1)-ARMA(3, 1). Obok niego mamy też fałszywe - wszystkie EGARCH: ARMA(0, 0), ARMA(1, 0), ARMA(2, 0), ARMA(3, 0). Zobaczmy które z nich zostaną wykluczone przez mcsTest, a które przyjęte do zbioru ufności. Zrobimy test w dwóch wariantach: blokowy i stacjonarny. Ustawienia: nblock = 20 oraz nboot = 1000. 

library("rugarch")

library("xts")

# 1. Definicje stałych parametrów

stale_3_1 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  ar2    = 0.5,

  ar3    = 0.7,

  ma1    = 0.4,

  gamma1 = 0.05

)

stale_0_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  gamma1 = 0.05

)

stale_1_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  gamma1 = 0.05

)

stale_2_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  ar2    = 0.5,

  gamma1 = 0.05

)

stale_3_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  ar2    = 0.5,

  ar3    = 0.7,

  gamma1 = 0.05

)

stale_0_1 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ma1    = 0.4,

  gamma1 = 0.05

)

# 2. Specyfikacje modeli EGARCH(1,1) z różnymi specyfikacjami średniej:

spec_3_1 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(3,1), include.mean = TRUE),

  distribution.model = "norm",

  fixed.pars = stale_3_1

)

spec_0_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(0,0), include.mean = TRUE),

  distribution.model = "norm",

  fixed.pars = stale_0_0

)

spec_1_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(1,0), include.mean = TRUE),

  distribution.model = "norm",

  fixed.pars = stale_1_0

)

spec_2_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(2,0), include.mean = TRUE),

  distribution.model = "norm",

  fixed.pars = stale_2_0

)

spec_3_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(3,0), include.mean = TRUE),

  distribution.model = "norm",

  fixed.pars = stale_3_0

)

spec_0_1 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(0,1), include.mean = TRUE),

  distribution.model = "norm",

  fixed.pars = stale_0_1

)

# 3. Symulacja 330 obserwacji według poprawnego modelu (ARMA(3,1))

set.seed(765)

symulacja <- ugarchpath(spec_3_1, n.sim = 330)

dane_symulowane <- as.numeric(fitted(symulacja))

# 4. Dopasowanie modeli do całego szeregu (in-sample)

dopasowanie_3_1 <- ugarchfit(spec = spec_3_1, data = dane_symulowane)

dopasowanie_0_0 <- ugarchfit(spec = spec_0_0, data = dane_symulowane)

dopasowanie_1_0 <- ugarchfit(spec = spec_1_0, data = dane_symulowane)

dopasowanie_2_0 <- ugarchfit(spec = spec_2_0, data = dane_symulowane)

dopasowanie_3_0 <- ugarchfit(spec = spec_3_0, data = dane_symulowane)

dopasowanie_0_1 <- ugarchfit(spec = spec_0_1, data = dane_symulowane)

# 5. Zapisanie wyników dopasowania w liście (dla przejrzystości)

dopasowania <- list(

  model_3_1 = dopasowanie_3_1,

  model_0_0 = dopasowanie_0_0,

  model_1_0 = dopasowanie_1_0,

  model_2_0 = dopasowanie_2_0,

  model_3_0 = dopasowanie_3_0,

  model_3_0 = dopasowanie_0_1

)

print("Dopasowane modele:")

print(names(dopasowania))

# 6. Definicja funkcji strat:

Funkcja_strat <- function(rzeczywiste, prognoza) {

  roznica <- rzeczywiste - prognoza

  strata <- ifelse(sign(rzeczywiste) != sign(prognoza), abs(roznica), 0)

  return(strata)

}

# 7. Obliczenie funkcji strat dla każdego modelu (in-sample)

strata_3_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_1))

strata_0_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_0))

strata_1_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_1_0))

strata_2_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_2_0))

strata_3_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_0))

strata_0_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_1))

# Łączenie strat w macierz – wiersze odpowiadają obserwacjom, kolumny modelom

straty <- cbind(

  model_3_1 = strata_3_1,

  model_0_0 = strata_0_0,

  model_1_0 = strata_1_0,

  model_2_0 = strata_2_0,

  model_3_0 = strata_3_0,

  model_3_0 = strata_0_1

)

# 8. Usunięcie wierszy składających się tylko z zer

straty <- straty[rowSums(straty != 0) > 0, ]

straty <- coredata(straty)

# 9. Test MCS

# 9.1. Wariant 1: stacjonarny

mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "stationary")

# 9.2. Wariant 1: blokowy

mcs_block <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "block")

# 10. Wyświetlenie wyników:

cat("\nWyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):\n")

print(mcs_stationary)

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 1 5

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.004
[5,] 0.697
[6,] 1.000

$excludedR
[1] 6 2 4 3

$includedSQ
[1] 1 5

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.003
[5,] 0.697
[6,] 1.000

$excludedSQ
[1] 6 2 4 3

Należy zauważyć, że mamy tu 2 typy testów: statystyka R (Range statistic) oraz SQ (Semi-Quadratic statistic). Pierwszy z nich skupia się na różnicach w maksymalnych stratach między modelami, a drugi na różnicach w kwadratach strat. Najlepiej, jeśli obydwa wskazują ten sam zbiór, bo wtedy dostajemy potwierdzenie, że zbiór ufności zawiera odpowiednie modele.

Jak widać obydwa testy rzeczywiście wskazują to samo: included 1 oraz 5, czyli zbiór ufności zawiera tylko dwa modele - pierwszy oraz piąty. Pozostałe modele zostają odrzucone (excluded), zarówno przez R jak SQ. Możemy powiedzieć, że test zadziałał poprawnie, mimo przyjęcia błędnego modelu ARMA(3, 0). Model ten jest podobny do prawidłowego ARMA(3, 1), a co więcej poprawny model ma największe p-value = 1 (piąty 0.697). O tym może jeszcze powiem. Otóż wydrukowana lista może być trochę myląca, bo kolejność nie pokazuje p-value w kolejności od modelu 1 do 6, tylko w kolejności od najbardziej istotnego do najmniej istotnego wyniku dotyczącego odrzucenia modelu ze zbioru ufności. Im mniejsze p-value, tym większe prawdopodobieństwo, że dany model należy wyrzucić ze zbioru. 

Żeby zawęzić zbiór ufności, należy zwiększyć próbę. Po zwiększeniu jej do 1600, dostałem:

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 1

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.047
[6,] 1.000

$excludedR
[1] 6 2 4 3 5

$includedSQ
[1] 1

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.047
[6,] 1.000

$excludedSQ
[1] 6 2 4 3 5

Przy tak dużej próbie test nie ma już wątpliwości, że jedyny prawidłowy model to ten pierwszy.

Bardzo podobne wyniki dostałem dla boot = "block" i nie ma sensu pokazywać ich. Podobnie jak dla "stationary" musiałem ustawić próbę  liczną 1600, aby uzyskać wykluczenie wszystkich fałszywych modeli. Co jednak zauważyłem, to to, że w przypadku "block" wykluczenie modelu 5 było ledwo, ledwo, na styk, bo p-value wyniosło 0.049 (dla stationary też, ale niższe, bo 0.047).  Stacjonarny test generalnei jest lepszą alternatywą.

Powyższe przykłady zakładały rozkład normalny reszt modelu. Dla rynków może być lepszy inny, np. skośny rozkład t-Studenta. W takiej sytuacji potrzebne są dwa dodatkowe parametry do stałych: skew oraz shape. Zamiast "norm" będzie "sstd". Przykładowo kod będzie taki:

# 1. Definicje stałych parametrów

stale_3_1 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  ar2    = 0.5,

  ar3    = 0.7,

  ma1    = 0.4,

  gamma1 = 0.05,

  skew = 1.5,

  shape = 5

)

stale_0_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  gamma1 = 0.05,

  skew = 1.5,

  shape = 5

)

stale_1_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  gamma1 = 0.05,

  skew = 1.5,

  shape = 5

)

stale_2_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  ar2    = 0.5,

  gamma1 = 0.05,

  skew = 1.5,

  shape = 5

)

stale_3_0 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ar1    = -1,

  ar2    = 0.5,

  ar3    = 0.7,

  gamma1 = 0.05,

  skew = 1.5,

  shape = 5

)

stale_0_1 <- list(

  omega  = 1,

  alpha1 = 0.3,

  beta1  = 0.2,

  mu     = 0.0001,

  ma1    = 0.4,

  gamma1 = 0.05,

  skew = 1.5,

  shape = 5

)

# 2. Specyfikacje modeli EGARCH(1,1) z różnymi specyfikacjami średniej:

spec_3_1 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(3,1), include.mean = TRUE),

  distribution.model = "sstd",

  fixed.pars = stale_3_1

)

spec_0_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(0,0), include.mean = TRUE),

  distribution.model = "sstd",

  fixed.pars = stale_0_0

)

spec_1_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(1,0), include.mean = TRUE),

  distribution.model = "sstd",

  fixed.pars = stale_1_0

)

spec_2_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(2,0), include.mean = TRUE),

  distribution.model = "sstd",

  fixed.pars = stale_2_0

)

spec_3_0 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(3,0), include.mean = TRUE),

  distribution.model = "sstd",

  fixed.pars = stale_3_0

)

spec_0_1 <- ugarchspec(

  variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

  mean.model     = list(armaOrder = c(0,1), include.mean = TRUE),

  distribution.model = "sstd",

  fixed.pars = stale_0_1

)

# 3. Symulacja 330 obserwacji według poprawnego modelu (ARMA(3,1))

set.seed(55)

symulacja <- ugarchpath(spec_3_1, n.sim = 330)

dane_symulowane <- as.numeric(fitted(symulacja))

# 4. Dopasowanie modeli do całego szeregu (in-sample)

dopasowanie_3_1 <- ugarchfit(spec = spec_3_1, data = dane_symulowane)

dopasowanie_0_0 <- ugarchfit(spec = spec_0_0, data = dane_symulowane)

dopasowanie_1_0 <- ugarchfit(spec = spec_1_0, data = dane_symulowane)

dopasowanie_2_0 <- ugarchfit(spec = spec_2_0, data = dane_symulowane)

dopasowanie_3_0 <- ugarchfit(spec = spec_3_0, data = dane_symulowane)

dopasowanie_0_1 <- ugarchfit(spec = spec_0_1, data = dane_symulowane)

# 5. Zapisanie wyników dopasowania w liście (dla przejrzystości)

dopasowania <- list(

  model_3_1 = dopasowanie_3_1,

  model_0_0 = dopasowanie_0_0,

  model_1_0 = dopasowanie_1_0,

  model_2_0 = dopasowanie_2_0,

  model_3_0 = dopasowanie_3_0,

  model_3_0 = dopasowanie_0_1

)

print("Dopasowane modele:")

print(names(dopasowania))

# 6. Definicja funkcji strat:

Funkcja_strat <- function(rzeczywiste, prognoza) {

  roznica <- rzeczywiste - prognoza

  strata <- ifelse(sign(rzeczywiste) != sign(prognoza), abs(roznica), 0)

  return(strata)

}

# 7. Obliczenie funkcji strat dla każdego modelu (in-sample)

strata_3_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_1))

strata_0_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_0))

strata_1_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_1_0))

strata_2_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_2_0))

strata_3_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_0))

strata_0_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_1))

# Łączenie strat w macierz – wiersze odpowiadają obserwacjom, kolumny modelom

straty <- cbind(

  model_3_1 = strata_3_1,

  model_0_0 = strata_0_0,

  model_1_0 = strata_1_0,

  model_2_0 = strata_2_0,

  model_3_0 = strata_3_0,

  model_3_0 = strata_0_1

)

# 8. Usunięcie wierszy składających się tylko z zer

straty <- straty[rowSums(straty != 0) > 0, ]

straty <- coredata(straty)

# 9.1. Wariant stacjonarny

mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "stationary")

# 10. Wyświetlenie wyników:

cat("\nWyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):\n")

print(mcs_stationary)

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 5 1

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.071
[6,] 1.000

$excludedR
[1] 6 2 4 3

$includedSQ
[1] 5 1

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.071
[6,] 1.000

$excludedSQ
[1] 6 2 4 3

Jak widać, ponownie można by podwyższyć p-value do 10%, żeby dostać tylko model 1. Tak czy inaczej test działa bardzo dobrze także dla innych rozkładów.

Oczywiście to bardzo wzorcowe przykłady i w praktyce będzie to trudniejsze, ale to jest dobry punkt wyjścia.

* Hansen, P., A. Lunde, and J. Nason (2011). The model confidence set. Econometrica 79, 453–497.

Krzyż śmierci na S&P500 - czy faktycznie ma znaczenie?

 Niektórzy zwrócili niedawno uwagę, że na S&P 500 występuje tzw. krzyż śmierci, czyli 50-dniowa średnia krocząca przecinająca od góry 200-dniową średnią kroczącą:

Na ten sygnał zwrócił m.in. portal xyz.pl , dodając jednak notkę, że analitycy firmy LPL Financial zaobserwowali, że historycznie w ciągu kolejnego miesiąca indeks spadał średnio niecały 1%, a w ciągu 3 miesięcy nawet rósł. Mowa jednak o medianie, a więc średnia i dominanta mogły się sporo różnić. Podanie samej mediany jest niewystarczające.

Sprawdźmy jak to rzeczywiście wyglądało dotychczas. Zaczniemy tak jak ta firma od 1950 r. Trzeba precyzyjnie określić warunki. Nie wystarczy, że SMA50 przetnie SMA200 z góry w dół, bo za chwilę może jeszcze zawrócić i średnie mogą się znowu dotknąć albo w ogóle SMA50 może przebić z dołu w górę. Jeśli taka sytuacja dzieje się w np. w ciągu tygodnia, to sygnału nie będzie. Tak więc określamy minimalne okno tego zdarzenia. Przyjmę cały tydzień, co oznacza 3 kroki do tyłu i 3 kroki do przodu od momentu przecięcia. Taki układ warunków jest już wystarczający, chociaż może nie być przekonujący, bo nawet po przecięciu się obie średnie mogą rosnąć. Tak na intuicję wydaje się, że  powinny spadać. Z drugiej strony już sama nazwa "krzyż" sugeruje dowolne przecięcie, np. rosnącej SM200 i spadkowej SMA50. Z powodu tej różnicy w interpretacji krzyża śmierci, zrobimy 2 sprawdziany:

1) Są warunki: (a) SMA50 <= SMA200, (b) SMA50[t-3] > SMA200[t-3], (c) SMA50[t+3] < SMA200[t+3] ,

2) To samo co p. (1) i dodatkowo warunek, że obie średnie spadają.

Ostatnia sprawa to sam okres przyszłej stopy zwrotu. Sprawdzimy 1-3 miesięczne zwroty do przodu. Za 1 miesiąc przyjmiemy 21 dni, czyli 21 sesji.

Ad 1) Kod w R: 

# 1. Pobieramy dane S&P500, ładujemy pakiety, np. xts.

# 2. Przekształcamy ceny na logarytmy

logCena <- log(sp500)

# 3. Obliczamy 50-dniową i 200-dniową średnią kroczącą

sma50  <- SMA(logCena, n = 50)

sma200 <- SMA(logCena, n = 200)

# 4. Ustalamy wartość przesunięcia (krok) dla warunków – tu krok = 3

krok <- 3

# Warunki dla sygnału

war1 <- sma50 <= sma200

war2 <- lag(sma50, k = krok) > lag(sma200, k = krok)           # SMA50[t-krok] > SMA200[t-krok]

war3 <- lag(sma50, k = -krok) < lag(sma200, k = -krok)                # SMA50[t+krok] < SMA200[t+krok]

war4 <- TRUE

# 5. Tworzymy wektor sygnału: 1, gdy wszystkie warunki są spełnione, NA w przeciwnym przypadku

sygnal <- ifelse(war1 & war2 & war3 & war4, 1, 0)

sygnal <- na.omit(sygnal)

datySygnalow <- index(sygnal)

# Przypisujemy daty wszystkich sygnałów

print("Daty sygnałów (pełny szereg):")

print(datySygnalow)

# 6. Teraz dodajemy warunek, aby wyłapać jedynie moment przejścia z 0 na 1.

# To pozwala zachować tylko pierwszy dzień pojawienia się sygnału.

sygnal_krok1 <- lag(sygnal, k = 1)  # przesunięcie danych o jeden dzień do tyłu (opóźnione)

sygnal_krok1[is.na(sygnal_krok1)] <- 0  # traktujemy pierwszy dzień jako 0

unikalneDatySygnalow <- which(sygnal_krok1 == 0 & sygnal == 1)

unikalneDatySygnalow <- datySygnalow[unikalneDatySygnalow] 

print("Unikalne daty sygnału (przejście z 0 na 1):")

print(unikalneDatySygnalow)

# 7. Obliczamy dzienne logarytmiczne stopy zwrotu

logZwrot <- na.omit(diff(logCena))

# 8. Używamy rollapply do sumowania logarytmicznych stóp zwrotu na horyzoncie 21 sesji,

# ale dnia następnego po sygnale 

oknoPrzyszlosci <- 21

logZwrot_suma <- lag(rollapply(logZwrot,

                           width = oknoPrzyszlosci,

                           FUN = sum,

                           align = "left"), -1)

# 9. Przekształcamy sumy logarytmicznych stóp zwrotu na zwykłe stopy zwrotu: exp(sum) - 1

miesZwrot <- exp(logZwrot_suma) - 1

# 10. Dopasowujemy daty unikalnych sygnałów do obliczonych miesięcznych stóp zwrotu

wyniki <- data.frame(Date = unikalneDatySygnalow,

                     MonthlyReturn = as.numeric(miesZwrot[unikalneDatySygnalow]))

wyniki <- na.omit(wyniki)

print("Wyniki (data unikalnego sygnału + miesięczna stopa zwrotu):")

print(wyniki)

mean(wyniki$MonthlyReturn)

# 11. Rysujemy histogram miesięcznych stóp zwrotu dla wyłapanych sygnałów

h <- hist(wyniki$MonthlyReturn,

     breaks = 7,

     col = "lightblue",

     border = "black",

     main = "Histogram miesięcznych stóp zwrotu",

     xlab = "Stopa zwrotu",

     ylab = "Liczba sygnałów",

     freq = NULL)

# 12. Przekształcamy liczebności na prawdopodobieństwa

czestosc <- h$counts / sum(h$counts)

# 13. Ustawiamy układ dwóch wykresów oraz marginesy

par(mfrow = c(2, 1), mar = c(2, 5, 2.2, 3))

slupkiZwrotow <- 100 * round(h$mids, 3)

# 14. Rysujemy wykres słupkowy prawdopodobieństwa 1-miesięcznych stóp zwrotu

b <- barplot(czestosc,

             names.arg = slupkiZwrotow,

             col       = "lightblue",

             border    = "black",

             main      = "Prawdopodobieństwo 1-miesięcznych stóp zwrotu",

             xlab      = "",

             ylab      = "Prawdopodobieństwo",

             yaxt      = "n",

             space     = 0,

             cex.main  = 0.9)

# 15. Dodajemy oś Y z wartościami prawdopodobieństwa

axis(2, at = round(czestosc, 2), las = 2, cex = 0.7)

# 16. Zaznaczamy dominantę prawdopodobieństwa linią pionową

abline(v  = b[which.max(czestosc)],

       col  = "red",

       lty  = 2,

       lwd  = 2)

# 17. Rysujemy wykres słupkowy skumulowanego prawdopodobieństwa

par(mar = c(4, 5, 1.5, 3))

barplot(cumsum(czestosc),

        names.arg = slupkiZwrotow,

        col       = "lightblue",

        border    = "black",

        main      = "Skumulowane prawdopodobieństwo 1-miesięcznych stóp zwrotu",

        xlab      = "",

        ylab      = "Prawdopodobieństwo",

        yaxt      = "n",

        space     = 0,

        cex.main  = 0.9)

# 18. Dodajemy oś Y z wartościami skumulowanego prawdopodobieństwa

axis(2, at = round(cumsum(czestosc), 2), las = 2, cex = 0.7)

# 19. Dodajemy opis osi X

mtext(text = "Stopa zwrotu, 1-miesięczna (%)", side = 1, padj = 3.5)

# 20. Zaznaczamy dominantę skumulowanego prawdopodobieństwa linią poziomą

abline(h    = cumsum(czestosc)[which.max(czestosc)],

       col  = "red",

       lty  = 2,

       lwd  = 2)

# 21. Zaznaczamy dominantę pojedynczego prawdopodobieństwa linią pionową

abline(v    = b[which.max(czestosc)],

       col  = "red",

       lty  = 2,

       lwd  = 2)

Liczba przypadków = 35.

Otrzymane statystyki:

Okazuje się, że dominuje ujemna stopa zwrotu z medianą -2,5%, ale jej częstość to niecałe 0,4, a  skumulowane prawdopodobieństwo empiryczne 0,57. Oznacza to, że na prawie 60% będzie spadać w kolejnym miesiącu.  Skośność jednak powoduje, że sporo też jest dodatnich stóp zwrotu wokół +2,5%.

Sprawdzamy 2-miesięczne zwroty, tj. oknoPrzyszlosci <- 42.

 W kolejnych dwóch miesiącach po krzyżu śmierci dominowały wzrosty z medianą 2,5%, ale znów skośność powoduje duże ryzyko spadków.

3-miesięczna stopa zwrotu (oknoPrzyszlosci <- 63):

Po 3 miesiącach dostajemy w zasadzie czystą przypadkowość - możliwe są zarówno wzrosty jak i spadki.

Ad 2) Cały kod będzie ten sam, jedynie war4 przyjmuje teraz postać

war4 <- lag(diff(sma50, lag = krok*2 + 1), k = -krok) < 0 & lag(diff(sma200, lag = krok*2 + 1), k = -krok) < 0  # oba spadają

Dostałem 20 przypadków.

a) miesięczne zwroty:

Dostajemy tu zbliżone wyniki do pierwszej definicji krzyża śmierci, z tym, że w tym przypadku prawdopodobieństwa spadków są wyższe niż poprzednio, bo aż 0,75.

b) 2-miesięczne zwroty:

Dominanta praktycznie nie istnieje. Można tylko powiedzieć, że zwroty przechylają się minimalnie na plus.

c) 3-miesięczne zwroty

Prawdopodobieństwo 3-miesięcznych wzrostów jest w tym wypadku obiektywnie większe niż spadków (chociaż okolice do -0% też nie jest takie małe). 

Porównując otrzymane statystyki ze wspomnianym na początku badaniem mogę powiedzieć, że się one pokrywają. Oni dostali średniomiesięcznie (medianę) -1%, ja dostałem dominantę -2,5%, a przy skumulowanym prawdopodobieństwie 0,57, a 0,57*2,5 = 1,43. Podobnie jak oni dostałem także dodatnią 3-miesięczną średnią. 

Naturalne staje się na koniec pytanie, czy bieżący krzyż śmierci spełnia w ogóle określone definicje? Jeśli spełnia drugą, to oczywiście spełnia też pierwszą (bo druga jest mocniejsza), więc sprawdzamy najpierw drugą. Wszystkie sygnały są zapisane w zmiennej unikalneDatySygnalow:

> unikalneDatySygnalow
 [1] "1953-05-11" "1957-09-26" "1960-02-15" "1962-05-08" "1965-07-22" "1968-02-27" "1969-06-23"
 [8] "1977-03-03" "1980-04-22" "1984-02-06" "1987-11-04" "1990-09-07" "1994-04-19" "2000-10-30"
[15] "2010-07-02" "2011-08-12" "2015-08-28" "2016-01-11" "2018-12-07" "2020-03-27" "2025-04-14"

Ostatnia pozycja to właśnie ostatni sygnał z 14 kwietnia. Czyli druga definicja zostaje spełniona. Stąd wnioskujemy, że kolejne 21 sesji będzie spadkowych, z "prawdopodobieństwem" 0,75. To znaczy, że do połowy maja należy spodziewać się generalnie spadków. Kolejne 2 miesiące są niejednoznaczne, z lekkim przechyleniem na plus, a 3 miesiące "powinny" być dodatnie. Trudno mi uwierzyć w to ostatnie, ale wszystkiego trzeba się spodziewać przy takim rozedrganiu, z jakim mamy do czynienia.

Co nam powie Equity Market Volatility Index?

Oczywiście na giełdzie nie ma nic pewnego, ale prognozuję, że kwiecień będzie spadkowy. Spójrzmy na Indeks Zmienności Rynku Akcji USA (US Equity Market Volatility Index) w porównaniu z S&P 500 (od stycznia 2005 do marca 2025):

Największe piki zmienności występowały w końcówce głębokich obsunięć. Mniejsze piki występowały również przy zakończeniu hossy. Wyjątkiem można nazwać rok 2022, ale też z perspektywy czasu nie były to jakieś gigantyczne spadki, co można tłumaczyć większym wpływem geopolityki niż gospodarki. Zupełnie inaczej sprawa się ma w przypadku roku 2020 - krótkie, ale duże załamanie, sztucznie wywołane przez przymusową izolację i zakazy handlu (tzw. lockdown). Pamiętam, jak nie mogłem uwierzyć, że rząd może sobie tak po prostu, z dnia na dzień niemal, zakazać prowadzić działalności gospodarczej. W sumie zobaczmy to z bliska:

Cofnijmy się kolejne kilka lat:

Znowu chwilowy wstrząs w 2011 r. nie był całkowicie naturalny, bo dotyczył państw UE, które musiały poczuć na własnej skórze, do czego prowadzi socjalizm utrzymywany przez długie lata przez frywolne zadłużanie się. Unia w swojej "dobroduszności" postanowiła uratować je, czyli zrefinansować ich niespłacalne długi, tzn. jeszcze je dodłużyć. Koszty tej operacji były niewidoczne z dnia na dzień, dlatego wielu mogło uznać, że Grecji, Włochom czy Hiszpanii upiekło się. I w zasadzie mieli rację, chociaż mogli nie zauważyć, że cała UE od tej chwili zaczęła powoli upadać.  Nasza giełda jest tego świetnym przykładem.

Chociaż jak powiedziałem zmienność 2011 nie była całkiem naturalna, to jednak była skutkiem kryzysu finansowego parę lat wcześniej - a ten miał już duży wpływ na gospodarkę i przez to zmienność była naturalna:

Wróćmy do teraźniejszości, mamy:

Wzrost niepewności, który tu obserwujemy, dopiero się zaczął. Jest pewne, że w kwietniu będzie dużo wyższy, ponieważ cła USA ogłosiły na początku kwietnia, a w kolejnych dniach zaczęły się pojawiać plotki, częściowo zdementowane, że cła będą zmniejszane w trakcie negocjacji. Takie biznesowe podejście, jakby się grało w pokera, zupełnie bez zwracania uwagi na otoczenie, to prawdziwy koszmar dla gospodarki. Światowej, ale w tym wypadku głównie USA, bo amerykańcy przedsiębiorcy nie mogą od niej uciec. A to właśnie zagraniczni inwestorzy będą uciekać z tak niestabilnego otoczenia (stąd dolar spada i będzie dalej spadał). 

Spadki muszą być potężne, bo nakładają dwa negatywne czynniki: cła oraz niepewność, czy za tydzień nie zostaną zniesione. Na efektywnym rynku, ale też czysto psychologicznie, nie będzie wiecznie reakcji na ewentualną zmianę decyzji Trumpa. Ta niepewność musi zostać zdyskontowana w postaci wyższej awersji do ryzyka. Psychologicznie natomiast powiemy, że ten sam bodziec będzie działał z czasem coraz słabiej.

Sama liniowa korelacja stopy zwrotu z EMV wynosi -33%, ale niewiele to mówi, bo ich zależność dużo lepiej pokazuje wykres punktowy:

Ta zależność jest ewidentnie nieliniowa i składa się z dwóch części. Pierwsza połowa zmienności to czyste ryzyko, które generuje albo zysk, albo stratę (zauważmy jak równo wyszło - ok. 11-12% na plus jak i na minus). Druga połowa zmienności jest ujemnie proporcjonalna do zmian procentowych. Dla EMV >= 40 dostajemy same ujemne stopy zwrotu i im większe EMV, tym mocniej ujemna jest zmiana indeksu. W sumie dostajemy wskaźnik mówiący, że jeśli dochodzi do poziomu 40 i przekracza go,  spadki są bardzo prawdopodobne.  

Ale właściwie co to jest ten EMV? Jego konstrukcja jest dość kontrowersyjna, ale wygląda na to, że poprawna. Otóż sam skrót ma dwa znaczenia. Pierwszy pochodzi od jego nazwy, czyli Equity Market Volatility. Drugi to akronim od economy, market, volatility. Każde z tych słów to hasła-skojarzenia,  stanowiące grupę różnych słów związanych z tym hasłem. Np. pod "market" jest m.in. equity, S&P, pod "volatility" jest risk, uncertainty itp. Te słowa-klucze szukane są następnie w prasie (gazetach) i zwyczajnie się je zlicza. Uzyskane liczby odpowiednio się skaluje i dostaje ten sposób wskaźnik dla każdego miesiąca. Ponieważ największe szczyty zmienności dostajemy blisko zakończenia bessy, to prawdopodobnie media napędzają wszechogarniający pesymizm "ulicy", która pod wpływem psychicznej presji najmocniej pozbywa się aktywów blisko dołka. 

Chociaż trudno w to uwierzyć, ale EMV - wskaźnik obliczony tylko w oparciu o zliczanie słów w prasie - ma korelację z VIX ok. 70%. Bo to by znaczyło, że bez danych rynkowych można poprawnie oszacować bieżącą zmienność rynku. Nie jest to może szokujące, ale mocno zaskakujące.

Może się wydawać, że EMV nie jest tak przydatny jak VIX, który podawany jest z dzienną częstością, a nie miesięczną. Tylko zobaczmy jedną rzecz. Porównajmy oba wskaźniki w tym samym okresie po dopasowaniu VIX do miesięcznej częstości:

Zoom od 2019:

W przypadku okresu pandemicznego oba wskaźniki rosły i spadały równolegle, natomiast ostatni wstrząs polityki celnej Trumpa był poprzedzony tak samo szybkim wzrostem EMV, ale już nie VIX. Różnica jest taka, że pierwsze wydarzenie było całkowicie nieprzewidywalne, natomiast drugie było mocno nagłaśniane w mediach. Przypuszczam, że różnica leży w tym, że informacje o cłach były z jednej strony stopniowo odkrywane, a z drugiej ciągle powtarzane, co spowodowało, że rynek zachował się trochę jak powoli gotowana żaba. VIX zaczął gwałtownie rosnąć dopiero w kwietniu, którego tutaj siłą rzeczy nie można było umieścić. Bo dziennie jest skok:

Wygląda więc na to, że EMV przewidział z miesięcznym wyprzedzeniem spadki amerykańskiej giełdy (ponieważ w marcu przekroczył poziom 40), podczas gdy VIX po prostu je odzwierciedla. Spójrzmy też na oba wskaźniki przy zakończeniu hossy 2007:

EMV praktycznie wskazał szczyt S&P500, który korespondował ze szczytem EMV (wartością min. 40). Tego nie można powiedzieć o VIX. Z kolei przy zakończeniu bessy taka rozbieżność nie występuje.