lab03.qmd

---
title: "Meus ativos"
format:
  html:
    code-fold: true
    monofont: "JetBrains Mono"
execute:
  freeze: true
---

Quiz para aquecer: <https://forms.gle/CcMvpyXdhaacQtPi8>

Nosso lab terá 2 partes:

- Primeiro, vamos fazer a descritiva e rodar os modelos GARCH
- Depois, vamos trabalhar com otimização do portfolio, CAPM e VaR

Vamos dividir em grupos. Quem estiver na frente, pode avançar.

Teremos 2 momentos para discutir os avanços.

```{r}
#| code-fold: show
#| message: false
library(fpp3)
library(rugarch)
```


```{r}
#| code-fold: show
#| message: false
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.

```{r}
#| code-fold: show
#| message: false
# 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
)

View(da)

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

Plotar

```{r}
#| fig.height: 10
#| fig.width: 8
da_tsibble |>
  autoplot(price_close, colour = "black") +
  facet_wrap(~ticker, scales = "free_y", ncol = 1)
```

```{r}
#| fig.height: 10
#| fig.width: 8
da_tsibble |>
  autoplot(ret_closing_prices, colour = "black") +
  facet_wrap(~ticker, scales = "free_y", ncol = 1)
```

Data mínima comum a todas as séries

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

data_corte
```

```{r}
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

```{r}
da_tsibble |>
  ACF(ret_closing_prices) |>
  autoplot()
```


```{r}
da_tsibble |>
  PACF(ret_closing_prices) |>
  autoplot()
```


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


```{r}
da_tsibble |>
  dplyr::mutate(ret2 = ret_closing_prices^2) |>
  ACF(ret2) |>
  autoplot()
```

```{r}
da_tsibble |>
  dplyr::mutate(ret2 = ret_closing_prices^2) |>
  PACF(ret2) |>
  autoplot()
```

Normalidade

```{r}
# histogram with geom_histogram of each ticker
da_train |>
  ggplot(aes(x = ret_closing_prices)) +
  geom_histogram(bins = 90) +
  facet_wrap(~ticker, ncol = 3)
```

```{r}
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()
```

Com outra distribuição

```{r}
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()
```

## Ajustando modelos garch

Função para ajustar um garch

```{r}
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
}
```

Testando para um ativo

```{r}
garch_individual(
  parms = list(
    p = 1, q = 1, m = 1, n = 1, dist = "std"
  ),
  ret = da_train |>
    dplyr::filter(ticker == "HGRE11.SA") |>
    pull(ret_closing_prices)
)
```

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

```{r}

### OMITIDO

```

Rodando as funções

```{r}
#| eval: false
melhores_por_ativo <- ativos |>
  purrr::set_names() |>
  purrr::map(melhor_garch, .progress = TRUE) |>
  dplyr::bind_rows(.id = "ticker")
```


```{r}
#| echo: false
melhores_por_ativo <- readr::read_rds("melhores_por_ativo.rds")
```


## Prever volatilidade um passo à frente

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

```{r}
prever_volatilidade <- function(parms, n_steps = 5) {
  usethis::ui_info("Prevendo volatilidade para {parms$ticker}...")

  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

```{r}
parametros_melhores <- melhores_por_ativo |>
  group_by(ticker) |>
  slice_head(n = 1) |>
  ungroup()

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

vol_futuro
```

## Comparar volatilidades entre os retornos selecionados

...

## Montagem 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/>

```{r}
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)
```


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

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

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

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

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

(wts <- wts/sum(wts))

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

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

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

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.~

```{r}

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, .progress = TRUE) |>
  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.

```{r}

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

```

```{r}
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"
  )
```

```{r}
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
  )
```

## VaR do portfolio

```{r}
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)))
```

## CAPM

```{r}
portfolio_returns <- da_train |>
  tidyquant::tq_portfolio(
    ticker,
    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)

all_returns <- market_returns |>
  dplyr::inner_join(portfolio_returns, "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_closing_prices, ibov) / var(ibov))
}

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

sum(betas * pesos_finais)
beta_geral

```

```{r}

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

capm_lm_individual
```