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

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

acf(SP500)

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

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

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

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

\(\color{deepskyblue}{\text{ARCH модели }}\)

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

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

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

\(\color{deepskyblue}{\text{GARCH модели }}\)

\(GARCH(p, q)\) модели (generalized ARCH) су уопштење \(ARCH\) модела облика: \[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,\] где је \(\alpha_0>0, \alpha_1, \dots, \alpha_p, \beta_1, \dots, \beta_q \ge 0\) и \(e_t \sim \mathcal{N}(0,1)\) је Гаусов бели шум.

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

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

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

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

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

Како дефиниција \(GARCH\) модела даје прогнозу за један корак унапред (јер се \(X_t\) рачуна помоћу \(X_{t−1}\)), изражавањем \(X_{t−1}\) преко ранијих вредности може се доћи до формуле за прогнозу за два корака унапред, итд.

Врло често се \(GARCH\) модели комбинују са \(ARMA\) или \(ARIMA\) моделом, којима моделујемо саму серују, а \(GARCH\) моделом серију резидуала. У том случају, \(GARCH\) модел нема утицај на прогнозу саме серије у неком будућем тренутку, али може утицати на дисперзију прогнозираних вредности.

\(\color{deepskyblue}{\text{Пример 1}}\)

library(tseries)
## Warning: package 'tseries' was built under R version 3.3.3
library(astsa)
## Warning: package 'astsa' was built under R version 3.3.3
plot(nyse)

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

Видимо да се волатилност знатно мења око 1000. дана, па зато користимо \(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

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

Погледајмо како изгледа прогноза \(GARCH\) модела.

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")

\(\color{deepskyblue}{\text{Пример 2}}\)

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

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

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

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
plot(stemp.garch.res)

acf(stemp.garch.res)

acf(stemp.garch.res^2)

\(\color{deepskyblue}{\text{Пример 3}}\)

Posmatramo vremensku seriju Daily Closing Prices of Major European Stock Indices, 1991–1998. Fitovaćemo GARCH(1,1) model.

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) 

\(\color{deepskyblue}{\text{Задатак}}\)

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