Séries temporais com Machine Learning — Parte 2
Uma breve comparação entre a abordagem usando a estatística clássica e o uso dos modelos de Machine Learning em quatro artigos
No primeiro artigo fizemos uma breve introdução do que é uma série temporal, citamos os modelos da estatística clássica (ARMA, ARIMA, SARIMA, VARIMA, etc.) e propusemos o uso de modelos de Machine Learning para esse tipo de aplicação.
Escolhemos um tema a partir de dados de uma competição do Kaggle e estabelecemos um baseline com o Facebook Prophet.
Melhorando a performance da predição
Primeiramente vamos entender um pouco melhor qual o comportamento das visualizações das páginas da NFL na Wiki, avaliaremos a distribuição e a média das visualizações, por fim a estacionariedade da série.
Abaixo, temos as médias mensal e para cada dia da semana das visualizações das páginas escolhidas.
O comportamento das médias mensais parece acompanhar o calendário da NFL, com um aumento à partir de setembro, quando começa temporada regular, e uma diminuição à partir de fevereiro, quando acontece o Super Bowl.
O comportamento das médias diárias tem maior valor as segundas-feira, o que pode ser explicado pelo maior número de jogos que acontece no domingo com o interesse das pessoas em acompanhar os resultados e notícias sobre os times.
Por fim, a distribuição das visualizações no formato de histograma. A distribuição concentra-se até 5.000 visualizações e apresenta uma cauda longa com valores indo de 10.000 a 30.000 visualizações.
Além do histograma com a distribuição das visualizações, o teste de normalidade de Kolgomorov-Smirnov ( KstestResult statistic=1.0, pvalue=0.0) revelou que a série é normal.
Com o entendimento de como as observações estão distribuídas, podemos construir alguns modelos simples para compreender mais profundamente a série e discutir suas implicações.
Mover, Suavizar, Avaliar
Uma abordagem para analisar esse comportamento é usar a média móvel. Ao suavizar a série temporal original conseguimos identificar qual a tendência.
Podemos fazer isso usando a função DataFrame.rolling(janela).mean() disponível no Pandas. Quanto maior a janela, mais suave a tendência. No caso de dados muito ruidosos, geralmente encontrados em finanças, esse procedimento pode ajudar a detectar padrões comuns.
Primeiro vamos testar uma janela de 4 dias.
E depois com uma janela de 7 dias.
Conforme esperado, quando maior a janela mais suave a média móvel. As primeiras observações em julho de 2015 coincidem com o final do período de treino dos times, a partir de setembro, quando começa a temporada regular, os valores aumentam, e tem um pico entre dezembro e janeiro, período de playoffs e SuperBowl, cai ao final do evento e volta a subir novamente em setembro de 2016.
Os intervalos de confiança ao redor da série observada, mostram que os dados tem alguns outliers, que não são capturados pela média móvel. Caso quiséssemos usar essa abordagem para prever essa série, teríamos problemas. Logo, é interessante uma abordagem mais robusta capaz de lidar com outliers.
Outro fator importante é que todas as observações tem o mesmo peso no cálculo da média. Adicionar pesos as observações pode nos ajudar a lidar com esse comportamento.
Uma abordagem tradicional é adicionar pesos iguais a todas as observações, nesse caso começarmos a ponderar todas as observações disponíveis enquanto diminuímos exponencialmente os pesos à medida que avançamos no tempo. Existe uma fórmula para suavização exponencial que nos ajudará com isso:
𝑦̂ 𝑡+1=𝛼⋅𝑦𝑡+(1−𝛼)⋅𝑦̂ 𝑡−1y^t+1=α⋅yt+(1−α)⋅y^t−1
Aqui, o valor do modelo é uma média ponderada entre o valor verdadeiro atual e os valores anteriores do modelo. O peso alpha é chamado de fator de suavização. Ele define a rapidez com que “esqueceremos” a última observação verdadeira disponível. Quanto menor o alpha, maior a influência das observações anteriores e mais suave a série.
Para maiores valores de alpha a série projetada tende a capturar melhor o comportamento da série observada. A medida que alpha diminui as médias se suavizam cada vez mais e a projeção se descola cada vez mais do observado.
O próximo passo será adicionar mais uma etapa de suavização da série.
A decomposição em série nos ajudará a obtermos dois componentes: intercepto (isto é, nível) ℓ e inclinação (isto é, tendência) 𝑏. Aplicaremos a mesma suavização exponencial à tendência, assumindo que a direção futura da série temporal mude dependendo das alterações ponderadas anteriores. Como resultado, obtemos o seguinte conjunto de funções:
ℓ𝑥=𝛼𝑦𝑥+(1−𝛼)(ℓ𝑥−1+𝑏𝑥−1)
𝑏𝑥=𝛽(ℓ𝑥−ℓ𝑥−1)+(1−𝛽)𝑏𝑥−1
𝑦̂ 𝑥+1=ℓ𝑥+𝑏𝑥
O primeiro descreve a interceptação, que, como antes, depende do valor atual da série. O segundo termo agora está dividido em valores anteriores do nível e da tendência. A segunda função descreve a tendência, que depende das mudanças de nível na etapa atual e do valor anterior da tendência. Nesse caso, o coeficiente 𝑏𝑒𝑡𝑎 é um peso para a suavização exponencial. A previsão final é a soma dos valores do modelo da interceptação e tendência.
Testamos algumas combinações de alpha e beta, e obtemos uma melhor predição com o alpha de 0.9 e beta de 0.02, com o resultado praticamente em cima da série observada em roxo. Dada a formulação da suavização exponencial dupla nossa série tem maior influência dos pontos anteriores em torno da tendência, dada por alpha, e menor relação com o nível anterior, a própria tendência, dado por beta.
Na sequência vamos discutir a suavização exponencial tripla e como escolher os parâmetros de forma a evitar que algumas combinações produzam resultados estranhos.
A idéia é adicionar um terceiro componente — a sazonalidade. Isso significa que não devemos usar esse método se não for esperado que nossa série temporal tenha sazonalidade. Os componentes sazonais no modelo explicam variações repetidas em torno da interceptação e tendência, e serão especificadas pela duração do evento, ou seja, pelo período após o qual as variações se repetem. Para cada observação na temporada, há um componente separado; por exemplo, se a duração da temporada for 7 dias (uma sazonalidade semanal), teremos 7 componentes sazonais, um para cada dia da semana.
Com isso, vamos escrever um novo sistema de equações:
ℓ𝑥=𝛼(𝑦𝑥−𝑠𝑥−𝐿)+(1−𝛼)(ℓ𝑥−1+𝑏𝑥−1)
𝑏𝑥=𝛽(ℓ𝑥−ℓ𝑥−1)+(1−𝛽)𝑏𝑥−1
𝑠𝑥=𝛾(𝑦𝑥−ℓ𝑥)+(1−𝛾)𝑠𝑥−𝐿
𝑦̂ 𝑥+𝑚=ℓ𝑥+𝑚𝑏𝑥+𝑠𝑥−𝐿+1+(𝑚−1)𝑚𝑜𝑑𝐿
A interceptação agora depende do valor atual da série menos qualquer componente sazonal correspondente. A tendência permanece inalterada, e o componente sazonal depende do valor atual da série menos a interceptação e do valor anterior do componente. Leve em consideração que o componente é suavizado por todas as estações disponíveis; por exemplo, se tivermos um componente de segunda-feira, será calculado apenas a média de outras segundas-feiras. Agora que temos o componente sazonal, podemos prever não apenas um ou dois passos à frente, mas também um futuro arbitrário de visualizações, o que é muito encorajador.
Abaixo o resultado do modelo triplo de suavização exponencial, também conhecido pelos sobrenomes de seus criadores, Charles Holt e seu aluno Peter Winters.
𝑦̂𝑚𝑎𝑥𝑥=ℓ𝑥−1+𝑏𝑥−1+𝑠𝑥−𝑇+𝑚⋅𝑑𝑡−𝑇
𝑦̂𝑚𝑖𝑛𝑥=ℓ𝑥−1+𝑏𝑥−1+𝑠𝑥−𝑇−𝑚−𝑑𝑡−𝑇
𝑑𝑡=𝛾∣𝑦𝑡−𝑦̂ 𝑡∣+(1−𝛾)𝑑𝑡−𝑇,
onde 𝑇 é a duração da temporada, 𝑑 é o desvio previsto. Outros parâmetros foram obtidos da tripla suavização exponencial.
Antes de começarmos a construir um modelo, vamos discutir como estimar os parâmetros do modelo automaticamente.
Não há nada incomum aqui; como sempre, temos que escolher uma função de perda adequada para a tarefa que nos dirá quão próximo o modelo se aproxima dos dados. Em seguida, usando a validação cruzada, avaliaremos nossa função de perda escolhida para os parâmetros de modelo fornecidos, calcularemos o gradiente, ajustamos os parâmetros do modelo e assim por diante, eventualmente descendo para o mínimo global.
Você pode estar se perguntando como fazer a validação cruzada para séries temporais porque as séries temporais possuem essa estrutura temporal e não é possível misturar valores aleatoriamente em uma dobra enquanto preserva essa estrutura. Com a randomização, todas as dependências de tempo entre as observações serão perdidas. É por isso que teremos que usar uma abordagem mais complicada para otimizar os parâmetros do modelo, a “validação cruzada em uma base contínua”.
A idéia é bastante simples — treinamos nosso modelo em um pequeno segmento da série cronológica desde o início até alguns 𝑡, fazemos previsões para as próximas etapas de 𝑡+𝑛 e calculamos um erro. Em seguida, expandimos nossa amostra de treinamento para o valor de 𝑡+𝑛, fazemos previsões de 𝑡+𝑛 até 𝑡+2∗𝑛 e continuamos movendo nosso segmento de teste da série temporal até atingirmos a última observação disponível. Como resultado, temos tantas dobras quanto 𝑛 entre a amostra de treinamento inicial e a última observação.
No modelo de Holt-Winters, assim como nos outros modelos de suavização exponencial, há uma restrição quanto ao tamanho dos parâmetros de suavização, cada um deles variando de 0 a 1. Portanto, para minimizar nossa função de perda, precisamos escolher um algoritmo que suporte restrições nos parâmetros do modelo. No nosso caso, usaremos o gradiente conjugado de Newton truncado.
Abaixo o gráfico com os resultados do modelo de Holt Winters para uma janela de predição de 60 dias. O MAPE calculado foi de 17,59%, um ganho de praticamente 10 p.p. em relação ao baseline.
O modelo conseguiu aproximar com sucesso as séries temporais iniciais, capturando a sazonalidade diária, a tendência geral de queda e até algumas anomalias, reagindo bem às mudanças na estrutura da série e depois retornando rapidamente aos valores normais, essencialmente “esquecendo” o passado. Esse recurso do modelo nos permite construir rapidamente sistemas de detecção de anomalias, mesmo para dados ruidosos da série, sem gastar muito tempo e dinheiro na preparação dos dados e no treinamento do modelo.
A abordagem econométrica
Antes de começarmos a modelar, devemos mencionar uma propriedade tão importante da série temporal: estacionaridade.
Se um processo é estacionário, isso significa que ele não altera suas propriedades estatísticas ao longo do tempo, ou seja, sua média e variância. A função de covariância não depende do tempo; deve depender apenas da distância entre as observações. Para mais informações clique aqui.
A estacionariedade é uma propriedade importante, pois nos permiti assumir que as propriedades estatísticas futuras não serão diferentes daquelas atualmente observadas. A maioria dos modelos de séries temporais, de uma maneira ou de outra, tenta prever essas propriedades (média ou variação, por exemplo). As previsões estariam erradas se a série original não estivesse estacionária. Infelizmente, a maioria das séries temporais que vemos fora dos livros didáticos não é estacionária, mas podemos (e devemos) mudar isso.
Então, para combater a não estacionariedade, precisamos conhecer nosso inimigo, por assim dizer.
O StatsModels nos fornece as função sm.tsa.seasonal_decompose para ajudar nessa tarefa. Esta é uma decomposição ingênua. Métodos mais sofisticados devem ser preferidos.
O modelo aditivo é Y [t] = T [t] + S [t] + e [t]
O modelo multiplicativo é Y [t] = T [t] * S [t] * e [t]
O componente sazonal é removido primeiro aplicando um filtro de convolução aos dados. A média dessa série suavizada para cada período é o componente sazonal retornado.
A decomposição da série revela uma tendência de aumento nas observações correspondentes aos meses de pós temporada (dezembro e janeiro), um decaimento em fevereiro, e uma retomada em setembro.
Vamos analisar a autocorrelação (ACF) e autocorrelação parcial (PACF) da série com a ajuda da função sm.tsa.stattools.adfuller.
O gráfico a seguir mostra a deterioração ACF com os intervalo de uma semana, seguindo a ocorrência dos jogos em temporada regular. As linhas azuis representam o intervalo de confiança dos valores estimados de autocorrelação. Se o valor estiver além do intervalo, consideraremos a correlação automática como significativa. Como observado, os valores do ACF decaem exponencialmente.
O PACF permite que a medida da correlação entre 𝑥𝑡 e 𝑥𝑡+ℎ seja medida diretamente, enquanto remove o efeito de correlação de qualquer ponto no tempo entre eles.
As séries com rolagem de 7 mostram-se são estacionárias; o teste de Dickey-Fuller rejeitou a hipótese nula de que uma raiz unitária está presente. Na verdade, podemos ver isso no próprio gráfico — não temos uma tendência visível; portanto, a média é constante e a variação é praticamente estável.
O modelo ARIMA
ARIMA é uma série que precisa ser diferenciada para tornar-se estacionária como uma série “integrada” (I).
As defasagens da série estacionarizada são chamadas de “auto-regressivas” que se referem aos termos (AR) e as defasagens dos erros de previsão são chamadas de “média móvel”, que se refere aos termos (MA). É basicamente usado para previsão
ARIMA é um modelo aleatório generalizado que é ajustado para eliminar toda a auto-correlação residual. É um modelo de suavização exponencial generalizada que pode incorporar tendências e sazonalidade a longo prazo.
O modelo de regressão estacionarizado usa atrasos das variáveis dependentes e / ou atrasos dos erros de previsão como regressores. Aqui, o modelo de previsão de séries temporais pode ser estacionarizado usando transformações como diferenciação, registro e deflação.
Com isso, podemos dizer que uma série temporal é “estacionária” se todas as propriedades estatísticas, como média, variância, auto-correlação etc., são constantes no tempo.
Construção de um modelo ARIMA:
- Estacionarizar a série, se necessário, por diferenças (e talvez também pelo uso de LOG, deflação, etc.).
- Estudar o padrão da auto-correlação e da auto-correlação parcial para determinar os LAGs ou atrasos da série estacionária, e/ou LAGs ou atrasos dos erros de predição que devem estar inclusos na equação de predição.
- Ajuste do modelo que é sugerido e checagem diagnóstica do seu resíduo, particularmente dos resíduos dos gráficos ACF e PACF, para ver se todos os coeficientes são significativos e todos os padrões estão explicados.
- Padrões que permanecem no ACF e PACF sugerem a necessidade de termos adicionais para AR e MA.
Terminologias em ARIMA
O modelo ARIMA pode ser (quase) completamente resumido por três números:
- p = o número de termos auto-regressivos
p é o número de termos auto-regressivos (parte AR). Permite incorporar o efeito de valores passados em nosso modelo. Intuitivamente, isso seria semelhante ao afirmar que é provável que esteja quente amanhã se estiver quente nos últimos três dias.
- d = o número de diferenças não sazonais
d é o número de diferenças não sazonais necessárias para a estacionariedade. Intuitivamente, isso seria semelhante ao afirmar que provavelmente haverá a mesma temperatura amanhã se a diferença de temperatura nos últimos três dias tiver sido muito pequena.
- q = o número de termos da média móvel
q é o número de erros de previsão atrasados na equação de previsão (parte MA). Isso nos permite definir o erro do nosso modelo como uma combinação linear dos valores de erro observados em momentos anteriores no passado.
Estes são os três números inteiros (p, d, q) usados para parametrizar os modelos ARIMA.
Agora que sabemos quais são os parâmetros do modelo, vamos testar algumas combinações deles à seguir:
p_values = [6, 7], no caso, o número de visualizações dessa semana será igual ao da semana anterior.
d_values = range(0, 2)
q_values = range(0, 4)
O objetivo é minimizar o MAPE de acordo com essa combinação de intervalos de parâmetros.
ARIMA(6, 0, 0) MAPE=21.869
ARIMA(6, 0, 1) MAPE=19.948
ARIMA(6, 0, 2) MAPE=19.093
ARIMA(6, 0, 3) MAPE=17.473
ARIMA(6, 1, 0) MAPE=17.988
ARIMA(6, 1, 1) MAPE=16.991
ARIMA(6, 1, 2) MAPE=17.007
ARIMA(6, 1, 3) MAPE=16.977
ARIMA(7, 0, 0) MAPE=17.639
ARIMA(7, 0, 1) MAPE=16.492
ARIMA(7, 0, 2) MAPE=16.497
ARIMA(7, 0, 3) MAPE=16.434
ARIMA(7, 1, 0) MAPE=16.973
ARIMA(7, 1, 1) MAPE=16.985
ARIMA(7, 1, 2) MAPE=16.970
ARIMA(7, 1, 3) MAPE=16.965
Best ARIMA(7, 0, 3) MAPE=16.434Abaixo o gráfico comparando a predição do modelo e os dados observados.
No final, obtivemos previsões bastante adequadas, em média, nosso modelo teve um MAPE de 16,43%, o que representa um ganho de 1.16 p.p. em relação ao modelo de Holt Winters, o que é muito, muito bom, mas os custos gerais de preparação dos dados, tornando a seleção de parâmetros estacionários e de força bruta para fazer o ajuste dos parâmetros, podem não valer essa precisão.
Sobrevivemos aos nove círculos do inferno para finalizar o modelo ARIMA. Nos próximos artigos vamos voltar ao Facebook Prophet e aos modelos de Machine Learning aplicados a esse problema.
Acesse o projeto completo no Github:
https://github.com/SouzaCadu/Time-Series
Referências
