Lab 3 (R): Séries Financeiras, GARCH e Portfólios

Mesma trilha do lab oficial, em R: fatos estilizados, GARCH, fronteira, VaR e CAPM

Objetivos do Lab

  • Aplicar download e tratamento de séries financeiras com yfR
  • Analisar fatos estilizados em fundos imobiliários brasileiros
  • Criar modelos GARCH (rugarch) e portfólios pela fronteira eficiente
  • Avaliar risco e exposição via VaR e betas (CAPM)
NotaQual versão usar?

Este é o espelho em R do laboratório, fiel ao script de referência. A versão oficial em Python está em Lab 3 oficial, e há ainda uma versão alternativa com blue-chips e fronteira por otimização. Os três cobrem o mesmo pipeline: escolha a linguagem com que você se sente mais à vontade.

DicaEstrutura do Lab
  • Bloco 3.1: séries financeiras e GARCH. Preços, log-retornos, fatos estilizados, grid de modelos GARCH, previsão de volatilidade.
  • Bloco 3.2: portfólio. Retorno e risco por ativo, fronteira eficiente por Monte Carlo, VaR via t de Student e betas/CAPM contra o Ibovespa.

O pool são quatro FIIs que o monitor já teve na carteira (HGRE11, BTLG11, HGRU11, VGIR11). A pergunta de negócio é honesta: foi um péssimo investimento ou apenas ruim? FIIs costumam ter beta baixo contra o Ibovespa, então a parte de CAPM rende boa discussão.

Setup compartilhado

Código 
library(fpp3)
library(rugarch)
library(yfR)
library(dplyr)
library(tidyr)
library(purrr)
library(tibble)
library(ggplot2)
library(patchwork)
library(timetk)
library(broom)
library(scales)
library(forcats)
Código 
start_date <- "2018-01-01"
# Fundos imobiliários que o monitor já teve na carteira.
ativos <- c("HGRE11.SA", "BTLG11.SA", "HGRU11.SA", "VGIR11.SA")

Trabalhamos tanto com os dados em formato de tibble quanto de tsibble.

Código 
da <- yfR::yf_get(
  ativos,
  first_date = start_date,
  type_return = "log",
  freq_data = "daily",
  do_complete_data = TRUE
)

da_tsibble <- da |>
  as_tsibble(key = ticker, index = ref_date, regular = FALSE)

# Data mínima comum a todas as séries (alguns FIIs listaram depois de 2018)
data_corte <- da |>
  dplyr::group_by(ticker) |>
  dplyr::filter(ref_date == min(ref_date)) |>
  dplyr::ungroup() |>
  with(max(ref_date))

da_train <- da |>
  dplyr::filter(ref_date > data_corte)

data_corte
[1] "2018-07-31"

Bloco 3.1: Séries Financeiras e GARCH

Parte 1: Preços e log-retornos

Código 
da_tsibble |>
  autoplot(price_close, colour = "black") +
  facet_wrap(~ticker, scales = "free_y", ncol = 1) +
  labs(x = NULL, y = "Preço de fechamento")

Código 
da_tsibble |>
  autoplot(ret_adjusted_prices, colour = "black") +
  facet_wrap(~ticker, scales = "free_y", ncol = 1) +
  labs(x = NULL, y = "Log-retorno")

Parte 2: Fatos estilizados

ACF e PACF dos retornos e dos retornos ao quadrado, mais a checagem de normalidade.

Código 
da_tsibble |>
  ACF(ret_adjusted_prices) |>
  autoplot() +
  labs(title = "ACF dos retornos")

Código 
da_tsibble |>
  dplyr::mutate(ret2 = ret_adjusted_prices^2) |>
  ACF(ret2) |>
  autoplot() +
  labs(title = "ACF dos retornos ao quadrado")

Código 
da_tsibble |>
  dplyr::mutate(ret2 = ret_adjusted_prices^2) |>
  PACF(ret2) |>
  autoplot() +
  labs(title = "PACF dos retornos ao quadrado")

  1. Os retornos simples têm pouca autocorrelação (mercado eficiente), mas os retornos ao quadrado têm autocorrelação forte e persistente: é o sinal de clusters de volatilidade.
  2. A ACF dos retornos ao quadrado mostra que a variância muda ao longo do tempo de forma previsível: períodos de alta volatilidade puxam mais alta volatilidade.
Código 
da_train |>
  ggplot(aes(x = ret_adjusted_prices)) +
  geom_histogram(bins = 90) +
  facet_wrap(~ticker, ncol = 2) +
  labs(x = "Log-retorno", y = NULL, title = "Distribuição dos retornos")

QQ-plot contra a normal e contra a t de Student (3 graus de liberdade).

Código 
da_train |>
  group_by(ticker) |>
  summarise(
    gg = list(
      ggplot(pick(everything()), aes(sample = ret_adjusted_prices)) +
        geom_qq() +
        geom_qq_line() +
        labs(title = cur_group()$ticker, x = NULL, y = NULL)
    )
  ) |>
  dplyr::pull(gg) |>
  patchwork::wrap_plots() +
  patchwork::plot_annotation(title = "QQ-plot vs. Normal")

Código 
da_train |>
  group_by(ticker) |>
  summarise(
    gg = list(
      ggplot(pick(everything()), aes(sample = ret_adjusted_prices)) +
        geom_qq(distribution = qt, dparams = list(df = 3)) +
        geom_qq_line(distribution = qt, dparams = list(df = 3)) +
        labs(title = cur_group()$ticker, x = NULL, y = NULL)
    )
  ) |>
  dplyr::pull(gg) |>
  patchwork::wrap_plots() +
  patchwork::plot_annotation(title = "QQ-plot vs. t de Student (df = 3)")

Parte 3: Modelos GARCH

Função que ajusta um GARCH com a especificação do rugarch. A média é um ARMA simples (p, q) e a variância um GARCH de ordem (m, n).

Código 
garch_individual <- function(parms, ret, prog = NULL) {
  if (!is.null(prog)) {
    prog()
  }
  garch_model <- ugarchspec(
    variance.model = list(
      model = "fGARCH",
      submodel = "GARCH",
      garchOrder = c(parms$m, parms$n)
    ),
    mean.model = list(
      armaOrder = c(parms$p, parms$q),
      include.mean = TRUE
    ),
    distribution.model = parms$dist
  )
  # às vezes não converge
  suppressWarnings({
    fit <- tryCatch(
      ugarchfit(garch_model, data = ret),
      error = function(e) NULL
    )
  })
  fit
}

Testando para um ativo:

Código 
library(rugarch)
fit_exemplo <- garch_individual(
  parms = list(p = 0, q = 0, m = 1, n = 1, dist = "std"),
  ret = da_train |>
    dplyr::filter(ticker == "HGRE11.SA") |>
    pull(ret_adjusted_prices)
)
fit_exemplo

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : fGARCH(1,1)
fGARCH Sub-Model    : GARCH
Mean Model  : ARFIMA(0,0,0)
Distribution    : std 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000067    0.000190  0.35295 0.724127
omega   0.000019    0.000005  3.56818 0.000359
alpha1  0.249339    0.056885  4.38322 0.000012
beta1   0.683382    0.059358 11.51298 0.000000
shape   3.360479    0.299568 11.21775 0.000000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000067    0.000223  0.30119 0.763270
omega   0.000019    0.000006  3.28707 0.001012
alpha1  0.249339    0.059356  4.20071 0.000027
beta1   0.683382    0.067374 10.14315 0.000000
shape   3.360479    0.280052 11.99948 0.000000

LogLikelihood : 6149.992 

Information Criteria
------------------------------------
                    
Akaike       -6.3090
Bayes        -6.2947
Shibata      -6.3090
Hannan-Quinn -6.3038

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.326 0.24959
Lag[2*(p+q)+(p+q)-1][2]     1.342 0.39949
Lag[4*(p+q)+(p+q)-1][5]     6.198 0.08077
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.1116  0.7383
Lag[2*(p+q)+(p+q)-1][5]    1.1823  0.8172
Lag[4*(p+q)+(p+q)-1][9]    5.3641  0.3779
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]   0.04307 0.500 2.000  0.8356
ARCH Lag[5]   1.70636 1.440 1.667  0.5398
ARCH Lag[7]   5.92820 2.315 1.543  0.1461

Nyblom stability test
------------------------------------
Joint Statistic:  1.9573
Individual Statistics:             
mu     0.4980
omega  0.7086
alpha1 0.1785
beta1  0.5146
shape  0.7183

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.28 1.47 1.88
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           1.4060 0.1599    
Negative Sign Bias  0.7681 0.4425    
Positive Sign Bias  1.1200 0.2629    
Joint Effect        2.3524 0.5025    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     21.22       0.3248
2    30     27.78       0.5297
3    40     46.21       0.1990
4    50     58.83       0.1588


Elapsed time : 0.6799998

Um grid de modelos por ativo

Em vez de chutar a ordem e a distribuição, varremos uma grade de (m, n, dist) e escolhemos pelo menor AIC. É o melhor_garch que ficava omitido no script de referência.

Código 
melhor_garch <- function(ativo, prog = NULL) {
  if (!is.null(prog)) {
    prog()
  }
  ret <- da_train |>
    dplyr::filter(ticker == ativo) |>
    dplyr::pull(ret_adjusted_prices)

  grid <- tidyr::expand_grid(
    p = 0,
    q = 0, # média (ARMA) simples
    m = 1:2,
    n = 1:2, # ordem do GARCH
    dist = c("norm", "std") # normal e t de Student
  )

  purrr::pmap(grid, function(p, q, m, n, dist) {
    parms <- list(p = p, q = q, m = m, n = n, dist = dist)
    fit <- garch_individual(parms, ret)
    ic <- tryCatch(rugarch::infocriteria(fit), error = function(e) NULL)
    if (is.null(ic)) {
      return(NULL)
    }
    tibble::tibble(
      p = p,
      q = q,
      m = m,
      n = n,
      dist = dist,
      aic = ic["Akaike", 1],
      bic = ic["Bayes", 1]
    )
  }) |>
    dplyr::bind_rows() |>
    dplyr::arrange(aic)
}
Código 
melhores_por_ativo <- ativos |>
  purrr::set_names() |>
  purrr::map(melhor_garch, .progress = TRUE) |>
  dplyr::bind_rows(.id = "ticker")

# Melhor modelo (menor AIC) por ativo
parametros_melhores <- melhores_por_ativo |>
  dplyr::group_by(ticker) |>
  dplyr::slice_head(n = 1) |>
  dplyr::ungroup()

parametros_melhores
# A tibble: 4 × 8
  ticker        p     q     m     n dist    aic   bic
  <chr>     <dbl> <dbl> <int> <int> <chr> <dbl> <dbl>
1 BTLG11.SA     0     0     1     2 std   -7.16 -7.15
2 HGRE11.SA     0     0     1     1 std   -6.31 -6.29
3 HGRU11.SA     0     0     1     2 std   -6.96 -6.94
4 VGIR11.SA     0     0     2     2 std   -7.26 -7.24

Previsão de volatilidade um passo à frente

Reajustamos o modelo escolhido de cada ativo e projetamos a volatilidade alguns passos à frente. Guardamos o shape (graus de liberdade da t) para reusar no VaR.

Código 
prever_volatilidade <- function(parms, n_steps = 5) {
  ret <- da_train |>
    dplyr::filter(ticker == parms$ticker) |>
    pull(ret_adjusted_prices)

  garch_model <- ugarchspec(
    variance.model = list(
      model = "fGARCH",
      submodel = "GARCH",
      garchOrder = c(parms$m, parms$n)
    ),
    mean.model = list(armaOrder = c(parms$p, parms$q), include.mean = TRUE),
    distribution.model = parms$dist
  )

  fit <- ugarchfit(garch_model, data = ret, out.sample = n_steps - 1)

  shape <- if (parms$dist == "std") {
    as.numeric(fit@fit$coef["shape"])
  } else {
    NA_real_
  }
  forecasts <- ugarchforecast(fit, n.ahead = n_steps)@forecast

  tibble::tibble(
    ticker = parms$ticker,
    passo = seq_len(n_steps),
    serie = as.numeric(forecasts$seriesFor),
    volatilidade = as.numeric(forecasts$sigmaFor),
    shape = shape
  )
}

vol_futuro <- parametros_melhores |>
  dplyr::group_split(ticker) |>
  purrr::map(\(x) prever_volatilidade(as.list(x), n_steps = 5)) |>
  dplyr::bind_rows()

vol_futuro
# A tibble: 20 × 5
   ticker    passo     serie volatilidade shape
   <chr>     <int>     <dbl>        <dbl> <dbl>
 1 BTLG11.SA     1 0.000217       0.00669  2.90
 2 BTLG11.SA     2 0.000217       0.00709  2.90
 3 BTLG11.SA     3 0.000217       0.00749  2.90
 4 BTLG11.SA     4 0.000217       0.00786  2.90
 5 BTLG11.SA     5 0.000217       0.00822  2.90
 6 HGRE11.SA     1 0.0000722      0.00910  3.37
 7 HGRE11.SA     2 0.0000722      0.00981  3.37
 8 HGRE11.SA     3 0.0000722      0.0104   3.37
 9 HGRE11.SA     4 0.0000722      0.0110   3.37
10 HGRE11.SA     5 0.0000722      0.0115   3.37
11 HGRU11.SA     1 0.0000759      0.00790  2.73
12 HGRU11.SA     2 0.0000759      0.00911  2.73
13 HGRU11.SA     3 0.0000759      0.00986  2.73
14 HGRU11.SA     4 0.0000759      0.0106   2.73
15 HGRU11.SA     5 0.0000759      0.0113   2.73
16 VGIR11.SA     1 0.000214       0.00842  3.63
17 VGIR11.SA     2 0.000214       0.00936  3.63
18 VGIR11.SA     3 0.000214      NA        3.63
19 VGIR11.SA     4 0.000214      NA        3.63
20 VGIR11.SA     5 0.000214      NA        3.63

Bloco 3.2: Portfólio

Parte 4: Retorno, risco e fronteira eficiente

Seguimos a receita do codingfinance: retorno médio por ativo, matriz de covariância (em escala diária) e milhares de carteiras aleatórias.

Código 
da_wide <- da_train |>
  dplyr::select(ref_date, name = ticker, value = ret_adjusted_prices) |>
  tidyr::pivot_wider()

da_xts <- da_wide |>
  timetk::tk_xts(select = -ref_date, date_var = ref_date)

mean_ret <- colMeans(da_xts, na.rm = TRUE)
cov_mat <- cov(da_xts, use = "complete.obs")
round(mean_ret, 5)
BTLG11.SA HGRE11.SA HGRU11.SA VGIR11.SA 
  0.00040   0.00021   0.00030   0.00030

Demonstrando em uma única carteira aleatória: pesos, retorno, risco e Sharpe.

Código 
set.seed(2)
wts <- runif(length(ativos))
wts <- wts / sum(wts)

port_returns <- sum(wts * mean_ret)
port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts))
sharpe_ratio <- port_returns / port_risk

tibble::tibble(
  ativo = ativos,
  peso = round(wts, 3)
)
# A tibble: 4 × 2
  ativo      peso
  <chr>     <dbl>
1 HGRE11.SA 0.114
2 BTLG11.SA 0.431
3 HGRU11.SA 0.352
4 VGIR11.SA 0.103
Código 
c(
  retorno = port_returns,
  risco = as.numeric(port_risk),
  sharpe = as.numeric(sharpe_ratio)
) |>
  round(5)
retorno   risco  sharpe 
0.00027 0.00891 0.03065

Agora 5000 carteiras aleatórias. As de menor risco e maior Sharpe são a mínima variância e a carteira de tangência.

Código 
sim_returns <- function(i) {
  wts <- rgamma(length(ativos), 1)
  wts <- wts / sum(wts)
  port_ret <- sum(wts * mean_ret)
  port_sd <- as.numeric(sqrt(t(wts) %*% (cov_mat %*% wts)))

  wts |>
    purrr::set_names(ativos) |>
    tibble::enframe() |>
    tidyr::pivot_wider() |>
    dplyr::mutate(
      return = port_ret,
      risk = port_sd,
      sharpe = port_ret / port_sd
    )
}

portfolio_values <- purrr::map(1:5000, sim_returns) |>
  dplyr::bind_rows(.id = "run")

min_var <- portfolio_values[which.min(portfolio_values$risk), ]
max_sr <- portfolio_values[which.max(portfolio_values$sharpe), ]
Código 
min_var |>
  tidyr::pivot_longer(all_of(ativos)) |>
  dplyr::mutate(name = forcats::fct_reorder(name, value)) |>
  ggplot(aes(name, value)) +
  geom_col(fill = "#3ACC9F") +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "Ativo", y = "Peso", title = "Carteira de mínima variância")

Código 
max_sr |>
  tidyr::pivot_longer(all_of(ativos)) |>
  dplyr::mutate(name = forcats::fct_reorder(name, value)) |>
  ggplot(aes(name, value)) +
  geom_col(fill = "#E50505") +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "Ativo", y = "Peso", title = "Carteira de tangência (máximo Sharpe)")

Código 
portfolio_values |>
  ggplot(aes(x = risk, y = return, color = sharpe)) +
  geom_point(size = 0.6) +
  scale_color_viridis_c() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(
    x = "Risco",
    y = "Retorno",
    title = "Fronteira eficiente por Monte Carlo (5000 carteiras)"
  ) +
  geom_point(
    aes(x = risk, y = return),
    data = min_var,
    color = "#3ACC9F",
    size = 3
  ) +
  geom_point(
    aes(x = risk, y = return),
    data = max_sr,
    color = "#E50505",
    size = 3
  )

Parte 5: VaR do portfólio

Usamos a carteira de mínima variância e o shape (graus de liberdade) mais conservador estimado pelos GARCHs. O VaR sai de um quantil da t de Student.

Código 
pesos_finais <- min_var |>
  dplyr::select(all_of(ativos)) |>
  as.numeric()

# retorno esperado um passo à frente, por ativo, na ordem de `ativos`
serie_1passo <- vol_futuro |>
  dplyr::group_by(ticker) |>
  dplyr::slice_head(n = 1) |>
  dplyr::ungroup() |>
  dplyr::arrange(match(ticker, ativos))

rt_final <- sum(serie_1passo$serie * pesos_finais)
st_dev_final <- as.numeric(sqrt(pesos_finais %*% cov_mat %*% pesos_finais))
nu <- min(vol_futuro$shape, na.rm = TRUE)
valor_t <- qt(0.95, nu)

VaR <- rt_final + valor_t * st_dev_final / sqrt(nu / (nu - 2))
c(rt_final = rt_final, st_dev = st_dev_final, nu = nu, VaR = VaR) |> round(5)
rt_final   st_dev       nu      VaR 
 0.00013  0.00749  2.72562  0.00962
NotaComo ler esse VaR

O VaR resume, num número só, “qual é uma perda plausível num dia ruim”, usando a cauda da t de Student (graus de liberdade baixos = caudas pesadas). O nu precisa ser maior que 2 para a variância existir. Na Aula 6 formalizamos VaR paramétrico, histórico e Monte Carlo.

Parte 6: CAPM e betas

Medimos a exposição ao mercado. O beta de cada ativo vem de cov / var contra o Ibovespa; o beta da carteira é a média ponderada dos betas, e checamos com uma regressão.

Código 
market_returns <- yfR::yf_get(
  "^BVSP",
  first_date = data_corte + 1,
  type_return = "log",
  freq_data = "daily",
  do_complete_data = TRUE
) |>
  dplyr::select(ref_date, ibov = ret_adjusted_prices)

# retornos da carteira = soma ponderada por data (sem tidyquant)
port_wide <- da_train |>
  dplyr::select(ref_date, ticker, ret = ret_adjusted_prices) |>
  tidyr::pivot_wider(names_from = ticker, values_from = ret) |>
  tidyr::drop_na()

port_wide$portfolio <- as.numeric(
  as.matrix(port_wide[, ativos]) %*% pesos_finais
)

all_returns <- market_returns |>
  dplyr::inner_join(
    dplyr::select(port_wide, ref_date, portfolio),
    "ref_date"
  ) |>
  tidyr::drop_na()

beta_geral <- with(all_returns, cov(portfolio, ibov) / var(ibov))

calcular_beta <- function(ativo) {
  da_train |>
    dplyr::filter(ticker == ativo) |>
    dplyr::inner_join(market_returns, "ref_date") |>
    tidyr::drop_na() |>
    with(cov(ret_adjusted_prices, ibov) / var(ibov))
}

betas <- purrr::map_dbl(ativos, calcular_beta) |>
  purrr::set_names(ativos)

round(betas, 3)
HGRE11.SA BTLG11.SA HGRU11.SA VGIR11.SA 
    0.223     0.188     0.213     0.160
Código 
c(beta_ponderado = sum(betas * pesos_finais), beta_regressao = beta_geral) |>
  round(3)
beta_ponderado beta_regressao 
         0.194          0.194
Código 
lm(portfolio ~ ibov, data = all_returns) |>
  broom::tidy()
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept) 0.000206  0.000174      1.18 2.37e- 1
2 ibov        0.194     0.0117       16.6  3.37e-58
ImportanteA resposta da pergunta de negócio

Os betas dos FIIs costumam ficar bem abaixo de 1: eles se movem pouco com o Ibovespa. Isso é exatamente o ponto de ter FIIs na carteira (descorrelação com a bolsa), e não um sinal de mau investimento por si só. “Péssimo ou apenas ruim” depende do retorno ajustado ao risco. Na Aula 7 formalizamos o CAPM e a leitura de beta, alpha e a reta de mercado.

Dinâmica: Monte sua Carteira

A dinâmica final é a mesma do lab oficial: o app interativo Monte sua Carteira (Shiny for Python). Em grupos, escolham ativos, mexam nos pesos e observem retorno, risco, Sharpe, VaR e o backtest contra o Ibovespa ao vivo. O passo a passo está no Lab 3 oficial e em labs/app-carteira/README.md.

De volta ao topo