Previsão de comportamento de ações com Machine Learning em Python

19 min readJul 18, 2022

Esse artigo é fruto do projeto de finalização do curso de Data Science & Machine Learning da TERA. Nosso grupo optou pelo tema “Mercado Financeiro” e aqui vamos mostrar nossos principais aprendizados e alguns dos resultados que chegamos.

Primeiramente, quero deixar bastante claro que esse artigo não vai te mostrar como ficar rico usando ML no mercado de ações. Porém, garanto que essa leitura vai trazer bastante conhecimento, e uns insights legais sobre Séries Temporais e como utiliza-las para fazer predições.

Contexto

Segundo dados da B3, o número de contas de pessoa física investidoras cresceu mais de 300% entre os anos de 2018 e 2021. Esse aumento é, além de espantoso, animador. Isso mostra que a população brasileira está amadurecendo seus investimentos, saindo de opções mais conservadoras e tradicionais, como renda fixa e poupança.

No entanto, essa evolução traz uma ressalva: será que os novos investidores sabem o que estão fazendo? Um estudo realizado por Bruno Giovannetti e Fernando Chague, da Fundação Getulio Vargas (FGV EESP), mostrou que 97% das pessoas que especulam na bolsa perdem dinheiro e quem ganha leva menos de R$ 300 por dia.

Dessa forma, é importante que sejam desenvolvidas ferramentas capazes de analisar os valores de ações e que tomem decisões rápidas e com o menor índice de perdas possível, fazendo com que os ganhos sejam maximizados, e reduzindo as possibilidades de perdas no mercado.

Tendo esse cenário em mente, para a compreensão e a previsão futura da evolução do mercado, a utilização e construção de modelos e softwares econométricos têm demandado maior atenção e investimentos.

Problema de Negócio

Sabendo-se que o mercado financeiro de ações pode se tornar altamente rentável para o investidor, porém ao mesmo tempo de alto risco, propõe-se neste projeto o desenvolvimento de um sistema inteligente, que maximize os ganhos financeiros no mercado de ações. Para o desenvolvimento preliminar deste trabalho, serão escolhidas algumas ações da bolsa de valores americana. Após a criação de um modelo robusto para algumas ações selecionadas, vamos ver se ele é replicável a outros ativos. Depois dessas etapas, podemos tentar melhorar as performances dos modelos e atacar outras frentes, como mercado de criptomoedas, por exemplo.

Questões a serem respondidas:

Conseguimos prever os próximos valores das ações, baseado no histórico de valores passados?
Qual o tempo máximo que conseguimos prever um preço com uma confiabilidade estatisticamente relevante?
Com uma estratégia automatizada, quanto dinheiro é possível ser ganho?

Etapas do Projeto

Apesar de o título deste artigo dar a entender que o foco do trabalho é a área de Machine Learning, vamos trabalhar todas as áreas do projeto, desde a obtenção dos dados, até o deploy do mesmo. As etapas serão as seguintes:

Extração e tratamento: nessa etapa utilizaremos a API do Yahoo Finance e bibliotecas como Pandas e Numpy, em Python
EDA (Exploratory Data Analysis): hipóteses são criadas, validadas ou rejeitadas, além da insights gerados através do entendimento dos dados
Modelos de ML: alguns modelos serão treinados e avaliados segundo as métricas adequadas
Teste do modelo: o melhor modelo será testado em dataset novo
Otimização de tempo e de processamento: tamanho do dataset será otimizado em relação ao processamento e tempo de execução
Otimização de hiperparâmetros: o melhor modelo será tunado em relação aos seus hiperparâmetros
Deploy: o app será colocado em produção, na cloud do Streamlit

Extração e tratamento dos dados

Como a ideia do projeto é utilizar dados de ações da bolsa de valores para recomendação de compra e venda, em um mesmo dia, precisamos de dados com variações horárias. Para isso, podemos utilizar de algumas API’s disponíveis no mercado. A escolhida foi a do Yahoo Finance, pela eficiência e facilidade de uso.

API Yahoo Finance

A API do Yahoo Finance é uma variedade de bibliotecas gratuitas com métodos para obter dados históricos e em tempo real para uma variedade de mercados, e produtos financeiros. Isso inclui dados de mercado, sobre cripto moedas, moedas regulares, ações e títulos, dados fundamentais e de opções e análises e notícias de mercado.

Além de livre, essa API é muito fácil de configurar e de utilizar, pois com apenas algumas linhas de código, conseguimos baixar os dados que precisamos diretamente ao nosso código contruído em Python.

Utilizando o Yahoo Finance para nosso caso

Como mencionado anteriormente, vamos começar nosso trabalho avaliando os dados da ação da Apple (AAPL). Para isso, vamos precisar importar a API do Yahoo Finance e baixar os dados que precisamos. Vamos baixar os dados horários de 01.01.2021 a 22.04.2022.

# importando a biblioteca da api da Yahoo Finance
import yfinance as yf# dados horários aapl
df_AAPL_hora = yf.download(tickers = 'AAPL',
                           start = '2021-01-01',
                           end = '2022-04-19',
                           interval = '1h',
                           ajusted = True)

Segundo a fonte dos dados, as features apresentadas possuem os seguintes significados:

Open: primeiro preço da unidade de tempo avaliada (hora, dia…)
High: maior preço durante a unidade de tempo avaliada (hora, dia…)
Low: menor preço durante a unidade de tempo avaliada (hora, dia…)
Close: preço de fechamento segundo a unidade de tempo avaliada (hora, dia…)
Adj Close: preço ajustado depois de splits distribuição de dividendos. Os dados são ajustados segundo os padrões estabelecidos por Center for Research in Security Prices (CRSP).

Avaliação inicial

Agora que temos os dados da AAPL, de hora em hora, vamos fazer as avaliações iniciais e verificar se há dados faltantes.

# verificando informações gerais de AAPL
print('Dataset AAPL HORA\n')
print(df_aapl_hora.info())Dataset AAPL HORA

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 2273 entries, 2021-01-04 14:30:00+00:00 to 2022-04-21 20:00:00+00:00
Data columns (total 6 columns):
 #   Column     Non-Null Count  Dtype  
---  ------     --------------  -----  
 0   Open       2273 non-null   float64
 1   High       2273 non-null   float64
 2   Low        2273 non-null   float64
 3   Close      2273 non-null   float64
 4   Adj Close  2273 non-null   float64
 5   Volume     2273 non-null   int64  
dtypes: float64(5), int64(1)
memory usage: 124.3 KB
None

Como podemos observar, não há dados faltantes e os tipos são todos numéricos. Agora, vamos avaliar as primeiras linhas do dataset.

# verificando as primeiras linhas dos datasets
df_aapl_hora.head()

Como os dados extraídos não possuem dados faltantes e já estãp aparentemente estruturados, podemos passar para a etapa de análise dos mesmos.

EDA

Agora já com os dados em mãos, vamos partir para as análises.

Para evitar dados errôneos devido a splits e distribuição de dividendos, vamos utilizar a feature Adj Close, que já considera essas questões.

O primeiro passo é visualizar a variação do valor da ação com o passar do tempo e a distribuição dos valores no período.

# plotando a variação de preço do ativo por hora
ax = df_aapl_hora['Adj Close'].plot(figsize=(12,5))
plt.title('Preço horário da AAPL entre 2021 e 2022')
plt.xlabel('')
plt.ylabel('US$')
plt.tick_params(rotation = 0)# plotando as distribuições dos preços, por DIA E HORA
plt.figure(figsize=(8,5))
ax = sns.distplot(df_aapl_hora['Adj Close'],
                  fit=norm,
                  kde=True)
                  ax.set_xlabel('Adj Close HORA')
ax.set_title('Distribuição dos PREÇOS de AAPL HORA')

Observa-se que os valores da ação da AAPL sofreram variações horárias e diárias, como esperado, e com tendência de alta no período avaliado.

Agora, vamos realizar as outras avaliações, específicas para alguns modelos que vamos utilizar. A começar pelos parâmetros do modelo ARIMA.

EDA para ARIMA

Os principais pontos para utilizar o modelo ARIMA são:
* Sazonalidade
* Série estacionária — definição do parâmetro d
* Média Móvel
* Outliers
* Função de Autocorrelação Parcial (p)
* Função de Autocorrelação (q)

Sazonalidade
Começando pela sazonalidade, vamos analisar se ela existe para esse ativo. Para isso, precisamos agrupar a série em valores em certa periodicidade. Assim, vamos definir em 40, correspondente a aproximadamente 1 semana.

# fazendo a decomposição mensal
decomposicao_aapl = seasonal_decompose(df_aapl_hora['Adj Close'],
model='additive', freq=40)
fig = decomposicao_aapl.plot()
fig.set_size_inches((12, 6))# Tight layout to realign things
fig.tight_layout()
plt.show()

Nota-se, pelo gráfico acima, que há componente de sazonalidade na nossa série temporal.
O resíduo fica maior no final do período analisado

Série Estacionária

Agora, vamos analisar se a série é estacionária, ou não. Para que seja, ela precisa ter a média e a variância constantes no tempo, além de ter correlação entre 2 períodos constantes.

Vamos utilizar o teste de Dickey Fuller para verificar se nossa série é estacionária ou não.
H0 -> Série não estacionária (p-value > 0.05)
H1 -> Série estacionária (p-value <= 0.05)

# rodando o método de Dickey Fuller
adfuller_aapl_hora = adfuller(df_aapl_hora['Adj Close'])
print('ADF Statistics AAPL HORA:', adfuller_aapl_hora[0])
print('p-value AAPL HORA:', adfuller_aapl_hora[1])ADF Statistics AAPL HORA: -1.2372893665712357 p-value AAPL HORA: 0.6573526804673867

O resultado estatístico mostrou que a série não é estacionária. (p-value > 0.05). Então, vamos utilizar a técnica de diferenciação para tornar a série estacionária, para podermos aplicar o modelo ARIMA

# Calculo do retorno percentual (1º diferenciação)
df_arima = df_aapl_hora.copy()
df_arima['return_percent'] = df_arima['Adj Close'].pct_change()
df_arima.dropna(inplace=True)# verificando as distribuições dos retornos por HORA
fig, axs = plt.subplots(1,2,figsize=(15,5))
axs[0].plot(df_arima['return_percent'])
axs[0].set_title('1º Diferenciação')
ax2 = sns.distplot(df_arima['return_percent'],
                            fit=norm,
                            kde=True,
                            ax=axs[1])
ax2.set_title('Distribuição dos retornos (1º Diferenciação)')

Aplicando novamente o teste de Dickey Fuller, verificamos que com apenas uma diferenciação (retorno percentual), a série se tornou estacionária. Assim, podemos definir o valor 1 para o parâmetro d do ARIMA.

Identificação de outliers

Para séries temporais, definimos os outliers em função das médias móveis e dos desvios padrão móveis também. Nesse caso, definimos os outliers aqueles que estão fora do intervalo de 3 desvios padrão da média móvel.

O resultado dessa avaliação está representado na imagem abaixo.

O tratamento dos outliers foi realizado com o método de winsorize da lib stats. O resultado é o seguinte:

Como podemos observar, após a utilização do método de winsorize, a variação dos valores ficou menor.

PACF — Função de Autocorrelação Parcial (p)

Definindo o parâmetro p do ARIMA.

# plotando PACF para df HORAS
fig, ax = plt.subplots(figsize=(12, 5))
plot_pacf(df_aapl_hora['return_percent'].dropna(), lags = 30, ax=ax)
ax.set_title('Autocorrelação Parcial PACF - AAPL Horas')
plt.show()

Como podemos observar, as lags 16, 21, 26 e 27 podem ser utilizados como valor de p.

ACF — Função de Autocorrelação (q)

Definindo o parâmetro q do ARIMA.

# plotando ACF para df HORAS
fig, ax = plt.subplots(figsize=(12, 5))
plot_acf(df_aapl_hora['return_percent'].dropna(), lags = 30, ax=ax)
ax.set_title('Autocorrelação ACF - AAPL Horas')
plt.show()

Como podemos observar, as lags 16, 21, 26 e 27 podem ser utilizados como valor de q.

EDA para outros modelos

Outras variáveis deverão ser criadas e avaliadas, de forma a complementarem os modelos de ML. Dessa forma, foram pesquisados alguns dos principais indicadores utilizados por especialistas em investimentos.

Indicador RSI

O indicador RSI compara os ganhos recentes com as perdas recentes e possui variação de 0 a 100. O calculo de RSI é:

Para realizar o cálculo, utilizamos a lib chamada talib e definimos 20 períodos para o cálculo.

RSI_PERIOD = 20 # definindo o período considerado para cálculo de RSI
df_var_exog['rsi'] = talib.RSI(df_var_exog['Adj Close'], RSI_PERIOD)  # criando a feature RSI# RSI < 30 -> 0 (comprar)
# RSI > 70 -> 1 (vender)
# Else -> 2 (nada)df_var_exog['rsi_indicator'] = df_var_exog['rsi'].apply(lambda x: 0 if x<30 else 1 if x>70 else 2) 
display(df_var_exog[['Adj Close','rsi','rsi_indicator']].tail(15))

No gráfico acima estão indicados a variação de RSI para o período, assim como a faixa correspondente à região neutra (sem indicação de compra ou venda).

Indicador de Bollinger

O Indicador de Bollinger mostra a volatilidade do ativo. Seu cálculo é realizado em função da média móvel e do desvio padrão móvel.
Se o indicador é maior que 1, é uma indicação de venda (1). Caso seja menor que 0, então é sinal de compra (0).

# calculando a média móvel e limites superior e inferiror
# limites com base em 2 desvios padrão 
up, mid, low = BBANDS(df_var_exog['Adj Close'], timeperiod=20, nbdevup=2, nbdevdn=2, matype=0)# criando features para a média e os limites
df_var_exog['upper'] = up
df_var_exog['mid'] = mid
df_var_exog['low'] = low
df_var_exog['bbp'] = (df_var_exog['Adj Close'] - df_var_exog['low'])/(df_var_exog['upper'] - df_var_exog['low'])
df_var_exog.dropna(inplace=True)# bbp < 0 -> 0 (comprar)
# RSI > 1 -> 1 (vender)
# Else -> 2 (nada)df_var_exog['bbp_indicator'] = df_var_exog['bbp'].apply(lambda x: 0 if x<0 else 1 if x>1 else 2) 
display(df_var_exog[['Adj Close','bbp','bbp_indicator']].tail(15))

A região em destaque (entre 0 e 1) é aquela em que não se deve nem comprar e nem vender (por indicação do Bollinger). Essa decisão é indicada pelo valor 2 em nosso dataset.

Suporte e resistência

Os indicadores de Suporte e Resistência também são bastante utilizados na avaliação técnica de ações.

Resumidamente, uma linha de suporte por exemplo caracteriza um fundo, onde a pressão compradora demonstrou-se maior que a pressão vendedora em algum (ou vários) momento no passado. Por outro lado, uma linha de resistência caracteriza um topo, onde a pressão vendedora foi maior que a compradora.

Suporte e Resistência podem ser definidos através de, no mínimo, 3 candles. Dessa forma, vamos criar features de suporte e resistência, baseadas em 5 candles consecutivos.

# definindo a função de resistencia
def is_resistance(df,i):
  resistance = (df['High'][i] > df['High'][i-1]
                and df['High'][i] > df['High'][i+1]
                and df['High'][i+1] > df['High'][i+2]
                and df['High'][i-1] > df['High'][i-2])
  return resistance# definindo a função de suporte
def is_support(df,i):
  support = (df['Low'][i] < df['Low'][i-1]
             and df['Low'][i] < df['Low'][i+1]
             and df['Low'][i+1] < df['Low'][i+2]
             and df['Low'][i-1] < df['Low'][i-2])
  return support# resistência verdadeiro -> 1 (vender)
# suporte verdadeiro -> 0 (comprar)
# outros (2)# criando feature com valores 2
df_var_exog['suport_resistencia'] = 2# definindo os valores 1 e 0
for i in range(2, df_var_exog.shape[0] - 2):
  if is_resistance(df_var_exog,i):
    df_var_exog['suport_resistencia'][i] = 1 # definindo 1 para resistência
  elif is_support(df_var_exog,i):
    df_var_exog['suport_resistencia'][i] = 0 # definindo 0 para suporte# contando os suportes e resistÊncias
df_var_exog.suport_resistencia.value_counts()

Após a aplicação das funções de suporte e resistência criadas, identificamos 149 pontos de suporte e 152 de resistência.

LTA e LTB

De forma geral, as variáveis lta e ltb devem nos dizer se, em determinado período de tempo, nossa série está subindo, caindo, ou andando horizontalmente.

As análises podem ter sensibilidades diferentes caso a tendência seja de alta ou de baixa, por isso essa variável pode ser interessante.

Para isso, vamos utilizar o método rolling com correlation do pandas. No nosso caso, definimos que correlações maiores que 0.5 são classificadas como tendência de alta, menos que -0.5 são tendência de baixa e o restante é a lateralidade.

df_var_exog2 = df_var_exog.reset_index()
df_var_exog['corr'] = (df_var_exog2['Adj Close'].rolling(20).corr(pd.Series(df_var_exog2.index))).tolist()
df_var_exog.dropna(inplace=True)def condition(x):
    if x<=-0.5:
        return -1
    elif x>-0.5 and x<0.5:
        return 0
    else:
        return 1
df_var_exog['corr_class'] = df_var_exog['corr'].apply(condition)# plotando os dois gráficos
sns.scatterplot(x = df_var_exog.index,y = df_var_exog['corr'])
sns.countplot(df_var_exog['corr_class'])

Podemos verificar que há um balanceamento entre as classes de tendências do valor da ação.

Conclusões da EDA

Nossa série de Adj Close da AAPL horaria apresenta uma série não estacionária
Como mostrado pelos dados agrupados mensalmente, a série possui sazonalidade
Há outliers na série (outliers considerados com média móvel +/- 3*sigma móvel)
Os outliers foram tratados com o método de Windorize
Com 1 diferenciação, a série se tornou estacionária (parâmetro d=1 do ARIMA)
Os resíduos ficam maiores no final do período analisado
O indicador RSI indicou 41 sinais de compra e 133 de venda
O indicador de Bollinger mostrou 137 sinais de compra e 123 de venda
Os indicadores de LTA e LTB mostram que há quantidades parecidas para indicação de alta, baixa e lateralidade

Modelos de Machine Learning

Nessa etapa, já entendemos como os dados se comportam, tivemos mais contato com o dataset e fizemos algum feature engineering, com conhecimentos adquiridos do mercado especialista.

Agora, portanto, vamos treinar e testar alguns modelos de ML, para verificar se podemos ter resultados expressivos na predição de valores das ações.

Assim, vamos ter um modelo como baseline. O escolhido foi o ARIMA, com os parâmetros mais básicos. Como métrica, utilizaremos o MAE (Mean Absolut Error).

Modelo Baseline — ARIMA

A primeira implementação do modelo ARIMA, foi com os parâmetros p=2, d=1, q=2, fazendo a previsão com apenas uma dataset de treinamento.

model = ARIMA(train_aapl_hora, order=(2,1,2))
model_fit = model.fit(disp=0)# verificando os resultados
step = len(test_aapl_hora)# fc (forecast), se ()
fc, se, conf = model_fit.forecast(step)fc = pd.Series(fc, index = test_aapl_hora[:step].index)
lower = pd.Series(conf[:,0], index = test_aapl_hora[:step].index)
upper = pd.Series(conf[:,1], index = test_aapl_hora[:step].index)plt.figure(figsize=(16,8))
plt.plot(train_aapl_hora[1000:], label = 'treino', color='skyblue')
plt.plot(test_aapl_hora[:step], label = 'teste')
plt.plot(fc, label='previsto')
plt.fill_between(lower.index, lower, upper, color = 'k', alpha = 0.1)
plt.title('Previsto vs atual')
plt.legend(loc='best')

Como observado no gráfico de previsões acima, a previsão feita dessa maneira ficou muito aquém de um resultado confiável. Assim, resolvemos fazer um treinamento a cada previsão, realimentando o dataset, o que diminuiria muito o impacto das previsões de longo prazo. O resultado foi o seguinte:

Para esse último teste, obtivemos um resultado visivelmente melhor, e com MAE=1.0450.

Modelos “tradicionais”

Partindo para modelos mais “tradicionais”, vamos testar uma metodologia parecida, porém com modelos que conhecemos melhor, como LGBM, Random Forest e Árvore de Decisão.

Para isso, precisamos construir features defasadas passadas com a seguinte função:

def constroi_features_defasadas(base,lista_features,defasagem_maxima):
    # Constrói features defasadas com base na base original
    # Copia a base
    base_cop = base.copy()
    for feat in lista_features:       
        for i in range(1,defasagem_maxima+1):
            base_cop[str(feat)+'_def_'+str(i)] = base_cop[feat].shift(i)
    return base_cop

Assim, ficamos com o seguinte dataset com 30 features defasadas.

Avaliando o resultado segundo a metrica MAE, encontramos que o melhor modelo foi o Random Forest.

modelos_candidatos = {'Árvore':DecisionTreeRegressor(max_depth=5),
                      'RandomForest':RandomForestRegressor(max_depth=5),
                      'LGBM':LGBMRegressor(max_depth=5)
                     }
avaliacao = {}
for nome,model in modelos_candidatos.items():
    # Avalia a crossvalidação
    score = cross_val_score(model, # Escolhendo o nosso modelo da vez
                            X_train, y_train, # Nossos dados, excluindo o teste
                            cv=TimeSeriesSplit(n_splits=5), # Validação cruzada temporal
                            scoring='neg_mean_absolute_error', # Usando a métrica MAE
                           n_jobs=8 # Número de processadores, para ser mais rápido
                           ).mean() # Tirando a média de todos os folds
    avaliacao[nome] = -score # -score e não score para tornar o número positivoavaliacao{'LGBM': 2.432275552975372,  'RandomForest': 2.246902277238342,  'Árvore': 2.3692681121851775}

Com o melhor modelo em mãos, podemos verificar qual o melhor número de features defasadas, com relação ao MAE.

O resultado sugere que o melhor número de features é 20.

Apesar de termos resultados aparentemente interessantes, não estávamos conseguindo chegar muito próximo (percentualmente) dos valores. Além disso, modelos como SARIMAX ou mesmo o ARIMA com os parâmetros otimizados se tornaram inviáveis para essa análise, uma vez que necessitavam de muito tempo para treinamento.

Assim, decidimos partir para uma outra questão.

Novo Insight

Após os resultados obtidos, tivemos um novo insight sobre nossos modelos: E se ao invés de tentarmos prever o valor exato, tentássemos identificar se o valor subiria ou desceria? Isso faz sentido do ponto de vista de negócio, uma vez que essa é uma das principais informações que um trader precisa para tomar suas decisões.

A partir daí, mudamos nosso foco, que era de modelos de regressão, para modelos de classificação.

Pre Processamento — Classificação

Agora, nosso target não será mais o valor da ação, e sim um binário [0,1], que indicam se o valor subiu ou desceu com relação ao valor da hora anterior.

Assim, criamos essa nossa target, com relação aos valores anteriores. Para isso, criamos a seguinte função:

def target(df):# criando feature com 1h de defasagem (pegando a linha de cima)
  df['def_1'] = df['Adj Close'].shift(1)
  # criando feature comparando valor atual com o defasado
  df['subt'] = df['Adj Close'] - df['def_1']# criando a target de subida ou descida do valor da ação
  #0 -> caiu (com relação ao anterior)
  #1 -> subiu (com relação ao anterior)
  #2 -> igual ao anteriordf['target'] = df['subt'].apply(lambda x: 0 if x<0 else 1 if x>0 else 2)return df

Além disso, vamos utilizar 20 features defasadas, como nos modelos de regressão, e algumas features exógenas, como rsi, bollinger, temporais, lta/ltb, e suporte/resistência.

É durante o pre-processamento quando separamos nosso dataset em treino e teste (lembrando de não fazer a separação aleatória).

def preprocessamento(base, corte_treino_teste,target):
  # separando a base da target
  X = base.drop(target, axis=1)
  y = base[target]  # cortando em treino e teste
  X_train = X[X.index<=corte_treino_teste]
  X_test = X[X.index>corte_treino_teste]  y_train = y[y.index<=corte_treino_teste]
  y_test = y[y.index>corte_treino_teste]  return X_train, X_test, y_train, y_test

Modelos Classificação

Os modelos a serem treinados nessa etapa serão: Árvores de Decisão, Random Forest e XGBoost. Primeiro testamos com o Cross Validation, para previsões com 1 hora de antecedência.

modelos_candidatos = {'Árvore':DecisionTreeClassifier(random_state=42),
                      'RandomForest':RandomForestClassifier(random_state=42),
                      'XGBoost':XGBClassifier(random_state=42)
                     }
def avaliacao_modelos(modelos,X_train,y_train):
  avaliacao = {}
  scores = []
  ts = TimeSeriesSplit(n_splits=5).split(X_train,y_train)
  for nome,model in modelos_candidatos.items():
      # Avalia a crossvalidação
      score = cross_val_score(model, # Escolhendo o nosso modelo da vez
                              X_train, y_train, # Nossos dados, excluindo o teste
                              cv=TimeSeriesSplit(n_splits=5).split(X_train,y_train), # Validação cruzada temporal
                              scoring='accuracy', # Usando a acurácia como métrica
                              ).mean() # Tirando a média de todos os foldsavaliacao[nome] = score
      scores.append(score)
  
  return avaliacao, scores

Para essa configuração, obtivemos os seguintes resultados de acurácia:

Agora, vamos verificar como ficam as previsões mais antecipadas, até 8h de antecedência. O resultado foi o seguinte para acurácia:

Tivemos o melhor resultado com previsão de 1 horas de antecedência, com acurácia de 0,744554, para o modelo XGBoost. Porém, nota-se que não há uma discrepância muito grande, e que talvez utilizar mais horas de antecedência seja possível.

Testando Modelo Campeão

Como nosso melhor modelo foi o XGBoost, vamos utiliza-lo para teste com valores novos.

base_cop = base.copy()base_cop = constroi_features_futuras(base_cop,'target',1)
base_cop.drop('target', axis=1, inplace=True)X1_train, X1_test, y1_train, y1_test = preprocessamento(base_cop, '2022-02-01', 'target_fut_1')# criando df de resultados
result = pd.DataFrame()
# fitando o modelo
xgb = XGBClassifier(random_state=42,max_depth=5)
xgb.fit(X1_train, y1_train)
y1_pred = xgb.predict(X1_test)
y1_proba = xgb.predict_proba(X1_test)
y1_proba = y1_proba[:, 1]# criando as colunas de resultados
result['y1_test'] = y1_test
result['y1_predict'] = y1_pred
result['proba'] = y1_probaprint(classification_report(y1_test, y1_pred))# criando a matriz de confusão
ax = sns.heatmap(confusion_matrix(y1_test, y1_pred),
                 annot=True,
                 annot_kws={"fontsize":10},
                 fmt = 'd',
                 cmap = 'Blues')

O resultado mostra bons resultados de acurácia e previsões para valores 0. Isso era mais esperado, por conta do desbalanceamento da nossa target (aproximadamente 3:1). O interessante desse resultado é a alta performance para identificação de valores 0.

Importância das features

Como possuímos muitas features, podemos fazer uma avaliação de quais são as mais importantes para nossa previsão. Faremos isso com o método feature_importances para o modelo XGBoost.

feature_importance = xgb.feature_importances_
sorted_idx = np.argsort(feature_importance)
fig = plt.figure(figsize=(12, 6))
plt.barh(range(len(sorted_idx)), feature_importance[sorted_idx], align='center')
plt.yticks(range(len(sorted_idx)), np.array(X1_test.columns)[sorted_idx])
plt.title('Feature Importance')

A partir dessa análise, observamos que features temporais são as menos importantes dentre as exógenas, enquanto suporte/resistência, lta/ltb e o indicador de Bollinger são os mais importantes.

Definindo tamanho do Treino

Como nosso treinamento deve ser realizado muitas vezes, se não a cada previsão, é importante que tenhamos um tamanho ótimo do dataset. Para isso, rodamos o modelo para diversos tamanhos de treino, da seguinte forma.

X_train, X_test, y_train, y_test = preprocessamento(base_cop, '2022-03-15', 'target_fut_1')
accuracys = []
f1s = []
tempos = []for pace in range(300, 1500, 50):
  print('teste com:',pace)
  X_tr = X_train[-pace:]
  y_tr = y_train[-pace:]
  X_t = X_test.copy()
  y_t = y_testy_preds = []
  y_preds_proba = []temp = []for i in range(len(y_test)):
      # medindo o tempo de treinamento para cada pace
      start = time.time()
      # fitando o modelo
      xgb = XGBClassifier(random_state=42,max_depth=5)
      xgb.fit(X_tr, y_tr)
      y_pred = xgb.predict(X_t[:1])
      y_pred = y_pred[0]
      y_proba = xgb.predict_proba(X_t[:1])
      y_proba = y_proba[:, 1]
      end = time.time()
      temp.append((end-start))y_preds.append(y_pred)
      y_preds_proba.append(y_proba)X_tr = pd.concat([X_tr, X_t[:1]])
      X_tr = X_tr[1:]
      y_tr = y_tr.append(pd.Series(y_t[0]))
      y_tr = y_tr[1:]X_t = X_t[1:]
      y_t = y_t[1:]accuracys.append(accuracy_score(y_test, y_preds))
  f1s.append(f1_score(y_test, y_preds, average='weighted'))
  tempos.append(np.mean(temp))df = pd.DataFrame({'accuracy': accuracys,
                   'f1_score': f1s,
                   'tempo(s)': tempos}, index=list(range(300, 1500, 50)))
df

Observamos que o tempo de treino e teste cresce linearmente, variando de 0.11s com 300 linhas no dataset de treino, para 0.5s com 1450 linhas.

Além disso, o dataset de treino com 600 linhas obteve melhores métricas para acurácia (0.78) e f1-score weighted (0.72).

Otimização de Hiperparâmetros

Agora, com o tamanho de dataset de treino e modelo ideais, podemos utilizar o GridSearchCV para otimizar os hiperparâmetros.

X_train, X_test, y_train, y_test = preprocessamento(base_cop, '2022-02-15', 'target_fut_1')# otimização para XGBoost undersampled
xgb = XGBClassifier(random_state=42)parametros_xgb = {
    'n_estimators': range(50,250,50),
    'max_depth': range(2,10,2),
    'gamma':[0.1, 0.5, 0.75]}
    
# Obtendo melhores parâmetros e melhor recall para XGBoost
grid_xgb = GridSearchCV(estimator = xgb,
                        param_grid = parametros_xgb,
                        cv=5,
                        scoring = 'f1_weighted')grid_xgb.fit(X_train[-600:], y_train[-600:]) # somente a melhor qtd encontrada
xgb_over_best_params = grid_xgb.best_params_
xgb_over_best_score = grid_xgb.best_score_print(grid_xgb.best_params_)
print(grid_xgb.best_score_)

Encontramos que os hiperparâmetros otimizados para nosso modelo são os seguintes: gamma = 0.1, max_depth = 8 e n_estimators=100.

Deploy do Modelo no Streamlit

O modelo final, com os hiperparâmetros otimizados foi colocado em produção na cloud do Streamlit, e pode ser encontrado no seguinte link: https://andrekuniyoshi-ter-streamlitmercado-financeiro-streamlit-6z1ehf.streamlitapp.com/.

O código do programa colocado em produção no Streamlit pode ser encontrado aqui.

Alguns pontos importantes sobre a plataforma gratuita e processamento do modelo:

A cada atualização da página, nosso modelo faz todo o processo, desde a coleta dos dados, até o treinamento do modelo
Isso já dificulta o processo, deixando-o mais lento
Porém, não fica dependente de métodos pagos da cloud, para requests programadas
nosso modelo rodou no Google Colab em 0.23s, enquanto em produção está demorando cerca de 4min.

Conclusões

Nosso produto final é capaz de fazer previsões para do valor da ação da AAPL, com acurácia de cerca de 75%. Isso ainda está longe do que esperávamos, uma vez que nossa base é desbalanceada em proporção semelhante à essa.

No entanto, no início do projeto partimos de acurácia de 50% e evoluímos incluindo novas features baseadas em estudos do funcionamento do mercado.

Além disso, fomos capazes de otimizar tanto o processamento do modelo, verificando tamanho de dataset ideal, como com relação aos hiperparâmetros.

Apesar de ainda estar longe do ideal, a ideia de mudar de modelos de regressão para classificação foi acertada, pois para a visão da área de negócio, esse último faz mais sentido e torna o modelo mais simples e robusto.

Recomendações para continuidade do trabalho

Para melhorar as métricas do modelo, deve-se identificar mais features que podem ajudar.
Utilizar previsões em outra escala, como a de minutos, por exemplo. Há plataformas em que se pode ganhar dinheiro caso adivinhe se a ação vai subir ou cair. Pode ser um mercado a ser explorado.
O treinamento não precisa ser realizado a cada requisição. Ele pode ser feito uma vez ao dia, por exemplo, o que garantiria uma performance em produção muito melhor.
Os dados não precisam ser coletados a cada requisição, caso tenhamos espaço alocado em alguma cloud. (indicaria isso, caso essa ideia virasse, de fato, um negócio)
O treshold utilizado para as análises foi o default (0.5). No entanto, uma análise mais profunda deve ser realizada, a fim de mostrar quais as faixas de treshold podemos confiar mais. Exemplo: só indicaria compras com probabilidade > 0.7 e vendas com probabilidade < 0.3.