bovespa_arch

Especificando as bibliotecas

library(yfR)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggthemes)
library(e1071)
library(FinTS)

## Carregando pacotes exigidos: zoo
## 
## Anexando pacote: 'zoo'
## 
## Os seguintes objetos são mascarados por 'package:base':
## 
##     as.Date, as.Date.numeric

library(WriteXLS)
library(xtable)
library(tbl2xts)
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## 
## Anexando pacote: 'forecast'
## 
## O seguinte objeto é mascarado por 'package:FinTS':
## 
##     Acf

library(tseries)
library(timeSeries)

## Carregando pacotes exigidos: timeDate
## 
## Anexando pacote: 'timeDate'
## 
## O seguinte objeto é mascarado por 'package:xtable':
## 
##     align
## 
## Os seguintes objetos são mascarados por 'package:e1071':
## 
##     kurtosis, skewness
## 
## 
## Anexando pacote: 'timeSeries'
## 
## O seguinte objeto é mascarado por 'package:zoo':
## 
##     time<-
## 
## O seguinte objeto é mascarado por 'package:dplyr':
## 
##     lag
## 
## Os seguintes objetos são mascarados por 'package:graphics':
## 
##     lines, points

library(xts)

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Anexando pacote: 'xts'
## 
## Os seguintes objetos são mascarados por 'package:dplyr':
## 
##     first, last

library(corrplot)

## corrplot 0.95 loaded

library(knitr)
library(kableExtra)

## 
## Anexando pacote: 'kableExtra'
## 
## O seguinte objeto é mascarado por 'package:dplyr':
## 
##     group_rows

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(rugarch)

## Carregando pacotes exigidos: parallel
## 
## Anexando pacote: 'rugarch'
## 
## O seguinte objeto é mascarado por 'package:purrr':
## 
##     reduce
## 
## O seguinte objeto é mascarado por 'package:stats':
## 
##     sigma

library(WriteXLS)
library(tbl2xts)

rm(list=ls())

tickers = c("^BVSP","EWZ","BRAX11.SA")
indices_dados <- yfR::yf_get(tickers,
                             first_date = '2023-01-01',
                             last_date = '2025-01-01',
                             type_return = "log",
                             freq_data = "daily")

## 

## ── Running yfR for 3 stocks | 2023-01-01 --> 2025-01-01 (731 days) ──

## 

## ℹ Downloading data for benchmark ticker ^GSPC

## ℹ (1/3) Fetching data for ^BVSP

## !    - not cached

## ✔    - cache saved successfully

## ✔    - got 499 valid rows (2023-01-02 --> 2024-12-30)

## ✔    - got 97% of valid prices -- Youre doing good!

## ℹ (2/3) Fetching data for BRAX11.SA

## !    - not cached

## ✔    - cache saved successfully

## ✔    - got 499 valid rows (2023-01-02 --> 2024-12-30)

## ✔    - got 97% of valid prices -- Good job agariel!

## ℹ (3/3) Fetching data for EWZ

## !    - not cached

## ✔    - cache saved successfully

## ✔    - got 502 valid rows (2023-01-03 --> 2024-12-31)

## ✔    - got 100% of valid prices -- Good job agariel!

## ℹ Binding price data

## 

## ── Diagnostics ─────────────────────────────────────────────────────────────────

## ✔ Returned dataframe with 1500 rows -- Well done agariel!

## ℹ Using 209.9 kB at /tmp/RtmpFv2MSr/yf_cache for 4 cache files

## ℹ Out of 3 requested tickers, you got 3 (100%)

visualizando os dados

str(indices_dados)

## tibble [1,500 × 11] (S3: tbl_df/tbl/data.frame)
##  $ ticker                : chr [1:1500] "^BVSP" "^BVSP" "^BVSP" "^BVSP" ...
##  $ ref_date              : Date[1:1500], format: "2023-01-02" "2023-01-03" ...
##  $ price_open            : num [1:1500] 109734 106377 104167 105336 107642 ...
##  $ price_high            : num [1:1500] 109734 106684 105627 107743 109433 ...
##  $ price_low             : num [1:1500] 105981 103852 103915 105333 107642 ...
##  $ price_close           : num [1:1500] 106376 104166 105334 107518 108836 ...
##  $ volume                : num [1:1500] 8130500 14466700 14451200 15512000 12626600 ...
##  $ price_adjusted        : num [1:1500] 106376 104166 105334 107518 108836 ...
##  $ ret_adjusted_prices   : num [1:1500] NA -0.021 0.0112 0.0205 0.0122 ...
##  $ ret_closing_prices    : num [1:1500] NA -0.021 0.0112 0.0205 0.0122 ...
##  $ cumret_adjusted_prices: num [1:1500] 1 0.979 0.99 1.011 1.023 ...
##  - attr(*, "df_control")= tibble [3 × 5] (S3: tbl_df/tbl/data.frame)
##   ..$ ticker              : chr [1:3] "^BVSP" "BRAX11.SA" "EWZ"
##   ..$ dl_status           : chr [1:3] "OK" "OK" "OK"
##   ..$ n_rows              : int [1:3] 499 499 502
##   ..$ perc_benchmark_dates: num [1:3] 0.966 0.966 1
##   ..$ threshold_decision  : chr [1:3] "KEEP" "KEEP" "KEEP"

summary(indices_dados)

##     ticker             ref_date            price_open       
##  Length:1500        Min.   :2023-01-02   Min.   :    22.56  
##  Class :character   1st Qu.:2023-07-03   1st Qu.:    32.20  
##  Mode  :character   Median :2024-01-03   Median :   103.61  
##                     Mean   :2023-12-31   Mean   : 40378.01  
##                     3rd Qu.:2024-07-03   3rd Qu.:114801.50  
##                     Max.   :2024-12-31   Max.   :137349.00  
##                                                             
##    price_high          price_low          price_close            volume        
##  Min.   :    22.60   Min.   :    22.26   Min.   :    22.40   Min.   :       0  
##  1st Qu.:    32.45   1st Qu.:    32.00   1st Qu.:    32.25   1st Qu.:    2433  
##  Median :   104.30   Median :   102.97   Median :   103.56   Median :10258000  
##  Mean   : 40646.25   Mean   : 40123.79   Mean   : 40383.24   Mean   :11207226  
##  3rd Qu.:115507.00   3rd Qu.:113980.75   3rd Qu.:114828.25   3rd Qu.:17767750  
##  Max.   :137469.00   Max.   :136664.00   Max.   :137344.00   Max.   :65477800  
##                                                                                
##  price_adjusted      ret_adjusted_prices  ret_closing_prices  
##  Min.   :    22.03   Min.   :-0.0689929   Min.   :-6.899e-02  
##  1st Qu.:    29.30   1st Qu.:-0.0063940   1st Qu.:-6.481e-03  
##  Median :   103.56   Median : 0.0000000   Median : 0.000e+00  
##  Mean   : 40382.37   Mean   : 0.0001658   Mean   : 7.431e-05  
##  3rd Qu.:114828.25   3rd Qu.: 0.0070353   3rd Qu.: 6.987e-03  
##  Max.   :137344.00   Max.   : 0.0516589   Max.   : 5.166e-02  
##                      NA's   :3            NA's   :3           
##  cumret_adjusted_prices
##  Min.   :0.9131        
##  1st Qu.:1.0917        
##  Median :1.1689        
##  Mean   :1.1496        
##  3rd Qu.:1.2056        
##  Max.   :1.3684        
## 

head(indices_dados)

Reestruturando o data frame. Filtrando somente o preco de fechamento, data e ticker.

indices_dados[is.na(indices_dados)] = 0 #lidando com valores NA

a = c(2,1,10,6) #data, ticker, retorno e preço de fechamento 
indices_dados = indices_dados[,a]
names(indices_dados) = c("Data","Ativo","Retorno","Preco")

Plotando os gráficos

Precos

p= ggplot(indices_dados) + geom_line(mapping = aes(x=Data,y=Preco),color='palegreen4') +
  facet_wrap(~Ativo)+
  labs(x = "",y='Preco Fechamento',title="Cotacao Diaria dos Ativos",
         subtitle = "Periodo: de 01/01/2023 a 01/01/2025", 
         caption = "Fonte: Yahoo Finance")+
  theme_grey()
p

Se quisermos comparar os ativos, visualizando-os no mesmo grafico

p1 = ggplot(indices_dados,aes(x=Data,y=Preco, color = factor(Ativo))) + geom_line()+
  labs(x = "",y='Preco Fechamento',title="Cotacao Diaria dos Ativos",
         subtitle = "Periodo: de 01/01/2014 a 10/01/2024", color= "Ativo",
         caption = "Fonte: B3") + theme_gray()
p1

Enquanto a bovespa é um índice de preços, os outros ativos são ETFs, e a escala é desproporcional. Usandos os log retornos é possível compará-los adequadamente.

Retornos diarios

p3 = ggplot(indices_dados) + geom_line(mapping = aes(x=Data,y=Retorno),color='tomato') +
  facet_wrap(~Ativo)+
  labs(x = "",y='Retornos',title="Retornos Diarios dos Ativos",
         subtitle = "Periodo: de 01/01/2014 a 10/01/2024", 
         caption = "Fonte: B3")+
  theme_grey()
p3

Retornos absolutos

p4 = ggplot(indices_dados) + geom_line(mapping = aes(x=Data,y=abs(Retorno)),color='orchid3') +
  facet_wrap(~Ativo)+
  labs(x = "",y='Retornos absolutos',title="Retornos Absolutos Diarios dos Ativos",
         subtitle = "Periodo: de 01/01/2014 a 10/01/2024", 
         caption = "Fonte: B3")+
  theme_grey()
p4

QQplot

p5=ggplot(indices_dados,aes(sample = Retorno)) + stat_qq()+ stat_qq_line() +
  facet_wrap(~Ativo)+
  labs(x = "Teorico",y='Amostra',title="QQplot",
         subtitle = "Periodo: de 01/01/2014 a 10/01/2024", 
         caption = "Fonte: B3")+
  theme_grey()
p5

Histograma

p6=ggplot(indices_dados) +
  geom_histogram(aes(x=Retorno, y = after_stat(density)),
                 color="grey60", fill="mintcream",linetype="solid",alpha = 0.8) +
  geom_density(aes(x = Retorno,y = after_stat(density)),color="black") +
  facet_wrap(~Ativo)+
  labs(x = "",y='Densidade',title="Histogramas",
         subtitle = "Periodo: de 01/01/2014 a 10/01/2024", 
         caption = "Fonte: B3")+
  theme_grey()
p6

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A fim de calcular as variáveis descritivas, vamos separar o Data Frame em que cada coluna é o log retorno de um ativo.

# Listar os ativos (atualizando com os ativos do seu código)
ativos = unique(indices_dados$Ativo)

# Criar um novo DataFrame com o log-retorno de cada ativo
atv.dados = indices_dados %>% 
  select(Data, Ativo, Retorno) %>% 
  spread(key = Ativo, value = Retorno)  # Espalhando as colunas para cada ativo

# Conferindo como ficou o DataFrame
head(atv.dados)

Valores mínimo, máximo, média, mediana e quartis.

summary(atv.dados)

##       Data                ^BVSP              BRAX11.SA         
##  Min.   :2023-01-02   Min.   :-0.0319903   Min.   :-0.0326752  
##  1st Qu.:2023-07-03   1st Qu.:-0.0054890   1st Qu.:-0.0050768  
##  Median :2024-01-02   Median :-0.0000580   Median : 0.0000939  
##  Mean   :2024-01-01   Mean   : 0.0002462   Mean   : 0.0002570  
##  3rd Qu.:2024-07-02   3rd Qu.: 0.0064691   3rd Qu.: 0.0060662  
##  Max.   :2024-12-31   Max.   : 0.0419842   Max.   : 0.0423149  
##                       NA's   :17           NA's   :17          
##       EWZ            
##  Min.   :-0.0689929  
##  1st Qu.:-0.0088478  
##  Median : 0.0000000  
##  Mean   :-0.0002787  
##  3rd Qu.: 0.0088597  
##  Max.   : 0.0516588  
##  NA's   :14

atv.dados[is.na(atv.dados)] = 0
summary(atv.dados)

##       Data                ^BVSP              BRAX11.SA         
##  Min.   :2023-01-02   Min.   :-0.0319903   Min.   :-0.0326752  
##  1st Qu.:2023-07-03   1st Qu.:-0.0052995   1st Qu.:-0.0048816  
##  Median :2024-01-02   Median : 0.0000000   Median : 0.0000000  
##  Mean   :2024-01-01   Mean   : 0.0002381   Mean   : 0.0002486  
##  3rd Qu.:2024-07-02   3rd Qu.: 0.0062688   3rd Qu.: 0.0058338  
##  Max.   :2024-12-31   Max.   : 0.0419842   Max.   : 0.0423149  
##       EWZ            
##  Min.   :-0.0689929  
##  1st Qu.:-0.0085956  
##  Median : 0.0000000  
##  Mean   :-0.0002711  
##  3rd Qu.: 0.0086657  
##  Max.   : 0.0516588

Desvio padrão, variância, curtose e assimetria.

resultados <- sapply(atv.dados, function(x) {
  if(is.numeric(x)) {
    list(
      desvio_padrao = sd(x, na.rm = TRUE),
      variancia = var(x, na.rm = TRUE),
      curtose = kurtosis(x, na.rm = TRUE),
      assimetria = skewness(x, na.rm = TRUE)
    )
  }
})
print(resultados)

## $Data
## NULL
## 
## $`^BVSP`
## $`^BVSP`$desvio_padrao
## [1] 0.009315028
## 
## $`^BVSP`$variancia
## [1] 8.676974e-05
## 
## $`^BVSP`$curtose
## [1] 0.8474032
## attr(,"method")
## [1] "excess"
## 
## $`^BVSP`$assimetria
## [1] 0.1076411
## attr(,"method")
## [1] "moment"
## 
## 
## $BRAX11.SA
## $BRAX11.SA$desvio_padrao
## [1] 0.009134101
## 
## $BRAX11.SA$variancia
## [1] 8.34318e-05
## 
## $BRAX11.SA$curtose
## [1] 1.154936
## attr(,"method")
## [1] "excess"
## 
## $BRAX11.SA$assimetria
## [1] -0.006789242
## attr(,"method")
## [1] "moment"
## 
## 
## $EWZ
## $EWZ$desvio_padrao
## [1] 0.01490216
## 
## $EWZ$variancia
## [1] 0.0002220743
## 
## $EWZ$curtose
## [1] 0.9773584
## attr(,"method")
## [1] "excess"
## 
## $EWZ$assimetria
## [1] -0.228032
## attr(,"method")
## [1] "moment"

ajustar_modelos_garch <- function(serie_retornos, nome_ativo) {
  modelos <- list(
    garch.norm = ugarchspec(
      variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
      mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
      distribution.model = "norm"
    ),
    garch.std = ugarchspec(
      variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
      mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
      distribution.model = "std"
    ),
    egarch.norm = ugarchspec(
      variance.model = list(model = "eGARCH", garchOrder = c(1,1)),
      mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
      distribution.model = "norm"
    ),
    egarch.std = ugarchspec(
      variance.model = list(model = "eGARCH", garchOrder = c(1,1)),
      mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
      distribution.model = "std"
    ),
    gjr.norm = ugarchspec(
      variance.model = list(model = "gjrGARCH", garchOrder = c(1,1)),
      mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
      distribution.model = "norm"
    ),
    gjr.std = ugarchspec(
      variance.model = list(model = "gjrGARCH", garchOrder = c(1,1)),
      mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
      distribution.model = "std"
    )
  )

  resultados <- list()
  for (nome_modelo in names(modelos)) {
    cat("Ajustando modelo:", nome_modelo, "para", nome_ativo, "\n")
    fit <- ugarchfit(spec = modelos[[nome_modelo]], data = serie_retornos)
    resultados[[nome_modelo]] <- fit
  }

  return(resultados)
}

modelo_bvsp <- ajustar_modelos_garch(atv.dados$`^BVSP`, "^BVSP")

## Ajustando modelo: garch.norm para ^BVSP 
## Ajustando modelo: garch.std para ^BVSP 
## Ajustando modelo: egarch.norm para ^BVSP 
## Ajustando modelo: egarch.std para ^BVSP 
## Ajustando modelo: gjr.norm para ^BVSP 
## Ajustando modelo: gjr.std para ^BVSP

modelo_ewz <- ajustar_modelos_garch(atv.dados$EWZ, "EWZ")

## Ajustando modelo: garch.norm para EWZ 
## Ajustando modelo: garch.std para EWZ 
## Ajustando modelo: egarch.norm para EWZ 
## Ajustando modelo: egarch.std para EWZ 
## Ajustando modelo: gjr.norm para EWZ 
## Ajustando modelo: gjr.std para EWZ

modelo_brax <- ajustar_modelos_garch(atv.dados$`BRAX11.SA`, "BRAX11.SA")

## Ajustando modelo: garch.norm para BRAX11.SA 
## Ajustando modelo: garch.std para BRAX11.SA 
## Ajustando modelo: egarch.norm para BRAX11.SA 
## Ajustando modelo: egarch.std para BRAX11.SA 
## Ajustando modelo: gjr.norm para BRAX11.SA 
## Ajustando modelo: gjr.std para BRAX11.SA

comparar_modelos <- function(modelos_fit) {
  sapply(modelos_fit, function(fit) {
    c(
      LogLikelihood = likelihood(fit),
      AIC = infocriteria(fit)[1],
      BIC = infocriteria(fit)[2],
      Persistence = persistence(fit)
    )
  }) %>% t() %>% round(4)
}

comparar_modelos(modelo_bvsp)

##             LogLikelihood     AIC     BIC Persistence
## garch.norm       1690.235 -6.5397 -6.5150      0.9926
## garch.std        1691.678 -6.5414 -6.5085      0.9922
## egarch.norm      1692.820 -6.5458 -6.5129      0.9955
## egarch.std       1693.484 -6.5445 -6.5034      0.9953
## gjr.norm         1694.087 -6.5507 -6.5178      0.9990
## gjr.std          1694.549 -6.5486 -6.5075      0.9990

comparar_modelos(modelo_ewz)

##             LogLikelihood     AIC     BIC Persistence
## garch.norm       1438.867 -5.5654 -5.5407      0.9926
## garch.std        1443.841 -5.5808 -5.5479      0.9977
## egarch.norm      1440.808 -5.5690 -5.5361      0.1757
## egarch.std       1445.276 -5.5825 -5.5413      0.1127
## gjr.norm         1438.838 -5.5614 -5.5285      0.9974
## gjr.std          1443.809 -5.5768 -5.5356      0.9985

comparar_modelos(modelo_brax)

##             LogLikelihood     AIC     BIC Persistence
## garch.norm       1701.410 -6.5830 -6.5583      0.9914
## garch.std        1704.556 -6.5913 -6.5584      0.9922
## egarch.norm      1703.455 -6.5870 -6.5541      0.9945
## egarch.std       1705.493 -6.5911 -6.5499      0.9950
## gjr.norm         1705.171 -6.5937 -6.5608      0.9988
## gjr.std          1706.752 -6.5959 -6.5548      0.9986

O melhor modelo para o Ibovespa

forecast_garch <- ugarchforecast(modelo_bvsp$gjr.std, n.ahead = 10)
plot(forecast_garch, which = 1)  # Série

plot(forecast_garch, which = 3)  # Sigma

O modelo prevê volatilidade elevada e estável nos dias subsequentes, com retornos médios próximos de zero.

O melhor modelo para o MSCI Brazil

forecast_garch <- ugarchforecast(modelo_ewz$egarch.std, n.ahead = 10)
plot(forecast_garch, which = 1)  # Série

plot(forecast_garch, which = 3)  # Sigma

O modelo prevê maior incerteza do que a modelagem para o ibovespa. O fato de ser Egarch captura a assimetria de picos negativos, e a incerteza pode estar relacionada à influência de eventos internacionais.

O melhor modelo para o índice IBrX-100

forecast_garch <- ugarchforecast(modelo_brax$gjr.std, n.ahead = 10)
plot(forecast_garch, which = 1)  # Série

plot(forecast_garch, which = 3)  # Sigma

Nos retornos previstos, as bandas de incerteza são semelhantes às do ibovespa, e o histórico de falta de picos indica volatilidade estável, pelo menos no curto prazo. A previsão da volatilidade condicional indica um valor menor que o dos outros índices, o que sugere um comportamento mais previsível. Estruturlmente, o IBrX-100 é mais diversificado, e por isso menos concentrado em empresas como a petrobras e vale. Isso reforça a ideia de que ele pode ser usado como índice de referência do mercado brasileiro.