Meus ativos

Quiz: https://forms.gle/CcMvpyXdhaacQtPi8

Code

library(fpp3)
library(rugarch)

Code

start_date <- '2018-01-01'
# esses são ativos de fundos imobiliários que eu ja tive
# e queria saber fiz um péssimo investimento
# ou apenas ruim.
ativos <- c(
  "HGRE11.SA", 
  "BTLG11.SA", 
  "HGRU11.SA", 
  "VGIR11.SA", 
  "MGFF11.SA"
) 

Vamos trabalhar tanto com os dados no formado de tibble quanto no formato de tsibble.

Code

# library(curl)
# has_internet_via_proxy <<- TRUE
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)

Plotar

Code

da_tsibble |> 
  autoplot(price_close, colour = "black") +
  facet_wrap(~ticker, scales = "free_y", ncol = 1)

Code

da_tsibble |> 
  autoplot(ret_closing_prices, colour = "black") +
  facet_wrap(~ticker, scales = "free_y", ncol = 1)

Warning: Removed 1 row(s) containing missing values (geom_path).

Data mínima comum a todas as séries

Code

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

data_corte

[1] "2018-07-31"

Code

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

Descritivas bacanas

• ACF/PACF dos retornos
• visualizar os retornos ao quadrado
• ACF/PACF dos retornos ao quadrado

Code

da_tsibble |> 
  ACF(ret_closing_prices) |> 
  autoplot()

Warning: Provided data has an irregular interval, results should be treated with
caution. Computing ACF by observation.

Code

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

Warning: Provided data has an irregular interval, results should be treated with
caution. Computing ACF by observation.

Code

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

Warning: Removed 1 row(s) containing missing values (geom_path).

Code

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

Warning: Provided data has an irregular interval, results should be treated with
caution. Computing ACF by observation.

Code

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

Warning: Provided data has an irregular interval, results should be treated with
caution. Computing ACF by observation.

Normalidade

Code

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

Com outra distribuição

Code

da_train |> 
  group_by(ticker) |> 
  summarise(
    gg = list(
      ggplot(cur_data(), 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()

Ajustando modelos garch

Função para ajustar um garch

Code

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
}

Função para ajustar uma grid de garchs e pegar as informações

Code

melhor_garch <- function(ativo, p = 0:3, q = 0:3, 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 = dist) |> 
    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}...")
  
  progressr::with_progress({
    prog <- progressr::progressor(nrow(grid))
    modelos <- grid |> 
      group_split(id) |> 
      purrr::map(garch_individual, ret = ret, prog)
  })
  
  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
  
}

Rodando as funções

Code

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

Prever volatilidade um passo à frente

Função que ajusta o modelo e faz as previsões

Code

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)
  
  # browser()
  
  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

Code

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

vol_futuro <- parametros_melhores |> 
  group_split(ticker) |> 
  purrr::map(prever_volatilidade, n_steps = 1) |> 
  dplyr::bind_rows()

vol_futuro

# A tibble: 5 × 4
  ticker          serie volatilidade shape
  <chr>           <dbl>        <dbl> <dbl>
1 BTLG11.SA  0.000292        0.00905  2.50
2 HGRE11.SA  0.000190        0.0144   3.11
3 HGRU11.SA -0.000130        0.0128   2.32
4 MGFF11.SA -0.00128         0.0130   3.65
5 VGIR11.SA -0.00000968      0.00649  2.90

Comparar volatilidades entre os retornos selecionados

…

Fazer um portfolio

Vamos prever a volatilidade um passo a frente para um portfolio gerado aleatoriamente

Code

set.seed(42)
pesos <- runif(5)
pesos <- pesos / sum(pesos)

portfolio_returns <- da_train |> 
  tidyquant::tq_portfolio(
    ticker, 
    ret_closing_prices, 
    weights = pesos,
    col_rename = "portfolio"
  )

Warning: `spread_()` was deprecated in tidyr 1.2.0.
ℹ Please use `spread()` instead.
ℹ The deprecated feature was likely used in the tidyquant package.
  Please report the issue at
  <https://github.com/business-science/tidyquant/issues>.

Warning in PerformanceAnalytics::Return.portfolio(., weights = weights, : NA's
detected: filling NA's with zeros

Code

# faz uma grid de hiperparâmetros 
grid <- expand_grid(p = 0:3, q = 0:3, m = 0:1, n = 0:1, dist = c("std", "norm")) |> 
  tibble::rownames_to_column("id")
  
# pega o retorno do ativo
ret <- portfolio_returns$portfolio

usethis::ui_info("Ajustando modelos para portfolio...")

ℹ Ajustando modelos para portfolio...

Code

future::plan(future::multisession, workers = 8)
progressr::with_progress({
  prog <- progressr::progressor(nrow(grid))
  modelos <- grid |> 
    group_split(id) |> 
    furrr::future_map(garch_individual, ret = ret, prog)
})

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]}"
))

ℹ Melhor modelo:
  p <- 2
  q <- 0
  m <- 1
  n <- 1

Code

garch_model = ugarchspec(
  variance.model = list(
    model = "fGARCH",
    submodel = "GARCH",
    garchOrder = c(melhores$m[1], melhores$n[1])
  ),
  mean.model = list(
    armaOrder = c(melhores$p[1], melhores$q[1]),
    include.mean = TRUE
  ),
  distribution.model = melhores$dist[1]
)

fit <- ugarchfit(garch_model, data = ret)

# browser()

if (melhores$dist[1] == "std") {
  shape <- as.numeric(fit@fit$coef["shape"])
} else {
  shape <- NA_real_
}

forecasts <- ugarchforecast(fit, n.ahead = 1)@forecast

tibble::tibble(
  serie = as.numeric(forecasts$seriesFor),
  volatilidade = as.numeric(forecasts$sigmaFor),
  shape = shape
)

# A tibble: 1 × 3
      serie volatilidade shape
      <dbl>        <dbl> <dbl>
1 -0.000867      0.00649  4.12

Otimização de portfolio

Reproduzindo código daqui:

https://www.codingfinance.com/post/2018-05-31-portfolio-opt-in-r/

Versão em python

https://www.codingfinance.com/post/2018-05-31-portfolio-opt-in-python/

Code

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)

Code

mean_ret <- colMeans(da_xts, na.rm = TRUE)
print(round(mean_ret, 5))

BTLG11.SA HGRE11.SA HGRU11.SA MGFF11.SA VGIR11.SA 
  0.00035   0.00007   0.00013  -0.00039   0.00000 

Next we will calculate the covariance matrix for all these stocks. We will NOT annualize it by multiplying by 252.

Code

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

          BTLG11.SA HGRE11.SA HGRU11.SA MGFF11.SA VGIR11.SA
BTLG11.SA  0.000144  0.000065  0.000055  0.000058  0.000040
HGRE11.SA  0.000065  0.000207  0.000079  0.000084  0.000068
HGRU11.SA  0.000055  0.000079  0.000148  0.000063  0.000056
MGFF11.SA  0.000058  0.000084  0.000063  0.000192  0.000060
VGIR11.SA  0.000040  0.000068  0.000056  0.000060  0.000143

Before we apply our methods to thousands of random portfolio, let us demonstrate the steps on a single portfolio.

To calculate the portfolio returns and risk (standard deviation) we will us need

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

Code

set.seed(2)
# Calculate the random weights
wts <- runif(n = length(ativos))

(wts <- wts/sum(wts))

[1] 0.07186944 0.27303447 0.22286964 0.06532697 0.36689948

Code

# Calculate the portfolio returns
(port_returns <- sum(wts * mean_ret))

[1] 4.950717e-05

Code

# Calculate the portfolio risk
(port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts)))

            [,1]
[1,] 0.009518973

Code

# Calculate the Sharpe Ratio
(sharpe_ratio <- port_returns/port_risk)

            [,1]
[1,] 0.005200893

We have everything we need to perform our optimization. All we need now is to run this code on 5000 random portfolios. For that we will use a for loop.

~Before we do that, we need to create empty vectors and matrix for storing our values.~

Code

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

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

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

Lets plot the weights of each portfolio. First with the minimum variance portfolio.

Code

min_var |> 
  pivot_longer(2:6) |> 
  mutate(name = forcats::fct_reorder(name, value)) |> 
  ggplot(aes(name, value)) +
  geom_col() +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Asset",
    y = "Weight",
    title = "Minimum variance portfolio weights"
  )

Code

max_sr |> 
  pivot_longer(2:6) |> 
  mutate(name = forcats::fct_reorder(name, value)) |> 
  ggplot(aes(name, value)) +
  geom_col() +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Asset",
    y = "Weight",
    title = "Tangency portfolio weights"
  )

Code

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'
  ) +
  geom_point(
    aes(x = risk, y = return), 
    data = max_sr, 
    color = 'red'
  )

VaR do portfolio

Code

pesos_finais <- min_var |> 
  dplyr::select(2:6) |> 
  as.numeric()


rt_final <- mean(vol_futuro$serie * pesos_finais)
st_dev_final <- sqrt(pesos_finais %*% cov_mat %*% pesos_finais)
nu <- min(vol_futuro$shape)
valor_t <- qt(.95, nu)

(VaR <- rt_final + valor_t * st_dev_final / sqrt(nu/(nu-2)))

            [,1]
[1,] 0.008717397

CAPM

Code

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

Warning in PerformanceAnalytics::Return.portfolio(., weights = weights, : NA's
detected: filling NA's with zeros

Code

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)

── Running yfR for 1 stocks | 2018-08-01 --> 2022-11-12 (1564 days) ──

ℹ Downloading data for benchmark ticker ^GSPC

ℹ (1/1) Fetching data for ^BVSP

!   - not cached

✔   - cache saved successfully

✔   - got 1061 valid rows (2018-08-01 --> 2022-11-11)

✔   - got 96% of valid prices -- You got it !

ℹ Binding price data

ℹ Completing data points (do_complete_data = TRUE)

── Diagnostics ─────────────────────────────────────────────────────────────────

✔ Returned dataframe with 1061 rows -- Well done !

ℹ Using 864.4 kB at C:\Users\julio\AppData\Local\Temp\Rtmp4yKLF3/yf_cache for 7 cache files

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

Code

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

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

[1] 0.2333843

Code

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(ativos, calcular_beta) |> 
  purrr::set_names(ativos)

sum(betas * pesos_finais)

[1] 0.2564028

Code

beta_geral

[1] 0.2333843

Code

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_tudo

[1] 0.2333843

Code

capm_lm_individual

HGRE11.SA BTLG11.SA HGRU11.SA VGIR11.SA MGFF11.SA 
0.2635135 0.2264689 0.2475435 0.1971465 0.2871617