O que é?
Este documento Markdown apresenta uma análise estatística sobre predição de valor do ETF SPY, que representa o índice S&P500, da bolsa americana. A ideia principal é baseada na estratégia de debit spreads utilizando opções.
- O objetivo é entender com qual frequência a valorização do ETF SPY é acima de 1 USD em 5 dias (1 semana de pregão).
- Isso porque a estratégia é comprar um debit spread de SPY com 1 USD de intervalo (por exemplo, $300-$301) e entender com que frequência essa valorização acontece.
- Caso essa valorização aconteça, existe um lucro interessante que pode ser obtido com esse debit spread.
- Caso não aconteça, o prejuízo é limitado ao valor de entrada no debit spread.
Após isso, vou apresentar modelos preditivos que possam aumentar a nossa chance de entrar em trades com maior chance de sucesso, ou seja, maiores chances de que a valorização de SPY seja acima de 1 USD para aqueles 5 dias.
Biblioteca quantmod e dados de abertura e fechamento de pregões
Para os dados, vamos utilizar a biblioteca quantmod para baixar os dados do YahooFinance com a cotação histórica de SPY.
library(quantmod)
library(tidyverse)
spy = getSymbols("SPY", from = "2015-01-01", to = "2021-10-09", env = NULL) %>%
data.frame()
knitr::kable(head(spy))| SPY.Open | SPY.High | SPY.Low | SPY.Close | SPY.Volume | SPY.Adjusted | |
|---|---|---|---|---|---|---|
| 2015-01-02 | 206.38 | 206.88 | 204.18 | 205.43 | 121465900 | 180.9762 |
| 2015-01-05 | 204.17 | 204.37 | 201.35 | 201.72 | 169632600 | 177.7079 |
| 2015-01-06 | 202.09 | 202.72 | 198.86 | 199.82 | 209151400 | 176.0341 |
| 2015-01-07 | 201.42 | 202.72 | 200.88 | 202.31 | 125346700 | 178.2277 |
| 2015-01-08 | 204.01 | 206.16 | 203.99 | 205.90 | 147217800 | 181.3903 |
| 2015-01-09 | 206.40 | 206.42 | 203.51 | 204.25 | 158567300 | 179.9367 |
Apenas o valor de abertura e fechamentos vão nos interessar nesse momento. Pois a ideia é comprar na abertura do 1o dia e verificar o valor de SPY no fechamento ao 5o dia após a compra.
spy = spy %>%
mutate(
ds = rownames(spy),
open = `SPY.Open`,
close = `SPY.Close`,
.keep = "none"
)
knitr::kable(head(spy))| ds | open | close |
|---|---|---|
| 2015-01-02 | 206.38 | 205.43 |
| 2015-01-05 | 204.17 | 201.72 |
| 2015-01-06 | 202.09 | 199.82 |
| 2015-01-07 | 201.42 | 202.31 |
| 2015-01-08 | 204.01 | 205.90 |
| 2015-01-09 | 206.40 | 204.25 |
Calculo dos diferences alvos: retorno em 5 e 10 dias
Em seguida vamos analisar os retornos em 5 dias, criando a nossa variável alvo. Por questões de teste, vamos também calcular os retornos para 10 dias e ver se há potencial também:
spy = spy %>%
arrange(ds) %>%
mutate(
target_price_5d = lead(close, 5),
target_return_5d = round(lead(close, 5) / open - 1,4),
price_diff_5d = lead(close, 5) - open,
target_price_10d = lead(close, 10),
target_return_10d = round(lead(close, 10) / open - 1, 4),
price_diff_10d = lead(close, 10) - open
)
knitr::kable(head(spy, 10))| ds | open | close | target_price_5d | target_return_5d | price_diff_5d | target_price_10d | target_return_10d | price_diff_10d |
|---|---|---|---|---|---|---|---|---|
| 2015-01-02 | 206.38 | 205.43 | 204.25 | -0.0103 | -2.130005 | 201.63 | -0.0230 | -4.750000 |
| 2015-01-05 | 204.17 | 201.72 | 202.65 | -0.0074 | -1.520004 | 202.06 | -0.0103 | -2.110000 |
| 2015-01-06 | 202.09 | 199.82 | 202.08 | 0.0000 | -0.009994 | 203.08 | 0.0049 | 0.990006 |
| 2015-01-07 | 201.42 | 202.31 | 200.86 | -0.0028 | -0.559997 | 206.10 | 0.0232 | 4.680008 |
| 2015-01-08 | 204.01 | 205.90 | 199.02 | -0.0245 | -4.989991 | 204.97 | 0.0047 | 0.960006 |
| 2015-01-09 | 206.40 | 204.25 | 201.63 | -0.0231 | -4.769989 | 205.45 | -0.0046 | -0.949997 |
| 2015-01-12 | 204.41 | 202.65 | 202.06 | -0.0115 | -2.350006 | 202.74 | -0.0082 | -1.669999 |
| 2015-01-13 | 204.12 | 202.08 | 203.08 | -0.0051 | -1.039993 | 200.14 | -0.0195 | -3.979996 |
| 2015-01-14 | 199.65 | 200.86 | 206.10 | 0.0323 | 6.450012 | 201.99 | 0.0117 | 2.340011 |
| 2015-01-15 | 201.63 | 199.02 | 204.97 | 0.0166 | 3.339996 | 199.45 | -0.0108 | -2.180008 |
library(ggplot2)
library(gridExtra)
library(skimr)
spy %>%
ggplot(aes(x = price_diff_5d)) +
geom_histogram(alpha = 0.5, fill = "darkblue") +
theme_light() +
labs(
title = "SPY Diferença de Preços entre abertura e fechamento após 5 pregões.",
x = "Diferença de Preços",
y = "Frequência",
caption = paste0("Período entre ", min(spy$ds), " e ", max(spy$ds), ".")
)
spy %>%
ggplot(aes(x = price_diff_10d)) +
geom_histogram(alpha = 0.5, fill = "darkblue") +
theme_light() +
labs(
title = "SPY Diferença de Preços entre abertura e fechamento após 10 pregões.",
x = "Diferença de Preços",
y = "Frequência",
caption = paste0("Período entre ", min(spy$ds), " e ", max(spy$ds), ".")
)
| Name | spy |
| Number of rows | 1705 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| ds | 0 | 1 | 10 | 10 | 0 | 1705 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| open | 0 | 1.00 | 275.78 | 66.18 | 182.34 | 216.08 | 267.96 | 303.47 | 453.32 | ▇▇▅▂▂ |
| close | 0 | 1.00 | 275.82 | 66.19 | 182.86 | 216.16 | 267.87 | 303.33 | 453.19 | ▇▇▃▂▂ |
| target_price_5d | 5 | 1.00 | 276.03 | 66.17 | 182.86 | 216.29 | 268.58 | 303.67 | 453.19 | ▇▇▃▂▂ |
| target_return_5d | 5 | 1.00 | 0.00 | 0.02 | -0.19 | -0.01 | 0.00 | 0.01 | 0.15 | ▁▁▇▃▁ |
| price_diff_5d | 5 | 1.00 | 0.72 | 6.53 | -56.87 | -1.63 | 1.02 | 3.91 | 33.46 | ▁▁▂▇▁ |
| target_price_10d | 10 | 0.99 | 276.25 | 66.14 | 182.86 | 216.41 | 269.02 | 304.12 | 453.19 | ▇▇▃▂▂ |
| target_return_10d | 10 | 0.99 | 0.01 | 0.03 | -0.22 | -0.01 | 0.01 | 0.02 | 0.17 | ▁▁▇▅▁ |
| price_diff_10d | 10 | 0.99 | 1.41 | 9.09 | -66.12 | -1.67 | 2.13 | 5.93 | 41.45 | ▁▁▂▇▁ |
Note que temos algumas informações interessantes aqui durante esse período:
- O returno médio de 5 dias foi de 0.72 USD.
- 50% das vezes (p50), o retorno de 5 dias foi de 1.02 USD - esta informação nos mostra que a estratégia de debit spread pode funcionar, uma vez que dentro esse período, 50% das vezes o retorno é acima de 1 USD.
- Para uma estratégia mais segura, pode-se utilizar o retorno de 10 dias, que foi de 2.13 USD para debit spreads de 1 USD.
Por curiosidade, vamos observar melhor os percentis para a diferença de preços de 10 dias:
## 0% 5% 10% 15% 20% 25%
## -66.1199950 -12.9079924 -7.4200012 -4.1930113 -2.7500000 -1.6699985
## 30% 35% 40% 45% 50% 55%
## -0.7540064 -0.0010009 0.8059994 1.4989949 2.1299900 2.7100015
## 60% 65% 70% 75% 80% 85%
## 3.3840086 4.2099950 4.9660002 5.9300085 6.8299924 7.8880016
## 90% 95% 100%
## 9.7559880 13.0779971 41.4500130O dado mostra um viés positivo dessa distribuição, indicando que em 55% das vezes, há um retorno superior a 1.50 USD. Como precisamos apenas de 1.00 USD, essa estratégia me parece bem promissora.
É possível predizer melhor quando o retorno será acima de 1 USD?
Nesta etapa vou desenvolver alguns modelos para verificar se é possível realizar a predição e melhorar nossas chances de realizar uma compra de um debit spread e ter lucro após 5 ou 10 dias. Podemos utilizar alguns métodos diferentes para isso que incluem análise de regressão, árvores de decisão e algoritmos de rede neural. No geral, o processo de modelagem será composto por:
- Criação de variáveis
- Divisão do conjunto de observações em treinamento e teste.
- Ajuste de modelo nos dados de treinamento, avaliação no conjunto de teste.
- Escolha do melhor modelo.
Criação de variáveis
A criação de variáveis será feita com base em dados passados, inclusive com o fechamento do dia anterior. Mas no geral, as variáveis são construções baseadas nos preços históricos, incluindo alguns indicadores conhecimentos de média móveis, e algumas outras variáveis que tirei da minha cabeça.
O preço do fechamento anterior será utilizado para evitar informação futura da abertura do dia atual. Em seguida, todas as variáveis serão criadas baseadas nesta nova coluna.
model.data = spy %>%
mutate(
# Saídas
target_return_5d = target_return_5d,
target_return_10d = target_return_10d,
# O preço do fechamento anterior será utilizado para evitar informação
# futura da abertura do dia atual:
fecha = lag(close, 1),
# Variáveis de retorno temporal
return_1d = fecha / lag(fecha, 1) - 1,
return_2d = fecha / lag(fecha, 2) - 1,
return_3d = fecha / lag(fecha, 3) - 1,
return_5d = fecha / lag(fecha, 5) - 1,
return_10d = fecha / lag(fecha, 10) - 1,
return_20d = fecha / lag(fecha, 20) - 1,
# Variação em relação à médias moveis
return_mean_3d = fecha / rollapply(fecha, 3, mean, align = "right", fill = NA) - 1,
return_mean_5d = fecha / rollapply(fecha, 5, mean, align = "right", fill = NA) - 1,
return_mean_10d = fecha / rollapply(fecha, 10, mean, align = "right", fill = NA) - 1,
return_mean_20d = fecha / rollapply(fecha, 20, mean, align = "right", fill = NA) - 1,
# Desvio padrão temporal
sdratio_5d = rollapply(fecha, 5, sd, align = "right", fill = NA) /
rollapply(fecha, 5, mean, align = "right", fill = NA),
sdratio_10d = rollapply(fecha, 10, sd, align = "right", fill = NA) /
rollapply(fecha, 10, mean, align = "right", fill = NA),
sdratio_20d = rollapply(fecha, 20, sd, align = "right", fill = NA) /
rollapply(fecha, 20, mean, align = "right", fill = NA),
# Razão de desvios
sdratio_10d_5d = sdratio_10d / sdratio_5d - 1,
sdratio_20d_10d = sdratio_20d / sdratio_10d - 1,
# Razão das máximas
return_max_3d = fecha / rollapply(fecha, 3, max, align = "right", fill = NA) - 1,
return_max_5d = fecha / rollapply(fecha, 5, max, align = "right", fill = NA) - 1,
return_max_10d = fecha / rollapply(fecha, 10, max, align = "right", fill = NA) - 1,
return_max_20d = fecha / rollapply(fecha, 20, max, align = "right", fill = NA) - 1,
# Razão das mínimas
return_min_3d = fecha / rollapply(fecha, 3, min, align = "right", fill = NA) - 1,
return_min_5d = fecha / rollapply(fecha, 5, min, align = "right", fill = NA) - 1,
return_min_10d = fecha / rollapply(fecha, 10, min, align = "right", fill = NA) - 1,
return_min_20d = fecha / rollapply(fecha, 20, min, align = "right", fill = NA) - 1,
.keep = "none" # Incluir apenas as variáveis mencionadas
) %>%
filter(complete.cases(.)) # Filtrar as linhas vazias
skim(model.data)| Name | model.data |
| Number of rows | 1674 |
| Number of columns | 26 |
| _______________________ | |
| Column type frequency: | |
| numeric | 26 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| target_return_5d | 0 | 1 | 0.00 | 0.02 | -0.19 | -0.01 | 0.00 | 0.01 | 0.15 | ▁▁▇▃▁ |
| target_return_10d | 0 | 1 | 0.01 | 0.03 | -0.22 | -0.01 | 0.01 | 0.02 | 0.17 | ▁▁▇▅▁ |
| fecha | 0 | 1 | 275.65 | 65.04 | 182.86 | 216.64 | 268.96 | 302.84 | 453.19 | ▇▇▃▂▂ |
| return_1d | 0 | 1 | 0.00 | 0.01 | -0.11 | 0.00 | 0.00 | 0.01 | 0.09 | ▁▁▇▁▁ |
| return_2d | 0 | 1 | 0.00 | 0.01 | -0.14 | 0.00 | 0.00 | 0.01 | 0.11 | ▁▁▇▂▁ |
| return_3d | 0 | 1 | 0.00 | 0.02 | -0.13 | -0.01 | 0.00 | 0.01 | 0.17 | ▁▂▇▁▁ |
| return_5d | 0 | 1 | 0.00 | 0.02 | -0.18 | -0.01 | 0.00 | 0.01 | 0.17 | ▁▁▇▁▁ |
| return_10d | 0 | 1 | 0.01 | 0.03 | -0.23 | -0.01 | 0.01 | 0.02 | 0.19 | ▁▁▇▃▁ |
| return_20d | 0 | 1 | 0.01 | 0.04 | -0.31 | -0.01 | 0.02 | 0.03 | 0.23 | ▁▁▇▇▁ |
| return_mean_3d | 0 | 1 | 0.00 | 0.01 | -0.08 | 0.00 | 0.00 | 0.00 | 0.05 | ▁▁▂▇▁ |
| return_mean_5d | 0 | 1 | 0.00 | 0.01 | -0.10 | 0.00 | 0.00 | 0.01 | 0.09 | ▁▁▇▁▁ |
| return_mean_10d | 0 | 1 | 0.00 | 0.02 | -0.15 | 0.00 | 0.00 | 0.01 | 0.07 | ▁▁▁▇▁ |
| return_mean_20d | 0 | 1 | 0.00 | 0.02 | -0.20 | 0.00 | 0.01 | 0.02 | 0.11 | ▁▁▁▇▁ |
| sdratio_5d | 0 | 1 | 0.01 | 0.01 | 0.00 | 0.00 | 0.01 | 0.01 | 0.07 | ▇▁▁▁▁ |
| sdratio_10d | 0 | 1 | 0.01 | 0.01 | 0.00 | 0.01 | 0.01 | 0.01 | 0.09 | ▇▁▁▁▁ |
| sdratio_20d | 0 | 1 | 0.01 | 0.01 | 0.00 | 0.01 | 0.01 | 0.02 | 0.11 | ▇▁▁▁▁ |
| sdratio_10d_5d | 0 | 1 | 0.65 | 1.06 | -0.32 | -0.02 | 0.31 | 0.92 | 12.78 | ▇▁▁▁▁ |
| sdratio_20d_10d | 0 | 1 | 0.61 | 0.91 | -0.30 | 0.01 | 0.34 | 0.88 | 7.98 | ▇▁▁▁▁ |
| return_max_3d | 0 | 1 | -0.01 | 0.01 | -0.14 | -0.01 | 0.00 | 0.00 | 0.00 | ▁▁▁▁▇ |
| return_max_5d | 0 | 1 | -0.01 | 0.01 | -0.17 | -0.01 | 0.00 | 0.00 | 0.00 | ▁▁▁▁▇ |
| return_max_10d | 0 | 1 | -0.01 | 0.02 | -0.23 | -0.01 | 0.00 | 0.00 | 0.00 | ▁▁▁▁▇ |
| return_max_20d | 0 | 1 | -0.02 | 0.03 | -0.29 | -0.02 | -0.01 | 0.00 | 0.00 | ▁▁▁▁▇ |
| return_min_3d | 0 | 1 | 0.01 | 0.01 | 0.00 | 0.00 | 0.00 | 0.01 | 0.11 | ▇▁▁▁▁ |
| return_min_5d | 0 | 1 | 0.01 | 0.01 | 0.00 | 0.00 | 0.01 | 0.02 | 0.17 | ▇▁▁▁▁ |
| return_min_10d | 0 | 1 | 0.02 | 0.02 | 0.00 | 0.01 | 0.02 | 0.03 | 0.17 | ▇▁▁▁▁ |
| return_min_20d | 0 | 1 | 0.03 | 0.03 | 0.00 | 0.01 | 0.03 | 0.05 | 0.29 | ▇▁▁▁▁ |
Este conjunto de variáveis já nos permite ter um início para trabalhar. Em seguida, precisamos dividir o conjunto de dados em treinamento e teste, para que possamos desenvolver o modelo.
library(caret)
library(Hmisc)
set.seed(100) # Pra reprodução dos resultados
# Criação da variável indicadora para treinamento e test
index = createDataPartition(
y = model.data$target_return_5d,
p = 0.7,
list = F
)
# Separação dos conjuntos de dados
data.train = model.data[index,]
data.test = model.data[-index,]
# Função para plotar o gráfico de resíduos - Será utilizado em todos os mdoelos.
residuals.plot = function(observed, predicted) {
res = data.frame(
observed,
predicted,
residual = observed - predicted
)
plot1 = res %>%
ggplot(aes(x = predicted, y = residual)) +
geom_point() +
geom_smooth() +
theme_light(base_size = 20) +
labs(
title = "Gráfico de Resíduos",
x = "Valor Predito",
y = "Resíduos (Observado - Estimado)"
)
plot2 = res %>%
ggplot(aes(x = residual)) +
geom_histogram(alpha = 0.5, fill = "darkblue") +
theme_light(base_size = 20) +
labs(
title = "Distribuição de Resíduos",
x = "Resíduos",
y = "Frequenia"
)
grid.arrange(plot1, plot2, nrow = 1)
}
# Função para taxa de acerto
# Esta função vai avaliar quantos % de acerto para predição acima de 1 USD.
hit.rate = function(price, observed, predicted) {
df = data.frame(
price,
observed,
predicted,
price_observed = price * observed,
price_predicted = price * predicted
)
df %>%
mutate(
price_observed_1usd = ifelse(price_observed > 1, "yes", "no"),
price_predicted_1usd = ifelse(price_predicted > 1, "yes", "no")
) %>%
filter(price_predicted_1usd == "yes") %>%
count(price_predicted_1usd, price_observed_1usd) %>%
mutate(
pct = n/sum(n)
)
}
# Simplificando cálculos para a matrix de confusão:
my.confusion.matrix = function(price, observed, predicted) {
df = data.frame(
price,
observed,
predicted,
price_observed = price * observed,
price_predicted = price * predicted
)
df.out = df %>%
mutate(
price_observed_1usd = ifelse(price_observed > 1, "yes", "no"),
price_predicted_1usd = ifelse(price_predicted > 1, "yes", "no")
)
caret::confusionMatrix(as.factor(df.out$price_predicted_1usd), as.factor(df.out$price_observed_1usd), positive = "yes")
}
# Função de perfil, para testar diferentes limites de predição:
profiling = function(price, observed, predicted) {
df = data.frame(
price,
observed,
predicted,
price_observed = price * observed,
price_predicted = price * predicted
)
df %>%
mutate(
price_observed_1usd = ifelse(price_observed > 1, "yes", "no"),
pred_intervals = cut2(price_predicted, g = 10),
) %>%
count(pred_intervals, price_observed_1usd) %>%
pivot_wider(
names_from = price_observed_1usd,
values_from = n,
names_prefix = "observed_"
) %>%
mutate(
hit_rate = observed_yes / (observed_yes + observed_no)
)
}Modelos para Returno de 5 dias
Vamos iniciar agora com os modelos para predição de retorno de 5 dias. Vamos fazer o de 10 dias posteriormente.
Regressão Linear Simples
Aqui temos o modelo de regressão linear simples como ponto de partida, utilizando procedimento stepwise para modelo parcimonioso.
mod.reg.null = lm(target_return_5d ~ 1, data = data.train[,-c(2,3)])
mod.reg.full = lm(target_return_5d ~ ., data = data.train[,-c(2,3)])
mod.reg.step = step(mod.reg.null, direction = "forward", scope = formula(mod.reg.full), trace = 0)
summary(mod.reg.step)##
## Call:
## lm(formula = target_return_5d ~ return_max_20d + return_max_10d +
## return_3d + return_2d + return_mean_20d + sdratio_20d + return_min_20d +
## return_min_3d + return_5d + return_min_5d + sdratio_10d +
## sdratio_10d_5d + return_10d, data = data.train[, -c(2, 3)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.151812 -0.008591 0.002218 0.012085 0.143570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0024818 0.0013831 1.794 0.073017 .
## return_max_20d -0.8344263 0.1423507 -5.862 5.96e-09 ***
## return_max_10d 0.7519354 0.1550565 4.849 1.41e-06 ***
## return_3d -0.2572269 0.0858848 -2.995 0.002802 **
## return_2d 0.2872759 0.0997717 2.879 0.004058 **
## return_mean_20d 0.3744827 0.1820480 2.057 0.039904 *
## sdratio_20d -1.1779313 0.3276572 -3.595 0.000338 ***
## return_min_20d 0.2245139 0.1098132 2.045 0.041129 *
## return_min_3d -0.4415899 0.1704695 -2.590 0.009706 **
## return_5d -0.2168892 0.0715576 -3.031 0.002492 **
## return_min_5d 0.0401672 0.1544889 0.260 0.794909
## sdratio_10d 0.7021622 0.2724562 2.577 0.010084 *
## sdratio_10d_5d -0.0016528 0.0008899 -1.857 0.063520 .
## return_10d -0.0967114 0.0649589 -1.489 0.136810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02278 on 1160 degrees of freedom
## Multiple R-squared: 0.07224, Adjusted R-squared: 0.06185
## F-statistic: 6.948 on 13 and 1160 DF, p-value: 4.394e-13
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 306 349
## yes 217 302
##
## Accuracy : 0.5179
## 95% CI : (0.4889, 0.5468)
## No Information Rate : 0.5545
## P-Value [Acc > NIR] : 0.9946
##
## Kappa : 0.0478
##
## Mcnemar's Test P-Value : 3.664e-08
##
## Sensitivity : 0.4639
## Specificity : 0.5851
## Pos Pred Value : 0.5819
## Neg Pred Value : 0.4672
## Prevalence : 0.5545
## Detection Rate : 0.2572
## Detection Prevalence : 0.4421
## Balanced Accuracy : 0.5245
##
## 'Positive' Class : yes
## | price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 217 | 0.4181118 |
| yes | yes | 302 | 0.5818882 |
| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-19.332,-0.670) | 48 | 70 | 0.5932203 |
| [ -0.670,-0.135) | 55 | 62 | 0.5299145 |
| [ -0.135, 0.272) | 63 | 55 | 0.4661017 |
| [ 0.272, 0.570) | 49 | 68 | 0.5811966 |
| [ 0.570, 0.839) | 59 | 58 | 0.4957265 |
| [ 0.839, 1.115) | 50 | 68 | 0.5762712 |
| [ 1.115, 1.456) | 51 | 66 | 0.5641026 |
| [ 1.456, 1.932) | 48 | 70 | 0.5932203 |
| [ 1.932, 3.019) | 53 | 64 | 0.5470085 |
| [ 3.019,35.544] | 47 | 70 | 0.5982906 |
Aplicação na base de teste a avaliação:
# Final results (testing)
residuals.plot(data.test$target_return_5d, predict(mod.reg.step, newdata = data.test))
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 129 143
## yes 98 130
##
## Accuracy : 0.518
## 95% CI : (0.4732, 0.5626)
## No Information Rate : 0.546
## P-Value [Acc > NIR] : 0.903480
##
## Kappa : 0.0437
##
## Mcnemar's Test P-Value : 0.004593
##
## Sensitivity : 0.4762
## Specificity : 0.5683
## Pos Pred Value : 0.5702
## Neg Pred Value : 0.4743
## Prevalence : 0.5460
## Detection Rate : 0.2600
## Detection Prevalence : 0.4560
## Balanced Accuracy : 0.5222
##
## 'Positive' Class : yes
## knitr::kable(hit.rate(438.66, data.test$target_return_5d, predict(mod.reg.step, newdata = data.test)))| price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 98 | 0.4298246 |
| yes | yes | 130 | 0.5701754 |
knitr::kable(profiling(438.66, data.test$target_return_5d, predict(mod.reg.step, newdata = data.test)))| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-9.878,-0.967) | 25 | 25 | 0.50 |
| [-0.967,-0.201) | 23 | 27 | 0.54 |
| [-0.201, 0.210) | 22 | 28 | 0.56 |
| [ 0.210, 0.501) | 23 | 27 | 0.54 |
| [ 0.501, 0.861) | 25 | 25 | 0.50 |
| [ 0.861, 1.200) | 26 | 24 | 0.48 |
| [ 1.200, 1.598) | 23 | 27 | 0.54 |
| [ 1.598, 2.229) | 26 | 24 | 0.48 |
| [ 2.229, 3.939) | 16 | 34 | 0.68 |
| [ 3.939,19.968] | 18 | 32 | 0.64 |
A conclusão que podemos tirar é que o modelo, meso com um baixo R2, tem algum poder preditivo, uma vez que o hit-rate foi acima de 50%. Houve estabilidade também, uma vez que os hit-rates são semelhantes entre treinamento e test:
- Hit-rate treinamento: 58.2%
- Hit-rate test: 57%
Esse modelo não considera interações entre variáveis, portanto, vamos verificar se é possível ter um modelo melhor com árvores de decisão:
Árvore de Decisão
library(rpart)
library(rpart.plot)
mod.tree = rpart(target_return_5d ~ ., data = data.train[,-c(2,3)])
rpart.plot(mod.tree)
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 497 547
## yes 26 104
##
## Accuracy : 0.5119
## 95% CI : (0.4829, 0.5409)
## No Information Rate : 0.5545
## P-Value [Acc > NIR] : 0.9985
##
## Kappa : 0.1002
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.15975
## Specificity : 0.95029
## Pos Pred Value : 0.80000
## Neg Pred Value : 0.47605
## Prevalence : 0.55451
## Detection Rate : 0.08859
## Detection Prevalence : 0.11073
## Balanced Accuracy : 0.55502
##
## 'Positive' Class : yes
## | price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 26 | 0.2 |
| yes | yes | 104 | 0.8 |
| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-23.75, 5.60) | 497 | 547 | 0.5239464 |
| 5.60 | 9 | 28 | 0.7567568 |
| [ 7.22,27.31] | 17 | 76 | 0.8172043 |
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 203 229
## yes 24 44
##
## Accuracy : 0.494
## 95% CI : (0.4493, 0.5387)
## No Information Rate : 0.546
## P-Value [Acc > NIR] : 0.9912
##
## Kappa : 0.0515
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.1612
## Specificity : 0.8943
## Pos Pred Value : 0.6471
## Neg Pred Value : 0.4699
## Prevalence : 0.5460
## Detection Rate : 0.0880
## Detection Prevalence : 0.1360
## Balanced Accuracy : 0.5277
##
## 'Positive' Class : yes
## | price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 24 | 0.3529412 |
| yes | yes | 44 | 0.6470588 |
| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-23.75, 5.60) | 203 | 229 | 0.5300926 |
| 5.60 | 6 | 13 | 0.6842105 |
| [ 7.22,27.31] | 18 | 31 | 0.6326531 |
Este modelo de árvore de decisão apresenta um resultado já mais interessante: * Hit rate treinamento: 80% * Hit rate teste: 64%
Apesar de um pouco instável devido a diferença nessas estatísticas, me parece que este modelo ainda sim tente a acertar um pouco mais do que aleatoriamente (64% vs 50.3%), com potencial para identificar entrada de trades.
Em seguida vamos testar um modelo de redes neurais, que é mais complexo e considera ainda mais interações.
Rede Neural
library(neuralnet)
mod.neural = neuralnet(target_return_5d ~ ., data = data.train[,-c(2,3)], hidden = 10)
plot(mod.neural)# Final Results (Training):
residuals.plot(data.train$target_return_5d, predict(mod.neural, data.train))
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 304 285
## yes 219 366
##
## Accuracy : 0.5707
## 95% CI : (0.5418, 0.5992)
## No Information Rate : 0.5545
## P-Value [Acc > NIR] : 0.138619
##
## Kappa : 0.1417
##
## Mcnemar's Test P-Value : 0.003788
##
## Sensitivity : 0.5622
## Specificity : 0.5813
## Pos Pred Value : 0.6256
## Neg Pred Value : 0.5161
## Prevalence : 0.5545
## Detection Rate : 0.3118
## Detection Prevalence : 0.4983
## Balanced Accuracy : 0.5717
##
## 'Positive' Class : yes
## | price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 219 | 0.374359 |
| yes | yes | 366 | 0.625641 |
| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-47.476,-2.988) | 75 | 43 | 0.3644068 |
| [ -2.988,-1.436) | 57 | 60 | 0.5128205 |
| [ -1.436,-0.474) | 58 | 60 | 0.5084746 |
| [ -0.474, 0.339) | 55 | 62 | 0.5299145 |
| [ 0.339, 0.989) | 59 | 58 | 0.4957265 |
| [ 0.989, 1.720) | 45 | 73 | 0.6186441 |
| [ 1.720, 2.516) | 53 | 64 | 0.5470085 |
| [ 2.516, 3.664) | 50 | 68 | 0.5762712 |
| [ 3.664, 5.571) | 44 | 73 | 0.6239316 |
| [ 5.571,40.025] | 27 | 90 | 0.7692308 |
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 134 119
## yes 93 154
##
## Accuracy : 0.576
## 95% CI : (0.5313, 0.6198)
## No Information Rate : 0.546
## P-Value [Acc > NIR] : 0.09619
##
## Kappa : 0.1529
##
## Mcnemar's Test P-Value : 0.08598
##
## Sensitivity : 0.5641
## Specificity : 0.5903
## Pos Pred Value : 0.6235
## Neg Pred Value : 0.5296
## Prevalence : 0.5460
## Detection Rate : 0.3080
## Detection Prevalence : 0.4940
## Balanced Accuracy : 0.5772
##
## 'Positive' Class : yes
## knitr::kable(hit.rate(438.66, data.test$target_return_5d, predict(mod.neural, newdata = data.test)))| price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 93 | 0.3765182 |
| yes | yes | 154 | 0.6234818 |
knitr::kable(profiling(438.66, data.test$target_return_5d, predict(mod.neural, newdata = data.test)))| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-29.520,-2.934) | 28 | 22 | 0.44 |
| [ -2.934,-1.595) | 25 | 25 | 0.50 |
| [ -1.595,-0.497) | 25 | 25 | 0.50 |
| [ -0.497, 0.182) | 30 | 20 | 0.40 |
| [ 0.182, 0.979) | 24 | 26 | 0.52 |
| [ 0.979, 1.855) | 18 | 32 | 0.64 |
| [ 1.855, 2.756) | 19 | 31 | 0.62 |
| [ 2.756, 4.185) | 26 | 24 | 0.48 |
| [ 4.185, 6.269) | 17 | 33 | 0.66 |
| [ 6.269,69.948] | 15 | 35 | 0.70 |
De forma interessante, a rede neural estimou um modelo estável, com um bom hit rate também:
- Hit-rate treinamento: 62.5%
- Hit-rate test: 62.3%
Desta forma, também é um modelo candidato a ser utilizado.
Juntar os modelos pra avaliação conjunta
Vamos observar as tabelas de predições e ver qual é o poder preditivo dos modelos em conjunto.
data.train %>%
mutate(
pred.lm = ifelse(438.66 * predict(mod.reg.step) > 1, "yes", "no"),
pred.rpart = ifelse(438.66 * predict(mod.tree) > 1, "yes", "no"),
pred.neural = ifelse(438.66 * predict(mod.neural, newdata = .) > 1, "yes","no"),
obs = ifelse(438.66 * target_return_5d > 1, "yes","no")
) %>%
count(pred.lm, pred.rpart, pred.neural, obs) %>%
pivot_wider(names_from = obs, values_from = n) %>%
mutate(
hit_rate = yes / (yes + no)
) %>%
knitr::kable()| pred.lm | pred.rpart | pred.neural | no | yes | hit_rate |
|---|---|---|---|---|---|
| no | no | no | 204 | 157 | 0.4349030 |
| no | no | yes | 90 | 152 | 0.6280992 |
| no | yes | no | 4 | 15 | 0.7894737 |
| no | yes | yes | 8 | 25 | 0.7575758 |
| yes | no | no | 93 | 105 | 0.5303030 |
| yes | no | yes | 110 | 133 | 0.5473251 |
| yes | yes | no | 3 | 8 | 0.7272727 |
| yes | yes | yes | 11 | 56 | 0.8358209 |
data.test %>%
mutate(
pred.lm = ifelse(438.66 * predict(mod.reg.step, newdata = data.test) > 1, "yes", "no"),
pred.rpart = ifelse(438.66 * predict(mod.tree, newdata = data.test) > 1, "yes", "no"),
pred.neural = ifelse(438.66 * predict(mod.neural, newdata = data.test) > 1, "yes","no"),
obs = ifelse(438.66 * target_return_5d > 1, "yes","no")
) %>%
count(pred.lm, pred.rpart, pred.neural, obs) %>%
pivot_wider(names_from = obs, values_from = n) %>%
mutate(
hit_rate = yes / (yes + no)
) %>%
knitr::kable()| pred.lm | pred.rpart | pred.neural | no | yes | hit_rate |
|---|---|---|---|---|---|
| no | no | no | 79 | 68 | 0.4625850 |
| no | no | yes | 44 | 63 | 0.5887850 |
| no | yes | no | 4 | 7 | 0.6363636 |
| no | yes | yes | 2 | 5 | 0.7142857 |
| yes | no | no | 44 | 41 | 0.4823529 |
| yes | no | yes | 36 | 57 | 0.6129032 |
| yes | yes | no | 7 | 3 | 0.3000000 |
| yes | yes | yes | 11 | 29 | 0.7250000 |
E em seguida, vamos ver se um modelo de regressão consegue ponderar estas predições:
blend.train = data.train %>%
mutate(
pred.lm = ifelse(438.66 * predict(mod.reg.step) > 1, "yes", "no"),
pred.rpart = ifelse(438.66 * predict(mod.tree) > 1, "yes", "no"),
pred.neural = ifelse(438.66 * predict(mod.neural, newdata = .) > 1, "yes","no"),
obs = ifelse(438.66 * target_return_5d > 1, "yes","no")
)
blend.lm = lm(target_return_5d ~ pred.lm * pred.rpart * pred.neural, data = blend.train)
summary(blend.lm)##
## Call:
## lm(formula = target_return_5d ~ pred.lm * pred.rpart * pred.neural,
## data = blend.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.181817 -0.008542 0.003095 0.012530 0.081582
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004683 0.001168 -4.010 6.44e-05
## pred.lmyes 0.004275 0.001962 2.179 0.02952
## pred.rpartyes 0.013930 0.005222 2.668 0.00774
## pred.neuralyes 0.010650 0.001843 5.778 9.69e-09
## pred.lmyes:pred.rpartyes 0.003796 0.008631 0.440 0.66020
## pred.lmyes:pred.neuralyes -0.008061 0.002812 -2.866 0.00423
## pred.rpartyes:pred.neuralyes -0.001979 0.006650 -0.298 0.76606
## pred.lmyes:pred.rpartyes:pred.neuralyes 0.009193 0.010041 0.916 0.36008
##
## (Intercept) ***
## pred.lmyes *
## pred.rpartyes **
## pred.neuralyes ***
## pred.lmyes:pred.rpartyes
## pred.lmyes:pred.neuralyes **
## pred.rpartyes:pred.neuralyes
## pred.lmyes:pred.rpartyes:pred.neuralyes
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02219 on 1166 degrees of freedom
## Multiple R-squared: 0.1155, Adjusted R-squared: 0.1102
## F-statistic: 21.76 on 7 and 1166 DF, p-value: < 2.2e-16
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 407 395
## yes 116 256
##
## Accuracy : 0.5647
## 95% CI : (0.5358, 0.5933)
## No Information Rate : 0.5545
## P-Value [Acc > NIR] : 0.2499
##
## Kappa : 0.1629
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.3932
## Specificity : 0.7782
## Pos Pred Value : 0.6882
## Neg Pred Value : 0.5075
## Prevalence : 0.5545
## Detection Rate : 0.2181
## Detection Prevalence : 0.3169
## Balanced Accuracy : 0.5857
##
## 'Positive' Class : yes
## | price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 116 | 0.311828 |
| yes | yes | 256 | 0.688172 |
| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| -2.054 | 204 | 157 | 0.4349030 |
| -0.179 | 93 | 105 | 0.5303030 |
| 0.957 | 110 | 133 | 0.5473251 |
| 2.617 | 90 | 152 | 0.6280992 |
| 4.056 | 4 | 15 | 0.7894737 |
| [ 7.597,11.897] | 22 | 89 | 0.8018018 |
O blend de modelos aparenta ter poder preditivo superior aos modelos individuais. Vamos checar agora se isso ocorre na base de teste também.
blend.test = data.test %>%
mutate(
pred.lm = ifelse(438.66 * predict(mod.reg.step, data.test) > 1, "yes", "no"),
pred.rpart = ifelse(438.66 * predict(mod.tree, data.test) > 1, "yes", "no"),
pred.neural = ifelse(438.66 * predict(mod.neural, newdata = data.test) > 1, "yes","no"),
obs = ifelse(438.66 * target_return_5d > 1, "yes","no")
)
residuals.plot(blend.test$target_return_5d, predict(blend.lm, newdata = blend.test))
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 159 166
## yes 68 107
##
## Accuracy : 0.532
## 95% CI : (0.4872, 0.5764)
## No Information Rate : 0.546
## P-Value [Acc > NIR] : 0.75
##
## Kappa : 0.0891
##
## Mcnemar's Test P-Value : 2.281e-10
##
## Sensitivity : 0.3919
## Specificity : 0.7004
## Pos Pred Value : 0.6114
## Neg Pred Value : 0.4892
## Prevalence : 0.5460
## Detection Rate : 0.2140
## Detection Prevalence : 0.3500
## Balanced Accuracy : 0.5462
##
## 'Positive' Class : yes
## knitr::kable(hit.rate(438.66, blend.test$target_return_5d, predict(blend.lm, newdata = blend.test)))| price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 68 | 0.3885714 |
| yes | yes | 107 | 0.6114286 |
knitr::kable(profiling(438.66, blend.test$target_return_5d, predict(blend.lm, newdata = blend.test)))| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| -2.054 | 79 | 68 | 0.4625850 |
| -0.179 | 44 | 41 | 0.4823529 |
| 0.957 | 36 | 57 | 0.6129032 |
| 2.617 | 44 | 63 | 0.5887850 |
| [ 4.056, 7.860) | 11 | 10 | 0.4761905 |
| [ 7.860,11.897] | 13 | 34 | 0.7234043 |
Treinando modelos usando caret
Em seguida, vamos treinar todos os modelos utilizando o pacote caret, para que a comparação entre modelos seja feita de forma eficiente e que possamos utilizar cross-validation de forma simples e rápida. Também, vamos incluir outros modelos como random forest e gradient boosting.
# Parâmetros para cross-validação
trainParams = trainControl(method = "repeatedcv", number = 5, repeats = 3)
set.seed(100)
# Modelo Linear
caret.lm = train(target_return_5d ~ .,
data = data.train[,-c(2,3)],
method = "lm",
trControl = trainParams)
# Árvore de Decisão
caret.dtree = train(target_return_5d ~ .,
data = data.train[,-c(2,3)],
method = "rpart",
trControl = trainParams)
# Rede Neural
caret.neuralnet = train(target_return_5d ~ .,
data = data.train[,-c(2,3)],
method = "neuralnet",
trControl = trainParams)
# Random Forest
caret.rf = train(target_return_5d ~ .,
data = data.train[,-c(2,3)],
method = "rf",
trControl = trainParams)
# Gradient Boosting Machine (gbm)
caret.gbm = train(target_return_5d ~ .,
data = data.train[,-c(2,3)],
method = "gbm",
trControl = trainParams,
verbose = F)
# Comparação de Modelos
results = resamples(
list(
linear_model = caret.lm,
decision_tree = caret.dtree,
neural_net = caret.neuralnet,
random_forest = caret.rf,
gradient_boosting = caret.gbm
)
)
dotplot(results, metric = "RMSE", main = "Comparação de Estimativas de Erro")
Podemos verificar que o modelo de GBM teve uma estimativa de erro um pouco menor, mas ainda sim um intervalo de confiança grande no geral, com grande overlap com modelos mais simples como o modelo de regressão linear e a árvore de decisão. Por fim, vamos utilizar nossas funções anteriores para ver como o resultado na base teste se mantém:
## Confusion Matrix and Statistics
##
## Reference
## Prediction no yes
## no 153 157
## yes 74 116
##
## Accuracy : 0.538
## 95% CI : (0.4932, 0.5824)
## No Information Rate : 0.546
## P-Value [Acc > NIR] : 0.6574
##
## Kappa : 0.096
##
## Mcnemar's Test P-Value : 6.845e-08
##
## Sensitivity : 0.4249
## Specificity : 0.6740
## Pos Pred Value : 0.6105
## Neg Pred Value : 0.4935
## Prevalence : 0.5460
## Detection Rate : 0.2320
## Detection Prevalence : 0.3800
## Balanced Accuracy : 0.5495
##
## 'Positive' Class : yes
## | price_predicted_1usd | price_observed_1usd | n | pct |
|---|---|---|---|
| yes | no | 74 | 0.3894737 |
| yes | yes | 116 | 0.6105263 |
knitr::kable(profiling(438.66, data.test$target_return_5d, predict(caret.gbm, newdata = data.test)))| pred_intervals | observed_no | observed_yes | hit_rate |
|---|---|---|---|
| [-9.361, 0.690) | 54 | 46 | 0.4600000 |
| [ 0.690, 0.983) | 99 | 110 | 0.5263158 |
| [ 0.983, 1.333) | 15 | 33 | 0.6875000 |
| [ 1.333, 1.499) | 19 | 24 | 0.5581395 |
| [ 1.499, 2.735) | 20 | 31 | 0.6078431 |
| [ 2.735,21.586] | 20 | 29 | 0.5918367 |
Considerações Finais
Este documento apresentou o procedimento geral de desenvolvimento de um modelo de machine learning para predição futura de preço de ações. No geral, os resultados não são muito positivos, o que se é esperado dado que esta tarefa é muito complexa e os mercados são bem otimizados. No entanto, ainda algumas coisas podem ser feitas ainda para explorar mais esta área:
- Desenvolver mais variáveis preditoras, baseados nos preços históricos, relação com outros índices e dados de mercado, etc.
- Remover eventos discrepantes do conjunto de dados - isso pode remover um pouco de viés dos modelos onde outliers podem alterar nossa capacidade preditiva, mas ao mesmo tempo também não nos permite prever momentos discrepantes.
- Testar mais modelos, outros algoritmos de machine learning, etc.