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 odjęliś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
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)
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 można w nim ustawić metodę "full". 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.75BIC 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 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))
premia_fit <- arfimafit(spec = premia_spec, data = premia, solver = "nlminb")
prognoza <- arfimaforecast(premia_fit, n.ahead = 2)
*----------------------------------*
* ARFIMA Model Forecast *
*----------------------------------*
Horizon: 2
Roll Steps: 0
Out of Sample: 0
0-roll forecast:
T+1 T+2
0.14876 0.02238
Na rok 2025 dostajemy premię niecałe 15% i na 2026 2,2%. Dodajemy nowe prognozy do wykresu:
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:
Podobny rezultat jak w przypadku premii nad obligacjami uzyskałem dla premii nad bonami:
Chow test
data: premia ~ 1
F = 0.22, p-value = 0.6acf(premia, 15)
pacf(premia, 15)
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 |
*----------------------------------*
* ARFIMA Model Fit *
*----------------------------------*
Mean Model : ARFIMA(4,0,5)
Distribution : std
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.161833 0.000000 2959700.74 0
ar1 0.151799 0.000009 17708.57 0
ar2 0.000000 NA NA NA
ar3 0.000000 NA NA NA
ar4 0.848201 0.000052 16219.51 0
ma1 -0.241419 0.000017 -13997.21 0
ma2 -0.672269 0.000053 -12747.99 0
ma3 0.327000 0.000022 14596.10 0
ma4 -1.166367 0.000080 -14596.10 0
ma5 0.063512 0.000001 68899.16 0
sigma 0.119640 0.000749 159.68 0
shape 100.000000 0.776030 128.86 0
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.161833 0.000002 66502.9198 0.000000
ar1 0.151799 0.000199 764.2664 0.000000
ar2 0.000000 NA NA NA
ar3 0.000000 NA NA NA
ar4 0.848201 0.002818 301.0231 0.000000
ma1 -0.241419 0.001761 -137.1293 0.000000
ma2 -0.672269 0.002135 -314.8811 0.000000
ma3 0.327000 0.001856 176.1863 0.000000
ma4 -1.166367 0.007626 -152.9473 0.000000
ma5 0.063512 0.000124 512.7090 0.000000
sigma 0.119640 0.037396 3.1992 0.001378
shape 100.000000 41.141275 2.4306 0.015072
LogLikelihood : 56.39
Information Criteria
------------------------------------
Akaike -1.15972
Bayes -0.86197
Shibata -1.18658
Hannan-Quinn -1.04034
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.3203 0.5714
Lag[2*(p+q)+(p+q)-1][26] 11.3668 0.9999
Lag[4*(p+q)+(p+q)-1][44] 21.5256 0.5907
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 1.925 0.1653
Lag[2*(p+q)+(p+q)-1][2] 2.722 0.1669
Lag[4*(p+q)+(p+q)-1][5] 3.404 0.3381
ARCH LM Tests
------------------------------------
Statistic DoF P-Value
ARCH Lag[2] 3.306 2 0.1915
ARCH Lag[5] 4.479 5 0.4827
ARCH Lag[10] 7.391 10 0.6881
Nyblom stability test
------------------------------------
Joint Statistic: 669.4
Individual Statistics:
mu 0.029974
ar1 0.017018
ar4 0.027913
ma1 0.009996
ma2 0.036438
ma3 0.056745
ma4 0.087832
ma5 0.021669
sigma 0.093384
shape 0.629655
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 2.29 2.54 3.05
Individual Statistic: 0.35 0.47 0.75Okazało się, że najlepszym modelem nie jest wcale AR(4), ale ARMA(4, 5), przy usuniętych ar2 i ar3.
Z poszczególnych współczynników jedynie shape, czyli kształt rozkładu reszt, jest trochę niestabilny. Wynika to po prostu z tego, że rozkład t-Studenta nie jest idealnym rozkładem tutaj, ale nie jest to wcale tak istotne. Ważne, że inne współczynniki są stabilne. Jednak całkowity model okazuje się niestabilny, bo Joint Statistic = 669. To może być jednak efekt tego kształtu rozkładu.
$includedR
[1] 2
$pvalsR
[,1]
[1,] 0
[2,] 1
$excludedR
[1] 1
$includedSQ
[1] 2
$pvalsSQ
[,1]
[1,] 0
[2,] 1
$excludedSQ
[1] 1
