Pacote PortfolioAnalytics

Para resolver estes problemas menos convencionais, podemos nos valer do pacote PortfolioAnalytics. Este pacote fornece soluções numéricas para problemas de portfólio com restrições complexas. Este tipo de generalidade das possibilidades de restrições e funções objetivas, distinguem o pacote de competidores como o pacote fPortfolio, que são mais indicados para problemas padrões de Markowitz.

Assim, o pacote PortfolioAnalytics pode resolver problemas genéricos do tipo

\[\min_w g(w) \] \[\text{sujeito a } h_1(w) \leq 0\] \[...\] \[h_q(w) \leq 0\]

Após ter o problema montado, o valor ótimo desta função é encontrado a partir de um algorítimo de Differential Evolution (DE).

Estatísticas descritivas

É sempre boa prática verificar a média e desvio padrão dos retornos, que representam o retorno médio por ativo e a volatidade do ativo. A tabela @reff(tab:media-desvio) mostram estas estatísticas descritivas.

tbl_retornos %>% 
  group_by(ativo) %>% 
  summarise(`Média` = mean(retorno),
            `Desvio-padrão` = sd(retorno)) %>% 
  knitr::kable(caption = "Média e Desvio-padrão dos ativos")

Já a figura 1 exibe o retorno acumulado para as empresas selecionadas.

tbl_retornos %>% 
  group_by(ativo) %>% 
  mutate(retorno_acumulado = cumsum(retorno)) %>% 
  ggplot(aes(x = data, y = retorno_acumulado, color = ativo)) +
  geom_line() + 
  geom_vline(xintercept = as.Date("2020-03-01"), linetype = 2, color = "red") +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "", y = "Retorno acumulado", 
       title = "Retorno acumulado por empresa selecionada, 2017-2020",
       subtitle = "Linha vermelha representa o início da crise de Covid-19") +
  facet_wrap(~ativo, scales = "free_y", ncol = 2) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle=45, size = 8))

Figure 1: Retorno acumulado por empresa selecionada, 2017-2021

Contudo, o pacote PortfolioAnalytics não gosta deste formato dos dados. O formato original, que podemos chamar de tidy, apresenta uma única coluna que indica o nome das empresas e uma única coluna que informa os retornos. Precisamos converter para o formato wider, em que cada empresa tem sua própria coluna de retornos.

A tabela 4 mostra o novo formato dos dados, onde cada empresa ocupa uma coluna. Também vamos criar um período de treinamento, em que o portfolio é otimizado (entre 2012 e 2020), para ser testado e comparado com um benchmarking durante 2021.

tbl_retornos_wide <- tbl_retornos %>% 
  pivot_wider(names_from = ativo, 
              values_from = retorno
              ) %>% 
  column_to_rownames("data") %>% 
  as.xts()

tbl_retornos_wide_training <- tbl_retornos %>% 
  filter(data < as.Date("2021-01-01")) %>% 
  pivot_wider(names_from = ativo, 
              values_from = retorno
              ) %>% 
  column_to_rownames("data") %>% 
  as.xts()

tbl_retornos_wide_training %>% 
  head() %>% 
  knitr::kable(caption = "Tabela de retornos no formato wide")

Este procedimento de treinamento e teste será feito manualmente na primeira parte do post, mas com a função optimize.portfolio.rebalancing podemos realizar um processo de maneira mais automática.

Portfólio com Variância Mínima Global

Começamos com um problema de otimização que busca obter os pesos que nos forneça a menor variâcia possível (a volatilidade do portfólio). Volatilidade é uma medida de risco, e portanto, quanto menor a volatilidade, menor o risco em um portfólio.

Dito isto, nosso objetivo se torna encontrar os pesos que minimizem nossa variância sujeito a duas restrições:

1. que a soma dos pesos dos ativos seja a unidade; e

2. que cada ativo represente no mínimo 5% do total investido e no máximo 60%.

A primeira restrição garante o investimento completo (full investiment), ou seja, que utilizaremos todos os recursos disponíveis. A segunda restrição garante que todos os ativos farão parte do portfólio, ao mesmo tempo que não permite que um ativo domine todo o investimento.

Este problema pode ser anunciado matematicamente da seguinte forma:

\[\min \sigma^2 = W^T \sum W\] \[\text{s.a.}\] \[\sum_{i=1}^N W_i = 1\] \[0.05 \leq W_i \leq 0.6\]

Técnicamente, como a função objetiva é quadratica e as restrições são lineares, utilizaremos um solucionador (solver) quadrático. O pacote PortfolioAnalytics utiliza o plugin ROI.plugin.quadprog para resolver este tipo de problema.

#### Portfólio com a variância mínima global
min_var_portfolio <- portfolio.spec(assets = tickers) %>% 
  add.constraint(type = "full_investiment") %>% 
  add.constraint(type = "box", min = 0.05, max = 0.6) %>% 
  add.objective(type = "risk", name = "var")

otimizacao_risco <- tbl_retornos_wide_training %>% 
  optimize.portfolio(
    portfolio = min_var_portfolio,
    optimize_method = "quadprog",
    trace = TRUE
    )
otimizacao_risco

## ***********************************
## PortfolioAnalytics Optimization
## ***********************************
## 
## Call:
## optimize.portfolio(R = ., portfolio = min_var_portfolio, optimize_method = "quadprog", 
##     trace = TRUE)
## 
## Optimal Weights:
## MGLU3.SA ABEV3.SA BBDC3.SA EMBR3.SA JBSS3.SA PETR3.SA VALE3.SA WEGE3.SA 
##   0.0500   0.3042   0.0500   0.1409   0.0789   0.0500   0.0765   0.2494 
## 
## Objective Measure:
##  StdDev 
## 0.05104

pesos_var_minima <- extractWeights(otimizacao_risco)

tibble(ativos = tickers, pesos = pesos_var_minima) %>% 
  mutate(soma_pesos = cumsum(pesos)) %>% 
  knitr::kable(caption = "Conjunto de pesos dos ativos que minimizam o risco do porfólio")

Portfólio de Retorno Máximo Esperado

Agora o objetivo é inverso. Queremos encontrar os pesos que produzem o portfólio com maior retorno esperado. O problema se torna maximizar o retorno \(\mu\) sujeito às mesmas duas restrições anteriores:

\[\max \mu^T W\] \[\text{s.a. } \sum_{i=1}^N W_i = 1\]

\[\text{ e } \epsilon_i \leq W_i \leq \delta_i\]

Como a função é linear, podemos utilizar um solver de programação linear. O PortfolioAnalytics fornece a plugin glpk como uma solução para problemas de otimização de programação linear.

portfolio_retorno_maximo <- portfolio.spec(assets = tickers) %>% 
  add.constraint(type = "full_investiment") %>% 
  add.constraint(type = "box", min = 0.05, max = 0.6) %>% 
  add.objective(type = "return", name = "mean")

otimizacao_retorno <- tbl_retornos_wide_training %>% 
  optimize.portfolio(
    portfolio = portfolio_retorno_maximo,
    optimize_method = "glpk",
    trace = TRUE
    )
otimizacao_retorno

## ***********************************
## PortfolioAnalytics Optimization
## ***********************************
## 
## Call:
## optimize.portfolio(R = ., portfolio = portfolio_retorno_maximo, 
##     optimize_method = "glpk", trace = TRUE)
## 
## Optimal Weights:
## MGLU3.SA ABEV3.SA BBDC3.SA EMBR3.SA JBSS3.SA PETR3.SA VALE3.SA WEGE3.SA 
##     0.60     0.05     0.05     0.05     0.05     0.05     0.05     0.10 
## 
## Objective Measure:
##    mean 
## 0.04271

pesos_retorno_maximo <- extractWeights(otimizacao_retorno)


tibble(ativos = tickers, pesos = pesos_retorno_maximo) %>% 
  mutate(soma_pesos = cumsum(pesos)) %>% 
  knitr::kable(caption = "Conjunto de pesos dos ativos que maximizam o retorno médio do porfólio")

Portfólio com Razão de Sharpe Máxima

Por fim, podemos encontrar um portfólio que maximiza a razão de Sharpe. Mas antes, uma breve explicação sobre esta medida de risco.

portfolio_sharpe <- portfolio.spec(assets = tickers) %>% 
  add.constraint(type = "full_investiment") %>% 
  add.constraint(type = "box", min = 0.05, max = 0.6) %>% 
  add.objective(type = "return", name = "mean") %>% 
  add.objective(type = "risk", name = "StdDev")

otimizacao_sharpe <- tbl_retornos_wide_training %>% 
  optimize.portfolio(
    portfolio=portfolio_sharpe,
    optimize_method="ROI",
    maxSR=TRUE, 
    trace=TRUE)

otimizacao_sharpe

## ***********************************
## PortfolioAnalytics Optimization
## ***********************************
## 
## Call:
## optimize.portfolio(R = ., portfolio = portfolio_sharpe, optimize_method = "ROI", 
##     trace = TRUE, maxSR = TRUE)
## 
## Optimal Weights:
## MGLU3.SA ABEV3.SA BBDC3.SA EMBR3.SA JBSS3.SA PETR3.SA VALE3.SA WEGE3.SA 
##   0.1369   0.0500   0.0500   0.0500   0.1050   0.0500   0.0500   0.5082 
## 
## Objective Measure:
## StdDev 
##  0.063 
## 
## 
##    mean 
## 0.02789

pesos_sharpe <- extractWeights(otimizacao_sharpe)

tibble(ativos = tickers, pesos = pesos_sharpe) %>% 
  mutate(soma_pesos = cumsum(pesos)) %>% 
  knitr::kable(caption = "Conjunto de pesos dos ativos que maximizam a Razão de Sharpe")

Construíndo o Portfólio

A figura 2 mostra os pesos de cada ativo em cada problema resolvido.

tibble(ativos = tickers, 
       pesos_sharpe = pesos_sharpe,
       pesos_retorno = pesos_retorno_maximo,
       pesos_risco = pesos_var_minima) %>% 
  pivot_longer(cols = "pesos_sharpe":"pesos_risco") %>% 
  ggplot(aes(x = name, y = value, fill = ativos)) +
  geom_col() + 
  scale_x_discrete(labels = c("Maximizar retorno", "Minimizar risco", "Maximizar Risco de Sharpe")) +
  labs(x = "Estratégia", y = "Peso no portfólio",
       title = "Seleção de carteiras de investimentos por estratégia",
       subtitle = "Empresas selecionadas, 2017-2020") +
  geom_text(aes(label = ativos),
            position = position_stack(vjust = .6),
            size = 2) + 
  theme(legend.position = "none")

Figure 2: Seleção de carteiras de investimentos por estratégia

Agora podemos construir uma série no tempo dos retornos do portfólio caso tivessemos investido exatamente como sugerido pelos três problemas definidos. Por fim, vamos compara-los a uma solução de benchmarking, em que investimos 12,5% em cada ativo.

Para construir esta nova série, vamos multiplicar cada retorno observado em 2021 pelo peso definido pelo otimizador. Assim, o retorno do portfólio em um mês \(t\) é dado por:

\[r_{\text{portfolio}, t} = \sum_{i=1}^5 W_i r_{i,t}\]

É fácil computar manualmente, ou utilizando a função Return.portfolio do mesmo pacote PortfolioAnalytics, mas como uma alternativa conveniente vamos utilizar o pacote tidyquant, que possui a função tq_portfolio(). Esta função utiliza os dados originais (no formato tidy) e o vetor de pesos produzidos pela otimização para construir um data.frame de retornos de portfólio. Abaixo fazemos o computo do retorno do portfólio para as três estratégias.

retorno_portfolio_retorno_maximo <- tbl_retornos %>% 
  filter(data > as.Date("2021-01-01")) %>% 
  tq_portfolio(
    assets_col = ativo,
    returns_col = retorno,
    weights = pesos_retorno_maximo,
    col_rename = "retorno_max_retorno",
    rebalance_on = "months"
  )

retorno_portfolio_risco_minimo <- tbl_retornos %>% 
  filter(data > as.Date("2021-01-01")) %>% 
  tq_portfolio(
    assets_col = ativo,
    returns_col = retorno,
    weights = pesos_var_minima,
    col_rename = "retorno_minimo_risco",
    rebalance_on = "months"
  )

retorno_portfolio_sharpe <- tbl_retornos %>% 
  filter(data > as.Date("2021-01-01")) %>% 
  tq_portfolio(
    assets_col = ativo,
    returns_col = retorno,
    weights = pesos_sharpe,
    col_rename = "retorno_max_sharpe",
    rebalance_on = "months"
  )

## Warning in check_weights(weights, assets_col, map, x): Sum of weights must be 1.

Vamos definir ainda um portfólio benchmarking de pesos iguais para os ativos. A intuição é que uma abordagem de pesos iguais não apresenta preferências para um ativo específico. A pergunta que faremos é: “os problemas de otimização de carteiras podem produzir resultados melhores do que simplesmente construir um portfólio de pesos iguais?”.

Como temos 8 empresas no nosso portfólio, vamos criar um portfólio naïve, onde o investidor aloca a mesma quantidade de recursos para as 8 empresas, ou seja, 12,5%.

peso_naive <- rep(1/8, 8)

retorno_portfolio_naive <- tbl_retornos %>% 
  filter(data > as.Date("2021-01-01")) %>% 
  tq_portfolio(
    assets_col = ativo,
    returns_col = retorno,
    weights = peso_naive,
    col_rename = "retorno_naive",
    rebalance_on = "months"
  )

A tabela 8 mostra os dados para as quatro estratégias empregadas.

retorno_portfolio <- retorno_portfolio_naive %>% 
  left_join(retorno_portfolio_retorno_maximo, by = "data") %>% 
  left_join(retorno_portfolio_risco_minimo, by = "data") %>% 
  left_join(retorno_portfolio_sharpe, by = "data")

retorno_portfolio %>% tail() %>% 
  knitr::kable(caption = "Retorno para os quatro portfólios")

Agora podemos visualizar o retorno acumulado no período caso o investidor tivesse, em janeiro de 2021, investido seus recursos empregando cada uma das estratégias.

A figura 3 mostra que a estratégia que produziu os maiores ganhos foi a de minimização do risco.

retorno_portfolio_acumulado <- retorno_portfolio %>%
  pivot_longer(cols = "retorno_naive":"retorno_max_sharpe",
               names_to = "estrategia",
               values_to = "retorno") %>% 
  mutate(estrategia = case_when(estrategia == "retorno_max_retorno" ~ "Maximizar Retorno Médio",
                                estrategia == "retorno_minimo_risco" ~ "Minimizar o risco",
                                estrategia == "retorno_naive" ~ "Alocação igualitária",
                                estrategia == "retorno_max_sharpe" ~ "Maximizar Razão de Sharpe")) %>% 
  group_by(estrategia) %>% 
  mutate(retorno_acumulado = cumsum(retorno)) %>% 
  ungroup()

retorno_portfolio_acumulado %>% 
  ggplot(aes(x = data, y = retorno_acumulado, color = estrategia)) + 
  geom_line(size = 1.2) + 
  scale_y_continuous(label = scales::percent) +
  labs(title = "Retorno acumulado por tipo de Portfólio, 2021", y = "Retorno acumulado do portfólio",
       x = "") + 
  facet_wrap(~estrategia) + 
  theme(legend.position = "none")

Figure 3: Retorno acumulado por tipo de Portfólio, 2017-2021

Por fim, algumas contas de padaria. A tabela 9 mostra que se um investidor em janeiro de 2020 tivesse alocado \(\text{R\$}\) 1000 dividos igualmente entre as 8 empresas (\(\text{R\$}\) 125 para cada), teria um retorno acumulado em dezembro de 2021 de \(\text{R\$}\) 1069. Este é um retorno semelhante ao de de um portfólio que tentasse minimizar os riscos (\(\text{R\$}\) 1148). No caso da estratégia de maximização de retornos, o investimento de \(\text{R\$}\) 1000 teria produzido \(\text{R\$}\) 339 reais. Maximizar a razão de Sharpe renderia \(\text{R\$}\) 941.

retorno_portfolio_acumulado %>% 
  select(-retorno) %>% 
  filter(data == "2021-12-31") %>% 
  mutate(retorno_1000_reais = 1000 + (1000 * retorno_acumulado)) %>% 
  knitr::kable(caption = "Retorno de R$ 1000 para cada tipo de portfólio")

Minimização de Risco com Rebalanceamento

Vamos repetir a otimização de minimização de risco, mas agora permitindo que o portfólio tenha seus pesos reequilibrados trimestralmente. Além dos parâmetros comuns à função optimize.portfolio(), a função optimize.portfolio.rebalance exige que informemos a periodicidade com que os pesos devem ser rebalanceados (rebalance_on = "quarters").

Devemos ainda definir um período de treinamento, com training_period = 60, e um valor para rolling_window, para que o rebalanceamento dos pesos utilize uma janela móvel de dados. Caso o usuário queira utilizar todos os dados para estimar os pesos trimestrais, pode fixar rolling_window = NULL.

otimizacao_risco_rebalance <- tbl_retornos_wide %>% 
  optimize.portfolio.rebalancing(
    optimize_method = "ROI",
    portfolio = min_var_portfolio,
    rebalance_on = "quarters",
    training_period = 60,
    rolling_window = 60
  )

## Warning: executing %dopar% sequentially: no parallel backend registered

otimizacao_risco_rebalance

## **************************************************
## PortfolioAnalytics Optimization with Rebalancing
## **************************************************
## 
## Call:
## optimize.portfolio.rebalancing(R = ., portfolio = min_var_portfolio, 
##     optimize_method = "ROI", rebalance_on = "quarters", training_period = 60, 
##     rolling_window = 60)
## 
## Number of rebalancing dates:  20 
## First rebalance date:
## [1] "2017-03-31"
## Last rebalance date:
## [1] "2021-12-31"
## 
## Annualized Portfolio Rebalancing Return:
## [1] 0.1815471
## 
## Annualized Portfolio Standard Deviation:
## [1] 0.221952

A saída da função mostra que foram feitos, durante o período de 2012-2021, um total de 20 rebalanceamentos, sendo o primeiro no primeiro trimestre de 2017. Como o treinamento exige 60 meses, precisamos acumular 60 meses de informações para realizar o primeiro rebalanceamento. No trimestre seguinte, os meses do primeiro trimestre de 2017 são incluídos no treinamento, e informações do primeiro trimestre de 2012 são excluídos.

Com a função chart.Weights(otimizacao_sharpe_rebalance) podemos visualizar como os pesos mudaram ao longo dos trimestres. Como ela não gera uma visualização das mais bonitas, vamos recorrer aos gráficos do ggplot2. Os pesos trimestrais podem ser visualizados na figura 4.

extractWeights(otimizacao_risco_rebalance) %>% 
  as.data.frame() %>% 
  rownames_to_column("data") %>% 
  mutate(data = as.Date(data)) %>% 
  pivot_longer(cols = "MGLU3.SA":"WEGE3.SA") %>% 
  ggplot(aes(x = data, y = value, fill = name)) +
  labs(fill = "", x = "", y = "Peso do ativo",
       title = "Peso dos ativos ao longo dos trimestres - Minimização do Risco") +
  geom_col() +
  scale_x_date(date_breaks = "6 months") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90))

Figure 4: Peso dos ativos ao longo dos trimestres

Interessante notar como a Ambev perde participação no portfólio ao longo do tempo. Movimento inverso é visto com a empresa WEG.

Maximização de Retorno com Rebalanceamento

A seguir, o mesmo exercício para diferentes configurações de portfólio. A figura 5 mostra os pesos para uma estratégia que maximiza retorno:

otimizacao_retorno_rebalance <- tbl_retornos_wide %>% 
  optimize.portfolio.rebalancing(
    optimize_method = "ROI",
    portfolio = portfolio_retorno_maximo,
    rebalance_on = "quarters",
    training_period = 60,
    rolling_window = 60
  )

extractWeights(otimizacao_retorno_rebalance) %>% 
  as.data.frame() %>% 
  rownames_to_column("data") %>% 
  mutate(data = as.Date(data)) %>% 
  pivot_longer(cols = "MGLU3.SA":"WEGE3.SA") %>% 
  ggplot(aes(x = data, y = value, fill = name)) +
  labs(fill = "", x = "", y = "Peso do ativo",
       title = "Peso dos ativos ao longo dos trimestres - Maximização do Retorno") +
  geom_col() +
  scale_x_date(date_breaks = "6 months") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90))

Figure 5: Peso dos ativos ao longo dos trimestres - Maximização do retorno

Nenhuma surpresa que uma estratégia que maximiza retorno deixa o portfólio exposto a Magazine Luiza.

Função Objetiva com Aversão ao Risco

Por fim, vamos tentar um portfólio diferente. Vamos estimar um portfólio que maximizar uma função utilidade quadrática. Neste caso, a função objetivo tem dois termos: um para o retorno médio do portfólio, e outro para a variância do portfólio, de modo muito semelhante ao Risco de Sharpe. Contudo, vamos incluir um parâmetro de aversão ao risco, \(\lambda > 0\) no termo de minimização de risco. Um investidor tomador de riscos teria um valor de \(\lambda\) baixo, enquanto que um investidor averso ao risco teria um \(\lambda\) elevado. Se \(\lambda = 0\), o termo de risco é ignorado e voltamos para o caso de maximização de retorno. Assim, \(\lambda\) representa um trade-off entre risco e retorno esperado.

A figura 6 mostra os pesos para um portfólio com \(\lambda=10\).

portfolio_quadratico <- portfolio.spec(assets = tickers) %>% 
  add.constraint(type = "full_investiment") %>% 
  add.constraint(type = "box", min = 0.05, max = 0.6) %>% 
  add.objective(type = "return", name = "mean") %>% 
  add.objective(type = "risk", name = "var", risk_aversion = 10)

otimizacao_quadratico_rebalance <- tbl_retornos_wide %>% 
  optimize.portfolio.rebalancing(
    optimize_method = "ROI",
    portfolio = portfolio_quadratico,
    rebalance_on = "quarters",
    training_period = 60,
    rolling_window = 60
  )

extractWeights(otimizacao_quadratico_rebalance) %>% 
  as.data.frame() %>% 
  rownames_to_column("data") %>% 
  mutate(data = as.Date(data)) %>% 
  pivot_longer(cols = "MGLU3.SA":"WEGE3.SA") %>% 
  ggplot(aes(x = data, y = value, fill = name)) +
  labs(fill = "", x = "", y = "Peso do ativo",
       title = "Peso dos ativos ao longo dos trimestres - Função Quadrática") +
  geom_col() +
  scale_x_date(date_breaks = "6 months") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90))

Figure 6: Peso dos ativos ao longo dos trimestres - Quadrática

Comparando a perfomance:

A figura 7 mostra o retorno acumulado entre 2017 (data do primeiro rebalanceamento) e 2021 para as três funções objetivas.

retorno_rebalance_1 <- tbl_retornos_wide %>% 
  Return.portfolio(weights = extractWeights(otimizacao_risco_rebalance)) %>% 
  as.data.frame() %>% 
  rename(min_risco = 1)

retorno_rebalance_2 <- tbl_retornos_wide %>% 
  Return.portfolio(weights = extractWeights(otimizacao_retorno_rebalance)) %>% 
  as.data.frame() %>% 
  rename(max_retorno = 1)

retorno_rebalance_3 <- tbl_retornos_wide %>% 
  Return.portfolio(weights = extractWeights(otimizacao_quadratico_rebalance)) %>% 
  as.data.frame() %>% 
  rename(quadratico = 1)


retorno_rebalance <- cbind(retorno_rebalance_1, retorno_rebalance_2, 
                           retorno_rebalance_3)

retorno_rebalance %>% 
  rownames_to_column("data") %>% 
  pivot_longer(cols = 2:4) %>% 
  mutate(data = as.Date(data)) %>% 
  group_by(name) %>% 
  mutate(retorno_acumulado = cumsum(value)) %>% 
  ggplot(aes(x = data, y = retorno_acumulado, color = name)) +
  geom_line(size = 1.5) +
  facet_wrap(~name) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "", y = "Retorno acumulado",
       title = "Retorno acumulado por portfólio") +
  theme(legend.position = "none")

Figure 7: Retorno acumulado por portfólio