Посматрамо податке са берзе, односно дневене вредности S&P500IndexaS\&P500Indexa који се рачуна користећи цене на берзи 500 великих корпорација. Временска серија у RR-у враћа 100ln(SPIt/SPIt1)100\ln(SPI_t/SPI_{t−1}), где је SPItSPI_t вредност S&P500IndexaS\&P500Indexa за дан tt. То је у ствари серија стопа приноса вредности индекса.

library(MASS)
data(SP500)
plot(SP500, type = 'l')

Са графика серије уочавамо да је дисперзија највећа у последњој трећини, а најмања у другој трећини серије.

acf(SP500)

Kада је дисперзија серије константна, тада је серија хомоскедастична. Уколико дисперзија није константна, већ је нека функција од времена tt, тада је серија хетероскедастична. У нашој серији S&P500S\&P500 примећујемо периоде повећане или смањене дисперзије, тзв. волатилности, односно дисперзија је корелисана у времену (зависи од претходних вредности дисперзије). Тада серију називамо условно хетероскедастичном.

Наведену корелисаност дисперзије можемо уочити на корелограму квадриране серије, то ће бити еквивалентно дисперзији, уколико је очекивање серије 0 (у супротном, најпре центрирамо серију око 0).

acf((SP500-mean(SP500))^2)

ARCH модели

За моделовање наведених временских серија користимо ARCH(p)ARCH(p) моделе, тј. Аutoregressive conditional heteroskedacity моделе са pp параметара, дефинисане на следећи начин: Xt=σtet,σt2=α0+i=1pαiXti2,X_t =\sigma_t\cdot e_t, \quad \sigma_t^2 = \alpha_0+\sum_{i=1}^p \alpha_i X_{t-i}^2, где је α0>0,α1,,αp0\alpha_0>0, \alpha_1, \dots, \alpha_p \ge 0, и etN(0,1)e_t \sim \mathcal{N}(0, 1) је Гаусов бели шум.

Посматрамо ARCH(1)ARCH(1) модел: Xt=σtet,σt=α0+α1Xt12.X_t=\sigma_t\cdot e_t, \quad \sigma_t = \alpha_0+\alpha_1 X_{t-1}^2. Очекивање и дисперзија су: E(Xt)=0,D(Xt)=α0+α1D(Xt1).E(X_t) = 0, \quad D(X_t) = \alpha_0 + \alpha_1D(X_{t-1}).

Уколико је серија стационарна знамо да важи D(Xt)=D(Xt1)=constD(X_t) = D(X_{t−1}) = const, тада је D(Xt)=α01α1D(X_t) = \frac{\alpha_0}{1-\alpha_1}. Пошто је дисперзија ненегативна намеће се додатни услов 0α1<10 \le \alpha_1 <1.

GARCH модели

GARCH(p,q)GARCH(p, q) модели (generalized ARCH) су уопштење ARCHARCH модела облика: Xt=σtet,σt2=α0+i=1pαiXti2+j=1qβjσtj2,X_t = \sigma_t \cdot e_t, \quad \sigma_t^2 = \alpha_0 + \sum_{i=1}^p \alpha_i X_{t-i}^2+ \sum_{j=1}^q \beta_j \sigma_{t-j}^2, где је α0>0,α1,,αp,β1,,βq0\alpha_0>0, \alpha_1, \dots, \alpha_p, \beta_1, \dots, \beta_q \ge 0 и etN(0,1)e_t \sim \mathcal{N}(0,1) је Гаусов бели шум.

Својства GARCHGARCH моделоване серије аналогна су својствима ARCHARCH моделоване серије. Квадрирана серија овог модела понаша се као ARMA(p,q)ARMA(p, q) модел.

Посматрајмо модел GARCH(1,1)GARCH(1, 1). Слично као и код ARCHARCH модела, добија се E(Xt)=0E(X_t) = 0 и D(Xt)=E(Xt2)=E(σt2)D(X_t) = E(X_t^2) = E(\sigma_t^2).

Услов стационарности модела задаје услов D(Xt)=const.D(X_t) = const., тј. E(Xt2)=const.E(X^2_t) = const., па је Е(Xt2)=α01α1β1Е(X_t^2) =\frac{\alpha_0}{1-\alpha_1-\beta_1}, што даје услов стационарности GARCHGARCH модела 0α1+β1<10 \le \alpha_1+\beta_1<1.

ARCHARCH и GARCHGARCH модели налазе велику примену у очекивању стопа приноса финансијских временских серија. На берзи у бурним периодима, велики скокови акција су често праћени великим скоковима дисперзијеи обрнуто. Ово указује на условну хетероскедастичност временске серије приноса.

Погодност GARCHGARCH модела утврђује се посматрањем серије резидуала тј. испитивањем да ли је та серија бели шум. Постоје и формални статистички тестови за проверавање корелисаности резидуала, као што је Ljung-Boks тест.

Како дефиниција GARCHGARCH модела даје прогнозу за један корак унапред, јер се XtX_t рачуна помоћу Xt1X_{t−1}. Изражавањем Xt1X_{t−1} преко ранијих вредности може се доћи до формуле за прогнозу за два корака унапред, итд.

Врло често се GARCHGARCH модели комбинују са ARMAARMA или ARIMAARIMA моделом, којима моделујемо саму серују, а GARCHGARCH моделом серију резидуала. У том случају, GARCHGARCH модел нема утицај на прогнозу саме серије у неком будућем тренутку, али може утицати на дисперзију прогнозираних вредности.

Пример 1

library(tseries)
## Warning: package 'tseries' was built under R version 4.0.4
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(astsa)
plot(nyse)

Посматрамо временску серију nyse из пакета astsa која представља стопе приноса акција на берзи у Њујорку у периоду од 2000 радних дана.

Видимо да се волатилност знатно мења око 1000. дана, па зато користимо GARCH(1,1)GARCH(1,1) модел.

model <- garch(nyse)
## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     8.753541e-05     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -8.362e+03
##      1    7 -8.377e+03  1.78e-03  2.29e-03  1.0e-04  1.9e+11  1.0e-05  2.20e+08
##      2    8 -8.380e+03  3.10e-04  6.10e-04  1.0e-04  2.0e+00  1.0e-05  1.30e+02
##      3    9 -8.381e+03  1.10e-04  1.39e-04  9.1e-05  2.0e+00  1.0e-05  1.21e+02
##      4   10 -8.381e+03  5.41e-06  5.07e-06  9.8e-05  2.0e+00  1.0e-05  1.22e+02
##      5   17 -8.421e+03  4.79e-03  6.96e-03  2.9e-01  2.0e+00  4.1e-02  1.22e+02
##      6   20 -8.473e+03  6.08e-03  9.77e-03  7.0e-01  2.0e+00  2.8e-01  5.42e+00
##      7   29 -8.478e+03  5.88e-04  1.43e-03  8.0e-06  6.4e+00  5.4e-06  7.76e-02
##      8   30 -8.478e+03  2.00e-05  1.62e-05  7.9e-06  2.0e+00  5.4e-06  1.96e-02
##      9   31 -8.478e+03  9.60e-07  8.80e-07  8.0e-06  2.0e+00  5.4e-06  2.03e-02
##     10   41 -8.529e+03  5.94e-03  6.45e-03  3.4e-01  5.7e-01  3.5e-01  2.01e-02
##     11   43 -8.530e+03  2.15e-04  1.86e-03  2.5e-02  1.9e+00  3.5e-02  1.05e-02
##     12   46 -8.547e+03  1.94e-03  1.69e-03  5.5e-02  8.7e-01  8.4e-02  5.98e-03
##     13   57 -8.547e+03  2.34e-05  3.88e-05  1.5e-07  4.8e+00  2.4e-07  7.69e-04
##     14   66 -8.548e+03  4.57e-05  6.60e-05  7.6e-03  1.7e+00  1.2e-02  6.74e-04
##     15   67 -8.548e+03  1.39e-05  9.92e-06  1.3e-03  0.0e+00  2.6e-03  9.92e-06
##     16   69 -8.548e+03  1.10e-05  7.88e-06  3.0e-03  0.0e+00  6.3e-03  7.88e-06
##     17   70 -8.548e+03  3.51e-06  3.43e-06  3.4e-03  1.0e-01  6.3e-03  3.45e-06
##     18   71 -8.548e+03  1.15e-07  1.88e-07  9.9e-04  0.0e+00  1.7e-03  1.88e-07
##     19   74 -8.548e+03  6.59e-10  5.05e-10  5.6e-06  1.8e+00  9.2e-06  2.16e-09
##     20   77 -8.548e+03  7.43e-11  5.21e-11  9.7e-07  1.9e+00  1.6e-06  1.86e-09
##     21   80 -8.548e+03  8.60e-14  4.93e-12  1.0e-07  2.0e+00  1.7e-07  1.81e-09
##     22   83 -8.548e+03  3.90e-13  2.79e-13  5.8e-09  2.0e+00  9.5e-09  1.83e-09
##     23   90 -8.548e+03 -3.04e-14  1.33e-17  1.6e-14  2.6e+01  2.5e-14  1.82e-09
## 
##  ***** FALSE CONVERGENCE *****
## 
##  FUNCTION    -8.547869e+03   RELDX        1.550e-14
##  FUNC. EVALS      90         GRAD. EVALS      23
##  PRELDF       1.330e-17      NPRELDF      1.820e-09
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    6.552055e-06     1.000e+00    -4.565e+00
##      2    1.117548e-01     1.000e+00     4.269e-01
##      3    8.086265e-01     1.000e+00     3.083e-01
summary(model)
## 
## Call:
## garch(x = nyse)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.66460 -0.45228  0.08785  0.59918  4.07508 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 6.552e-06   6.761e-07    9.691   <2e-16 ***
## a1 1.118e-01   4.056e-03   27.554   <2e-16 ***
## b1 8.086e-01   1.292e-02   62.566   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 3983.9, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 1.5874, df = 1, p-value = 0.2077

Жарк-Беров тест проверава да ли су вредности временске серије нормално расподељене на основу вредности коефицијената спљоштености и симетричности. Како је мала pp-вредност мала, одбацујемо хипотезу. Љунг-Боксов тест проверава корелисаност елемената временске серије. Прихватамо хипотезу о некорелисаности временске серије.

Погледајмо како изгледа прогноза GARCHGARCH модела.

pred<- predict(model)
plot(800:1000, nyse[800:1000], type="l", xlab="Time",ylab="NYSE Returns")
lines(pred[,1], col="blue", lty="dashed") # Primetimo da povećanje ili smanjenje vrednosti disperzije kasni za jedan dan
lines(pred[,2], col="blue", lty="dashed")

Пример 2

Посматрамо временску серију која садржи температуру на јужној хемисфери и периоду од 1850-2007. године.

www <- "stemp.txt"
stemp <- scan(www)
stemp.ts<-ts(stemp,start=1850,frequency=12)
plot(stemp.ts)

Моделоваћемо ARIMAARIMA моделом. Правимо функцију која бира најбољи ARIMAARIMA модел

get.best.arima <- function(x.ts, maxord = c(1,1,1,1,1,1))
{
  best.aic <- 1e8
  n <- length(x.ts)
  for (p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
    for (P in 0:maxord[4]) for(D in 0:maxord[5]) for(Q in 0:maxord[6])
    {
      fit <- arima(x.ts, order = c(p,d,q),
                   seas = list(order = c(P,D,Q),
                               frequency(x.ts)), method = "CSS")
      fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef)
      if (fit.aic < best.aic)
      {best.aic <- fit.aic
      best.fit <- fit
      best.model <- c(p,d,q,P,D,Q)
      }
    }
  list(best.aic, best.fit, best.model)
}

stemp.best <- get.best.arima(stemp.ts, maxord = rep(2,6))
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1

## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1

## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1

## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1

## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1

## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
stemp.best[[3]]
## [1] 1 1 2 2 0 1
stemp.arima <- arima(stemp.ts, order = c(1,1,2),
                                    seas = list(order = c(1,0,1), 12))
stemp.res <- resid(stemp.arima)
plot(stemp.res)

acf(stemp.res) # korelogram serije reziduala 

acf(stemp.res^2) # korelogram serije kvadriranih reziduala

Видимо да су квадрирани резидуали позитивно корелисани за већину корака, што значи да постоји волатилност у серији.

stemp.garch <- garch(stemp.res, trace = F)
t(confint(stemp.garch))
##                  a0         a1        b1
## 2.5 %  1.064206e-05 0.03299211 0.9249312
## 97.5 % 1.485815e-04 0.06525721 0.9630785
stemp.garch.res <- resid(stemp.garch)[-1] #-1 je uklonio prvu vrednost u seriji reziduala jer je NA zbog sigma=0
acf(stemp.garch.res)

acf(stemp.garch.res^2)

Пример 3

Посматраћемо временску серију Daily Closing Prices of Major European Stock Indices, 1991–1998. Фитоваћемо у GARCH(1,1)GARCH(1,1) модел.

data(EuStockMarkets)  
dax <- diff(log(EuStockMarkets))[,"DAX"]
dax.garch <- garch(dax)
## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     9.549651e-05     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -7.584e+03
##      1    8 -7.585e+03  1.45e-05  2.60e-05  1.4e-05  1.0e+11  1.4e-06  1.35e+06
##      2    9 -7.585e+03  1.88e-07  1.97e-07  1.3e-05  2.0e+00  1.4e-06  1.50e+00
##      3   18 -7.589e+03  6.22e-04  1.10e-03  3.5e-01  2.0e+00  5.5e-02  1.50e+00
##      4   21 -7.601e+03  1.58e-03  1.81e-03  6.2e-01  1.9e+00  2.2e-01  3.07e-01
##      5   23 -7.634e+03  4.22e-03  3.55e-03  4.3e-01  9.6e-01  4.4e-01  3.06e-02
##      6   25 -7.646e+03  1.61e-03  1.85e-03  2.9e-02  2.0e+00  4.4e-02  5.43e-02
##      7   27 -7.646e+03  3.82e-05  5.23e-04  1.3e-02  2.0e+00  2.0e-02  1.46e-02
##      8   28 -7.648e+03  1.86e-04  1.46e-04  6.5e-03  2.0e+00  9.9e-03  1.54e-03
##      9   29 -7.648e+03  3.12e-05  4.83e-05  6.4e-03  2.0e+00  9.9e-03  3.34e-03
##     10   30 -7.648e+03  1.39e-05  6.31e-05  6.2e-03  1.9e+00  9.9e-03  1.86e-03
##     11   31 -7.650e+03  2.70e-04  3.24e-04  6.0e-03  1.9e+00  9.9e-03  4.99e-03
##     12   34 -7.656e+03  8.42e-04  8.57e-04  2.2e-02  1.7e-01  3.9e-02  2.22e-03
##     13   36 -7.661e+03  6.12e-04  6.40e-04  1.9e-02  4.2e-01  3.9e-02  2.09e-03
##     14   38 -7.665e+03  4.87e-04  8.63e-04  4.9e-02  4.1e-01  9.6e-02  9.69e-04
##     15   48 -7.666e+03  1.02e-04  1.86e-04  1.9e-07  4.5e+00  3.5e-07  3.94e-04
##     16   49 -7.666e+03  1.12e-07  1.01e-07  1.9e-07  2.0e+00  3.5e-07  6.22e-05
##     17   57 -7.666e+03  1.60e-05  2.70e-05  2.0e-03  9.3e-01  3.7e-03  6.10e-05
##     18   59 -7.666e+03  5.23e-06  7.01e-06  3.7e-03  3.9e-01  8.0e-03  7.77e-06
##     19   60 -7.666e+03  4.08e-08  3.74e-08  1.4e-04  0.0e+00  3.1e-04  3.74e-08
##     20   61 -7.666e+03  2.31e-09  8.57e-10  8.6e-06  0.0e+00  2.0e-05  8.57e-10
##     21   62 -7.666e+03  5.35e-11  2.25e-13  7.6e-07  0.0e+00  1.6e-06  2.25e-13
##     22   63 -7.666e+03  1.81e-12  7.06e-16  1.7e-08  0.0e+00  3.4e-08  7.06e-16
##     23   64 -7.666e+03  6.98e-14  1.69e-17  1.0e-09  0.0e+00  2.4e-09  1.69e-17
##     24   65 -7.666e+03 -1.16e-14  1.76e-20  1.9e-10  0.0e+00  4.0e-10  1.76e-20
## 
##  ***** X- AND RELATIVE FUNCTION CONVERGENCE *****
## 
##  FUNCTION    -7.665775e+03   RELDX        1.874e-10
##  FUNC. EVALS      65         GRAD. EVALS      24
##  PRELDF       1.760e-20      NPRELDF      1.760e-20
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    4.639289e-06     1.000e+00    -2.337e-02
##      2    6.832875e-02     1.000e+00    -8.294e-07
##      3    8.890666e-01     1.000e+00    -2.230e-06
summary(dax.garch)       
## 
## Call:
## garch(x = dax)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -12.18398  -0.47968   0.04949   0.65746   4.48048 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 4.639e-06   7.560e-07    6.137 8.42e-10 ***
## a1 6.833e-02   1.125e-02    6.073 1.25e-09 ***
## b1 8.891e-01   1.652e-02   53.817  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 12947, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.13566, df = 1, p-value = 0.7126
plot(dax.garch) 

Задатак

  • Симулирати 1000 вредности временске серије која је дата моделом аt=wthtа_t =w_t\sqrt{h_t}, где је ht=0.1+0.4at1+0.2ht1h_t=0.1+0.4 a_{t-1} + 0.2h_{t-1}, wtw_t Гаусов бели шум.
  • Нацртати корелограм симулиране серије и корелограм квадриране серије.