Introdução

Escolher um ativo de uma série histórica de sete anos, extrair o retorno diário e realizar a análise estatística dos dados. Primeiramente é realizada a análise dos outliers, verificação de liquidez, assimetria, curtose, teste de normalidade e correlação. Na segunda parte será realizada a estimação da volatilidade do ativo através da volatilidade histórica, média móvel simples variando o tamanho da janela, EWMA (exponentially weighted moving average) realizando variações com o parâmetro de decaimento e modelo GARCH. Desta forma, serão abordados os tópicos da disciplina.

Preparação do Ambiente R

Implementação

Determinação de alguns parâmetros que serão utilizados na análise

Usaremos uma série de 7 anos, nesse intervalo de tempo será possível capturar dias de estresses na economia como o Joesley Day (17/05/2017) , greve dos caminhoneiros (21/05/2018) e pandemia do corona vírus em março de 2020.

anos <- 7 #####
first.date <- Sys.Date() - 365*anos
last.date <- Sys.Date()

Seleção do ativo que será analizado

O ativo escolhido para realizar o a aplicação dos modelos é ELET3, a Eletrobras é a maior e mais relevante concessionaria de energia da América Latina, na qual representa 37% da capacidade de geração instalada do Brasil e 45% de rede de transmissão nacional. Os múltiplos da empresa, vem apresentando uma melhora assim como a redução de alavancagem de aproximadamente 5 anos para cá dado a reformulação da empresa. A Eletrobras é uma sociedade de economia mista e de capital aberto sob controle acionário do Governo Federal, criada em 1962 para coordenar as empresas do setor elétrico. Atua na geração, transmissão e distribuição do sistema elétrico brasileiro. Composta por Eletrobras Chesf, Eletronorte, Eletronuclear, CGT Eletrosul, Furnas, Eletropar e Cepel. Em 2018, o presidente Michel Temer enviou ao Congresso Nacional o projeto de lei que dispõe sobre a privatização da Eletrobras. Recentemente fortes rumores de mercado levam em consideração a privatização da empresa, podendo assim destravar valores significativos e assim trazer um aumento no preço do papel.

ATIVO_input <-'ELET3.SA'

E agora podemos fazer o dowload dos que nos interessam.

# 4 - FAZER DOWNLOAD DOS DADOS DE COTACAO DE FECHAMENTO
  #Baixa o ativo
ativo_ticker.out <- getSymbols(ATIVO_input, src = "yahoo", from = first.date, to = last.date, auto.assign = FALSE)
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
## Warning: ELET3.SA contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.

É possível usar comandos específicos para entender melhor os dados que foram baixados. Assim, concluímos que nossa série é composta pelo preço do ativo na abertura (ELET3.SA.Open), maior e menor preço intraday (ELET3.SA.High, ELET3.SA.Low), preço de fechamento (ELET3.SA.Close) , volume negociado no dia (ELET3.SA.Volume) e o preço de ajuste (ELET3.SA.Adjusted). Nesta publicação usaremos o preço de ajusto, devido a mudança no preço no ativo em caso de dividendos, juros, agrupamentos, desdobramentos e outras mudanças que impactam na rentabilidade atrelada a precificação do papel.

ativo_ticker.out %>% glimpse()
## An 'xts' object on 2014-01-13/2021-01-08 containing:
##   Data: num [1:1737, 1:6] 5,81 5,8 5,74 5,66 5,46 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:6] "ELET3.SA.Open" "ELET3.SA.High" "ELET3.SA.Low" "ELET3.SA.Close" ...
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
## List of 2
##  $ src    : chr "yahoo"
##  $ updated: POSIXct[1:1], format: "2021-01-11 21:49:04"
names(ativo_ticker.out)
## [1] "ELET3.SA.Open"     "ELET3.SA.High"     "ELET3.SA.Low"     
## [4] "ELET3.SA.Close"    "ELET3.SA.Volume"   "ELET3.SA.Adjusted"
DT::datatable(head(ativo_ticker.out))

O pacote getSymbols gera um gráfico Candlesticks, no qual é possível ter a sensibilidade da variação do preço do ativo intraday e o volume negociado no dia.

chartSeries(ativo_ticker.out) 
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both ',', which could be confusing

## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both ',', which could be confusing

Transformação dos dados em séries de tempo

Esta etapa é importante para que se possa aplicar as técnicas de séries temporais aos dados de interesse. Além de transformar em séries de tempo, os dados não disponíveis são retirados. Conforme mencionado anteriormente, serão extraídos apenas o preço de ajuste e são removidas as datas não possuem informações devido a feriados.

ATIVO <- na.omit(ativo_ticker.out[,6])

Gráfico dos preços

O primeiro exercício que pode ser realizado é a construção dos gráficos da série temporal dos preços do ativo selecionado. Abaixo é feito o gráfico mais simples que pode ser feito usando o R.

  highchart(type = "stock") %>% 
    hc_add_series(ATIVO$ELET3.SA.Adjusted, type = "line",color = "darkblue",name = "Preço") %>% 
    hc_title(text = paste0("Preço de ajuste ",ATIVO_input)) %>% 
        hc_yAxis(opposite = FALSE) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

Podemos verificar na série uma tendencia de crescimento, com períodos de alta volatilidade do mercado. Há alguns fatos estilizados sobre séries temporais de ativos de financeiros que podem ser apontados.

  1. A série de preços apresenta características de ser um passeio aleatório, sem reversão à média. Esta é uma característica de um Passeio Aleatório, cuja forma mais simples pode ser escrita como:

\[ P_t = P_{t-1} + \epsilon_t, \] na qual \(\epsilon_t \sim D(0, \sigma^2), t = 1, \ldots, T\) e \(Cov(P_t, P_{t-s}) = \sigma_s\), $s = 1, 2 , $ Assim, o preço no instante de preço t é igual ao preço no instante t-1 somada a um $_t $.

Para que esta série tenha reversão à média, deveria haver um \(\beta < 1\), módulo. Veja, por exemplo:

\[ P_t = \beta\cdot P_{t-1} + \epsilon_t, |\beta| < 1 \]

\[ P_t = \beta \cdot ( \beta\cdot P_{t-2} + \epsilon_{t-2}), |\beta| < 1 \] \[ P_t = \beta^2\cdot P_{t-2} + \beta\cdot \epsilon_{t-1} + \epsilon_{t-2}, |\beta| < 1 \] Assim, podemos encontrar $_t $ através da diferença de preços de t e t-1.

\[ P_t - P_{t-1} = \epsilon_t.\] O comportamento da série se mantém, ao utilizamos o recurso matemático de transformar a série de preços da ações em logaritmo. Assim, quando plotamos o gráfico \(\log P_t\) verificamos a características de um passeio aleatório, sem reversão à média.

\[ \log P_t - \log P_{t-1} = \varepsilon_t,\] \[ \varepsilon_t \sim D(0, \sigma_t) \]

    highchart(type = "stock") %>% 
    hc_add_series(log(ATIVO$ELET3.SA.Adjusted), type = "line",color = "orange",name = "LogPreço") %>% 
    hc_title(text = paste0("Log Preço de ajuste ",ATIVO_input)) %>% 
         hc_yAxis(opposite = FALSE, labels = list(format = "{value}%")) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

CALCULAR RETORNOS

Retorno Simples

Há duas formas de obter a diferença. Podemos obtê-la usando o comando head(100*(ATIVO - lag(ATIVO, 1))/lag(ATIVO, 1)) ou usando o comando diff. Matematicamente,

\[ R_t = \frac{P_t} {P_{t-1}} - 1 ,\] \[ R_t = \frac{P_t-P_{t-1}} {P_{t-1}} .\]

rATIVO = na.omit(100* diff(ATIVO, lag = 1, differences = 1, arithmetic = TRUE, log = FALSE, na.pad = TRUE)/stats::lag(ATIVO, 1))
head(rATIVO)
##            ELET3.SA.Adjusted
## 2014-01-14             -0,69
## 2014-01-15             -1,57
## 2014-01-16             -4,07
## 2014-01-17             -0,18
## 2014-01-20              0,37
## 2014-01-21             -0,92

O gráfico abaixo permite observar outro fato estilizado.

  1. Há formação de clusters de volatilidade.
    highchart(type = "stock") %>% 
    hc_add_series((rATIVO$ELET3.SA.Adjusted), type = "line",color = "darkred",name="Retorno Simples") %>% 
    hc_title(text = paste0("Retorno Simples Preço de ajuste ",ATIVO_input)) %>% 
        hc_yAxis(opposite = FALSE) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

Retorno logaritmo

O uso do log retorno ou retorno contínuo, possui vantagens em relação ao retorno discreto (ou simples).São elas: 1. Se considerarmos múltiplos períodos podemos somar o log retorno de cada instante de tempo t; 2. Existem vantagens estatística de trabalhar com log retorno.

\[ r_t = \log{P_t}-log {P_{t-1}} ,\]

  rlogATIVO = na.omit(100* diff(log(ATIVO$ELET3.SA.Adjusted), lag=1))

Os clusters de volatilidade também são observados com os log retornos.

    highchart(type = "stock") %>% 
    hc_add_series((rlogATIVO$ELET3.SA.Adjusted), type = "line",color = "red",name = "LogRetorno") %>% 
    hc_title(text = paste0("Log Retorno Preço de ajuste ",ATIVO_input)) %>% 
        hc_yAxis(opposite = FALSE, labels = list(format = "{value}%")) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

Retorno logaritmo ao quadrado

O gráfico dos quadrados dos retornos logaritmos (simplesmente retornos a partir de agora) também evidencia a existência de clusters de volatilidade. O uso dessa visualização do quadrado do retorno isola o efeito positivo ou negativo do retorno no período, assim podemos verificar o cluster e a alta volatilidade graficamente.

    highchart(type = "stock") %>% 
    hc_add_series((rlogATIVO$ELET3.SA.Adjusted)^2,type= "bar",color = "blue",name = "Preço^2") %>% 
    hc_title(text = paste0("Log Retorno ao Quadrado ",ATIVO_input))

Além dos gráficos, estatísticas resumo ajudam a compreender outros fatos estilizados de séries de retornos financeiros. Para esta parte, o conteúdo apresentado em Torres-Reyna pode ser muito útil.

A tabela abaixo apresenta estatísticas descritivas da série de interesse.
Verificamos, que a série possui média de 0,11, desvio padrão de 3,7, assimetria 0,71 e curtose 12,25. Podemos verificar que não apresenta a característica de uma distribuição normal, dado que a normal possui uma curtose próxima de 3 e uma assimetria 0. Quando uma distribuição possui uma curtose maior que 3 verificamos um peso maior na cauda, que é um comportamento dominante no mercado financeiro.

Para observar melhor a distribuição serão apresentados: o histograma; box-plot dos log retornos e qq-plot. Através do histograma e distribuição do log retorno dos preços podemos verificar o peso na cauda o que difere da distribuição normal. A comparação com a distribuição normal é verificada no qqplot na qual os valores extremos não seguem o comportamento de uma distribuição gaussiana.

Medida Estimativa
Percentil 1% -8,93
Percentil 5% -5,05
Percentil 10% -3,69
Média 0,11
Mediana 0
Percentil 90% 4,06
Percentil 95% 5,83
Percentil 99% 9,76
Desvio-padrão 3,7
Variância 13,68
Assimetria 0,71
Curtose 12,25

É possível obter um contjunto de estatísticas descritivas usando apenas um comando:

basicStats(rlogATIVO)
##             ELET3.SA.Adjusted
## nobs                 1733,000
## NAs                     0,000
## Minimum               -23,534
## Maximum                40,076
## 1. Quartile            -1,802
## 3. Quartile             1,893
## Mean                    0,115
## Median                  0,000
## Sum                   198,809
## SE Mean                 0,089
## LCL Mean               -0,060
## UCL Mean                0,289
## Variance               13,682
## Stdev                   3,699
## Skewness                0,707
## Kurtosis               12,233
summary(rlogATIVO)
##      Index            ELET3.SA.Adjusted
##  Min.   :2014-01-14   Min.   :-24      
##  1st Qu.:2015-10-13   1st Qu.: -2      
##  Median :2017-07-12   Median :  0      
##  Mean   :2017-07-10   Mean   :  0      
##  3rd Qu.:2019-04-09   3rd Qu.:  2      
##  Max.   :2021-01-08   Max.   : 40
chart.Drawdown(rlogATIVO$ELET3.SA.Adjusted)

chart.Histogram(rlogATIVO$ELET3,methods = c("add.density", "add.normal"),ylim = c(0,0.25))

chart.Histogram(rlogATIVO$ELET3,methods = c("add.density", "add.rug"),ylim = c(0,0.25))

t_logret <- as_tibble(rlogATIVO) %>% arrange((ELET3.SA.Adjusted))
## Warning: Calling `as_tibble()` on a vector is discouraged, because the behavior is likely to change in the future. Use `tibble::enframe(name = NULL)` instead.
## This warning is displayed once per session.
t_logret %>% 
  ggplot(aes(x = ELET3.SA.Adjusted/100)) +
  geom_histogram(alpha = 0.25, binwidth = .01,)

  ggtitle("Log Retorno Histogram")
## $title
## [1] "Log Retorno Histogram"
## 
## attr(,"class")
## [1] "labels"
ggplot(rlogATIVO) + 
  geom_density(aes(x = ELET3.SA.Adjusted),color='black') +
  ggtitle("Dist. Log Retorno")

qqnorm(rlogATIVO)
qqline(rlogATIVO)

box_plot_df <- as_tibble(rlogATIVO) %>% mutate(NOME = ATIVO_input)
 


bp_df <- box_plot_df %>%
  ggplot( aes(x=NOME, y=ELET3.SA.Adjusted, fill=NOME)) +
    geom_boxplot() +
   
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("BoxPlotLog Retorno") +
    xlab("")
ggplotly(bp_df)

Análise de autocorrelação Serial

Autocorrelação serial do preço de fechamento do ativo

A Função de autocorrelação serial simples dos preços (ACF) captura a dinâmica linear, assim podemos verificar que a série sugere uma forte correlação do preço no instante t possui uma forte correlação linear no instante t-1. O gráfico da função de autocorrelação serial simples apresenta comportamento típico de uma acf-simples de séries não estacionária. Há grande persistência, a acf-simples apresenta decaimento muito lento.

acf(ATIVO)
title(ATIVO_input, line = -1)

Por outro lado, o gráfico da função de autocorrelação serial parcial apresenta um valor significativo para a primeira defasagem (muito próximo a um) e não significativo para dos demais. Este comportamento é típico de um passeio aleatório.

  acf(ATIVO, type="partial")
  title(ATIVO_input, line = -1)

Autocorrelação serial do retorno do ativo

Por outro lado, os gráficos de autocorrelação simples e parcial dos retornos não apresentam nenhuma estrutura de dependência temporal.

acf(rlogATIVO)
title(paste("ACF simples do log retorno do", ATIVO_input), line = +1)

  ###### Rafael falou isto!
  acf(rlogATIVO, type="partial")
  title(paste("ACF parcial do log retorno do", ATIVO_input), line = +1)

Testes de Raiz unitária

TESTE DE RAIZ UNITARIA

O teste da raiz unitária tem por objetivo verificar a estacionaridade nas séries realizando um estudo da existência de alguma raiz dos operadora dentro do círculo unitário. O teste é dado por:

\[ H_0 = Existe\; pelo\; menos\; uma\; raiz\; dentro\; do\; circulo \; unitário\] \[ H_1 = Não\; Existem\; raizes\; dentro\; do\; circulo\; unitário\]

Nesta publicação abordaremos o deste de Dickey-Fuller. Assim, para verificar a raiz unitário de um \(AR(p)\) temos que: \[ x_t=c_t +\beta x_{t-1}+\sum _{i=1}^{p-1}\phi\Delta_{t-i}+\epsilon_t \] Assim para verificar a estacionaridade, $H_0 : =1 $ e $H_1 : |t|)
## z.lag.1 0,000284 0,000846 0,34 0,737
## z.diff.lag -0,043886 0,024051 -1,82 0,068 . ## — ## Signif. codes: 0 ‘’ 0,001 ’’ 0,01 ’’ 0,05 ‘.’ 0,1 ‘’ 1 ## ## Residual standard error: 0,71 on 1730 degrees of freedom ## Multiple R-squared: 0,00196, Adjusted R-squared: 0,000807 ## F-statistic: 1,7 on 2 and 1730 DF, p-value: 0,183 ## ## ## Value of test-statistic is: 0,34 ## ## Critical values for test statistics: ## 1pct 5pct 10pct ## tau1 -2,6 -1,9 -1,6 ```

  summary(ur.df(ATIVO, type = "drift",  selectlags = "BIC"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4,848 -0,246 -0,029  0,224  6,100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0,04546    0,03162    1,44    0,151  
## z.lag.1     -0,00161    0,00156   -1,03    0,304  
## z.diff.lag  -0,04324    0,02405   -1,80    0,072 .
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Residual standard error: 0,71 on 1729 degrees of freedom
## Multiple R-squared:  0,00256,    Adjusted R-squared:  0,0014 
## F-statistic: 2,22 on 2 and 1729 DF,  p-value: 0,109
## 
## 
## Value of test-statistic is: -1 1,1 
## 
## Critical values for test statistics: 
##      1pct 5pct 10pct
## tau2 -3,4 -2,9  -2,6
## phi1  6,4  4,6   3,8
  summary(ur.df(ATIVO, type = "trend", selectlags = "BIC"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4,902 -0,257 -0,018  0,226  6,051 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  8,26e-03   3,42e-02    0,24   0,8090   
## z.lag.1     -1,05e-02   3,50e-03   -2,99   0,0028 **
## tt           2,17e-04   7,67e-05    2,83   0,0047 **
## z.diff.lag  -3,91e-02   2,40e-02   -1,63   0,1040   
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Residual standard error: 0,71 on 1728 degrees of freedom
## Multiple R-squared:  0,00715,    Adjusted R-squared:  0,00543 
## F-statistic: 4,15 on 3 and 1728 DF,  p-value: 0,00612
## 
## 
## Value of test-statistic is: -3 3,4 4,5 
## 
## Critical values for test statistics: 
##      1pct 5pct 10pct
## tau3 -4,0 -3,4  -3,1
## phi2  6,1  4,7   4,0
## phi3  8,3  6,2   5,3

SERIE COM RETORNO LOGARITMO

O teste de Dickey-Fuller será aplicado agora aos retornos logaritmos do ativo ELET3. Verificamos que atraves do p-value: <2e-16 nos tres modelos de regressão (sem intercepto e sem tendencia, com intercepto e com intercepto e tendencia) que rejeitamos \(H_0\), ou seja a série é estacionária.

  summary(ur.df(rlogATIVO, type = "none", selectlags = "BIC"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23,56  -1,81   0,00   1,90  40,10 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1,00736    0,03412  -29,52   <2e-16 ***
## z.diff.lag  0,00108    0,02405    0,04     0,96    
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Residual standard error: 3,7 on 1729 degrees of freedom
## Multiple R-squared:  0,503,  Adjusted R-squared:  0,503 
## F-statistic:  875 on 2 and 1729 DF,  p-value: <2e-16
## 
## 
## Value of test-statistic is: -30 
## 
## Critical values for test statistics: 
##      1pct 5pct 10pct
## tau1 -2,6 -1,9  -1,6
  summary(ur.df(rlogATIVO, type = "drift", selectlags = "BIC"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23,68  -1,93  -0,12   1,78  39,98 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0,11723    0,08908    1,32     0,19    
## z.lag.1     -1,00932    0,03414  -29,56   <2e-16 ***
## z.diff.lag   0,00206    0,02406    0,09     0,93    
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Residual standard error: 3,7 on 1728 degrees of freedom
## Multiple R-squared:  0,504,  Adjusted R-squared:  0,503 
## F-statistic:  876 on 2 and 1728 DF,  p-value: <2e-16
## 
## 
## Value of test-statistic is: -30 437 
## 
## Critical values for test statistics: 
##      1pct 5pct 10pct
## tau2 -3,4 -2,9  -2,6
## phi1  6,4  4,6   3,8
  summary(ur.df(rlogATIVO, type = "trend", selectlags = "BIC"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23,69  -1,94  -0,13   1,79  39,98 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1,32e-01   1,78e-01    0,74     0,46    
## z.lag.1     -1,01e+00   3,42e-02  -29,55   <2e-16 ***
## tt          -1,72e-05   1,78e-04   -0,10     0,92    
## z.diff.lag   2,06e-03   2,41e-02    0,09     0,93    
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Residual standard error: 3,7 on 1727 degrees of freedom
## Multiple R-squared:  0,504,  Adjusted R-squared:  0,503 
## F-statistic:  584 on 3 and 1727 DF,  p-value: <2e-16
## 
## 
## Value of test-statistic is: -30 291 437 
## 
## Critical values for test statistics: 
##      1pct 5pct 10pct
## tau3 -4,0 -3,4  -3,1
## phi2  6,1  4,7   4,0
## phi3  8,3  6,2   5,3

Estimando a volatilidade

Volatilidade histórica

A volatilidade histórica é estimada por meio do desvio-padrão dos retornos dos ativos. Na operação abaixo foi feita uma operação neutra da adição apenas para que a nova variável gerada pudesse ter o formato de dados mais conveniente para plotagem. Ou seja, somar rlogATIVO*0 ao desvio padrão é uma forma de passar para VOL_historica a estrutura de dados presente em rlogATIVO. Assim, a volatidade histórica será igual ao da tabela apresentada anteriormente, dado que estamos calculando o desvio padrão do período de 7 anos.

VOL_historica <- sd(rlogATIVO)+0*rlogATIVO
  str(rlogATIVO)
## An 'xts' object on 2014-01-14/2021-01-08 containing:
##   Data: num [1:1733, 1] -0,694 -1,58 -4,156 -0,185 0,369 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "ELET3.SA.Adjusted"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
## List of 3
##  $ src      : chr "yahoo"
##  $ updated  : POSIXct[1:1], format: "2021-01-11 21:49:04"
##  $ na.action: 'omit' int 1
##   ..- attr(*, "index")= num 1,39e+09

O Gráfico abaixo evidencia que a volatilidade história, apesar se ser uma média que pode ser calculada facilmente, não é capaz de capturar corretamente o fato estilizado de que séries de tempo têm clusters de volatilidade. Neste caso, em cenário de calmaria do ativo estaríamos sobrestimando e em cenário de estresse estaríamos subestimando o risco do papel.

  colnames(VOL_historica) <- "VOL_HIST"
  rlog_estima <-   cbind(rlogATIVO,VOL_historica)
  
  highchart(type = "stock") %>% 
    hc_add_series((rlog_estima$ELET3.SA.Adjusted), type = "line",color = "red",name = "LogRetorno") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_HIST), type = "line",color = "blue",name = "Vol Historica") %>% 
    hc_title(text = paste0("Log Retorno Preço de ajuste - Vol Historica ",ATIVO_input)) %>% 
         hc_yAxis(opposite = FALSE, labels = list(format = "{value}%")) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

Modelo Media móvel

Com o objetivo de capturar melhor o movimento em determinado período da série, utilizamos uma janela móvel no cálculo da volatilidade. Neste publicação vamos utilizar a prox do mercado mensal (21 dias), trimestral (63 dias) e anual (252).

rlogATIVO_square <- rlogATIVO^2
VOL_MediaMovelSimples_21 <- rlogATIVO_square*0+sqrt(movavg(rlogATIVO_square, 21, "s"))
VOL_MediaMovelSimples_63 <- rlogATIVO_square*0+sqrt(movavg(rlogATIVO_square, 63, "s"))
VOL_MediaMovelSimples_252 <- rlogATIVO_square*0+sqrt(movavg(rlogATIVO_square, 252, "s"))

Podemos verificar que quando utilizamos um período melhor conseguimos mensurar melhor que os demais períodos, esse comportamento é explicado porque todos os retornos possuem pesos iguais na volatilidade. Ou seja, uma volatidade baixa num passado distante tem impacto relevante na série temporal. Explicado pelo mesmo fato, quando utilizamos períodos maiores a volatidade em tempo de calmaria a demora na saída de um dado da amostra está aumentando a volatidade de um período.

  colnames(VOL_MediaMovelSimples_21) <- "VOL_MM21"
  colnames(VOL_MediaMovelSimples_63) <- "VOL_MM63"
  colnames(VOL_MediaMovelSimples_252) <- "VOL_MM252"
  rlog_estima <-   cbind(rlog_estima,VOL_MediaMovelSimples_21,VOL_MediaMovelSimples_63,VOL_MediaMovelSimples_252)
  
  highchart(type = "stock") %>% 
    hc_add_series((rlog_estima$ELET3.SA.Adjusted), type = "line",color = "red",name = "LogRetorno") %>% 
    #hc_add_series(qnorm(0.05)*(rlog_estima$VOL_HIST), type = "line",color = "blue",name = "LogRetorno") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_MM21), type = "line",color = "yellow",name = "Vol Média Móvel 21") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_MM63  ), type = "line",color = "lightgreen",name = "Vol Média Móvel 100") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_MM252), type = "line",color = "darkgreen",name = "Vol Média Móvel 252") %>%
    hc_title(text = paste0("Log Retorno Preço de ajuste - Vol Media Movel ",ATIVO_input)) %>% 
     hc_yAxis(opposite = FALSE, labels = list(format = "{value}%")) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

EWMA

O EWMA dá a opção de colocar pesos maiores nas observações atuais, assim podemos aprimorar matematicamente o modelo de volatilidade histórica. Dado por:

\[ \sigma_t^2= (1-\lambda)r^2_{t-1} + \lambda \sigma _{t-1}^2. \] Assim, podemos variar o \(\lambda\), com o objetivo de dividir o peso nos retornos e volatilidades colocando peso nas observações recentes.

VOL_EWMA_50 <- rlogATIVO_square*0+sqrt(EWMAvol(rlogATIVO, lambda = 0.50)[1]$Sigma.t)

VOL_EWMA_94 <- rlogATIVO_square*0+sqrt(EWMAvol(rlogATIVO, lambda = 0.94)[1]$Sigma.t)
VOL_EWMA_96 <- rlogATIVO_square*0+sqrt(EWMAvol(rlogATIVO)[1]$Sigma.t)
VOL_EWMA_98 <- rlogATIVO_square*0+sqrt(EWMAvol(rlogATIVO, lambda = 0.98)[1]$Sigma.t)
      ### acima novamente aquele truque de usar pela série original multiplicada por zero para carregar o formato
  colnames(VOL_EWMA_50) <- "VOL_EWMA_50"
  colnames(VOL_EWMA_94) <- "VOL_EWMA_94"
  colnames(VOL_EWMA_96) <- "VOL_EWMA_96"
  colnames(VOL_EWMA_98) <- "VOL_EWMA_98"
  rlog_estima <-   cbind(rlog_estima,VOL_EWMA_50,VOL_EWMA_94,VOL_EWMA_96,VOL_EWMA_98)
  
  highchart(type = "stock") %>% 
    hc_add_series((rlog_estima$ELET3.SA.Adjusted), type = "line",color = "red",name = "LogRetorno") %>% 
    #hc_add_series(qnorm(0.05)*(rlog_estima$VOL_HIST), type = "line",color = "blue") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_50), type = "line",color = "orange",name = "EWMA 50%") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_94), type = "line",color = "pink",name = "EWMA 94%") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_96), type = "line",color = "darkblue",name = "EWMA 96%") %>%
     hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_98), type = "line",color = "black",name = "EWMA 98%") %>%
    hc_title(text = paste0("Log Retorno Preço de ajuste - Vol EWMA ",ATIVO_input)) %>% 
         hc_yAxis(opposite = FALSE, labels = list(format = "{value}%")) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)

Modelagem GARCH

O modelo GARCH(p,q) Heretocedasticidade Condicional Autoregressivo Geral pode ser escrito como:

\[ \sigma^2_t=\alpha_0+\sum_{i=1}^{p}\epsilon_{t-i}^2\alpha _i +\sum_{j=1}^{q}\beta_j\sigma^2_{t-j} \] Em finanças o modelo mais utilizado é o GARCH(1,1) \[ \sigma^2_t=\alpha_0+\alpha_1\cdot \epsilon_ {t-1}^2 +\beta_1 \cdot \sigma^2_{t-1} \]

Análise das autocorrelações dos quadrados dos retornos

Através da ACF verificamos um decaimento lento nos lags, assim podemos realizar teste e variação do modelo GARCH como GARCH(2,1), além disso iremos variar a distribuição em t-student e normal.

acf(rlogATIVO^2)
title(ATIVO_input, line = -1)

acf(rlogATIVO^2, type="partial")
title(ATIVO_input, line = -1)

Modelo GARCH

gfit.std <- garchFit(data=rlogATIVO, cond.dist="std")
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          std
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          1733
##  Recursion Init:            mci
##  Series Scale:              3,7
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                   U      V params includes
##     mu     -3,1e-01   0,31  0,031     TRUE
##     omega   1,0e-06 100,00  0,100     TRUE
##     alpha1  1,0e-08   1,00  0,100     TRUE
##     gamma1 -1,0e+00   1,00  0,100    FALSE
##     beta1   1,0e-08   1,00  0,800     TRUE
##     delta   0,0e+00   2,00  2,000    FALSE
##     skew    1,0e-01  10,00  1,000    FALSE
##     shape   1,0e+00  10,00  4,000     TRUE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1  shape 
##      1      2      3      5      8 
##  Persistence:                  0,9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     2222.4950: 0.0310144 0.100000 0.100000 0.800000  4.00000
##   1:     2222.2649: 0.0310116 0.101425 0.102861 0.802022  4.00014
##   2:     2222.1345: 0.0310051 0.0984407 0.104940 0.801008  4.00032
##   3:     2222.0612: 0.0309870 0.0972734 0.111957 0.803461  4.00092
##   4:     2221.8744: 0.0309363 0.0910026 0.114962 0.801351  4.00228
##   5:     2221.8016: 0.0307729 0.0906850 0.117796 0.802814  4.00668
##   6:     2221.3339: 0.0243558 0.0896743 0.0959237 0.815827  4.17521
##   7:     2220.7961: 0.0178701 0.0814594 0.101990 0.809193  4.34251
##   8:     2220.3926: 0.0147130 0.107943 0.119072 0.773189  4.41535
##   9:     2220.1971: 0.0138014 0.0964635 0.113101 0.790224  4.41198
##  10:     2220.1738: 0.0138002 0.0954703 0.112927 0.789668  4.41197
##  11:     2220.1623: 0.0137967 0.0950192 0.113881 0.790114  4.41199
##  12:     2220.1550: 0.0137326 0.0943972 0.114057 0.789855  4.41273
##  13:     2220.1461: 0.0135994 0.0945398 0.114591 0.790002  4.41430
##  14:     2219.9855: 0.00635351 0.100302 0.118730 0.777432  4.49935
##  15:     2219.8322: -0.000743501 0.0818110 0.109963 0.806853  4.40782
##  16:     2219.8116: -0.00157855 0.0862256 0.113247 0.801243  4.35211
##  17:     2219.8024: -0.00109467 0.0874646 0.113472 0.799004  4.41244
##  18:     2219.8015: -0.00107508 0.0880124 0.114074 0.797554  4.41132
##  19:     2219.8015: -0.00109396 0.0877962 0.113938 0.797999  4.41008
##  20:     2219.8015: -0.00109354 0.0877926 0.113933 0.798008  4.41011
## 
## Final Estimate of the Negative LLH:
##  LLH:  4487    norm LLH:  2,6 
##     mu  omega alpha1  beta1  shape 
## -0,004  1,201  0,114  0,798  4,410 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##            mu   omega   alpha1  beta1  shape
## mu     -216,1     5,4      1,5     51    2,1
## omega     5,4  -114,0   -824,6  -1239  -23,8
## alpha1    1,5  -824,6  -9624,6 -10581 -223,1
## beta1    51,4 -1239,3 -10581,0 -14593 -283,8
## shape     2,1   -23,8   -223,1   -284  -10,1
## attr(,"time")
## Time difference of 0,1 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0,68 secs
gfit.norm <- garchFit(data=rlogATIVO, cond.dist="norm")
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          1733
##  Recursion Init:            mci
##  Series Scale:              3,7
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                   U      V params includes
##     mu     -3,1e-01   0,31  0,031     TRUE
##     omega   1,0e-06 100,00  0,100     TRUE
##     alpha1  1,0e-08   1,00  0,100     TRUE
##     gamma1 -1,0e+00   1,00  0,100    FALSE
##     beta1   1,0e-08   1,00  0,800     TRUE
##     delta   0,0e+00   2,00  2,000    FALSE
##     skew    1,0e-01  10,00  1,000    FALSE
##     shape   1,0e+00  10,00  4,000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1 
##      1      2      3      5 
##  Persistence:                  0,9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     2381.7395: 0.0310144 0.100000 0.100000 0.800000
##   1:     2381.4699: 0.0310135 0.109520 0.0998004 0.802775
##   2:     2379.8672: 0.0310123 0.110608 0.0924809 0.796171
##   3:     2379.0774: 0.0310081 0.127949 0.0848068 0.790351
##   4:     2378.2876: 0.0310013 0.132156 0.0874684 0.771151
##   5:     2377.2734: 0.0309940 0.141624 0.0992132 0.758273
##   6:     2376.8054: 0.0309786 0.149995 0.102384 0.740577
##   7:     2375.6493: 0.0308847 0.172760 0.117259 0.711850
##   8:     2375.0952: 0.0304367 0.186835 0.129206 0.679843
##   9:     2374.8651: 0.0302295 0.196998 0.121909 0.679229
##  10:     2374.7789: 0.0298356 0.197237 0.125631 0.674091
##  11:     2374.7729: 0.0295888 0.198913 0.130203 0.673247
##  12:     2374.7103: 0.0294685 0.198830 0.128923 0.670956
##  13:     2374.6642: 0.0292447 0.204078 0.128671 0.668118
##  14:     2374.5894: 0.0286767 0.204773 0.128516 0.664235
##  15:     2374.4823: 0.0273939 0.206930 0.135927 0.657689
##  16:     2374.1953: 0.0193124 0.233940 0.140592 0.629709
##  17:     2373.9402: 0.0113119 0.255950 0.159641 0.584528
##  18:     2373.9229: 0.0150736 0.253838 0.157245 0.589977
##  19:     2373.9215: 0.0137775 0.250940 0.156041 0.594294
##  20:     2373.9211: 0.0140674 0.252048 0.156671 0.592509
##  21:     2373.9211: 0.0140714 0.252042 0.156698 0.592508
##  22:     2373.9211: 0.0140697 0.252039 0.156710 0.592505
## 
## Final Estimate of the Negative LLH:
##  LLH:  4641    norm LLH:  2,7 
##     mu  omega alpha1  beta1 
##  0,052  3,448  0,157  0,593 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##            mu  omega alpha1 beta1
## mu     -150,8    2,7    -60    19
## omega     2,7  -39,4   -280  -464
## alpha1  -59,8 -280,1  -3851 -3924
## beta1    19,5 -463,7  -3924 -5850
## attr(,"time")
## Time difference of 0,032 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0,13 secs
summary(gfit.std)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(data = rlogATIVO, cond.dist = "std") 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000002f09b818>
##  [data = rlogATIVO]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##         mu       omega      alpha1       beta1       shape  
## -0,0040449   1,2011825   0,1139328   0,7980078   4,4101069  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     -0,004045    0,068234   -0,059  0,95273    
## omega   1,201182    0,410453    2,926  0,00343 ** 
## alpha1  0,113933    0,027981    4,072 4,66e-05 ***
## beta1   0,798008    0,049686   16,061  < 2e-16 ***
## shape   4,410107    0,477043    9,245  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Log Likelihood:
##  -4487    normalized:  -2,6 
## 
## Description:
##  Mon Jan 11 20:49:22 2021 by user: Livia Rodrigues 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value
##  Jarque-Bera Test   R    Chi^2  36148     0      
##  Shapiro-Wilk Test  R    W      0,91      0      
##  Ljung-Box Test     R    Q(10)  5,1       0,88   
##  Ljung-Box Test     R    Q(15)  8,8       0,89   
##  Ljung-Box Test     R    Q(20)  10        0,96   
##  Ljung-Box Test     R^2  Q(10)  0,71      1      
##  Ljung-Box Test     R^2  Q(15)  1,1       1      
##  Ljung-Box Test     R^2  Q(20)  1,3       1      
##  LM Arch Test       R    TR^2   0,99      1      
## 
## Information Criterion Statistics:
##  AIC  BIC  SIC HQIC 
##  5,2  5,2  5,2  5,2
summary(gfit.norm)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(data = rlogATIVO, cond.dist = "norm") 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000002c7762d8>
##  [data = rlogATIVO]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##       mu     omega    alpha1     beta1  
## 0,052043  3,448404  0,156710  0,592505  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu       0,05204     0,08271    0,629    0,529    
## omega    3,44840     0,73274    4,706 2,52e-06 ***
## alpha1   0,15671     0,03458    4,532 5,83e-06 ***
## beta1    0,59251     0,07446    7,957 1,78e-15 ***
## ---
## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
## 
## Log Likelihood:
##  -4641    normalized:  -2,7 
## 
## Description:
##  Mon Jan 11 20:49:22 2021 by user: Livia Rodrigues 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value
##  Jarque-Bera Test   R    Chi^2  20221     0      
##  Shapiro-Wilk Test  R    W      0,92      0      
##  Ljung-Box Test     R    Q(10)  6,5       0,77   
##  Ljung-Box Test     R    Q(15)  11        0,78   
##  Ljung-Box Test     R    Q(20)  12        0,91   
##  Ljung-Box Test     R^2  Q(10)  0,39      1      
##  Ljung-Box Test     R^2  Q(15)  0,78      1      
##  Ljung-Box Test     R^2  Q(20)  1         1      
##  LM Arch Test       R    TR^2   0,54      1      
## 
## Information Criterion Statistics:
##  AIC  BIC  SIC HQIC 
##  5,4  5,4  5,4  5,4
gfit.std_garch21 <- (garchFit(formula = ~ garch(2, 1),data=rlogATIVO, cond.dist="std")) 
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(2, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               2 1
##  Max GARCH Order:           2
##  Maximum Order:             2
##  Conditional Dist:          std
##  h.start:                   3
##  llh.start:                 1
##  Length of Series:          1733
##  Recursion Init:            mci
##  Series Scale:              3,7
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                   U      V params includes
##     mu     -3,1e-01   0,31  0,031     TRUE
##     omega   1,0e-06 100,00  0,100     TRUE
##     alpha1  1,0e-08   1,00  0,050     TRUE
##     alpha2  1,0e-08   1,00  0,050     TRUE
##     gamma1 -1,0e+00   1,00  0,100    FALSE
##     gamma2 -1,0e+00   1,00  0,100    FALSE
##     beta1   1,0e-08   1,00  0,800     TRUE
##     delta   0,0e+00   2,00  2,000    FALSE
##     skew    1,0e-01  10,00  1,000    FALSE
##     shape   1,0e+00  10,00  4,000     TRUE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1 alpha2  beta1  shape 
##      1      2      3      4      7     10 
##  Persistence:                  0,9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     2222.8070: 0.0310144 0.100000 0.0500000 0.0500000 0.800000  4.00000
##   1:     2222.4938: 0.0310117 0.100823 0.0534403 0.0524567 0.801216  4.00013
##   2:     2222.2807: 0.0310076 0.0974376 0.0549417 0.0524897 0.798703  4.00021
##   3:     2222.0100: 0.0309952 0.0961080 0.0628055 0.0565108 0.798400  4.00063
##   4:     2221.7881: 0.0309670 0.0921207 0.0675178 0.0547119 0.792276  4.00135
##   5:     2221.5805: 0.0309188 0.0967252 0.0737099 0.0523801 0.788972  4.00273
##   6:     2221.5023: 0.0308368 0.0996138 0.0777206 0.0486644 0.783492  4.00495
##   7:     2221.4384: 0.0306698 0.100387 0.0800762 0.0522623 0.780275  4.00962
##   8:     2221.3971: 0.0304785 0.101446 0.0791874 0.0536399 0.776751  4.01471
##   9:     2221.3566: 0.0302784 0.104188 0.0789611 0.0533119 0.775084  4.02001
##  10:     2220.9562: 0.0250851 0.0831223 0.0992615 0.0254797 0.802558  4.15258
##  11:     2220.7209: 0.0197008 0.0862359 0.0945800 0.0247030 0.790081  4.28552
##  12:     2220.0791: 0.0143268 0.103067 0.0810663 0.0705785 0.751131  4.40337
##  13:     2219.8564: 0.0126074 0.109568 0.0991528 0.0349430 0.760994  4.37661
##  14:     2219.7898: 0.0126062 0.108045 0.0982015 0.0342373 0.759776  4.37659
##  15:     2219.7830: 0.0125398 0.108164 0.0982071 0.0347154 0.760134  4.37708
##  16:     2219.7738: 0.0124045 0.107567 0.0976543 0.0350474 0.760131  4.37810
##  17:     2219.6173: 0.00563138 0.104880 0.0879892 0.0354498 0.769427  4.42947
##  18:     2219.5808: -0.00122113 0.117731 0.101005 0.0456592 0.736458  4.44272
##  19:     2219.5348: -0.000917866 0.110527 0.0974422 0.0461932 0.747988  4.39403
##  20:     2219.5001: -0.000474187 0.110530 0.0905469 0.0469619 0.752219  4.39107
##  21:     2219.4970: -0.000213925 0.109556 0.0900920 0.0468560 0.753768  4.40540
##  22:     2219.4963: 7.15144e-05 0.108447 0.0894917 0.0468781 0.755789  4.40215
##  23:     2219.4963: 0.000129021 0.108354 0.0894093 0.0468456 0.755768  4.40689
##  24:     2219.4962: 7.57132e-05 0.108461 0.0894921 0.0468406 0.755677  4.40502
##  25:     2219.4962: 7.90654e-05 0.108456 0.0894897 0.0468433 0.755687  4.40490
## 
## Final Estimate of the Negative LLH:
##  LLH:  4486    norm LLH:  2,6 
##      mu   omega  alpha1  alpha2   beta1   shape 
## 0,00029 1,48390 0,08949 0,04684 0,75569 4,40490 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##            mu omega alpha1 alpha2 beta1  shape
## mu     -216,9     4    -10     22    36    2,1
## omega     4,0   -78   -543   -568  -844  -19,7
## alpha1  -10,0  -543  -6572  -6276 -6955 -180,8
## alpha2   21,6  -568  -6276  -6800 -7375 -186,9
## beta1    36,5  -844  -6955  -7375 -9880 -234,8
## shape     2,1   -20   -181   -187  -235  -10,2
## attr(,"time")
## Time difference of 0,11 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0,31 secs

Segundo o Critério de Informação Bayesiano (BIC), o modelo mais adequado entre os dois testados é aquele que usa a distribução t-student, pois apresentou o menor valor segundo este critério. O modelo GARCH(2,1) com distribuicao t-student apresentou o mesmo BIC, todavia iremos adotar o modelo GARCH(1,1)

coef(gfit.std)
##     mu  omega alpha1  beta1  shape 
## -0,004  1,201  0,114  0,798  4,410
# plot in-sample volatility estimates
length(gfit.std@sigma.t)
## [1] 1733
 str(gfit.std@sigma.t)
##  num [1:1733] 3,7 3,49 3,35 3,48 3,3 ...
vol_garch.std = rlogATIVO*0 + gfit.std@sigma.t
vol_garch.norm = rlogATIVO*0 + gfit.norm@sigma.t
vol_garch.21 =  rlogATIVO*0 +gfit.std_garch21@sigma.t
 

  colnames(vol_garch.std) <- "GARCH_STD"
  colnames(vol_garch.norm) <-  "GARCH_NORM"
  colnames(vol_garch.21) <- "GARCH_STD_21"
  rlog_estima <-   cbind(rlog_estima,vol_garch.std,vol_garch.norm,vol_garch.21)
  
  highchart(type = "stock") %>% 
    hc_add_series((rlog_estima$ELET3.SA.Adjusted), type = "line",color = "red",name = "LogRetorno") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_HIST), type = "line",color = "blue",name = "Vol Historica") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_50), type = "line",color = "orange",name = "EWMA 50%") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_94), type = "line",color = "pink",name = "EWMA 94%") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_96), type = "line",color = "darkblue",name = "EWMA 96%") %>%
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_EWMA_98), type = "line",color = "black",name = "EWMA 98%") %>%
    hc_add_series(qnorm(0.05)*(rlog_estima$GARCH_STD), type = "line",name = "GARCH(1,1) NORMAL") %>%
    hc_add_series(qnorm(0.05)*(rlog_estima$GARCH_NORM), type = "line",name = "GARCH(1,1) STUDENT") %>%
    hc_add_series(qnorm(0.05)*(rlog_estima$GARCH_STD_21), type = "line",name = "GARCH(2,1) STUDENT") %>%
        hc_add_series(qnorm(0.05)*(rlog_estima$VOL_MM21), type = "line",name = "Vol Média Móvel 21") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_MM63  ), type = "line",name = "Vol Média Móvel 100") %>% 
    hc_add_series(qnorm(0.05)*(rlog_estima$VOL_MM252), type = "line",name = "Vol Média Móvel 252") %>%
    hc_title(text = paste0("Log Retorno Preço de ajuste - Vol Comparativo ",ATIVO_input)) %>% 
         hc_yAxis(opposite = FALSE, labels = list(format = "{value}%")) %>% 
      hc_exporting(enabled = TRUE) %>%
      hc_legend(enabled = TRUE)
# Breve resumo do modelo
gfit.std@description
## [1] "Mon Jan 11 20:49:22 2021 by user: Livia Rodrigues"

Conclusão

O objetivo do projeto é analisar os métodos e resultados dos modelos de estimação de volatilidade do ativo ELET3. Mostramos graficamente a comparação do histórico das perdas estimadas pelo VaR com os retornos observados na série, com o objetivo de avaliar se os modelos de VaR testados estão ajustados e qual deles possibilitam uma melhor estimativa para as perdas. Cada modelo apresentando possui uma limitação e a escolha entre eles pode ser feita usando backtesting, por exemplo, o de Kupiec. Podemos verificar que utilizar o desvio padrão (volatilidade) do período, não captura as mudanças de volatilidade dentre de um determinado período. Para contornar a situação da volatilidade histórica, utilizamos o modelo de médias móveis com janelas diferentes de 21 dias, 63 dias e 252 dias. Esse modelo depende da janela adotada, assim podendo superestimar a volatilidade em um período de calmaria ou subestimar num período de estresse. Uma variação do modelo de médias móveis é o EWMA (exponentially weighted moving average), no qual podemos definir um parâmetro de peso podendo colocar maior peso nas observações recente e diminuir a significância das observações passadas. Um ponto importante desse modelo é que utilizamos peso no retorno, assim, um retorno positivo pode penalizar a estimação de volatilidade.

Todos os modelos anteriores não possuem fundamento estatístico teórico, são simplificações matemáticas para calcular a volatilidade. Assim, desenvolvemos o modelo GARCH (heterocedasticidade Condicional Autorregressiva Geral) e o GARCH(1,1) com uma t-student ajustou bem os dados do preço do papel, mesmo fazendo variações GARC(2,1).

  1. UERJ↩