Trabalho Final - Finantial Analytics

Ativos

Vamos primeiro carregar os ativos selecionados

start_date <- '2020-01-01'

ativos <- c(
  "SIFY",
  "INFY",
  "HAPPSTMNDS.NS",
  "AFFLE.BO",
  "ZENSARTECH.BO",
  "BCG.BO"
)

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)

1. Plotar gráficos dos preços e retornos e tecer comentários sobre possível heterocedasticidade condicional.

Gráfico dos preços e retornos:

Vemos que os preços das ações,apresentaram uma tendência de crescimento e depois decrescimento em algum instante do tempo, com períodos de maior e menor volatilidade.

Também observamos que a nuvem de pontos de retorno se concentra em torno do zero, indicando que a média dos retoros é nula.

Funções de Autocorrelação dos retornos:

#Data mínima comum a todas as séries

data_corte <- da |>
  dplyr::group_by(ticker) |>
  dplyr::filter(ref_date == min(ref_date)) |>
  dplyr::ungroup() |>
  with(max(ref_date))

data_corte

## [1] "2020-09-17"

# "2020-09-17"

#Colocando todos os tickers na mesma data inicial
da_train <- da |>
  dplyr::filter(ref_date > data_corte)

#ACF e PACF dos retornos
da_tsibble |>
  ACF(ret_closing_prices) |>
  autoplot()

da_tsibble |>
  PACF(ret_closing_prices) |>
  autoplot()

Pelo comportamento dos preços vemos que a ACF e PACF de AFFLE, BCG e HAPPSTMNDS apresentam autocorrelação até os dois primeiros lags, decaindo para zero depois.

Normalidade dos retornos

da_train |>
  group_by(ticker) |>
  summarise(
    gg = list(
      ggplot(pick(everything()), aes(sample = ret_closing_prices)) +
        geom_qq() +
        geom_qq_line() +
        labs(title = cur_group())
      )
  ) |>
  dplyr::pull(gg) |>
  patchwork::wrap_plots()

Vemos que os retornos não possuem distruibuição normal, já que suas caudas são mais pesadas que a normal; exceto para BCG, que possui caudas mais leves.

Abaixo vamos analisar se os retornos se ajustam melhor a distribuição t-Student

da_train |>
  group_by(ticker) |>
  summarise(
    gg = list(
      ggplot(pick(everything()), aes(sample = ret_closing_prices)) +
        geom_qq(distribution = qt, dparams = list(df = 3)) +
        geom_qq_line(distribution = qt, dparams = list(df = 3)) +
        labs(title = cur_group())
    )
  ) |>
  dplyr::pull(gg) |>
  patchwork::wrap_plots()

Pelo Q-Plot, vemos que os retornos se ajustaram melhor a distribuição t-Student.

Vamos agora analisar os retornos ao quadrado, que é uma proxy da volatilidade.

Retornos ao quadrado:

da_tsibble |>
  dplyr::mutate(ret2 = ret_closing_prices^2) |>
  autoplot(ret2, colour = "black") +
  facet_wrap(~ticker, ncol = 1)

Dentre as séries financeiras, observamos que SIFY é a que apresenta clusters com picos mais elevados de volatilidade. Tal comportamento também é observado nas outras séries de forma mais discreta, sendo um indício de presença de heterocedasticidade condicional.

Funções de Autocorrelação dos Retornos ao Quadrado:

da_tsibble |>
  dplyr::mutate(ret2 = ret_closing_prices^2) |>
  ACF(ret2) |>
  autoplot()

da_tsibble |>
  dplyr::mutate(ret2 = ret_closing_prices^2) |>
  PACF(ret2) |>
  autoplot()

Tanto na ACF quanto na PACF vemos forte presença de autocorrelação. Desse modo, entendemos que as séries financeiras em questão apresentam alguns os fatos estilizados esperados (média zero e distribuição não normal) e heterocedasticidade condicional.

2. Ajustar modelos de vola lidade univariados e escolher o mais adequado (pode usar critérios de informação e/ou validação cruzada);

Primeiro criamos a função que ajusta o garch

garch_individual <- function(parms, ret, prog = NULL) {
  if (!is.null(prog)) prog()
  # daria para adicionar mais hiperparametros!!!
  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
  )
  # as vezes ele nao converge
  suppressWarnings({
    fit <- ugarchfit(garch_model, data = ret)
  })
  fit
}

Depois criamos um grid para encontrar a melhor especificação de cada GARCH:

# Função para ajustar uma grid de garchs e pegar as informações
melhor_garch <- function(ativo, p = 0:2, q = 0:2, m = 0:2, n = 0:2, dist = "std") {
  
  # faz uma grid de hiperparâmetros
  grid <- expand_grid(p = p, q = q, m = m, n = n, dist = "std") |>
    tibble::rownames_to_column("id")
  
  # pega o retorno do ativo
  ret <- da_train |>
    dplyr::filter(ticker == ativo) |>
    pull(ret_closing_prices)
  
  usethis::ui_info("Ajustando modelos para {ativo}...")
  
  modelos <- grid |>
    group_split(id) |>
    purrr::map(purrr::possibly(garch_individual), ret = ret, .progress = TRUE) %>% 
    purrr::compact()
  
  safe_info <- purrr::possibly(infocriteria, tibble::tibble())
  
  suppressWarnings({
    informacao <- purrr::map(modelos, safe_info) |>
      purrr::map(tibble::as_tibble, rownames = "criteria") |>
      dplyr::bind_rows(.id = "id")
  })
  
  melhores <- informacao |>
    dplyr::inner_join(grid, "id") |>
    tidyr::pivot_wider(names_from = criteria, values_from = V1) |>
    janitor::clean_names() |>
    dplyr::arrange(akaike)
  
  usethis::ui_info(c(
    "Melhor modelo:",
    "p <- {melhores$p[1]}",
    "q <- {melhores$q[1]}",
    "m <- {melhores$m[1]}",
    "n <- {melhores$n[1]}"
  ))
  
  melhores
  
}

Agora vamos rodar as possíveis combinações GARCH para cada ativo:

SIFY:

best_model_sify <- melhor_garch("SIFY") #p,q,m,n 2022

INFY:

best_model_infy <- melhor_garch("INFY")

HAPPSTMNDS.NS:

best_model_happ <-  melhor_garch("HAPPSTMNDS.NS")

AFFLE.BO:

best_model_affle <- melhor_garch("AFFLE.BO")

ZENSARTECH.BO:

best_model_zensartech <- melhor_garch("ZENSARTECH.BO")

Escolhendo as melhores especificações baseando-se no critério AIC

melhores_por_ativo <- ativos |>
  purrr::set_names() |>
  purrr::map(melhor_garch, .progress = TRUE) |>
  dplyr::bind_rows(.id = "ticker")

melhores_por_ativo %>% 
  arrange(akaike) %>% 
  distinct(ticker, .keep_all = TRUE)

# readr::write_rds(melhores_por_ativo, "melhores_por_ativo.rds")
melhores_por_ativo <- readr::read_rds("melhores_por_ativo.rds")

Notamos que para AFFLE o melhor modelo preditivo é um GARCH(0,0), ou seja, mesmo tendo autocorrelação nos retornos quadráticos, ele não apresentou heterocedasticidade condicional no ajuste otimizado do modelo. Por esse motivo, não faremos previsão de volatilidade para esse ativo.

3. Prever a volatilidade um passo à frente usando o modelo selecionado no item anterior

prever_volatilidade <- function(parms, n_steps = 5) {
  
  ret <- da_train |>
    dplyr::filter(ticker == parms$ticker) |>
    pull(ret_closing_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)
  
  if (parms$dist == "std") {
    shape <- as.numeric(fit@fit$coef["shape"])
  } else {
    shape <- NA_real_
  }
 
  forecasts <- ugarchforecast(fit, n.ahead = n_steps)@forecast
  tibble::tibble(
    ticker = parms$ticker,
    serie = as.numeric(forecasts$seriesFor),
    volatilidade = as.numeric(forecasts$sigmaFor),
    shape = shape
  )
}

# Ajustando modelos finais e prevendo volatilidade futura

parametros_melhores <- melhores_por_ativo |>
  group_by(ticker) |>
  slice_head(n = 1) |>
  ungroup()

vol_futuro <- parametros_melhores %>%
  filter(ticker != "AFFLE.BO") %>%  # exclusao do ativo Zensartech que apresentou 0 volatilidade 
  group_split(ticker) %>%
  #purrr::map(~ head(.x, 6)) %>%          # Apply head(6) after splitting
  purrr::map(prever_volatilidade, n_steps = 1, .progress = TRUE) %>%
  dplyr::bind_rows()

vol_futuro

## # A tibble: 5 × 4
##   ticker            serie volatilidade shape
##   <chr>             <dbl>        <dbl> <dbl>
## 1 BCG.BO        -0.00527        0.0409 15.9 
## 2 HAPPSTMNDS.NS -0.00112        0.0149  3.98
## 3 INFY           0.000395       0.0154  6.91
## 4 SIFY          -0.00662        0.0412  4.36
## 5 ZENSARTECH.BO -0.000141       0.0237  3.78

4. Comparar volatilidades entre os retornos selecionados (quais são maiores e menores, relacionando com algum storytelling);

## Comparar volatilidades entre os retornos selecionados

da_wide <- da_train |>
  dplyr::select(ref_date, name = ticker, value = ret_closing_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)
print(round(mean_ret, 6))

##      AFFLE.BO        BCG.BO HAPPSTMNDS.NS          INFY          SIFY 
##      0.000806      0.001972      0.001026      0.000322      0.000438 
## ZENSARTECH.BO 
##      0.001322

# Matriz de Covariancia sem anualização

cov_mat <-cov(da_xts, use = "complete.obs")
print(round(cov_mat,6))

##               AFFLE.BO   BCG.BO HAPPSTMNDS.NS     INFY      SIFY ZENSARTECH.BO
## AFFLE.BO      0.000572 0.000110      0.000140 0.000063  0.000060      0.000185
## BCG.BO        0.000110 0.001894      0.000196 0.000004  0.000009      0.000118
## HAPPSTMNDS.NS 0.000140 0.000196      0.000568 0.000078 -0.000027      0.000202
## INFY          0.000063 0.000004      0.000078 0.000308  0.000196      0.000074
## SIFY          0.000060 0.000009     -0.000027 0.000196  0.002777      0.000011
## ZENSARTECH.BO 0.000185 0.000118      0.000202 0.000074  0.000011      0.000768

# Matriz de correlação 
cor_mat <-cor(da_xts, use = "complete.obs")
print(round(cor_mat,6))

##               AFFLE.BO   BCG.BO HAPPSTMNDS.NS     INFY      SIFY ZENSARTECH.BO
## AFFLE.BO      1.000000 0.106136      0.245489 0.150661  0.047736      0.278590
## BCG.BO        0.106136 1.000000      0.188464 0.004894  0.003990      0.097995
## HAPPSTMNDS.NS 0.245489 0.188464      1.000000 0.185462 -0.021461      0.305925
## INFY          0.150661 0.004894      0.185462 1.000000  0.212511      0.151886
## SIFY          0.047736 0.003990     -0.021461 0.212511  1.000000      0.007628
## ZENSARTECH.BO 0.278590 0.097995      0.305925 0.151886  0.007628      1.000000

Como o oitavo exercício pede para refazer as questões de 5 a 7, vamos começar por ele

8. Fazer uma otimização dos pesos do seu portifólio inicial e repetir as etapas 5 e 7, comparando o seu portifólio inicial com otimizado.

Para calcular risco e retorno do portfolio precisamos calcular:

• Mean assets returns
• Portfolio weights
• Covariance matrix of all assets
• Random weights

Vamos começar calculando os random weights

set.seed(2)

wts <- runif(n = length(ativos))

(wts <- wts/sum(wts))

## [1] 0.05258389 0.19976799 0.16306447 0.04779703 0.26844513 0.26834149

Portfolio returns:

(port_returns <- sum(wts * mean_ret))

## [1] 0.001091481

Portfolio risk:

(port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts)))

##            [,1]
## [1,] 0.02029622

Sharpe Ratio

(sharpe_ratio <- port_returns/port_risk)

##            [,1]
## [1,] 0.05377757

Vamos agora fazer a otimização:

# criação de matriz vazia para armazenar valores 
sim_returns <- function(i) {
  wts <- runif(length(ativos))
  wts <- wts / sum(wts)
  port_ret <- sum(wts * mean_ret)
  port_sd <- as.numeric(sqrt(t(wts) %*% (cov_mat %*% wts)))
  sr <- port_ret / port_sd
  
  wts |>
    purrr::set_names(ativos) |>
    tibble::enframe() |>
    tidyr::pivot_wider() |>
    dplyr::mutate(
      return = port_ret,
      risk = port_sd,
      sharpe = sr
    )
}
# Loop para 5 mil vezes 
portfolio_values <- purrr::map(1:5000, sim_returns, .progress = TRUE) |>
  bind_rows(.id = "run")

min_var <- portfolio_values[which.min(portfolio_values$risk),]
max_sr <- portfolio_values[which.max(portfolio_values$sharpe),]

Plotando os pesos começando pelo portfolio de minima variancia:

min_var |> 
  pivot_longer(2:7) |> 
  mutate(name = forcats::fct_reorder(name, value)) |> 
  ggplot(aes(name, value, fill = name)) +  
  geom_col() +
  scale_fill_brewer(palette = "Set3") +  # Using a predefined color palette
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Asset",
    y = "Weight",
    title = "Minimum variance portfolio weights"
  ) +
  theme_minimal()

# Tangente do Portfolio 
max_sr %>%
  pivot_longer(2:7) %>%
  mutate(name = forcats::fct_reorder(name, value)) %>%
  ggplot(aes(name, value, fill = name)) +  
  geom_col() +
  scale_fill_brewer(palette = "Set2") +  
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Asset",
    y = "Weight",
    title = "Tangency portfolio weights"
  ) +
  theme_minimal() 

Portfolio Otimizado, Fronteira Eficiente:

portfolio_values |>
  ggplot(aes(x = risk, y = return, color = sharpe)) +
  geom_point() +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(
    x = 'Risk',
    y = 'Returns',
    title = "Portfolio Optimization & Efficient Frontier"
  ) +
  geom_point(
    aes(x = risk, y = return),
    data = min_var,
    color = 'red',
    size = 3
  ) +
  geom_point(
    aes(x = risk, y = return),
    data = max_sr,
    color = 'orange',
    size = 3
  )

Agora, vamos calcular o beta do portfolio assumindo as premissas do modelo CAPM:

pesos_finais <- max_sr %>%
  dplyr::select(2:7) %>%  # Selecionando pesos para os 6 ativos
  as.numeric()

# print(length(unique_assets))
print(length(pesos_finais))

## [1] 6

portfolio_returns <- da_train %>%
  tidyquant::tq_portfolio(
    assets_col = ticker,
    returns_col = ret_closing_prices,
    weights = pesos_finais,
    col_rename = "portfolio"
  )

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_closing_prices)

### Retorno
all_returns <- market_returns |> 
    dplyr::inner_join(portfolio_returns, by = "ref_date", copy = TRUE) |> 
    tidyr::drop_na()

Estimar o Beta do portfolio levando em consideração o Indice Bovespa

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

## [1] 0.09999931

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

betas <- purrr::map_dbl(unique(da_train$ticker), calcular_beta) |> 
  purrr::set_names(unique(da_train$ticker))

sum(betas * pesos_finais) # Beta médio ponderado do portfolio, leva em consideração o peso de cada ativo

## [1] 0.1243562

beta_geral # O beta geral sugere que o portfolio tem uma volatilidade muito inferior ao mercado.

## [1] 0.09999931

Beta para cada Ativo

CAPM do Portfolio e do Ativo

capm_lm_tudo <- lm(portfolio ~ ibov, data = all_returns) |>
  broom::tidy() |>
  dplyr::filter(term == "ibov") |>
  with(estimate)

capm_lm_individual <- purrr::map_dbl(ativos, \(ativo) {
  da_model <- da_train |>
    dplyr::filter(ticker == ativo) |>
    dplyr::inner_join(market_returns, "ref_date")
  lm(ret_closing_prices ~ ibov, data = da_model) |>
    broom::tidy() |>
    dplyr::filter(term == "ibov") |>
    dplyr::pull(estimate)
}) |>
  purrr::set_names(ativos)

capm_lm_individual

##          SIFY          INFY HAPPSTMNDS.NS      AFFLE.BO ZENSARTECH.BO 
##    0.72735806    0.29774791    0.02926001    0.08351927    0.11264343 
##        BCG.BO 
##    0.08435989

capm_lm_tudo

## [1] 0.09999931

betas_capm <- purrr::set_names(capm_lm_individual, unique(all_returns$ticker))
print(betas_capm)

## [1] 0.72735806 0.29774791 0.02926001 0.08351927 0.11264343 0.08435989

# Plot do Beta do Portfolio 
data_beta <- tibble::tibble(
  Ativo = c(names(capm_lm_individual), "Portfólio Total"),
  Beta = c(capm_lm_individual, capm_lm_tudo)
)

ggplot(data_beta, aes(x = Ativo, y = Beta, fill = Ativo)) +
  geom_col() +
  geom_hline(yintercept = capm_lm_tudo, linetype = "dashed", color = "red") +
  labs(title = "Comparação dos Betas CAPM", x = "Ativo", y = "Beta") +
  theme_minimal() +
  theme(legend.position = "none")