Посматрамо податке са берзе, односно дневене вредности 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)
За моделованје наведених временских серија користимо \(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\).
\(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\) модел нема утицај на прогнозу саме серије у неком будућем тренутку, али може утицати на дисперзију прогнозираних вредности.
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")
Посматрамо временску серију која садржи температуру на јужној хемисфери и периоду од 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)
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)