Poprawka do analizy premii za ryzyko S&P 500

Wykonana ostatnio analiza premii za ryzyko S&P 500 wskazywała, że roczne zmiany można modelować za pomocą ARMA(7, 4). Dziś widzę, że zawierała 3 problemy. Po pierwsze metoda "full" w funkcji autoarfima (pakiet rugarch dla R), którą wówczas uważałem za bardziej precyzyjną, lepszą, wcale taka nie jest. Poprzez symulacje przetestowałem ją i okazało się, że "partial" dużo lepiej znajduje prawidłowe rzędy. Wydawało się, że skoro "full" jest tak potwornie wolna, bo sprawdza wszystkie możliwe kombinacje (a więc np. 1 rzędu nie ma, ale jest drugi i trzeci, albo drugiego nie ma, ale jest pierwszy i trzeci, albo pierwszego i drugiego nie ma, ale jest trzeci itd.), to będzie najdokładniejsza. Chociaż na obrazkach rzeczywiście dopasowanie było niezłe, to jednak wiadomo, że w próbie zawsze może być zbyt przypadkowe. 

Drugi problem to założenie rozkładu t-studenta w resztach modelu. Właściwie dlaczego założyłem taki? Stopy zwrotu są zazwyczaj dalekie od normalności, a rozkład Studenta uwzględnia grubsze ogony. Proces MA można zapisać jako AR wielu rzędów z przeszłości (zob. ten wpis), a skoro AR to właśnie stopa zwrotu, tylko opóźniona ze stałymi współczynnikiami, to znaczy, że proces MA ma ten sam rozkład - zakładając stacjonarność stopy zwrotu. A skoro (stacjonarny) MA ma ten sam rozkład, to znaczy, że składnik losowy też ma ten sam rozkład (bo MA jedynie mnoży jego wartość z przeszłości przez stałą). W sumie więc, jeśli stopa zwrotu ma rozkład Studenta, to reszty w ARMA też muszą taki posiadać.

Tylko że okazało się, że nie można odrzucić hipotezy o gaussowskości ostatnich 80 rocznych stóp zwrotu lub premii:

> shapiro.test(as.numeric(stopa_zwrotu))

	Shapiro-Wilk normality test

data:  as.numeric(stopa_zwrotu)
W = 0.9804, p-value = 0.257

Nawet ostatnie 120 danych (czyli od 1904 r.) nie zmieniło tego wyniku:

	Shapiro-Wilk normality test

data:  as.numeric(stopa_zwrotu)
W = 0.9854, p-value = 0.221

 

Nienormaność pojawia się dopiero przy większych częstościach. Bardzo możliwe, że przy większej liczbie rocznych obserwacji należy odrzucić rozkład normalny, ale dla danych, które analizowałem, nie można tego zrobić.

Po trzecie porównałem autoarfima z auto.arima z pakietu forecast, zrobiłem różne symulacje ARMA i czasami auto.arima wskazywała prawidłowy model, czasami autoarfima. Stąd doszedłem do wniosku, że należy wykonywać testy obiema funkcjami. W dodatku najlepiej stosować zarówno kryterium nie tylko BIC, ale też AIC.

Premia za ryzyko (od samych akcji)

autoarfima:

klaster <- makeCluster(detectCores() - 1)

premia_model_arma <- autoarfima(data = premia, ar.max = 10, ma.max = 10, cluster = klaster, solver = "nlminb")

stopCluster(klaster)

Wszystkie kryteria wskazały ten sam model:

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.077126    0.021203   3.6376 0.000275
sigma  0.189641    0.014992  12.6491 0.000000

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.077126    0.019321   3.9917 0.000066
sigma  0.189641    0.017623  10.7608 0.000000

LogLikelihood : 19.4946 

Information Criteria
------------------------------------
                     
Akaike       -0.43737
Bayes        -0.37781
Shibata      -0.43857
Hannan-Quinn -0.41349

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.7322  0.3922
Lag[2*(p+q)+(p+q)-1][2]    1.5233  0.3557
Lag[4*(p+q)+(p+q)-1][5]    4.0949  0.2425

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.392  0.2381
Lag[2*(p+q)+(p+q)-1][2]     1.749  0.3081
Lag[4*(p+q)+(p+q)-1][5]     2.110  0.5928


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      2.253   2  0.3241
ARCH Lag[5]      2.403   5  0.7911
ARCH Lag[10]     5.979  10  0.8171

Nyblom stability test
------------------------------------
Joint Statistic:  0.2461
Individual Statistics:            
mu    0.1522
sigma 0.0835

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.61 0.749 1.07
Individual Statistic:	 0.35 0.47 0.75

Model jest dobrze dopasowany: stabilny, nie występują autokorelacje w resztach ani ARCH. Czyli nie żadna ARMA, ale zwykła średnia. Średnioroczna premia za ryzyko wynosi 7,7% +/- 19% (nie mylmy tu błędu standardowego reszt, który wynosi 1,9% - on wskazuje jak średnia przypadkowo się odchyla, a nie sama stopa zwrotu).

To samo auto.arima:

premia_model_arma <- auto.arima(premia, max.p=10, max.q=10, stepwise = FALSE,  approximation = FALSE)

Argumenty stepwise i approximation ustawiam na FALSE, bo zwiększają precyzję algorytmu. Niezależnie od kryterium dostaniemy model:

Series: premia 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
       mean
      0.077
s.e.  0.021

sigma^2 = 0.0364:  log likelihood = 19.49
AIC=-34.99   AICc=-34.83   BIC=-30.23

Obydwie funkcje - na obydwu kryteriach BIC i  AIC wskazały ten sam model - zwykłą średnią. To jest silny argument za odrzuceniem poprzedniego modelu ARMA.

Można pokazać to na wykresie:

Premia za ryzyko od akcji i obligacji skarbowych

Obligacje skarbowe są obarczone głównie ryzykiem stopy procentowej. Aby je uwzględnić metodą ad hoc, odejmujemy od stopy zwrotu z akcji rentowność bonów skarbowych. Nie jest to prawdziwa premia, bo nie uwzględnia zależności między akcjami a obligacjami, ale daje pewien obraz. 

Test na normalność:

	Shapiro-Wilk normality test

data:  premia
W = 0.9845, p-value = 0.447

Czyli też zachowuje się normalnie.

Analogicznie jak poprzednio używamy autoarfima i auto.arima.

autoarfima:

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.089739    0.019332   4.6421 0.000003
sigma  0.172907    0.013670  12.6491 0.000000

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.089739    0.017760   5.0529        0
sigma  0.172907    0.012238  14.1292        0

LogLikelihood : 26.885 

Information Criteria
------------------------------------
                     
Akaike       -0.62212
Bayes        -0.56257
Shibata      -0.62333
Hannan-Quinn -0.59825

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4299  0.5121
Lag[2*(p+q)+(p+q)-1][2]    1.4963  0.3619
Lag[4*(p+q)+(p+q)-1][5]    4.3507  0.2134

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.9464  0.3306
Lag[2*(p+q)+(p+q)-1][2]    3.1373  0.1288
Lag[4*(p+q)+(p+q)-1][5]    4.8670  0.1640


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      5.497   2 0.06403
ARCH Lag[5]      5.442   5 0.36440
ARCH Lag[10]    10.107  10 0.43112

Nyblom stability test
------------------------------------
Joint Statistic:  0.1892
Individual Statistics:            
mu    0.1121
sigma 0.0580

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.61 0.749 1.07
Individual Statistic:	 0.35 0.47 0.75

auto.arima:

Series: premia 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
       mean
      0.090
s.e.  0.019

sigma^2 = 0.0303:  log likelihood = 26.88
AIC=-49.77   AICc=-49.61   BIC=-45.01

Ponownie obie funkcje wskazały ten sam model zwykłej średniej: "pełna premia" wynosi 9% +/- 17,3%. taki wynik uzyskałem niezależnie od użytego kryterium. 

I wykres:

Chociaż brak zależności w rocznych stopach pogarsza sytuację prognostyczną, to z punktu widzenia inwestora może się ona nawet poprawić. Po pierwsze ułatwia mu wycenę - nie musi uwzględniać zmienności stopy dyskontowej. Po drugie zwiększa szansę, że wycena będzie szybko zbiegać do jego wyceny, a nie od niej uciekać. Na ten moment, nawet bez wyceny, możemy spekulować, że premia w 2025 spadnie, choć powinna pozostać dodania.