Inicialmente, vamos importar as bibliotecas necessárias.

library(tseries)
library(timeSeries)
library(quantmod)
library(fGarch)
library(sn)

Vamos utilizar o getSymbols para obter os dados do Ibovespa, de janeiro de 2007 até a data presente.

ibov = getSymbols('^BVSP', src='yahoo', from= '2007-01-01', auto.assign = F)[,4]
colnames(ibov)= c('ibov')
ibov = ibov[is.na(ibov)==F]   # elimina os na's da amostra

Para ver as seis primeiras observações da amostra usarei o comando head.

head(ibov)
##             ibov
## 2007-01-02 45383
## 2007-01-03 44445
## 2007-01-04 44020
## 2007-01-05 42245
## 2007-01-08 42830
## 2007-01-09 42007
plot(ibov)

Agora vou calcular os retornos logaritmos do Ibov.

ret = diff(log(ibov))
colnames(ret) = c('ret')
ret = ret[is.na(ret)==F]  # Drop na to work

O próximo passo é plotar os retornos e a autocorrelação desses retornos.

par(mfrow=c(2,1))
plot(ret, main='Evolução dos retornos do IBOV')
acf(ret)   # There isn't autocorrelation in returns. So isn't possible predict ret

Note que os retornos não são autocorrelacionados. Logo não há uma estrutura de dependência entre eles. Sendo assim não é possivel prever tais retornos. Tente você mesmos estimar um Modelo ARIMA com esses dados. Você perceberá que os coeficientes estimados não são significativos.

Agora vamos fazer o mesmo para o quadrado dos retornos.

par(mfrow=c(2,1))
plot(ret^2, main= 'Evolução do quadrado dos retornos do IBOV')
acf(ret^2)   # In squares of returns there is autocorrelation

As estatíscas descritivas do IBOV e seus retornos podem ser calculadas utilizando a função basicStats da biblioteca fGarch.

basicStats(ret)
##                     ret
## nobs        3277.000000
## NAs            0.000000
## Minimum       -0.159930
## Maximum        0.136766
## 1. Quartile   -0.008486
## 3. Quartile    0.009549
## Mean           0.000159
## Median         0.000679
## Sum            0.520295
## SE Mean        0.000318
## LCL Mean      -0.000465
## UCL Mean       0.000783
## Variance       0.000332
## Stdev          0.018222
## Skewness      -0.475159
## Kurtosis       9.841399
basicStats(ibov)
##                     ibov
## nobs        3.278000e+03
## NAs         0.000000e+00
## Minimum     2.943500e+04
## Maximum     1.195280e+05
## 1. Quartile 5.233175e+04
## 3. Quartile 6.828175e+04
## Mean        6.299761e+04
## Median      5.944500e+04
## Sum         2.065062e+08
## SE Mean     2.857887e+02
## LCL Mean    6.243727e+04
## UCL Mean    6.355795e+04
## Variance    2.677313e+08
## Stdev       1.636250e+04
## Skewness    1.221181e+00
## Kurtosis    1.387047e+00

Agora vamos estimar o modelo GARCH(1, 1). Usualmente estimamos vários modelos GARCH de ordens diferentes, calculamos o critério de informação de Akaike ou o critério de informação Bayesiano, e o modelo com o maior desses indicadores é considerado o mais adequado. No entanto, já sabemos que para séries financeiras o GARCH(1, 1) resolve bem o problema.

O modelo GARCH(p, q) pode ser escrito como:

\[ \sigma^2_t = \omega + \sum_{i=1}^p \alpha_i \epsilon^2_{t-i} + \sum_{j=1}^q \beta_j \sigma^2_{t-j}\]

Nessa primeira estimação, vou utilizar a biblioteca fGarch.

garch1 = garchFit(formula = ~garch(1,1),data = ret)
garch1@fit$matcoef
##            Estimate   Std. Error   t value     Pr(>|t|)
## mu     5.417886e-04 2.461217e-04  2.201304 2.771454e-02
## omega  7.180395e-06 1.515319e-06  4.738537 2.152672e-06
## alpha1 8.566411e-02 9.960498e-03  8.600384 0.000000e+00
## beta1  8.895676e-01 1.296377e-02 68.619516 0.000000e+00

Agora que temos um objeto chamado garch1 vamos extrair desse objeto a volatilidade estimada.

data = index(ret)
vol = garch1@sigma.t
vol = as.xts(vol, order.by = data)
par(mfrow=c(2, 2))
plot(ibov, main='Evolution of IBOV')
plot(ret, main='Evolution of IBOV returns')
plot(ret^2, main='Evolution of square \n of IBOV returns')
plot(vol, main= 'Volatility of IBOV \n returns by GARCH(1,1)')

Note que a volatilidade estimada pelo GARCH é bastante elevada já a partir de janeiro de 2020. Tal volatilidade é maior que da crise de 2008. Esse fato pode ser creditado a crise do COVID-19.

Agora vamos estimar o GARCH(1,1) utilizando a biblioteca rugarch.

library(rugarch)
spec = ugarchspec()
spec1 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), 
                       mean.model=list(armaOrder=c(0,0), include.mean=TRUE),  
                       distribution.model="norm")

garch2 = ugarchfit(spec = spec1, data= ret)
garch2@fit$matcoef
##            Estimate   Std. Error    t value     Pr(>|t|)
## mu     5.414678e-04 2.458425e-04   2.202499 2.763010e-02
## omega  7.137185e-06 1.280061e-06   5.575659 2.465949e-08
## alpha1 8.552554e-02 5.697098e-03  15.012125 0.000000e+00
## beta1  8.898728e-01 7.473327e-03 119.073185 0.000000e+00

Use o ugarchforecast para prever 5 passos a frente.

pred2 = ugarchforecast(fitORspec = garch2, n.ahead = 5)
pred2 = pred2@forecast[['sigmaFor']]
pred2
##     2020-04-07
## T+1 0.05060747
## T+2 0.05005243
## T+3 0.04950505
## T+4 0.04896524
## T+5 0.04843291

Assumimos que os retornos do IBOV tem distribuição normal para estimar o modelo GARCH(1, 1). No entanto, sabemos que isso não é verdade. Utizando o teste de Kolmogorov-Smirnov cuja hipótese nula é de normalidade na série, nota-se que tal hipótese é rejeitada.

ksnormTest(ret)
## 
## Title:
##  One-sample Kolmogorov-Smirnov test
## 
## Test Results:
##   STATISTIC:
##     D: 0.4728
##   P VALUE:
##     Alternative Two-Sided: < 2.2e-16 
##     Alternative      Less: < 2.2e-16 
##     Alternative   Greater: < 2.2e-16 
## 
## Description:
##  Wed Apr 08 14:40:55 2020 by user: user

Olhando para o histograma e para o QQ plot abaixo, nota-se que a distribuição dos retornos do IBOV possui caldas mais pesadas do que a distribuição normal teórica.

par(mfrow=c(2,1))
qqnorm(ret)
qqline(ret, col=2)
hist(ret, breaks = 35, col = 'lightgreen', main = 'Histograma dos retornos \n do IBOV')

Agora vamos supor que os dados dos retornos do IBOV tem distribuição skew normal, que tem um peso maior nas caudas, e estimar o modelo.

spec2 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), 
                       mean.model=list(armaOrder=c(0,0), include.mean=TRUE),  
                       distribution.model="snorm")

garch3 = ugarchfit(spec = spec2, data= ret)
garch3@fit$matcoef
##            Estimate   Std. Error    t value     Pr(>|t|)
## mu     4.424236e-04 2.456837e-04   1.800785 0.0717367099
## omega  6.402579e-06 1.813073e-06   3.531341 0.0004134578
## alpha1 8.222655e-02 5.697258e-03  14.432653 0.0000000000
## beta1  8.954728e-01 6.970409e-03 128.467755 0.0000000000
## skew   9.041227e-01 1.966559e-02  45.974871 0.0000000000

Fazendo a previsão para 5 pass a frente.

pred3 = ugarchforecast(fitORspec = garch3, n.ahead = 5)
pred3 
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 5
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=2020-04-07]:
##        Series   Sigma
## T+1 0.0004424 0.05145
## T+2 0.0004424 0.05094
## T+3 0.0004424 0.05043
## T+4 0.0004424 0.04993
## T+5 0.0004424 0.04943

Também posso estimar o modelo GARCH utilizando a biblioteca arch do Python.

import matplotlib.pyplot as plt
from arch import arch_model

am = arch_model(r.ret, lags=1, vol='Garch')

print(f'\033[1;033m {am.fit() } ')
## Iteration:      1,   Func. Count:      6,   Neg. LLF: -9027.629664690852
## Iteration:      2,   Func. Count:     20,   Neg. LLF: -9029.000722060973
## Iteration:      3,   Func. Count:     31,   Neg. LLF: -9029.0969488853
## Iteration:      4,   Func. Count:     46,   Neg. LLF: -9029.160078156841
## Iteration:      5,   Func. Count:     60,   Neg. LLF: -9029.160126808416
## Iteration:      6,   Func. Count:     74,   Neg. LLF: -9029.160209442904
## Iteration:      7,   Func. Count:     87,   Neg. LLF: -9029.160212299781
## Optimization terminated successfully.    (Exit mode 0)
##             Current function value: -9029.160213633248
##             Iterations: 9
##             Function evaluations: 98
##             Gradient evaluations: 7
##                       Constant Mean - GARCH Model Results                      
## ==============================================================================
## Dep. Variable:                      y   R-squared:                      -0.000
## Mean Model:             Constant Mean   Adj. R-squared:                 -0.000
## Vol Model:                      GARCH   Log-Likelihood:                9029.16
## Distribution:                  Normal   AIC:                          -18050.3
## Method:            Maximum Likelihood   BIC:                          -18025.9
##                                         No. Observations:                 3277
## Date:                Wed, Apr 08 2020   Df Residuals:                     3273
## Time:                        14:40:58   Df Model:                            4
##                                  Mean Model                                 
## ============================================================================
##                  coef    std err          t      P>|t|      95.0% Conf. Int.
## ----------------------------------------------------------------------------
## mu         5.6235e-04  4.356e-05     12.909  4.005e-38 [4.770e-04,6.477e-04]
##                               Volatility Model                              
## ============================================================================
##                  coef    std err          t      P>|t|      95.0% Conf. Int.
## ----------------------------------------------------------------------------
## omega      6.7764e-06  4.160e-10  1.629e+04      0.000 [6.776e-06,6.777e-06]
## alpha[1]       0.1000  1.482e-02      6.746  1.515e-11   [7.095e-02,  0.129]
## beta[1]        0.8800  1.197e-02     73.497      0.000     [  0.857,  0.903]
## ============================================================================
## 
## Covariance estimator: robust 
## 
## C:\Users\user\Anaconda3\lib\site-packages\arch\univariate\base.py:260: DataScaleWarning: y is poorly scaled, which may affect convergence of the optimizer when
## estimating the model parameters. The scale of y is 0.0003319. Parameter
## estimation work better when this value is between 1 and 1000. The recommended
## rescaling is 100 * y.
## 
## This warning can be disabled by either rescaling y before initializing the
## model or by setting rescale=False.
## 
##   DataScaleWarning)