1 BIBLIOTECAS

library(magrittr)
library(fst)
library(TTR)
library(caret)
library(FSelector)

2 FUNÇÕES AUXILIARES

Para um elemento na posição \(i\) de x, a função retorna o valor da função f aplicada no subconjunto x[(i-n+1):i] de x.

movMeas <- function(x,f,n,...){
  N <- length(x)
  
  out <- rep(NA,n-1)
  for (i in n:N) {
    v <- matrix(x[(i-n+1):i],nrow = 1)
    out <- c(out,
             apply(v,1,f,...))
  }
  out
}

A função retorna o vetor x com a substituição de cada valor NA pelo valor imediatamente anterior. Caso o valor anterior seja NA, não há alteração.

attrNA <- function(x){
  for (i in 2:length(x)) {
    if (is.na(x[i])) {
      x[i] <- x[i-1]
    }
  }
  return(x)
}

A função retorna uma matriz que condensa os preços de todos os ativos, onde cada linha corresponde a uma data e cada coluna, a um ativo.

genMatPrice <- function(prices_list){
  num_tick <- length(prices_list)
  tickers <- names(prices_list)
  
  # Sincronização ----
  datas <- prices_list[[1]][,1] %>% 
    as.character
  for (i in 2:num_tick) {
    datas <- c(datas,
               prices_list[[i]][,1] %>% 
                 as.character) %>% 
      unique
  }
  datas <- datas %>% 
    as.Date %>% 
    sort %>% 
    as.character
  
  # Merging ----
  mat_prices_list <- NULL
  for (i in 1:num_tick) {
    mat_prices_list <- cbind(mat_prices_list,
                             ifelse(datas %in% (prices_list[[i]][,1] %>% 
                                                  as.character),
                                    prices_list[[i]][datas,2],
                                    NA))
  }
  
  rownames(mat_prices_list) <- datas
  colnames(mat_prices_list) <- tickers
  
  # Atribuição de valores faltantes ----
  mat_prices_list <- mat_prices_list %>% apply(2,attrNA)
  mat_prices_list <- mat_prices_list %>% na.omit
  
  datas <- rownames(mat_prices_list)
  
  mat_prices_list
}

A função retorna uma matriz com os rendimentos (variações de preço) dentro do spread setado. Cada linha corresponde a uma data e cada coluna, a um ativo. O retorno é logarítmico:

\[ rcc = \ln \left( \frac{P[D+s]}{P[D+0]} \right) \]

Na fórmula acima, \(s\) é o spread e \(P[x]\) é o preço do ativo no dia \(x\).

genMatRend <- function(mat_prices_list,spread_DC){
  datas <- rownames(mat_prices_list)
  
  # Matriz de Rendimentos ----
  fcs <- datas %>% 
    as.Date %>% 
    diff %>% 
    mean %>% 
    as.numeric
  
  n <- round(spread_DC/fcs)
  
  rend <- mat_prices_list %>%
    ROC(n=n)
  
  rend
}

Sendo \(\widetilde{X}\) o vetor x normalizado, a função retorna o resultado de:

\[ \widetilde{X} = \frac{x - \bar{x}}{s_x}, \] onde \(\bar{x}\) é a média de x e \(s_x\) é o desvio-padrão de x.

normalize <- function(x) (x - mean(x))/sd(x)

O processamento ocorre nos moldes de normalize.

last_norm <- function(x) (tail(x,1) - mean(x))/sd(x)

São selecionados todos os argumentos com ganho de informação superior ao ganho médio dos preditores. Para conceitos sobre ganho de informação e entropia: https://victorzhou.com/blog/information-gain/.

featSel <- function(df){
  feat_imp <- information.gain(outcome ~ .,df) %>%
    apply(2,sort,decreasing=T)
  
  mi_imp <- quantile(feat_imp,0.4,names = F)
  pred_sel <- rownames(feat_imp)[feat_imp[,1]>=mi_imp]
  
  pred_sel
}

3 SETUP

O preço utilizado será o de fechamento e o spread para avaliação dos retornos será de 180 dias corridos.

type_price <- "close"
spread_DC <- 180

4 CARREGAMENTO DOS DADOS

Os preços dos ativos foram carregados para uma pasta chamada fst no diretório de trabalho. O script de carregamento usado é o que segue:

rm(list = ls())

setwd("~/fst")

library(magrittr)
library(fst)
library(quantmod)

acoes <- c("abev3","b3sa3","bbas3","bbdc4","bbse3","bova11","brfs3","brkm5",
           "brml3","ccro3","ciel3","cmig4","cogn3","csan3","csna3","cyre3",
           "ecor3","egie3","elet3","embr3","eqtl3","ggbr4","goau4","hype3",
           "itsa4","itub4","jbss3","lame4","lren3","mrfg3","petr4",
           "qual3","radl3","rail3","rent3","sbsp3","taee11","usim5",
           "vale3","vivt3","vvar3","wege3","yduq3",
           "brap4","enbr3","flry3","klbn4","mult3","leve3","sanb11",
           "ugpa3")

tickers <- acoes %>% casefold(upper = T) %>% paste0(".SA")

for (i in 1:length(tickers)) {
  X <- getSymbols(tickers[i],
                  to = "2020-12-31",
                  auto.assign = F)
  X <- X %>% na.omit
  datas <- index(X)
  
  X <- as.data.frame(X)
  names(X) <- c("open","high","low","close","volume","close.adj")
  X$date <- datas
  
  X <- X[X$volume>0,]
  
  if (nrow(X)>=3000) {
    write_fst(X,path = tickers[i],compress = 50)
    cat(tickers[i]);cat(" - ");cat(datas %>% as.character %>% tail(1));cat(" // ")
  }
}

A data-limite para obtenção dos dados foi setada como 31/12/2020. Foram excluídas as instância com volume 0, além de excluir os ativos com menos de 3000 pregões válidos.

Os dados do IBovespa seguiram uma dinâmica parecida, mas foram salvos diretamente no diretório de trabalho:

rm(list = ls())

library(magrittr)
library(fst)
library(quantmod)

X <- getSymbols("^BVSP",auto.assign = F)
X <- X %>% na.omit
datas <- index(X)

X <- as.data.frame(X)
names(X) <- c("open","high","low","close","volume","close.adj")
X$date <- datas

X <- X[X$volume>0,]

write_fst(X,path = "IBOVESPA",compress = 50)
cat("IBOVESPA");cat(" - ");cat(datas %>% as.character %>% tail(1));cat(" // ")
aim_price <- list()

precos <- read_fst("IBOVESPA")

aim_price[["IBOVESPA"]] <- precos[,c("date",type_price)]
rownames(aim_price[["IBOVESPA"]]) <- aim_price[["IBOVESPA"]]$date

tickers <- list.files("fst/")
for (ticker in tickers) {
  precos <- read_fst(paste0("fst/",ticker))
  aim_price[[ticker]] <- precos[,c("date",type_price)]
  rownames(aim_price[[ticker]]) <- aim_price[[ticker]]$date
}

4.1 Matriz de Rendimentos

De posse dos preços, é possível calcular a matriz de rendimentos, que é a matéria-prima da análise:

price_stocks <- genMatPrice(aim_price)
rend_stocks <- genMatRend(price_stocks,spread_DC)

N <- nrow(price_stocks)

fcs <- price_stocks %>% 
  rownames %>% 
  as.Date %>% 
  diff %>% 
  mean %>% 
  as.numeric

spread <- round(spread_DC/fcs)
>  Spread, em pregões:
>  122

A variável fcs é o tempo médio, em dias corridos, entre dois pregões consecutivos. Ela permite a conversão entre o spread em dias corridos para o spread em pregões.

5 FEATURES

Os preditores serão calculados a partir dos preços e rendimentos dos ativos, e serão divididos em 5 grupos: Indicadores Técnicos, Mercado, Ativo, Artificiais de Preço e Artificiais de Retorno. Para exemplificar, o desenvolvimento será feito para a PETR4, mas é o mesmo para todos os outros ativos.

5.1 Preditores de Indicadores Técnicos

Serão usados os seguintes períodos para o MACD: 12, para a média móvel rápida; 26, para a média móvel lenta; e 9, para o sinal.

Será usado o RSI de 14 períodos.

Serão usadas duas MM, uma de período igual ao spread em pregões e outra com período igual à metade do spread em pregões.

Será usado o CTI de 20 períodos.

ticker <- "PETR4.SA"
P <- price_stocks[,ticker]

macd <- MACD(P)
rsi <- RSI(P)
mmr <- SMA(P,round(spread/2))
mml <- SMA(P,spread)
cti <- CTI(P)

colnames(mmr) <- "mmr"
colnames(mml) <- "mml"

5.2 Preditores de Mercado

Serão usados o retorno do IBOV, relativo ao spread, e a sua volatilidade dos últimos \(2 \times spread\) pregões. A volatilidade será quantificada pelo desvio-padrão:

roc_ibov <- rend_stocks[,"IBOVESPA"]
volat_ibov <- roc_ibov %>% movMeas(f=sd,n=2*spread)

roc_ibov <- matrix(roc_ibov,ncol=1)
volat_ibov <- matrix(volat_ibov,ncol=1)

colnames(roc_ibov) <- "roc_ibov"
colnames(volat_ibov) <- "volat_ibov"

5.3 Preditores do Ativo

Serão os mesmos conceitos dos preditores de mercado:

roc_acao <- rend_stocks[,ticker]
volat_acao <- roc_acao %>% movMeas(f=sd,n=2*spread)

roc_acao <- matrix(roc_acao,ncol=1)
volat_acao <- matrix(volat_acao,ncol=1)

colnames(roc_acao) <- "roc_acao"
colnames(volat_acao) <- "volat_acao"

5.4 Preditores Artificiais de Preço

Serão as primeira e segunda derivadas discretas do preço do ativo e o preço normalizado com os dados (média e desvio-padrão) dos últimos spread pregões:

dP <- c(NA,P %>% diff) %>% matrix(ncol = 1)
d2P <- c(NA,NA,P %>% diff(differences = 2)) %>% matrix(ncol = 1)

zP <- P %>% movMeas(f=last_norm,n=spread) %>% matrix(ncol=1)

colnames(dP) <- "dP"
colnames(d2P) <- "d2P"
colnames(zP) <- "zP"

5.5 Preditores Artificiais de Retorno

Além das derivadas primeira e segunda e do valor normalizado aplicados ao retorno, como no tipo anterior, serão adicionadas o valor do retorno ao quadrado e as médias móveis dos retornos com períodos spread e spread/2:

drcc <- c(NA,roc_acao[,1] %>% diff) %>% matrix(ncol = 1)
d2rcc <- c(NA,NA,roc_acao[,1] %>% diff(differences = 2)) %>% matrix(ncol = 1)

rcc2 <- roc_acao^2

mmr_rcc <- roc_acao %>% SMA(n=round(spread/2)) %>% matrix(ncol=1)
mml_rcc <- roc_acao %>% SMA(n=spread) %>% matrix(ncol=1)

zrcc <- roc_acao[,1] %>% movMeas(f=last_norm,n=spread) %>% matrix(ncol=1)

colnames(drcc) <- "drcc"
colnames(d2rcc) <- "d2rcc"
colnames(rcc2) <- "rcc2"
colnames(mmr_rcc) <- "mmr_rcc"
colnames(mml_rcc) <- "mml_rcc"
colnames(zrcc) <- "zrcc"

6 MONTAGEM DO DATASET

Todas os preditores deverão ser defasados em spread linhas do target (os retornos do ativo), para simular a situação de previsão em D+spread:

df0 <- cbind(macd,rsi,mmr,mml,cti,
             roc_ibov,volat_ibov,
             roc_acao,volat_acao,
             dP,d2P,zP,
             drcc,d2rcc,rcc2,mmr_rcc,mml_rcc,zrcc)

df <- data.frame(outcome = roc_acao[(1+spread):N],
                 df0[1:(N-spread),]) %>% na.omit

6.1 Análise de Correlação

A fim de evitar processamento desnecessário mantendo preditores altamente correlacionados, é interessante analisar a matriz de correlação do dataset df:

mat_cor <- cor(df[,-1])
diag(mat_cor) <- 0
high_cor <- which(mat_cor>=0.9,arr.ind = T)
names(df[,-1])[high_cor] %>% matrix(ncol=2)
>       [,1]     [,2]    
>  [1,] "signal" "macd"  
>  [2,] "macd"   "signal"
>  [3,] "mml"    "mmr"   
>  [4,] "mmr"    "mml"
mat_cor[high_cor]
>  [1] 0.9571002 0.9571002 0.9735118 0.9735118
cor(df$outcome,df[,c("macd","signal","mmr","mml")])
>             macd     signal        mmr        mml
>  [1,] -0.1993699 -0.2092566 -0.3280258 -0.2698379

Nota-se altíssima correlação linear entre o MACD e o seu sinal, bem como entre as médias móveis de preços. Serão eliminadas a média móvel lenta e o MACD, devido a suas menores correlações absolutas com o target:

df$macd <- df$mml <- NULL

7 SIMULAÇÃO

As datas constantes na matriz de rendimento são (relativas a D+0):

datas <- df %>% 
  rownames %>% 
  as.Date
summary(datas)
>          Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
>  "2009-12-30" "2012-08-11" "2015-03-30" "2015-03-30" "2017-11-08" "2020-07-03"

Será considerado um intervalo de 5 anos de simulações:

data_final <- tail(datas,1)
dif <- (datas - (data_final - 5*365)) %>% as.numeric %>% abs
ini <- which(dif == min(dif))
data_inicial <- datas[ini]
data_inicial
>  [1] "2015-07-06"

A simulação iniciará em 06/07/2015, equivalente à linha 1364 do dataset (aproximadamente 52% das instâncias disponíveis serão usadas para o primeiro treinamento).

O procedimento de validação será time slice de janela móvel, onde um novo modelo será gerado com a adição da última instância à base de treino. Além disso, a seleção de features será realizada também cada iteração, simulando o que aconteceria em uma operação real.

Como o algoritmo exemplificado aqui será o kNN, será feita também, a cada iteração, uma normalização dos dados:

n <- nrow(df)
ini_teste <- ini

t0 <- Sys.time()

outpred <- NULL
for (i in ini_teste:n) {
  # normalização
  modelo_norm <- df[1:(i-1),-1] %>% preProcess
  
  df_treino <- cbind(outcome=df$outcome[1:(i-1)],
                     predict(modelo_norm,df[1:(i-1),-1]))
  
  # feature selection
  pred_sel <- featSel(df_treino)
  df1 <- df_treino[,c("outcome",pred_sel)]
  
  # modelagem
  modelo_knn <- train(outcome ~ .,
                      df1,
                      method = "knn",
                      tuneGrid = expand.grid(k=5),
                      metric = "MAE",
                      trControl = trainControl(method = "repeatedcv",
                                               repeats = 3,
                                               number = 10,
                                               summaryFunction = defaultSummary))
  # predição
  outpred <- c(outpred,
               predict(modelo_knn,
                       predict(modelo_norm,df[i,-1])))
}

t1 <- Sys.time()

deltaT <- as.numeric(t1) - as.numeric(t0)
if (deltaT < 60) {
  cat("Tempo de processamento (s): ")
  cat(deltaT)
}else{
  if (deltaT < 3600) {
    cat("Tempo de processamento (min): ")
    cat(deltaT/60)
  }else{
    cat("Tempo de processamento (h): ")
    cat(deltaT/3600)
  }
}
>  Tempo de processamento (min): 19.82468

7.1 Avaliação

Num primeiro momento, vamos avaliar o \(R^2\) e o MAE dos resultados obtidos:

r2 <- R2(outpred,df$outcome[ini_teste:n])

residuo <- outpred-df$outcome[ini_teste:n]
mae <- residuo %>% abs %>% mean
plot(residuo %>% density)

cat("R²: ")
>  R²:
cat(r2)
>  0.9278306
cat("MAE: ")
>  MAE:
cat(mae)
>  0.05352392