Como seria uma carteira com as “melhores” ações?

Um simples exemplo de otimização de carteira de ações.

A carteira de Markowitz

A teoria moderna da carteira (MPT) afirma que os investidores são avessos ao risco e, dado um nível de risco, escolherão as carteiras que mais oferecem retorno. Para isso, precisamos otimizar os portfólios.

Para realizar a otimização, precisaremos:

  1. Baixar os dados de preços dos ativos,
  2. Calcular os retornos médios para o período desejado e
  3. Atribuir pesos aos ativos para, em seguida, usá-los para construir uma fronteira eficiente.

Ações que mais subiram em setembro de 2020

Carregamento dos pacotes

library(tidyquant) # Para obter os dados
library(timetk)    # Para manipular as séries temporais
library(tidyverse) # Nosso velho amigo

Vamos obter os dados de preços das ações.

acoes <- c('AZUL4.SA', 'GGBR4.SA', 'PCAR3.SA', 'QUAL3.SA', 'RENT3.SA')

inicio <- '2020-09-01'

fim <- '2020-09-30'

# tq(get() carrega dados no formato 'tibble'
precos <- tq_get(acoes,
                 from = inicio,
                 to = fim,
                 get = 'stock.prices')

precos

Cálculo dos retornos (logarítimicos - contínuos).

Nota: Tente tq_transmute_fun_options() para conhecer as opções. São muitas!

Para a função periodReturn() veja este link.

ret_log_tidy <- precos %>%
    group_by(symbol) %>%
    tq_transmute(select = adjusted,
                 mutate_fun = periodReturn,  
                 period = 'daily', # retornos diários
                 col_rename = 'ret',
                 type = 'log')

ret_log_tidy

Os dados estão em formato tibble, ou long. Usaremos a função tidyr::pivot_wider() para convertê-los em um formato wide. Também vamos convertê-los em um objeto de série temporal usando a função tk_xts() e calcular as médias dos retornos diários.

ret_log_xts <- ret_log_tidy %>%
    pivot_wider(
      names_from = symbol, 
      values_from = ret) %>%
    tk_xts(date_var = date)

ret_log_xts
##                AZUL4.SA     GGBR4.SA     PCAR3.SA     QUAL3.SA     RENT3.SA
## 2020-09-01  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
## 2020-09-02  0.009989182 -0.010555402 -0.005877878  0.027709567  0.012814297
## 2020-09-03  0.036074398 -0.020935241 -0.026567586 -0.021026901 -0.007076899
## 2020-09-04  0.024702779  0.027482738  0.007934094  0.073949274  0.001183070
## 2020-09-08  0.065707013 -0.011613260 -0.012404770  0.013748934  0.057059676
## 2020-09-09 -0.001524429  0.041776122 -0.005455696  0.002321607  0.020995328
## 2020-09-10 -0.004587126 -0.033683012  0.138046545 -0.033601831 -0.055262662
## 2020-09-11  0.018223739  0.003018037  0.038897341 -0.033529421 -0.025255286
## 2020-09-14  0.060934404  0.018412947 -0.003917282  0.027214696  0.044314171
## 2020-09-15 -0.007818057  0.056089510 -0.000676784 -0.005748103 -0.000189333
## 2020-09-16  0.037119072 -0.012195228  0.016520299  0.003634247 -0.002654027
## 2020-09-17 -0.006206916  0.011728993  0.015465341 -0.013390399  0.005112142
## 2020-09-18 -0.024865543 -0.022642522 -0.012406120 -0.020740591 -0.004732506
## 2020-09-21 -0.081225475 -0.017810450 -0.010010282  0.006857845 -0.036522258
## 2020-09-22  0.010331066  0.006293846  0.008415489 -0.007170757  0.018137902
## 2020-09-23  0.018106870 -0.003384007 -0.040722338 -0.022785725  0.130797734
## 2020-09-24  0.003731348  0.001451889  0.036055764  0.048424309 -0.016758309
## 2020-09-25 -0.002610443 -0.008254676 -0.010613213  0.013931139  0.017266939
## 2020-09-28 -0.031483729  0.001461737 -0.043758249 -0.029606854 -0.040645544
## 2020-09-29 -0.080203283 -0.016196759 -0.003533965  0.001547721  0.002468389
autoplot(ret_log_xts) +
  ggtitle('Retornos diários (log)') +
  xlab('dia')

ret_medio <- colMeans(ret_log_xts)
print(round(ret_medio, 5))
## AZUL4.SA GGBR4.SA PCAR3.SA QUAL3.SA RENT3.SA 
##  0.00222  0.00052  0.00427  0.00159  0.00605

Em seguida, vamos calcular a matriz de covariância para todas essas ações. Vamos anualizar os valores multiplicando por 252 (dias úteis no ano).

Passo-a-passo

Antes de aplicarmos nossos métodos a milhares de portfólio aleatórios, vamos demonstrar os passos em um único portfólio.

Para calcular os retornos e riscos da carteira (medidos pelo desvio-padrão) precisaremos calcular:

  • Retornos médios de ativos
  • Pesos do portfólio
  • A matriz de covariância de todos os ativos
  • Os pesos aleatórios

Vamos calcular a matriz de covariância e criar os pesos aleatórios.

mat_cov <- cov(ret_log_xts) * 252
AZUL4.SA GGBR4.SA PCAR3.SA QUAL3.SA RENT3.SA
AZUL4.SA 0.3563967 0.0348715 0.0110948 0.0538902 0.1635244
GGBR4.SA 0.0348715 0.1234847 -0.0420649 0.0447650 0.0403520
PCAR3.SA 0.0110948 -0.0420649 0.3586800 -0.0210146 -0.1821512
QUAL3.SA 0.0538902 0.0447650 -0.0210146 0.1924187 0.0323459
RENT3.SA 0.1635244 0.0403520 -0.1821512 0.0323459 0.3965526 
pesos <- runif(n = length(acoes))
print(pesos)
## [1] 0.7209039 0.8757732 0.7609823 0.8861246 0.4564810
print(sum(pesos))
## [1] 3.700265

Vamos tornar a soma dos pesos igual a 1

pesos <- pesos/sum(pesos)
print(pesos)
## [1] 0.1948249 0.2366785 0.2056562 0.2394760 0.1233644
print(sum(pesos))
## [1] 1

O retorno anualizado da carteira (252 dias úteis)

\[ \begin{aligned} & \bar r_{diário} = \sum_i^ n \, (peso_i \,\times \, \bar r_i) \\ & r_{carteira} = (\bar r_{diário} \, + \, 1)^{252} \, - \, 1 \end{aligned} \]

ret_carteira <- (sum(pesos * ret_medio) + 1)^252 - 1
ret_carteira * 100
## [1] 90.50473

O retorno anualizado também pode ser aproximado considerando o retorno diário (logarítimico) como um retorno contínuo, utilizando a fórmula:

\[r_{carteira} = e^ { \; (\bar r_{diário} \, \times \, {252})} \; - \, 1\]

ret_carteira_continuo <- exp(sum(pesos * ret_medio) * 252) - 1
ret_carteira_continuo * 100
## [1] 90.66194

Em seguida, vamos calcular o risco da carteira (desvio-padrão).

Um pouco de álgebra linear não faz mal a ninguém.

\[\sigma_{carteira} = (p^t \, \Sigma \, p)^{1/2}\]

risco_carteira <- sqrt(t(pesos) %*% (mat_cov %*% pesos))

# Apresentação mais elegante no console utilizando o pacote glue.
glue::glue("Risco da carteira: {round(risco_carteira*100, 2)}", "%")
## Risco da carteira: 25.22%

Utilizando o índice de Sharpe

rf <- 0 # Renda Fixa

indice_sharpe <- (ret_carteira - rf) / risco_carteira

glue::glue("Índice Sharpe da carteira: {round(indice_sharpe*100, 2)}", "%")
## Índice Sharpe da carteira: 358.84%

Temos tudo o que precisamos para realizar a otimização. O que precisamos fazer agora é executar este código em 5000 carteiras aleatórias. Para isso usaremos um loop.

Antes de fazermos isso, precisamos criar vetores vazios e uma matriz para armazenar nossos valores.

num_cart <- 5000

# Cria uma matriz para armazenar os pesos

todos_pesos <- matrix(nrow = num_cart,
                      ncol = length(acoes))

# Cria um vetor vazio para armazenar os retornos da carteira

ret_carteira <- vector('numeric', length = num_cart)

# Cria um vetor vazio para armazenar os desvios-padrão da carteira 

risco_carteira <- vector('numeric', length = num_cart)

# Cria um vetor vazio para armazenar os índices de Sharpe da carteira 

indice_sharpe <- vector('numeric', length = num_cart)

Um loop 5000 vezes:

id <- diag(length(acoes))

for (i in seq_along(ret_carteira)) {
    
    if (i <= length(acoes)) # carteiras com apenas uma ação
      pesos <- id[i,]
    else  
      pesos <- runif(length(acoes))
    
    pesos <- pesos/sum(pesos)
   
    # Guarda os pesos em uma matriz
    todos_pesos[i,] <- pesos
    
    # Retornos das carteiras
    
    port_ret <- sum(pesos * ret_medio)
    port_ret <- ((port_ret + 1)^252) - 1
    
    # Salvando os retornos das carteiras
    ret_carteira[i] <- port_ret
    
    # Cria e salva os riscos das carteiras
    port_sd <- sqrt(t(pesos) %*% (mat_cov  %*% pesos))
    risco_carteira[i] <- port_sd
    
    # Cria e salva os índices de Sharpe das carteiras
    
    i_sharpe <- (port_ret - rf) / port_sd
    indice_sharpe[i] <- i_sharpe
    
}

Vamos armazenar os valores, converter a matriz em um ´tibble`, alterar os nomes das variáveis e combinar os valores

valores_carteira <- tibble(Retorno = ret_carteira,
                           Risco = risco_carteira,
                           IndiceSharpe = indice_sharpe)


todos_pesos <- tk_tbl(todos_pesos)

colnames(todos_pesos) <- colnames(ret_log_xts)

valores_carteira <- tk_tbl(cbind(todos_pesos, valores_carteira))

valores_carteira

Temos os pesos de cada ativo com os riscos e os retornos, juntamente com o índice Sharpe de cada carteira.

Em seguida, vamos olhar para as carteiras que mais importam:

  • A carteira de variância mínima
  • A carteira de tangência (a carteira com maior índice sharpe)
min_var <- valores_carteira[which.min(valores_carteira$Risco),]
max_i_sharpe <- valores_carteira[which.max(valores_carteira$IndiceSharpe),]

Vamos plotar os pesos de cada carteira. Primeiro com a carteira de variância mínima.

Como podemos observar, a carteira de variância mínima quase não tem alocação para as ações da Azul. A maior parte da carteira é investida em ações da Gerdau.

Em seguida, vamos olhar para a carteira de tangência ou a carteira com o maior índice Sharpe.

Que diferença! Localiza e Pão de Açúcar agora dominam completamente a carteira e Gerdau tem uma pequena participação.

Vamos calcular novamente o retorno e o risco anual de cada ação individualmente:

nomes <- c('Azul', 'Gerdau', 'Localiza', 'PAO', 'Qualicorp')
ativos <- data.frame(Acoes = nomes, 
                     Retorno = (colMeans(ret_log_xts, na.rm = TRUE)+ 1)^252 - 1, 
                     Risco = sapply(ret_log_xts, sd, na.rm = TRUE) * sqrt(252))
ativos

A fronteira eficiente

Segundo a teoria moderna de porfolio, desenvolvida por Markowitz em 1952, podemos definir o problema de seleção de uma carteira como o seguinte problema de otimização quadrática:

\[ \begin{aligned} \underset {p} {min} &= p^T \; \Sigma \; p \\ s.a. \\ p^T \hat \mu &= \bar r \\ p^T \; 1 &= 1 \end{aligned} \]

A expressão indica que minimizamos o risco induzido pela variância-covariância \(\small \sigma^2 = p^T \; \Sigma \; p\), onde a matriz \(\small \Sigma\) é uma estimativa da covariância entre os ativos da carteira. O vetor p denota as participações (pesos) dos investimentos individuais sujeitos à condição \(\small p^T 1=1\), ou seja, que o capital disponível está totalmente investido. O retorno esperado ou retorno-alvo \(\small \bar r\) é expresso pela condição \(\small p^T \hat \mu = \bar r,\) onde o vetor n-dimensional \(\small \hat \mu\) estima a média esperada dos ativos. O portfólio com vendas curtas (short-selling) ilimitadas podem ser resolvidos analiticamente. No entanto, se os pesos não podem ser negativos, o que significa a impossibilidade de vendas curtas, então a otimização tem que ser feita numericamente. A estrutura do problema é quadrática e, portanto, podemos usar um otimizador quadrático para calcular os pesos da carteira.

Consideremos o seguinte problema padrão:

\[ \begin{aligned} \underset {p} {min} &= p^T \; \Sigma \; p \\ s.a. \\ A \; p & \leq b \\ p^T 1 &= 1 \end{aligned} \]

Pode-se mostrar que, se \(\small \Sigma\) é uma matriz definida positiva, o problema do portfólio de Markowitz é um problema de otimização convexa. Como tal, as soluções ótimas locais também são soluções ótimas globais. O pacote quadprog fornece a função solve.QP(), que provê uma interface uma subrotina FORTRAN. Esta subrotina implementa o método dual de Goldfarb e Idnani para resolver problemas quadráticos de programação da forma \(\small min(-d^T b + 1/2 \; b^TD \; b)\) com as restrições \(\small A^T b \geq b_0\).

Utilizaremos essa função para a otimização do portfólio.

Duas funções auxiliares

pesosCarteira <- function(retornosAtivos, retornoAlvo)
{
## Argumentos:
# retornosAtivos - conjunto de dados dos retornos dos ativos
# retornoAlvo - o retorno-alvo da carteira 
  
##  A função solve.QP() do pacote quadprog implementa o método dual de Goldfarb e Idnani (1982, 1983) 
##  para a solução do problema de otimização quadrática na forma 
##  min(-d'b + 1/2 b' Db) com as restrições A'T b >= b0. 
    
## Para detalhes, veja D. Goldfarb and A. Idnani (1983). "A numerically stable dual method for solving strictly convex #quadratic programs". Mathematical Programming, 27, 1–33.    
    
## A solução aqui são os pesos que minimizam o risco para o retorno em 'retornoAlvo'    
    
    if(!require("quadprog")) install.packages("quadprog")
    suppressMessages(suppressWarnings(library(quadprog)))
    
    nAtivos  <-  ncol(retornosAtivos)
    portfolio <- solve.QP(
        Dmat <- cov(retornosAtivos),                        # matriz D
        dvec <- rep(0, times = nAtivos),                    # vetor  d
        Amat <- t(rbind(retorno = colMeans(retornosAtivos), # matriz A de restrições
                        orcamento = rep(1, nAtivos), 
                        longa = diag(nAtivos))),
        bvec <- c(retorno = retornoAlvo,                    # vetor  b0
                  orcamento = 1,
                  longa = rep(0, times = nAtivos)),
        meq = 2)                                            # as primeiro meq restrições são igualdades
    
    pesos  <-  portfolio$solution # vetor contendo a solução do problema
    pesos
}

fronteiraCarteira <- function(retornosAtivos, nPontos = 40)
{
    # Quantidade de ativos
    nAtivos <- ncol(retornosAtivos)
    # Retornos-alvo
    mu <- colMeans(retornosAtivos)
    retornoAlvo <- seq(min(mu), max(mu), length = nPontos)
    # Pesos ótimos
    pesos <- rep(0, nAtivos)
    pesos[which.min(mu)] <- 1
    for (i in 2:(nPontos-1)) {
        novosPesos <- pesosCarteira(retornosAtivos, retornoAlvo[i])
        pesos <- rbind(pesos, novosPesos)
    }
    novosPesos <- rep(0, nAtivos)
    novosPesos[which.max(mu)] <- 1
    pesos <- rbind(pesos, novosPesos)
    pesos <- round(pesos, 4)
    colnames(pesos) <- colnames(retornosAtivos)
    rownames(pesos) <- 1:nPontos
    
    # Valor do retorno
    pesos
}

Visualizando todas as carteiras e a fronteira eficiente

retornosAtivos <- ret_log_xts 
names(retornosAtivos)
## [1] "AZUL4.SA" "GGBR4.SA" "PCAR3.SA" "QUAL3.SA" "RENT3.SA"
pesos_front <- fronteiraCarteira(retornosAtivos, nPontos = 200)
as.data.frame(pesos_front)
mu  <-  colMeans(retornosAtivos)
retornoAlvos  <-  seq(min(mu), max(mu), length = nrow(pesos_front))

riscosAlvo  <-  NULL
for (i in 1:nrow(pesos_front)) {
    novoRiscoAlvo  <-  sqrt(pesos_front[i, ] %*% 
                            cov(retornosAtivos) %*%
                            pesos_front[i, ])
    riscosAlvo  <-  c(riscosAlvo, novoRiscoAlvo)
}

fronteiraEficiente <- data.frame(risco = riscosAlvo, retorno = retornoAlvos)

fronteiraEficiente
# Anualizaremos os retornos e os riscos das carteiras na fronteira eficiente
fronteiraEficiente$retorno <- (fronteiraEficiente$retorno+1)^252 - 1
fronteiraEficiente$risco <- fronteiraEficiente$risco * sqrt(252)
library(RColorBrewer)

p <- ggplot() +
      geom_point(data = valores_carteira, aes(x = Risco, y = Retorno, color = IndiceSharpe)) +
      #theme_bw() +
      scale_color_gradient(low="blue", high="yellow") +
      scale_y_continuous(labels = scales::percent) +
      scale_x_continuous(labels = scales::percent, limits = c(.15, .65)) +
      geom_point(data = min_var, aes(x = Risco, y = Retorno),  color = 'red', size = 3.5) +
      geom_point(data = max_i_sharpe, aes(x = Risco, y = Retorno),  color = 'red', size = 3.5) +
      geom_point(data = valores_carteira[1:5,], aes(x = Risco, y = Retorno),  color = 'red', size = 3.5) +
      geom_point(data = fronteiraEficiente, aes(x = risco, y = retorno),  color = 'black', size = 1.0) +
      annotate('text', x = max_i_sharpe$Risco-.05, y = max_i_sharpe$Retorno, label = "Max-Sharpe", size = 3.5) +
      annotate('text', x = min_var$Risco-.035, y = min_var$Retorno, label = "Min-Var", size = 3.5) +
      annotate('text', x = ativos$Risco-.035, y = ativos$Retorno+.021, label = nomes, size = 3.5) +
      labs(x = 'Risco Anualizado',
           y = 'Retorno Anualizado',
           title = "Fronteira Eficiente") 
p

No gráfico podemos observar todas as quase 5000 carteiras, a fronteira eficiente, a carteira de mínima variância e aquela com o melhor índice Sharpe. As ações individuais também são apresentadas no espaço risco-retorno. Como mencionado acima, um investidor avesso ao risco exigirá um maior retorno para um determinado nível de risco. Em outras palavras, tentará obter carteiras que residem na fronteira eficiente.

Referência

  • “Basic R for Finance”, RMetrics Association & Finance Online Publishing, Capítulo 32