Tag

arima Archives - Análise Macro

{tidyverts}: séries temporais no R

By | Data Science

{tidyverts} é uma família de pacotes de R criada para ser a próxima geração de ferramentas para modelagem e previsão de séries temporais, substituindo o famoso pacote {forecast}.

Utilizando uma interface simples e integrada com os pacotes do tidyverse, é possível construir uma ampla gama de modelos de previsão univariados e multivariados: ARIMA, VAR, suavização exponencial via espaço de estado (ETS), modelo linear (TSLM), autorregressivo (AR), passeio aleatório (RW), autoregressão de rede neural (NNETAR), Prophet, etc.

Neste exercício, daremos uma breve introdução aplicada a modelagem e previsão de séries temporais usando os pacotes do {tidyverts}.

Conhecendo os pacotes

Um breve resumo do que a família de pacotes do {tidyverts} tem a oferecer:

fable

  • Coleção de modelos univariados e multivariados de previsão
  • Modelagem de séries temporais em formato "tidy"
  • Especificação de modelos utiliza terminologia de fórmula (y ~ x)

fabletools

  • Extensões e ferramentas para construção de modelos
  • Combinação de modelos, previsão hierárquica e extração de resultados
  • Obtenção de medidas de acurácia e visualização de dados

feasts

  • Decomposição de séries temporais
  • Extração e visualização de componentes de séries temporais
  • Análise de autocorrelação, testes de raiz unitária, etc.

tsibble

  • Estrutura de dados tidy para séries temporais no R
  • Funções para tratamento de dados
  • Objeto orientado aos dados e a modelos, integrado ao tidyverse

Fluxo de trabalho para previsão

Com o tidyverts o processo de construir um modelo de previsão pode ser dividido em poucos passos:

Com esse esquema em mente, vamos ilustrar esse processo com um exercício prático e didático: construir um modelo de previsão para a taxa de crescimento do PIB brasileiro.

Pacotes

Para reproduzir o exercício a seguir você precisará dos seguintes pacotes:

Dados tidy

Utilizaremos o dataset global_economy armazenado como um objeto tsibble, trazendo variáveis econômicas em frequência anual para diversos países. Nosso interesse é a série da taxa de crescimento do PIB brasileiro:

Visualização de dados

Visualização é uma etapa essencial para entender os dados, o que permite identificar padrões e modelos apropriados. No nosso exemplo, criamos um gráfico de linha para plotar a série do PIB brasileiro usando a função autoplot():

Podemos também plotar os correlogramas ACF e PACF para identificar o processo estocástico da série, obtendo alguns modelos candidatos:

Especificação do modelo

Existem muitos modelos de séries temporais diferentes que podem ser usados para previsão, e especificar um modelo apropriado para os dados é essencial para produzir previsões.

Os modelos no framework do fable são especificados usando funções com nomenclatura abreviada do nome do modelo (por exemplo, ARIMA(), AR(), VAR(), etc.), cada uma usando uma interface de fórmula (y ~ x). As variáveis de resposta são especificadas à esquerda da fórmula e a estrutura do modelo é escrita à direita.

Por exemplo, um modelo ARIMA(1,0,2) para a taxa de crescimento do PIB pode ser especificado com: ARIMA(Growth ~ pdq(1, 0, 2)).

Neste caso, a variável resposta é Growth e está sendo modelada usando a estrutura de um modelo ARMA(1, 2) especificada na função especial pdq().

Existem diversas funções especiais para definir a estrutura do modelo e em ambos os lados da fórmula pode ser aplicado transformações. Consulte detalhes da documentação do fable.

Estimar o modelo

Identificado um modelo (ou mais) apropriado, podemos em seguida fazer a estimação usando a função model()1.Neste exemplo, estimaremos os seguintes modelos: ARIMA(1,0,2), ARIMA(1,0,0), ARIMA(0,0,2), o algoritmo de seleção automatizada do auto ARIMA criado pelo prof. Rob Hyndman e um passeio aleatório.

Diagnóstico do modelo

O objeto resultante é uma "tabela de modelo" ou mable, com a saída de cada modelo em cada coluna:

Para obter os critérios de informação use a função glance():

Os critérios de informação indicam que, dos modelos estimados, o modelo automatizado ARIMA(1,1,1) apresentou o menor valor de AICc - seguido pelos demais identificados pelos correlogramas ACF e PACF. Com a função gg_tsresiduals() podemos verificar o comportamento dos resíduos deste modelo, indicando que os resíduos se comportam como ruído branco:

Um teste de autocorrelação (Ljung Box) retorna um p-valor grande, também indicando que os resíduos são ruído branco:

Também pode ser interessante visualizar o ajuste do modelo. Utilize a função augment() para obter os valores estimados:

Previsão

Com o modelo escolhido, previsões podem ser geradas com a função forecast() indicando um horizonte de escolha.

Perceba que os pontos de previsão médios gerados são bastante similares a um processo de passeio aleatório (equivalente a um ARIMA(0,1,0)). O trabalho adicional de especificar termos AR e MA trouxe pouca diferença para os pontos de previsão neste exemplo, apesar de ser perceptível que os intervalos de confiança do modelo auto ARIMA são mais estreitos do que de um passeio aleatório.

Além disso, a previsão fora da amostra gerada ficou bastante aquém dos dados reais para a taxa de crescimento do PIB brasileiro observados no horizonte em questão, configurando apenas um exercício didático.

Saiba mais

Estes são apenas alguns dos recursos e ferramentas disponíveis na família de pacotes do tidyverts. Para uma referência aprofundada, confira o livro Forecasting: Principles and Practice, 3rd Edition, de Hyndman e Athanasopoulos (2021).

Aprenda com maior profundidade sobre sobre os temas abordados nos cursos de Séries Temporais e de Previsão Macro. Faça parte do Clube AM para acesso completo aos códigos de R e Python deste e de outros exercícios.

Confira outros exercícios aplicados com pacotes do tidyverts:

 


[1] A função suporta estimação dos modelos com computação paralela usando o pacote future, veja detalhes na documentação e este post para saber mais sobre o tema.

 

 

Modelagem e previsão com o {modeltime}

By | Data Science

O pacote {modeltime} é um framework para modelagem e previsão de séries temporais no R, com a possibilidade de integração com o ecossistema do {tidymodels} para uso de técnicas de machine learning. O pacote promete ser uma interface veloz e moderna que possibilita prever séries temporais em grande escala (i.e. 10+ mil séries).

Utilizando uma interface simples e integrada com outros pacotes populares na linguagem, é possível construir uma ampla gama de modelos de previsão univariados e multivariados, tradicionais ou "modernos", dentre eles:

  • ARIMA
  • Suavização exponencial (ETS)
  • Regressão linear (TSLM)
  • Autoregressão de rede neural (NNETAR)
  • Algoritmo Prophet
  • Support Vector Machines (SVM)
  • Random Forest
  • Boosted trees

Neste exercício, daremos uma breve introdução aplicada exemplificando o fluxo de trabalho do {modeltime} para modelagem e previsão, além de fornecer nossa opinião sobre o mesmo. Para se aprofundar nos detalhes de cada tópico confira os cursos de Séries Temporais e Modelos Preditivos da Análise Macro.

Fluxo de trabalho

O fluxo de trabalho do pacote {modeltime} para criar um modelo de previsão de séries temporais pode ser resumido em 5 etapas:

  1. Colete e prepare os dados no ambiente do R;
  2. Separe amostras de treino/teste para avaliar modelos;
  3. Crie e estime múltiplos modelos de previsão;
  4. Avalie os modelos (resíduos, acurácia no teste, etc.);
  5. Reestime os modelos na amostra completa e gere previsões.

Com esse fluxo em mente, vamos agora realizar um exercício prático e didático: construir um modelo de previsão para a taxa de inflação brasileira, medida pelo IPCA (% a.m.) do IBGE. O objetivo é apenas demonstrar as ferramentas ( {modeltime}) de maneira prática no R, não estamos preocupados em criar o melhor modelo possível e nem em avaliar especificações profundamente. Para tal dê uma olhada nos cursos da Análise Macro.

Para reproduzir o exercício você precisará dos seguintes pacotes de R:


Passo 1: dados

Nesse exercício vamos criar modelos de previsão para o IPCA (% a.m.), portanto começamos por coletar e tratar os dados de maneira online diretamente do Sidra/IBGE (veja um tutorial neste link). Nossa amostra parte de 1999 e vai até o último dado disponível.

Abaixo criamos um gráfico de linha simples da série temporal do IPCA (% a.m.), do qual podemos observar uma flutuação em torno de 0,5% e um padrão sazonal, além de outliers (veja mais sobre análise exploratória neste link):

Passo 2: separar amostras

Com o objetivo de avaliar os modelos que criaremos em termos de acurácia da previsão pontual, vamos criar duas amostras da série temporal do IPCA: a de treino é usada para estimar um modelo especificado e gerar previsões com o mesmo e a de teste é usada para comparar as previsões geradas com os dados observados (verdadeiros). Nesse exercício vamos fazer essa separação em 95% para treino e o restante para teste:

Passo 3: modelagem

Como exemplo, vamos criar dois modelos univariados simples: família (auto) ARIMA e um passeio aleatório como benchmark. No {modeltime} a modelagem é feita em 3 etapas:

  1. Especificar o modelo (usando funções como arima_reg() e outras);
  2. Definir o pacote a ser usado (alguns modelos são implementados em mais de um pacote, e você tem flexibilidade em escolher qual usar com a função set_engine());
  3. Estimar o modelo (usando a função genérica fit(y ~ date, data = df), onde date não é um regressor, é apenas para uso interno do {modeltime}).

Note que, por esse fluxo de trabalho, o {modeltime} traz a vantagem de integrar em um único lugar diversos outros pacotes/modelos (consulte a documentação).

Os modelos ficam salvos como objetos do tipo lista no R e ao imprimi-los temos os principais resultados:

Passo 4: avaliar modelos

Uma vez que os modelos tenham sido estimados, vamos então investigar brevemente suas especificações: focaremos apenas em testar a normalidade e autocorrelação dos resíduos dos modelos estimados e em calcular as métricas MAE e RMSE dos pontos de previsão pseudo fora da amostra (comparando com a amostra de teste).

Para tais tarefas de diagnóstico de modelos o {modeltime} disponibiliza algumas funções úteis:

  • modeltime_table() adiciona um ou mais modelos para criar uma "tabela de modelos";
  • modeltime_calibrate() calcula os resíduos dos modelos;
  • modeltime_residuals_test() aplica testes de normalidade/autocorrelação;
  • modeltime_forecast() gera previsões;
  • plot_modeltime_forecast() visualiza as previsões;
  • modeltime_accuracy() calcula métricas de acurácia.

Vamos começar aplicando os testes de normalidade (Shapiro-Wilk) e autocorrelação (Box-Pierce e Ljung-Box) dos resíduos:

O resultado é uma tabela com os p-valores de cada teste, que indicam que os resíduos do modelo auto ARIMA, ao nível de 5% de significância, não são normalmente distribuídos e não apresentam autocorrelação. Para dados econômicos esses resultados são, em grande parte, comuns, de modo que iremos ignorar aqui o problema de normalidade e seguir adiante.

Agora vamos verificar a performance preditiva dos modelos nesse esquema pseudo fora da amostra. Para tal, primeiro juntamos os modelos em uma tabela de modelos e geramos n previsões fora da amostra, onde n é o número de observações da amostra de teste. Abaixo plotamos o gráfico das previsões geradas:

Pela avaliação visual os modelos estão muito aquém do desejado, mas esse exercício é apenas didático e, conforme mencionado, não estamos interessados em criar um "super" modelo que desempenhe melhor que as expectativas do Focus/BCB.

Em seguida, usamos a tabela de modelos para comparar os pontos de previsões gerados com os valores observados da amostra de teste, computado as métricas MAE e RMSE:

Como esperado, o modelo ARIMA apresenta um RMSE menor em relação ao benchmark, para o horizonte e amostra utilizados. Portanto, se fossemos decidir um modelo final entre os dois, a escolha seria pelo ARIMA. Note que em exercícios mais robustos o mais interessante é que você faça essa avaliação usando validação cruzada por horizontes de previsão (1, 2, 3, ..., meses a frente).

Passo 5: previsão

Com os modelos treinados e avaliados, o último passo é reestimar o(s) modelo(s) usando a amostra completa de dados para gerar previsões (efetivamente) fora da amostra. Use a função modeltime_refit() para tal:

E assim terminamos o fluxo de trabalho para modelagem e previsão com o {modeltime}. O que você achou do pacote?

Comentários

Como qualquer pacote/software, o {modeltime} tem suas vantagens e desvantagens, aqui vou pontuar a minha visão.

  • Principais vantagens:
    • Interface simples para iniciantes;
    • Integra com o {tidymodels};
    • Ecossistema ativo de desenvolvimento.
  • Principais desvantagens:
    • Interface fora do comum dos pacotes econométricos (i.e. fit(y ~ date, data = df));
    • Documentação orientada a negócios/sem background estatístico;
    • Tem cerca de 3x mais dependências (pacotes) do que concorrentes como o {fabletools};
    • A comunidade de usuários é pequena (quando comparado pelo nº de perguntas no StackOverflow).

Em suma, se você é um iniciante o pacote {modeltime} deve ser ótimo para dar os primeiros passos no mundo de séries temporais e modelos de previsão, mas se você perguntar a uma pessoa mais experiente se esse é o pacote que ela prefere, provavelmente a resposta será não. Conforme você avança na área haverá a necessidade de ir além e essa lacuna provavelmente será preenchida por outro pacote. Contudo, se sua intenção não é mergulhar profundamente em séries temporais, basta instalar o pacote para começar usá-lo!

Saiba mais

Caso o assunto tenha despertado sua curiosidade, você pode se aprofundar nos tópicos com os cursos de Séries Temporais e Modelos Preditivos da Análise Macro.

Como exemplo, no curso de Modelos Preditivos utilizamos modelos de machine learning para a previsão do IPCA. Os resultados são sumarizados nessa dashboard (clique para acessar), em uma abordagem automatizada:

Para ter acesso aos códigos completos deste exercício e vários outros, faça parte do Clube AM.

 

 

Modelos ARIMA na abordagem bayesiana

By | Data Science

A boa e velha econometria nos proporciona métodos que são ótimos para relacionar variáveis e construir previsões, sendo quase sempre o primeiro ponto de partida para praticantes, além de ser a base de grande parte dos modelos "famosos" da atualidade. Como exemplo, em problemas de previsão de séries temporais a família de modelos ARIMA, popularizada por Box e Jenkins (1970), costuma ser a primeira opção que vem à mente.

Uma das suposições destes modelos é a de que os erros — ou seja, a diferença entre o que é observado e o que é estimado — seguem uma distribuição normal, o que no mundo real nem sempre é factível com os dados que temos em mãos. E então, o que fazer nessa situação? Devemos descartar os modelos ARIMA e usar outra classe de modelos? Ou devemos descartar os dados, pois, claramente, "eles não se encaixam no modelo"?

Bom, nada disso! Em primeiro lugar, eu afirmaria que os modelos ARIMA ainda podem ser úteis. Em segundo lugar, neste exercício apresentarei uma visão diferente no que se refere ao método de estimação do modelo, trazendo uma breve introdução à abordagem bayesiana. Por fim, exemplificarei com um modelo ARIMA simples as diferenças entre previsões geradas pelas duas visões, isto é, a frequentista versus a bayesiana.

O foco aqui é ser objetivo, apresentando brevemente os conceitos básicos, as vantagens e diferenças entre as abordagens, em uma tentativa bem humilde de fechar o gap entre as duas visões em termos introdutórios. Para mais detalhes e aprofundamento nos tópicos confira os cursos de Séries Temporais e Estatística Bayesiana.

ARIMA

Para não complicar, vamos considerar um modelo autoregressivo (AR) simples neste exercício. Como este modelo se parece?

yt = c + ϕ1yt-1 + ϕ2 yt-2 + + ϕpyt-p + εt

onde:

c → uma constante;

ϕ1, ⋯, ϕp → parâmetros desconhecidos do modelo;

εt → termo de erro que segue uma distribuição normal, εt N(μ, σ2).

O modelo é autoregressivo pois, como pode ser observado, se trata de uma regressão da variável contra ela mesma (defasagens). Esse modelo é comumente denominado com o acrônimo AR(p), onde p é a ordem de defasagens da variável de interesse yt e é um caso especial de um modelo ARIMA(p, 0, 0)1.

Ok, mas visualmente, como um processo AR(p) se parece?

Podemos visualizar uma série temporal que segue um processo, por exemplo, de AR(1) com uma simulação simples (estou usando o R para isso), conforme abaixo:

Estimador MQO

Você deve ter notado que valores para a constante e para o primeiro parâmetro do AR(1)e ϕ1 , "caíram do céu". Na vida real, fora da simulação, devemos estimar esses valores, eles são desconhecidos. Portanto, como os parâmetros do modelo são estimados?

Existem algumas maneiras de estimar os parâmetros do modelo, o padrão comum é através do estimador de mínimos quadrados ordinários (MQO), onde o objetivo é encontrar valores ótimos para os parâmetros e ϕp do modelo minimizando a soma dos quadrados dos resíduos, sob diversas premissas. Ao usar essa abordagem, os parâmetros do modelo autoregressivo são fixos: cada observação t da série temporal y é associada com um mesmo valor estimado do parâmetro ϕp.

Isso é factível? E se você acredita que o valor ϕp não deveria ser fixo? Como poderíamos incorporar essa crença ao modelo? É aqui que entra a abordagem bayesiana! Mas antes, vamos dar uma olhada em dados reais para tentar entender a importância e uma motivação adjacente do assunto.

Exemplo: IPCA

Como exemplo, abaixo apresento visualizações gráficas da série do IPCA sazonalmente ajustado, uma medida de inflação divulgada pelo IBGE, expresso em variação % ao mês.

Note que pela distribuição dos dados o IPCA apresenta uma média de aproximadamente 0,5% com uma cauda pesada à direita. Isso é característico do que chamam de outliers ou valores extremos e pode ser um problema em métodos econométricos. Além disso, a série apresenta uma autocorrelação significativa na primeira defasagem, o que pode ser um sinal para modelagem ARIMA, se supormos que tal classe de modelos univariados é suficiente em explicar a inflação.

Dada essa simples inspeção visual dos dados o que frequentistas e bayesianos poderiam tirar de conclusão? Um frequentista diria, dentre outras coisas, que a variância da série é a mesma para qualquer janela amostral, enquanto que um bayesiano diria que a variância pode variar.

Nessa mesma linha, um frequentista poderia estimar um modelo autoregressivo por MQO para explicar a inflação, mas ao fazê-lo estaria encontrando um valor fixo associado a defasagens da inflação, o que é equivalente a dizer que a maldita inércia inflacionária — um fenômeno econômico — não se alterou em todos esses mais de 20 anos de Plano Real. Isso é plausível? Certamente não é uma pergunta trivial de ser respondida, mas exercícios prévios do colega Vitor Wilher sugerem que a inércia inflacionária vem mudando recentemente.

Ou seja, são nesses pontos, para citar apenas alguns, que as abordagens divergem. Em outras palavras, frequentistas e bayesianos fazem diferentes suposições sobre os dados e sobre os parâmetros de um modelo em questão.

Portanto, voltando à questão, como a modelagem bayesiana incorpora sua visão sobre os parâmetros ao modelo?

A abordagem bayesiana

Na abordagem bayesiana, ao invés de tentarmos encontrar um valor ótimo para os parâmetros do modelo, nesse caso um AR(p), começaríamos por assumir um distribuição de probabilidade a priori para os mesmos e, então, usaríamos o teorema de Bayes. Ou seja, em nosso exemplo o que muda é que você assume uma distribuição a priori para os valores desconhecidos de c, ϕ1,,ϕp, σ2, antes mesmo dos dados serem observados.

Mas de onde essa distribuição a priori vem? Bom, há muito o que se discutir sobre esse ponto, mas em resumo:

  • Estudos prévios ou trabalhos publicados;
  • Intuição do analista/pesquisador;
  • Conhecimento da área/expertise;
  • Conveniência;
  • Métodos não paramétricos e outras fontes de dados.

Você pode pensar nessa distribuição a priori como uma crença/conhecimento prévio ou, no caso de você possuir estimativas prévias dos parâmetros do modelo, essas estimativas prévias se tornam a distribuição a priori.

Voltando ao exemplo do IPCA, supondo que queiramos estimar um modelo AR(1), yt = c + ϕ1yt-1 + εt, como definiríamos distribuições a priori para os parâmetros? Vamos começar admitindo que não sabemos muito sobre os parâmetros, apenas sabemos que para um AR(1) ser estacionário o parâmetro ϕ1 precisar ser menor do que 1 em módulo2, então vamos assumir distribuições genéricas:

Como é um primeiro modelo, vamos ser bem vagos sobre as escolhas. Partindo de um entendimento básico sobre o modelo, sabemos que o parâmetro  para caracterizar um processo estacionário e que a variância dos erros deve ser positiva. Dessa forma utilizaremos, respectivamente, a distribuição uniforme e a gamma inversa e, quando não, a normal, assumindo os valores abaixo:

Essas são as nossas priors, como é referido no jargão bayesiano. Visualmente, como essas distribuições se parecem? Vamos plotar um gráfico (novamente, estou usando o R aqui):

A interpretação destas distribuições a priori é de que não estamos informando em nosso modelo AR(1) bayesiano sobre quais poderiam ser os valores dos parâmetros. A estimativa final, pelo teorema de Bayes, dependerá dos dados observados e das priors assumidas. Vamos definir isso brevemente:

O que essa equação nos diz, em termos simples? A priori é a probabilidade de algo acontecer antes de incluirmos a probabilidade dos dados (verossimilhança), e a posteriori é a probabilidade após a incorporação dos dados. Ou seja, o teorema de Bayes nos fornece um método para atualizar nossa crença sobre parâmetros desconhecidos para cada informação nova (observações), condicionando a probabilidade aos dados observados.

Para estimar um modelo AR(p) na abordagem bayesiana será utilizado, geralmente, um método chamado Monte Carlo Markov Chain (MCMC). Apesar do nome apavorante, é um procedimento relativamente simples. O objetivo é "caminhar" aleatoriamente por um espaço de parâmetros para encontrar as estimativas e, assim, desenhar uma estimativa da posteriori, a nossa distribuição de principal interesse.

Comparação

Agora vamos comparar as estimativas de um modelo AR(1) para o IPCA geradas pela visão frequentista vs. bayesiana3. Indo direto ao ponto, abaixo plotamos os resultados:

Aqui eu omito a distribuição a posteriori, que representaria a incerteza da estimativa bayesiana, exibindo apenas a a estimativa média dessa distribuição.

Note que as estimativas pontuais dentro da amostra são praticamente idênticas (eu separei em dois gráficos para melhor visualizar cada linha), refletindo semelhanças nos parâmetros estimados da abordagem bayesiana e da frequentista. O valor estimado para a constante é aproximadamente 0,20 e 0,19, respectivamente, enquanto que o é 0,6 e 0,6. Por quê? A razão é decorrente do fato que a estimativa pontual do método frequentista, por MQO, é, em verdade, a mesma coisa que maximizar a probabilidade a posteriori de uma regressão linear bayesiana, que no caso representa a moda da distribuição.

Ok, mas comparamos somente estimativas pontuais um passo a frente dentro da amostra. O que acontece se compararmos ambas as abordagens em termos de previsão fora da amostra? No gráfico acima as últimas 12 observações do IPCA foram omitidas dos modelos, justamente para ter um período pseudo fora da amostra que nos permita gerar previsões e compará-las com os dados observados, podendo dessa forma avaliar a performance preditiva de cada modelo para dados que os mesmos ainda não conhecem.

Gerando previsões pseudo fora da amostra 12 períodos (meses) a frente, temos os seguintes resultados:

Comentários

Apesar de ambas as previsões serem visivelmente ruins — acredito eu que não devido ao método de estimação em si, mas sim por conta da especificação/modelo simples para os dados e pela janela amostral (choques) "desfavorável" —, podemos notar uma grande diferença nas previsões pontuais e nos intervalos. A abordagem frequentista simplesmente previu o retorno a média da série (o que nesse caso é menos pior), enquanto que a abordagem bayesiana previu uma deflação (o que não é a realidade brasileira no momento).

Pela simplicidade da especificação do modelo AR em questão não podemos esperar bons resultados, por qualquer método de estimação considerado, e uma comparação mais exigente e justa implementaria a técnica de validação cruzada, o que eu não fiz neste exercício. Além disso, questões adjacentes como a inspeção dos resíduos e dos parâmetros são importantes para diagnosticar o modelo, as quais também deixamos de lado aqui.

Em suma, estimamos um modelo AR(1) para o IPCA por dois métodos diferentes que fundamentalmente fazem a mesma coisa: eles usam a defasagem da inflação para explicá-la. A primeira abordagem foi um modelo autoregressivo clássico, usando o método dos mínimos quadrados ordinários (MQO). O segundo foi um mesmo modelo autoregressivo, mas com estimação bayesiana, explorando a ideia de usar distribuições a priori para obter estimativas pontuais de distribuições a posteriori. Apesar disso, os resultados não foram satisfatórios. Ou seja, há trabalho ainda a se fazer, além de outros modelos a se testar.

Se o assunto te interessou e se você quiser saber mais sobre os temas abordados, dê uma olhada nos cursos de Séries Temporais e Estatística Bayesiana. Para ter acesso aos códigos deste exercício faça parte do Clube AM.

 


1,2 Veja mais em Hyndman e Athanasopoulos (2021).

3 Aqui eu relaxei a restrição de distribuição uniforme para o parâmetro  por recomendação dos desenvolvedores do pacote utilizado e por falha de inicialização da estimação. A priori desse parâmetro passou então para .

 

Referências

Box, G. E. P., & Jenkins, G. M. (1970). Time series analysis: Forecasting and control. Holden-Day.

Durbin, J., & Koopman, S. J. (2012). Time series analysis by state space methods (Vol. 38). OUP Oxford.

Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3. Accessed on 2022-07-20.

Modelagem e Previsão em 6 Passos: fluxo de trabalho com séries temporais

By | Data Science

Produzir modelos preditivos para séries temporais é uma tarefa difícil e bastante valorizada no mercado de trabalho. São diversos os procedimentos que um profissional dessa área precisa sempre ter no radar e, se generalizarmos isso em um fluxo, a distância entre a idealização do modelo e sua implementação pode ser mais curta e menos árdua. Neste texto descrevemos resumidamente o dia a dia de um profissional que trabalha com previsões e mostramos um exemplo de rotinas em R!

Fluxo de trabalho

O processo de desenvolver modelos preditivos para séries temporais pode ser generalizado e dividido em 6 etapas:

Fonte: Hyndman and Athanasopoulos (2021)

Este fluxo, em resumo, descreve como é o dia a dia de trabalho de modelagem e previsão em algumas etapas organizadas. São elas:

Preparar os dados

O primeiro passo é preparar os dados no formato correto para uso em modelos. Esse processo pode envolver a coleta e importação de dados, a identificação de valores ausentes, agregação/sumarização/ajustes/filtros nas séries e outras tarefas de pré-processamento (veja um post sobre coleta de dados aqui). O objetivo será, em geral, construir uma matriz de variáveis onde cada variável terá sua própria coluna, cada linha formará uma observação e os valores são armazenados nas células. Esse é o princípio de dados "tidy". A facilidade que o pacote {tsibble} e outros do {tidyverse} oferecem simplifica consideravelmente esta etapa usando o R.

Vale dizer que muitos modelos têm pressupostos e requisitos diferentes, sendo que você precisará levar isso em consideração ao preparar os dados. Alguns exigem que a série seja estacionária, outros exigem que não haja valores ausentes. Dessa forma, você precisará conhecer bem os seus dados enquanto os prepara e a análise exploratória é outra etapa que caminhará lado a lado.

Visualizar os dados

A visualização é uma etapa essencial para a compreensão dos dados. Suas séries temporais apresentam tendência? Possuem sazonalidade? Há quebras ou observações extremas nas séries (outliers)? Observar seus dados, através de uma análise exploratória (veja um post sobre aqui no blog) permite identificar estes padrões comuns e, posteriormente, especificar um modelo apropriado. Essa etapa pode andar em conjunto com a etapa de preparação dos dados, de modo que após entender os dados você talvez precise voltar a um passo anterior e aplicar uma transformação, por exemplo, nas séries. Os pacotes {ggplot2}, {feasts}, {fabletools} e outros oferecem formas elegantes e práticas de visualizar seus dados.

Definir o modelo

Existem muitos modelos de séries temporais diferentes que podem ser usados para previsão. Especificar um modelo apropriado para os dados é essencial para produzir previsões. Nesse sentido, além de conhecer os seus dados você precisará ter conhecimento sobre os modelos que pretende trabalhar. Papers publicados são, em geral, boas fontes para buscar esse conhecimento técnico.

Grande parte dos modelos econométricos é descrita em Hyndman and Athanasopoulos (2021) e podem ser implementados pelo pacote {fable} através de funções específicas como ARIMA(), TSLM(), VAR(), etc., cada uma usando uma interface de fórmula (y ~ x). As variáveis de resposta são especificadas à esquerda da fórmula e a estrutura do modelo, que pode variar, é escrita à direita. Para outros modelos fora do escopo do {fable} a lógica e sintaxe é semelhante, quase sempre utilizando uma interface de fórmula.

Estimar o modelo

Uma vez que um modelo apropriado é definido, passamos à estimação do modelo com os dados. Uma ou mais especificações de modelo podem ser estimadas usando a função model() do pacote {fabletools} — framework que contempla grande parte dos modelos econométricos; para outros modelos há pacotes e funções específicas, mas com sintaxe semelhante.

Avaliar a performance

Após termos um modelo estimado, é importante verificar o desempenho dele nos dados. Existem várias ferramentas de diagnóstico disponíveis para verificar o ajuste do modelo, assim como medidas de acurácia que permitem comparar um modelo com outro; o RMSE é a métrica mais comumente utilizada para a maioria dos problemas de previsão. Conforme o diagnóstico do modelo, em seguida possivelmente sejam necessárias readequações, seja na especificação ou até mesmo nos dados utilizados. Em outras palavras, o fluxo de trabalho não é simplesmente um amontoado de procedimentos a serem implementados sequencialmente, mas sim um processo de descobrimento que envolve, na vida real, sucessivas tentativas e erros.

Além disso, se o interesse é previsão, há técnicas como a validação cruzada que auxiliam na tomada de decisão entre mais de um modelo (veja um post sobre isso aqui no blog). É sempre preferível ter mais de um modelo "candidato" potencialmente usado para fazer previsões, além de modelos básicos para simples comparação.

Realizar previsões

Com um modelo especificado, estimado e diagnosticado, é hora de produzir as previsões fora da amostra. Para alguns modelos você poderá simplesmente chamar uma função, como a forecast() ou predict() no R, especificando o número de períodos (horizonte de previsão) que deseja obter previsões; para outros você precisará prover uma tabela com os valores futuros das variáveis regressoras utilizadas no modelo, que servirá para produzir as previsões da variável de interesse, ou seja, você precisará de cenários.

Exemplo no R

Com esse esquema em mente, vamos ilustrar o fluxo de trabalho com um exercício prático e didático: construir um modelo de previsão para a taxa de crescimento anual do PIB brasileiro.

Para reproduzir o exercício a seguir você precisará dos seguintes pacotes:

Preparar os dados

Por conveniência, utilizaremos o dataset global_economy armazenado como um objeto tsibble, trazendo variáveis econômicas em frequência anual para diversos países. Nosso interesse é a série da taxa de crescimento do PIB brasileiro:

Visualizar os dados

Visualização é uma etapa essencial para entender os dados, o que permite identificar padrões e modelos apropriados. No nosso exemplo, primeiro criamos um gráfico de linha para plotar a série do PIB brasileiro usando a função autoplot():

Podemos também plotar os correlogramas ACF e PACF para identificar o processo estocástico da série, obtendo alguns modelos candidatos:

Definir o modelo

Existem muitos modelos de séries temporais diferentes que podem ser usados para previsão, e especificar um modelo apropriado para os dados é essencial para produzir previsões. Nesse exercício didático focaremos em modelos univariados simples para explicar o PIB brasileiro.

Se formos analisar somente pelos correlogramas ACF e PACF, podemos identificar alguns modelos possivelmente candidatos: ARIMA(1,0,2), ARIMA(1,0,0) e ARIMA(0,0,2). Contudo, se nossa leitura estiver incorreta, podemos contar ainda com a tecnologia ao nosso favor utilizando a seleção automatizada da estrutura do modelo. Isso é possível graças ao algoritmo de seleção automatizada da especificação do ARIMA criado pelo prof. Rob Hyndman.

Por fim, além destes possíveis modelos candidatos, também estimaremos um modelo de benchmark na forma de um passeio aleatório, apenas para ter uma base de comparação de previsões dos modelos.

Estimar o modelo

Identificado um modelo (ou mais) apropriado, podemos em seguida fazer a estimação usando a função model()1. Existem diversas funções especiais para definir a estrutura do modelo e em ambos os lados da fórmula podem ser aplicadas transformações. Nesse caso definiremos apenas a estrutura básica do ARIMA; a função ARIMA() também define automaticamente a estrutura sazonal, mas você pode desabilitar isso (consulte detalhes da documentação do {fable}).

O objeto resultante é uma "tabela de modelo" ou mable, com a saída de cada modelo em cada coluna.

Diagnóstico do modelo

Para obter os critérios de informação use a função glance():

Os critérios de informação indicam que, dos modelos estimados, o modelo automatizado ARIMA(1,1,1) apresentou o menor valor de AICc — seguido pelos demais identificados pelos correlogramas ACF e PACF. Com a função gg_tsresiduals() podemos verificar o comportamento dos resíduos deste modelo, indicando que os resíduos se comportam como ruído branco:

Um teste de autocorrelação (Ljung Box) retorna um p-valor grande, também indicando que os resíduos são ruído branco:

Também pode ser interessante visualizar o ajuste do modelo. Utilize a função augment() para obter os valores estimados:

Outros diagnósticos e avaliações podem ser realizados. Se você estiver especialmente interessado em previsão, considere implementar o método de validação cruzada, como explicado neste post.

Realizar previsões

Com o modelo escolhido, previsões podem ser geradas com a função forecast() indicando um horizonte de escolha.

Perceba que os pontos de previsão médios gerados são bastante similares a um processo de passeio aleatório (equivalente a um ARIMA(0,1,0)). O trabalho adicional de especificar termos AR e MA trouxe pouca diferença para os pontos de previsão neste exemplo, apesar de ser perceptível que os intervalos de confiança do modelo auto ARIMA são mais estreitos do que de um passeio aleatório.

Além disso, a previsão fora da amostra gerada ficou bastante aquém dos dados reais para a taxa de crescimento do PIB brasileiro observados no horizonte em questão, configurando apenas um exercício didático.

Saiba mais

Este é apenas um exercício simples e didático demonstrando um fluxo de trabalho de modelagem e previsão. Existem diversos tópicos que podem ser aprofundados e se você quiser saber mais confira o curso de Modelos Preditivos da Análise Macro.

Confira também outros exercícios aplicados com pacotes do {tidyverts}:

Referências

Hyndman, R. J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3. Accessed on 2022-04-01.

 


1A função suporta estimação dos modelos com computação paralela usando o pacote {future}, veja detalhes na documentação e este post para saber mais sobre o tema.

Validação Cruzada de Séries Temporais

By | Data Science

Quando estimamos um modelo para previsão de séries temporais precisamos avaliar sua performance preditiva. Separar os dados em uma amostra de treino e outra de teste é a forma mais simples de verificar como o modelo performa em dados pseudo fora da amostra, através do cálculo do erro de previsão. Entretanto, a escolha das amostras é arbitrária, o que traz a tona um método mais sofisticado: a validação cruzada. Neste exercício abriremos a caixa preta por trás do assunto mostrando como implementar a técnica no R.

Neste método há uma série de amostras de teste, cada uma consistindo em h observações, ou seja, são os períodos usados para gerar previsões a partir do modelo. A amostra de treino correspondente, usada para estimação do modelo, consiste apenas de observações que ocorreram antes das observações que formam a amostra de teste. Assim, nenhuma observação futura pode ser usada na construção da previsão.

O diagrama a seguir ilustra um esquema de validação cruzada de séries temporais, com uma série de amostras de treino e de teste, para quando se deseja avaliar a performance preditiva do modelo em h = 1 períodos à frente, onde as observações azuis formam as amostras de treino e as observações laranja formam as amostras de teste.

Fonte: Hyndman e Athanasopoulos (2021).

Operacionalização

Dado esse esquema de separação de uma série de amostras de treino e de teste de uma série temporal, o método de validação cruzada pode ser resumido como um procedimento recursivo:

  1. Para cada amostra de treino estime o modelo;
  2. Utilize o modelo estimado para gerar previsões h períodos à frente;
  3. Calcule o erro de previsão com a amostra de teste e a saída da etapa 2.

Ao final, obtém-se a métrica de acurácia (ME, RMSE, MAE, etc.) do modelo pela média de cada iteração.

Esse método é relativamente simples, mas bastante útil para validar a performance de modelos independentemente da amostra de treino/teste escolhida, trazendo uma "visão global" do modelo. A validação cruzada é largamente utilizada por praticantes de machine learning, além de ser uma prática comum em papers que tratam de previsão de séries temporais. Vale pontuar que existem variações do método, que deixaremos para explorar em uma outra oportunidade.

Implementação no R

Mostraremos duas formas de implementar validação cruzada de séries temporais no R:

  1. Através da família de pacotes {tidyverts};
  2. Abrindo a caixa preta e escrevendo o código manualmente.

Em ambos os casos utilizaremos modelos univariados simples, apenas para demonstrar o método, aplicados ao principal índice de preços da economia brasileira, o IPCA, com a finalidade de avaliar as previsões 12 períodos à frente. A primeira amostra de validação cruzada conterá 180 observações e será acrescentado 1 observação sucessivamente até o total de observações excluindo as últimas 12.

Antes de começar você deve ter disponível os seguintes pacotes no seu ambiente de R:

O código abaixo prepara os dados de exemplo. Para saber mais sobre coleta de dados macroeconômicos confira este post no blog da Análise Macro.

Exemplo com {tidyverts}

O primeiro exemplo de implementação de validação cruzada utilizará os pacotes da família {tidyverts}, que é um ótimo framework para trabalhar com séries temporais no R (confira este post sobre o assunto). Sendo assim, o primeiro passo que precisamos fazer é criar as amostras de treino de validação cruzada (observações azuis na ilustração). Isso pode ser feito conforme abaixo (note que removemos as últimas 12 observações, para poder calcular o erro de previsão da última amostra corretamente):

Agora temos uma série de amostras, partindo de 180 observações e acrescentando 1 sucessivamente, que utilizaremos para estimar os modelos de exemplo e gerar previsões. Isso é feito conforme o código abaixo. Note que criamos uma coluna para salvar a informação do horizonte temporal da previsão (isso é opcional, sendo útil apenas quando você quer analisar a acurácia por horizonte preditivo).

Por fim, podemos calcular as métricas de acurácia média considerando todas as amostras e previsões geradas, conforme abaixo:

Como resultado podemos gerar a visualização abaixo, onde é possível analisar os modelos em termos de acurácia (RMSE) por horizonte preditivo:

Exemplo manual

Conforme visto, os procedimentos adotados para implementar o método de validação cruzada de séries temporais são extremamente simples utilizando os pacotes do {tidyverts}. Entretanto, nem sempre teremos à disposição, dentro deste ou outros frameworks, todos os modelos de interesse para implementação rápida através de um pacote. Na vida real, é comum que surjam novos modelos ou que os pacotes implementem e disponibilizem apenas uma parcela do que existe por aí. Nestes casos, será seu trabalho escrever um código que implemente a validação cruzada com tal modelo.

Para exemplificar, abaixo escrevo uma função que implementa o esquema de validação cruzada de forma semelhante ao que acabamos de ver. O modelo utilizado é, por conveniência, o ARIMA automatizado e a definições de amostras são idênticas ao exemplo anterior. Note que o código fica significativamente mais complexo, principalmente para iniciantes, mas é um destino natural de todo profissional que deseje se aventurar por estes caminhos.

A função implementa o que foi destacado previamente, referente ao método recursivo da validação cruzada, retornando uma lista com uma tabela do RMSE por horizonte preditivo e outra com as previsões geradas a cada iteração. Fique à vontade para adaptar ou melhorar a função conforme suas necessidades.

Com a função criada, basta utilizá-la junto aos dados para obter os resultados:

Por fim, plotamos os resultados desta abordagem manual com a abordagem anterior. Note que o RMSE calculado difere levemente. Não investiguei o motivo, mas acredito que seja decorrente de alguma parametrização diferente entre o algoritmo de forecast::auto.arima() e fable::ARIMA().

Saiba mais

Este exercício simples é um ponto de partida para você começar a implementar a validação cruzada em seus modelos de previsão de séries temporais. Caso deseje se aprofundar, confira o curso de Modelos Preditivos da Análise Macro, onde implementamos diversos modelos para séries temporais da economia brasileira.

Referências
Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3. Accessed on <2022-05-12>.

Receba diretamente em seu e-mail gratuitamente nossas promoções especiais
e conteúdos exclusivos sobre Análise de Dados!

Assinar Gratuitamente