Jak wykorzystać model GARCH w trejdingu?

Model GARCH jako taki nie służy do prognozowania cen czy stóp zwrotu, ale ich zmienności. Można go jednak rozumieć w szerszy sposób, tj. ARIMA-GARCH. Ten artykuł jest kontynuacją wcześniejszych: Dlaczego oddzielne badanie autokorelacji jest błędem? oraz Modelowanie ryzyka z rugarch (język R). Omówię te dwa sposoby z punktu widzenia trejdingu. 

Przygotowanie danych

Na początek ładuję dwa pakiety: xts oraz rugarch

# Przedwstęp

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

  install.packages("xts")

}

library("xts")

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

  install.packages("rugarch")

}

library("rugarch")

Analizę wykonam na WIG20, plik .csv - dzienne dane - pobrane ze stooq.pl. Przekształcam je w tygodniowe i ucinam do 200 ostatnich notowań (ostatnie to 13.09.2024). Przekształcam w stopy zwrotu.

# Wstęp

nazwaPliku = "wig20_d.csv"

Obs <- 200

plik <- read.csv(nazwaPliku)

plik$Data <- as.Date(plik$Data)

tabela <- tail(to.weekly(as.xts(plik), indexAt="endof", OHLC=FALSE), Obs)

cena <- tabela$Zamkniecie

stopaZwrotu <- diff.xts(log(cena))

tabela$stopaZwrotu <- stopaZwrotu

tabela <- na.omit(tabela)

stopaZwrotu <- tabela$stopaZwrotu

1. ARIMA-GARCH

W szerszym ujęciu GARCH dzieli się na dwie części: najpierw szukamy optymalnego ARIMA za pomocą funkcji autoarfima(). Ustawiam max 10 rzędów zarówno dla części AR, jak i MA.  Na początek zakładam rozkład skośno-normalny głównie z powodu szybkości obliczania jego parametrów:  

# 1. część ARMA

arimaDopas <- autoarfima(stopaZwrotu, ar.max = 10, ma.max = 10, criterion = "HQIC",

                        distribution.model = "snorm", 

                        method="partial")

                       

Po wpisaniu zmiennej dostajemy statystyki

arimaDopas

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(5,0,3)
Distribution	: snorm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error    t value Pr(>|t|)
mu     0.002644    0.000001    2376.47        0
ar1   -0.616537    0.000031  -19828.09        0
ar2    1.054657    0.000034   31037.94        0
ar3    0.533467    0.000028   19356.08        0
ar4   -0.173088    0.000010  -17629.16        0
ar5    0.053956    0.000006    8470.41        0
ma1    1.038622    0.000009  118359.33        0
ma2   -1.083726    0.000008 -139217.06        0
ma3   -1.131477    0.000060  -18831.84        0
sigma  0.024462    0.000096     254.15        0
skew   0.833062    0.002360     352.95        0

Robust Standard Errors:
       Estimate  Std. Error    t value Pr(>|t|)
mu     0.002644    0.000000   6484.164        0
ar1   -0.616537    0.000119  -5183.800        0
ar2    1.054657    0.000090  11723.967        0
ar3    0.533467    0.000040  13248.164        0
ar4   -0.173088    0.000041  -4179.572        0
ar5    0.053956    0.000019   2915.708        0
ma1    1.038622    0.000358   2897.527        0
ma2   -1.083726    0.000012 -88962.282        0
ma3   -1.131477    0.000060 -18854.153        0
sigma  0.024462    0.000949     25.785        0
skew   0.833062    0.003603    231.216        0

LogLikelihood : 456.9 

Information Criteria
------------------------------------
                    
Akaike       -4.4819
Bayes        -4.2999
Shibata      -4.4876
Hannan-Quinn -4.4082

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                       2.763 0.09649
Lag[2*(p+q)+(p+q)-1][23]    11.182 0.91984
Lag[4*(p+q)+(p+q)-1][39]    17.964 0.71904

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      2.091  0.1482
Lag[2*(p+q)+(p+q)-1][2]     2.238  0.2262
Lag[4*(p+q)+(p+q)-1][5]     4.670  0.1816


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      2.347   2  0.3093
ARCH Lag[5]      9.088   5  0.1056
ARCH Lag[10]    13.467  10  0.1987

Nyblom stability test
------------------------------------
Joint Statistic:  7.683
Individual Statistics:             
mu    0.00450
ar1   0.21211
ar2   0.20711
ar3   0.19583
ar4   0.09643
ar5   0.12618
ma1   0.19549
ma2   0.19475
ma3   0.16916
sigma 0.27608
skew  0.10861

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 2.49 2.75 3.27
Individual Statistic:	 0.35 0.47 0.75

Najlepszy to ARMA(5, 3). Zapisujemy optymalne rzędy:

rzadAR <- arimaDopas$rank.matrix[1,1]

rzadMA <- arimaDopas$rank.matrix[1,2]

Model jest niestabilny (Joint Statistics > 2,75). Prawdopodobnie wybraliśmy zły rozkład. Zmieńmy na skośny uogólniony rozkład błędów, tj. "sged". Ponieważ wybór innego rozkładu niż normalny i skośno-normalny wydłuża znacząco cały proces obliczeń, dodamy parametr cluster, który wykorzystuje pakiet parallel do równoczesnego przetwarzania danych przez pozostałe rdzenie procesora (na koniec wyłączamy tę opcję); to przyspieszy proces: 

klaster <- makeCluster(detectCores() - 1)

arimaDopas <- autoarfima(stopaZwrotu, ar.max = 10, ma.max = 10, criterion = "HQIC",

                        distribution.model = "sged", cluster = klaster,

                        method="partial")

stopCluster(klaster)

arimaDopas

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(5,0,3)
Distribution	: sged 

Optimal Parameters
------------------------------------
       Estimate  Std. Error      t value Pr(>|t|)
mu     0.002994    0.000001    2333.1750        0
ar1   -0.735343    0.000231   -3185.0574        0
ar2    1.065479    0.000335    3185.0574        0
ar3    0.715202    0.000225    3185.0574        0
ar4   -0.164068    0.000052   -3185.0574        0
ar5   -0.012163    0.000005   -2435.7292        0
ma1    1.043899    0.000029   36290.9045        0
ma2   -1.073710    0.000050  -21428.3317        0
ma3   -1.130274    0.000008 -139344.0095        0
sigma  0.024216    0.001212      19.9784        0
skew   0.872133    0.075043      11.6218        0
shape  1.711239    0.172880       9.8984        0

Robust Standard Errors:
       Estimate  Std. Error     t value Pr(>|t|)
mu     0.002994    0.000005    556.1676 0.000000
ar1   -0.735343    0.021895    -33.5850 0.000000
ar2    1.065479    0.053812     19.8000 0.000000
ar3    0.715202    0.020475     34.9306 0.000000
ar4   -0.164068    0.001227   -133.6667 0.000000
ar5   -0.012163    0.000004  -2972.8839 0.000000
ma1    1.043899    0.000706   1479.4296 0.000000
ma2   -1.073710    0.001064  -1009.1810 0.000000
ma3   -1.130274    0.000055 -20500.1353 0.000000
sigma  0.024216    0.001543     15.6918 0.000000
skew   0.872133    0.112199      7.7731 0.000000
shape  1.711239    0.705025      2.4272 0.015216

LogLikelihood : 459.4 

Information Criteria
------------------------------------
                    
Akaike       -4.4965
Bayes        -4.2979
Shibata      -4.5033
Hannan-Quinn -4.4161

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                    0.002713  0.9585
Lag[2*(p+q)+(p+q)-1][23]  8.977574  1.0000
Lag[4*(p+q)+(p+q)-1][39] 16.731900  0.8385

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.5398  0.4625
Lag[2*(p+q)+(p+q)-1][2]    0.6712  0.6184
Lag[4*(p+q)+(p+q)-1][5]    3.1718  0.3766


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]     0.8306   2  0.6601
ARCH Lag[5]     7.7465   5  0.1708
ARCH Lag[10]   12.4933  10  0.2534

Nyblom stability test
------------------------------------
Joint Statistic:  1.744
Individual Statistics:             
mu    0.05618
ar1   0.03768
ar2   0.03560
ar3   0.03446
ar4   0.03642
ar5   0.03377
ma1   0.06555
ma2   0.06355
ma3   0.07824
sigma 0.24067
skew  0.25045
shape 0.27590

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 2.69 2.96 3.51
Individual Statistic:	 0.35 0.47 0.75

Rzeczywiście model staje się stabilny, co oznacza, że problemem był rozkład. Nadal najlepszy jest ARMA(5, 3). Zaskakujące, że diagnostyka nie wykrywa ani autokorelacji w kwadratach reszt, ani efektu ARCH.

2. GARCH

Skoro diagnostyka nie wychwyciła występowania ARCH, można by pominąć ten krok i zająć się od razu prognozą. Jednak zastosowanie modelu wariancji pozwoli sprawdzić również efekt Arch-In-Mean. W tym wpisie pokazałem, że dzienne stopy zwrotu mogą być modelowane za pomocą TGARCH-ARCH-In-Mean, a miesięczne za pomocą EGARCH-ARCH-In-Mean. Sugeruje to, że tygodniowe będą czymś pomiędzy. Dla przypomnienia  model ARCH-in-Mean oznacza, że stopa zwrotu jest zależna od wariancji warunkowej. Wówczas jednak rzędy ARMA były pominięte. Co by się stało, gdybyśmy je dodali? Oczywiście możemy to zrobić poprzez dodanie rzędów ARMA(5, 3), ale wcale nie jest pewne, że będzie to model najlepszy. Dlatego warto zrobić własną funkcję szukającą najlepszego modelu ARMA-GARCH. Funkcja wykorzysta otrzymane z autoarfimy optymalne rzędy jako maksymalne rzędy przeszukiwania, tzn. będziemy sprawdzać rzędy od 0 do rzadAR oraz od 0 do rzadMA. Co więcej będzie uwzględniał zarówno Arch-In-Mean, jak i bez niego. Kod dla TGARCH ma postać:

# 2. część GARCH

kryteriaInfor <- list()

for (k in 0:1) {

  for (i in 0:rzadAR) {

    for (j in 0:rzadMA) {

      mojGarch <- ugarchspec(variance.model = list(model = "fGARCH", submodel = "TGARCH",

                                                   garchOrder = c(1, 1)), 

                             mean.model = list(armaOrder = c(i, j), archm = as.logical(k)), 

                             distribution.model = "sged")

      dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu)

      

      nr <- paste0("ARMA(", i, ", ", j, ") - ARCH In Mean")

      # usuwamy ARCH in mean dla k = 1

      nr <- ifelse(k == 1, nr, sub(" - ARCH In Mean", "", nr))

      kryteriaInfor[[nr]] <- infocriteria(dopasGarch)

    }

  }

}

kryteriaDf <- do.call(cbind, kryteriaInfor)

colnames(kryteriaDf) <- names(kryteriaInfor)

minKryteria <- apply(kryteriaDf, 1, function(wiersz) {

  

  wartoscMin = min(wiersz)

  kolumnaMin = colnames(kryteriaDf)[which.min(wiersz)]

  return(c(Minimum = wartoscMin, Kolumna = kolumnaMin))

})

minKryteriaDf <- as.data.frame(t(minKryteria))

minKryteriaDf

Otrzymany wynik:

                       Minimum    Kolumna
Akaike       -4.45083487203567 ARMA(5, 3)
Bayes         -4.2045670547213 ARMA(3, 3)
Shibata      -4.46117146205014 ARMA(5, 3)
Hannan-Quinn -4.35036611102375 ARMA(5, 3)

 

Oznacza to, że wynik autoarfimy koresponduje z TGARCH. Dla EGARCH kod jest oczywiście taki sam poza fragmentem:

list(model = "fGARCH", submodel = "TGARCH"

który zamienia się w:

list(model = "eGARCH", submodel = NULL 

I z niego otrzymałem:

                       Minimum    Kolumna
Akaike       -4.43263086741263 ARMA(2, 2)
Bayes        -4.25058889217661 ARMA(2, 2)
Shibata      -4.43832575443847 ARMA(2, 2)
Hannan-Quinn  -4.3589537760039 ARMA(2, 2)

Sytuacja okazuje się niejednoznaczna. Z punktu widzenia AIC TGARCH jest lepszy. Z punktu widzenia BIC i HQIC lepszy jest EGARCH. Jeśli chodzi o SHC, czyli kryterium Shibata, to ja to pomijam, po pierwsze niemal zawsze wskazuje ten sam kierunek co AIC, a po drugie jest słabo opisany w literaturze.

Należy zauważyć, że oba modele nie zawierają elementu Arch-In-Mean. Można powiedzieć, że ARMA go zastąpiła.

Otrzymane modele wymagają diagnostyki. Zacznijmy od ARMA(5, 3)-TGARCH(1, 1):

mojGarch <- ugarchspec(variance.model = list(model = "fGARCH", submodel = "TGARCH",

                                             garchOrder = c(1, 1)), 

                       mean.model = list(armaOrder = c(rzadAR, rzadMA)), distribution.model = "sged")

dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu)

dopasGarch

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: fGARCH(1,1)
fGARCH Sub-Model	: TGARCH
Mean Model	: ARFIMA(5,0,3)
Distribution	: sged 

Optimal Parameters
------------------------------------
        Estimate  Std. Error    t value Pr(>|t|)
mu      0.001222    0.000004   287.4520 0.000000
ar1    -0.646631    0.000094 -6851.4518 0.000000
ar2     1.106298    0.000084 13106.6297 0.000000
ar3     0.627234    0.000090  6953.5872 0.000000
ar4    -0.102406    0.000040 -2533.7961 0.000000
ar5     0.099523    0.000039  2578.2384 0.000000
ma1     1.035966    0.000078 13332.0845 0.000000
ma2    -1.079756    0.000171 -6301.9226 0.000000
ma3    -1.128152    0.000256 -4412.9976 0.000000
omega   0.002007    0.001332     1.5071 0.131781
alpha1  0.054355    0.024212     2.2450 0.024768
beta1   0.876710    0.053657    16.3391 0.000000
eta11   0.996087    0.508698     1.9581 0.050217
skew    0.864175    0.062759    13.7697 0.000000
shape   1.601437    0.251582     6.3655 0.000000

Robust Standard Errors:
        Estimate  Std. Error      t value Pr(>|t|)
mu      0.001222    0.000018    66.728331  0.00000
ar1    -0.646631    0.004215  -153.419278  0.00000
ar2     1.106298    0.007908   139.895269  0.00000
ar3     0.627234    0.004506   139.206023  0.00000
ar4    -0.102406    0.000244  -419.427330  0.00000
ar5     0.099523    0.000060  1667.478136  0.00000
ma1     1.035966    0.000476  2177.437996  0.00000
ma2    -1.079756    0.000322 -3357.818162  0.00000
ma3    -1.128152    0.004830  -233.579659  0.00000
omega   0.002007    0.006530     0.307381  0.75855
alpha1  0.054355    0.686095     0.079224  0.93685
beta1   0.876710    0.741582     1.182216  0.23712
eta11   0.996087    8.183445     0.121720  0.90312
skew    0.864175    0.042192    20.481868  0.00000
shape   1.601437    0.680752     2.352454  0.01865

LogLikelihood : 457.9 

Information Criteria
------------------------------------
                    
Akaike       -4.4508
Bayes        -4.2026
Shibata      -4.4612
Hannan-Quinn -4.3504

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                      0.3174  0.5732
Lag[2*(p+q)+(p+q)-1][23]   10.3414  0.9983
Lag[4*(p+q)+(p+q)-1][39]   18.5933  0.6475
d.o.f=8
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                         statistic p-value
Lag[1]                  0.00004328  0.9948
Lag[2*(p+q)+(p+q)-1][5] 3.39042514  0.3404
Lag[4*(p+q)+(p+q)-1][9] 5.51638020  0.3576
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]     4.431 0.500 2.000 0.03530
ARCH Lag[5]     5.177 1.440 1.667 0.09349
ARCH Lag[7]     6.140 2.315 1.543 0.13209

Nyblom stability test
------------------------------------
Joint Statistic:  3.741
Individual Statistics:              
mu     0.01584
ar1    0.01865
ar2    0.02030
ar3    0.01897
ar4    0.02063
ar5    0.01946
ma1    0.01911
ma2    0.01778
ma3    0.01841
omega  0.06211
alpha1 0.04989
beta1  0.08077
eta11  0.31789
skew   0.07723
shape  0.04888

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 3.26 3.54 4.07
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.9408 0.3480    
Negative Sign Bias  0.3877 0.6987    
Positive Sign Bias  0.2705 0.7871    
Joint Effect        1.0407 0.7914    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     18.99       0.4575
2    30     34.82       0.2105
3    40     44.82       0.2409
4    50     40.95       0.7865

Model jest lekko niestabilny (przez to parametry (G)ARCH nieistotne). Teraz ARMA(2, 2)-EGARCH(1,1).

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: eGARCH(1,1)
Mean Model	: ARFIMA(2,0,2)
Distribution	: sged 

Optimal Parameters
------------------------------------
        Estimate  Std. Error     t value Pr(>|t|)
mu      0.000872    0.000007  1.3068e+02  0.0e+00
ar1    -1.695986    0.000000 -1.6489e+07  0.0e+00
ar2    -0.735113    0.000000 -9.7011e+06  0.0e+00
ma1     1.874230    0.000648  2.8906e+03  0.0e+00
ma2     0.868271    0.001679  5.1706e+02  0.0e+00
omega  -0.384768    0.000038 -1.0045e+04  0.0e+00
alpha1 -0.189127    0.046551 -4.0628e+00  4.8e-05
beta1   0.947281    0.001075  8.8104e+02  0.0e+00
gamma1 -0.075395    0.000542 -1.3914e+02  0.0e+00
skew    0.836575    0.112010  7.4688e+00  0.0e+00
shape   2.646523    0.666196  3.9726e+00  7.1e-05

Robust Standard Errors:
        Estimate  Std. Error     t value Pr(>|t|)
mu      0.000872    0.000542  1.6082e+00 0.107795
ar1    -1.695986    0.000013 -1.3200e+05 0.000000
ar2    -0.735113    0.000008 -8.7846e+04 0.000000
ma1     1.874230    0.091483  2.0487e+01 0.000000
ma2     0.868271    0.315734  2.7500e+00 0.005959
omega  -0.384768    0.107277 -3.5867e+00 0.000335
alpha1 -0.189127    4.249978 -4.4501e-02 0.964505
beta1   0.947281    0.106032  8.9339e+00 0.000000
gamma1 -0.075395    0.025790 -2.9234e+00 0.003462
skew    0.836575    0.494241  1.6926e+00 0.090523
shape   2.646523   61.861149  4.2782e-02 0.965876

LogLikelihood : 452.0468 

Information Criteria
------------------------------------
                    
Akaike       -4.4326
Bayes        -4.2506
Shibata      -4.4383
Hannan-Quinn -4.3590

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                       1.428  0.2320
Lag[2*(p+q)+(p+q)-1][11]     3.101  1.0000
Lag[4*(p+q)+(p+q)-1][19]     9.156  0.6119
d.o.f=4
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.1386  0.7097
Lag[2*(p+q)+(p+q)-1][5]    0.8697  0.8884
Lag[4*(p+q)+(p+q)-1][9]    2.8003  0.7914
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.2651 0.500 2.000  0.6066
ARCH Lag[5]    0.8789 1.440 1.667  0.7692
ARCH Lag[7]    2.6093 2.315 1.543  0.5901

Nyblom stability test
------------------------------------
Joint Statistic:  2.4751
Individual Statistics:              
mu     0.03485
ar1    0.03350
ar2    0.03469
ma1    0.04462
ma2    0.03174
omega  0.15840
alpha1 0.14496
beta1  0.16233
gamma1 0.20550
skew   0.09238
shape  0.09576

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 2.49 2.75 3.27
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.6633 0.5080    
Negative Sign Bias  0.1653 0.8689    
Positive Sign Bias  0.3918 0.6956    
Joint Effect        0.5216 0.9141    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     15.57       0.6856
2    30     39.04       0.1009
3    40     39.99       0.4258
4    50     54.52       0.2728

EGARCH jest dużo bardziej stabilny od TGARCH, ale efekt GARCH tylko częściowo występuje (alpha1 nieistotna, beta1 istotna). W tym miejscu powinniśmy zadać sobie pytanie, czy w ogóle nie należy pominąć GARCH i nie wybrać samego modelu ARMA(5, 3), który został uznany za najlepszy przez autoarfimę, a jej diagnostyka w ogóle nie wykryła autokorelacji w wariancjach składnika losowego. Dla ARMA(5, 3) dostaliśmy wartości kryteriów

Akaike       -4.4965
Bayes        -4.2979
Shibata      -4.5033
Hannan-Quinn -4.4161

Dla ARMA(5, 3)-TGARCH(1,1) dostaliśmy

Akaike       -4.4508
Bayes        -4.2026
Shibata      -4.4612
Hannan-Quinn -4.3504

Dla ARMA(2, 2)-EGARCH(1, 1):

Akaike       -4.4326
Bayes        -4.2506
Shibata      -4.4383
Hannan-Quinn -4.3590

Wydawałoby się, że należy po prostu zastosować zwykły ARMA(5, 3). Jednakże takie porównanie 1:1 nie wydaje się rozsądne. W przypadku modelu ARIMA prognozujemy jedynie stopy zwrotu, podczas gdy w GARCH (czy precyzyjniej ARIMA-GARCH) nie tylko stopy zwrotu, ale i odchylenie od średniej.

W sumie należy oddzielić dokładność prognozy stopy zwrotu od prognozy jej odchylenia. Zwykła korelacja nie nadaje się, bo skupia się na kierunku zmian, a nie ich znaku. Gdyby nasza prognoza miała zawsze taki sam znak jak prawdziwa stopa zwrotu, ale różny kierunek (tzn. gdy prawdziwa zmiana rośnie, to prognoza spada), to dostalibyśmy ujemną korelację, a więc zupełnie nie ma to sensu. Dlatego potrzebna jest miara badająca bardziej wspólność znaków niż kierunku. Oczywiście najbardziej naturalnym podejściem jest obliczenie częstości poprawności prognozy i porównanie jej z teoretycznym prawdopodobieństwem co najmniej tylu sukcesów w rozkładzie dwumianowym - tak jak pokazałem w tym artykule. Test dwumianowy wskaże czy liczba sukcesów jest istotna statystycznie, ale nie poda już  faktycznej skuteczności - tj. prawdopodobieństwa po uwzględnieniu przypadkowych odchyleń. 

Taką miarą jest zaś tzw. dokładność kierunkowa. rugarch oferuje specjalny test DAC (ang. Directional Accuracy Test - test kierunkowej dokładności), który ją oblicza i testuje jej istotność. Zastosowanie DAC jest proste. jego konstrukcja jest następująca:

DACTest(forecast, actual, test="PT")

Symbol "PT" oznacza test Pesaran-Timmermann. Zamiast tego można wpisać "AG" - test Anatolyev-Gerko. Oba dają bardzo zbliżone rezultaty.

Znajdujemy więc prognozy w każdym modelu. Dla ARMA(5, 3):

# prognoza w próbie

stopaZwrotuDopas <- fitted(dopasGarch)

# Testujemy zdolność prognostyczną modelu:

DACTest(forecast = stopaZwrotuDopas, actual = as.numeric(stopaZwrotu), test="PT")

$Test
[1] "Pesaran and Timmermann"

$Stat
[1] 5.390263

$p.value
[1] 3.51773e-08

$H0
[1] "Independently Distributed"

$Decision
[1] "Reject  H0"

$DirAcc
[1] 0.6884422

Prognoza zdaje test na 68,8% i jest istotna.

Następnie sprawdzamy to samo dla ARIMA(5, 3)-TGARCH(1, 1) oraz ARIMA(2, 2)-EGARCH(1, 1). I tak test DAC wskazał dla TGARCH dokładność kierunkową = 0.693, a dla EGARCH 0.658. Wynika z tego, że TGARCH jest najlepszy spośród wszystkich trzech modeli. Dostaliśmy dowód, że dodanie do ARMA modelu GARCH może poprawić prognozę samej części ARMA, choć mogłoby się wydawać, że są to części niezależne.

Przed utworzeniem wykresu z prognozą, warto jeszcze spojrzeć na dopasowanie reszt TGARCH do rozkładu SGED:

plot(dopasGarch, which=8)

Wygląda w miarę dobrze.

Możemy robić prognozę: dwa tygodnie do przodu, pokażemy 30 ostatnich obserwacji na wykresie:

prognoza_plus_wykres <- function(stopaZwrotu, dopasGarch, ileDoPrzodu, Ost) {

  # warunkowe odchylenie std w próbie

  ryzyko <- sigma(dopasGarch)

  

  # prognoza poza próbą

  prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = ileDoPrzodu)

  prognozaStopaZwrotu <- fitted(prognozaGarch)

  prognozaRyzyko <- sigma(prognozaGarch)

  

  # zamiana na ts i łączenie z prognozą dla wykresu

  stopaZwrotuWykres <- tail(as.numeric(stopaZwrotu), Ost)

  stopaZwrotuDopasWykres <- tail(as.numeric(stopaZwrotuDopas), Ost)

  ryzykoWykres <- tail(as.numeric(ryzyko), Ost)

  stopaZwrotuWykres <- c(stopaZwrotuWykres, rep(NA, ileDoPrzodu))

  stopaZwrotuDopasWykres <- c(stopaZwrotuDopasWykres, rep(NA, ileDoPrzodu))

  prognozaStopaZwrotuWykres <- c(rep(NA, Ost-1), stopaZwrotuWykres[Ost], 

                                 as.numeric(prognozaStopaZwrotu))

  

  # Wyłuskaj nazwę modelu wariancji

  modelWariancji <- dopasGarch@model$modeldesc$vmodel

  podmodelWariancji <- dopasGarch@model$modeldesc$vsubmodel

  modelWariancji <- ifelse(!is.null(podmodelWariancji), podmodelWariancji, modelWariancji)

  # Wyłuskaj nazwę modelu średniej

  modelSredniej <- paste0("ARMA(", dopasGarch@model$modelinc["ar"], ",", dopasGarch@model$modelinc["ma"], ")")

  # zmienne do wykresu

  ryzykoWykres <- c(ryzykoWykres, rep(NA, ileDoPrzodu))

  prognozaRyzykoWykres <- c(rep(NA, Ost-1), ryzykoWykres[Ost], as.numeric(prognozaRyzyko))

  okres <- index(stopaZwrotu)

  okresPrognoza <- seq.Date(from=max(okres), by="week", length.out=ileDoPrzodu+1)[-1]

  okresWykres <- c(tail(okres, Ost), okresPrognoza)

  # wykres

  zakresX <- c(min(okresWykres), max(okresWykres))

  zakresY <- c(min(na.omit(-ryzykoWykres))*3, max(na.omit(ryzykoWykres))*3)

  plot(x=okresWykres, y=stopaZwrotuWykres, main="stopa zwrotu WIG20", xaxt="n", xlim=zakresX, ylim=zakresY, lwd=2, type="l", xlab="")

  axis(side=1, at=okresWykres, labels=okresWykres, las=2, cex.axis=0.7)

  lines(x=okresWykres, y=stopaZwrotuDopasWykres, col="lightblue", lwd=2)

  lines(x=okresWykres, y=prognozaStopaZwrotuWykres, lwd=2, col="green")

  lines(x=okresWykres, y=ryzykoWykres, lwd=2, col="red")

  lines(x=okresWykres, y=-ryzykoWykres, lwd=2, col="red")

  lines(x=okresWykres, y=prognozaRyzykoWykres, lwd=2, col="brown")

  lines(x=okresWykres, y=-prognozaRyzykoWykres, lwd=2, col="brown")

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

  abline(v=okresWykres, col="grey")

  legend(x="topleft", legend=c("stopa zwrotu", 

                               paste("modelowana stopa zwrotu -",  modelSredniej), 

                               paste("prognoza stopy zwrotu -",  modelSredniej),

                               paste("modelowane ryzyko -",  modelWariancji), 

                               paste("prognoza ryzyka -",  modelWariancji)), 

         bg = "transparent", cex = 0.8, 

         lwd=c(2,2,2,2,2), col=c("black",                                                                                                        "lightblue",

                                 "green",

                                 "red",

                                 "brown"))

  

  legend("bottomleft", legend=paste("korelacja stopa zwrotu / prognoza stopy zwrotu:", korelacja), 

         bty="o", bg="white", box.lwd=1, text.font=2)

}

prognoza_plus_wykres(stopaZwrotu, dopasGarch, ileDoPrzodu = 2, Ost = 30)

Tak więc dostaliśmy prognozę dodatnią na ten tydzień i ujemną na następny.

Pozostaje kwestia stabilności. Przypomnę, że ARMA(5,3)-TGARCH(1,1) jest niestabilny, ale sam ARMA(5, 3) był stabilny. Logiczne więc, że cała niestabilność bierze się z części modelu wariancji. To dobra wiadomość dla prognozowania stopy zwrotu. Z kolei ARMA(2, 2)-EGARCH(1, 1) jest stabilny. Wobec tego możemy użyć tego ostatniego do prognozy wariancji. Dla niego dostałem:

Zauważmy, że prognozowane odchylenia (ryzyko) lepiej się tu dopasowują do faktycznych niż przy zastosowaniu TGARCH. Jednak końcówka i prognoza odchyleń - w zasadzie są identyczne. Natomiast nie sugerowałbym się już prognozą stóp zwrotu, chociaż w tym wypadku kierunek jest ten sam.

W ogólnym przypadku zmienność, a więc i sam GARCH (bez ARMA) można wykorzystać w trejdingu na różne sposoby. Jeżeli występuje efekt Arch-In-Mean, to wzrost ryzyka będzie powiązany z oczekiwaną stopą zwrotu, co pokazałem w tym art. Jeżeli Arch-In-Mean nie występuje, tak jak tutaj, to nadal możemy wykorzystać zmienność ryzyka. W przypadku EGARCH rosnące ryzyko oznacza, że w kolejnych okresach nadal będzie ono duże, a co za tym idzie możemy spodziewać się większych zmian. W takim razie, przykładowo, czekamy na odpowiednio duży spadek, aby w tym momencie wejść na rynek. Ze względu na większą zmienność, rośnie szansa, że nastąpi odbicie. Sprzedawać można, gdy stopa zbliża się do krawędzi prognozy odchylenia. Obecnie ryzyko raczej pozostaje podwyższone. W połączeniu z dodatnią prognozą stopy spodziewałbym się większego impulsu w górę, a za tydzień w dół - oczywiście na ten moment. W przypadku spadającego ryzyka można dłużej utrzymywać stratną pozycję i czekać na zajęcie pozycji przeciwnej - gdy stopa przekroczy odchylenie - w dół - pozycja długa, w górę - pozycja krótka, mówiąc umownie.