rm(list=ls())
library(tidyverse)
library(forecast)
library(sidrar)
library(zoo)LDE-Séries Temporais
Análise de Séries Temporais com ARIMA no R
Primeiro, vamos acionar os pacotes que serão utilizados neste estudo. Vejam a utilidade de cada pacote pesquisando na internet.
Séries temporais são dados organizados de forma sequencial no tempo. São dados muito comuns em Economia.
Vamos pegar, como exemplo, dados sobre a Pesquisa Industrial Mensal - Produção Física (PIM-PF) do IBGE, Tabela 8888, no Sidra.
pim <- get_sidra(api="/t/8888/n1/all/v/11601/p/all/c544/129314/d/v11601%201")Agora vamos limpar nossa série. Você pode usar uma LLM para explicar os detalhes do código.
pim_pf <- pim %>% # Vamos trabalhar com pipes, para canalizar o que queremos do df salvo.
select(Mes = `Mês (Código)`, # Usamos essa função para selecionar as colunas de nosso interesse.
Valor = Valor) %>%
mutate(Mes = as.numeric(Mes), # Estamos transformando as colunas selecionadas.
Valor = as.numeric(Valor)) %>% # Queremos que elas represente números.
arrange(Mes) # Essa função é usada para ordenar os dados.A série começa em janeiro de 2002. Vamos montar agora o objeto série temporal (ts) propriamente dito.
pim_ts <- ts(pim_pf$Valor,
start = c(2002, 01), # Queremos que a série começe nessa data.
frequency = 12) # Frequência, no caso, 12 meses em um ano.
pim_ts Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2002 0.0 2.2 -0.3 0.2 0.0 0.3 0.4 -0.9 0.4 0.4 -0.1 -0.3
2003 0.4 -0.4 -0.2 -0.2 -1.5 0.6 -0.9 1.4 2.8 0.5 1.1 -1.5
2004 2.0 1.5 0.6 0.8 0.9 0.4 0.7 1.1 0.4 -0.6 0.0 0.6
2005 0.5 -1.4 1.0 0.0 3.5 -2.1 -0.6 0.1 -0.7 -0.5 0.9 1.7
2006 0.4 0.8 -0.5 0.5 0.7 -1.2 0.8 -0.1 -1.2 0.9 1.3 0.6
2007 -0.2 1.3 1.3 1.0 0.2 0.6 -0.3 1.4 -1.2 2.8 -0.9 0.3
2008 1.7 -0.3 0.2 0.1 0.8 1.3 0.5 -1.1 1.1 -2.7 -6.1 -11.1
2009 2.1 2.3 0.5 1.3 1.2 1.4 1.3 1.6 1.4 2.2 1.3 0.4
2010 1.5 -0.2 2.3 0.4 -0.3 -0.4 -0.6 -0.2 0.2 -0.3 0.7 0.5
2011 0.3 1.2 0.3 -2.6 3.1 -1.8 -0.2 -0.8 -1.6 -1.1 0.5 2.6
2012 -4.2 0.8 -0.4 -0.1 0.5 0.2 1.4 2.2 -2.0 1.5 -1.3 0.6
2013 1.0 -2.9 2.0 1.6 -0.2 0.4 -0.6 0.3 0.5 -1.1 0.5 -3.5
2014 1.8 0.4 -0.4 -0.8 -1.0 -2.8 2.5 0.4 0.2 0.1 -0.9 -2.8
2015 0.5 -0.9 -0.4 -1.9 -1.1 -0.8 -1.8 0.4 -2.3 -0.1 -1.7 -1.5
2016 0.9 -1.1 0.6 0.0 0.4 -1.4 0.7 -1.9 1.1 -1.4 1.0 1.7
2017 0.8 1.2 -1.6 0.0 0.4 1.0 0.2 0.1 0.3 1.0 0.4 2.7
2018 -2.3 -0.6 1.4 0.4 -10.9 12.4 -0.2 -0.9 -2.9 1.4 -0.1 -0.2
2019 -0.5 0.3 -0.6 0.5 0.0 -0.9 -0.5 1.0 0.0 2.3 -2.4 -1.4
2020 1.4 0.6 -7.6 -19.7 7.8 10.0 9.3 2.9 2.7 1.6 0.4 -0.5
2021 0.7 -1.5 -2.1 -1.8 1.1 -0.4 -1.2 -0.1 -0.3 -0.4 0.0 1.5
2022 -0.1 0.9 -0.8 0.6 -0.5 -0.1 1.2 -1.6 -1.5 1.5 0.0 0.0
2023 0.0 -0.1 0.8 -0.4 0.6 -0.4 -0.5 0.4 -0.4 0.4 1.1 1.6
2024 -1.2 0.3 0.4 -0.3 -1.2 4.5 -1.4 -0.3 0.9 0.0 -0.7 -0.4
2025 0.2 0.0 1.7 -0.6 -0.5 0.1 -0.1 0.8 Gráfico da série com o ggplot2
Para desenhar o gráfico da série temporal usando o pacote ggplot precisamos da função autoplot. Mas para que isso funcione, precisamos do pacote forecast, desenvolvido para se trabalhar séries temporais.
autoplot(pim_ts, colour = "red", size = 1.2) + # Essa função no ggplot permite fazer gráficos de objetos como séries temporais.
ggtitle("Produção Física - 2002-2025") +
ylab("Variação Percentual") +
xlab("Mês/Ano") +
theme_minimal()
Analisando a série
Vamos começar decompondo a série em seus elementos, principais, a saber, tendência (trend), componente sazonal (seazonal) e ruído (random). Normalmente, a série temporal é a soma desses três componentes.
decomposicao <- decompose(pim_ts)
plot(decomposicao)
Vamos usar o pacote monthplot para visualizar a sazonalidade da série de maneira mais nítida. Esse pacote calcula a médias de cada mês.
monthplot(pim_ts, col = "darkgreen", lwd = 2,
main = "Sazonalidade da PIM-PF", ylab = "Variação percentual")
Pelo gráfico, notamos que a menor média é em abril de cada ano e a maior em julho.
Previsão (modelo ARIMA)
ARIMA é a sigla para:
AutoRegressive Integrated Moving Average, ou seja, modelo autorregresivo, integrado e com média móvel.
É um modelo usado para fazer previsões sobre valores futuros de uma série. O modelo é composto de três fatores:
O fator autorregressivo, AR(p).
O fator de integração, I(d).
O fator da média móvel, MA(q).
Os parâmetros p, d & q representam as dimensões de cada parte do modelo ARIMA. Vamos procurar entender o que são essas partes.
Componente autorregressivo, AR(p)
A autorregressão quer dizer que procuramos entender como uma série se comporta a partir de seus valores passados. Em símbolos,
\(y_t = f(y_{t-1})\)
Se a série é explicada, com duas defasagens, ou seja, \(p = 2\), o componente autorregressivo do modelo pode ser escrito da seguinte maneira:
\(y_t = \phi_1y_{t-1} + \phi_2y_{t-2} + \epsilon_t\)
Onde, c é um parâmetro constante que dá conta dos fatores que afetam a série e não tem nada a ver com seus componentes autorregressivos. Seria o intercepto de uma equação linear.
Os parâmetros \(\phi\) revelam como os componentes defasados da série a afetam. Por exemplo, digamos que temos os seguintes valores para cada item da série:
\(\phi_1 = 0,7\)
\(\phi_2 = 0,2\)
\(y_{t-1} = 4\)
\(y_{t-2} = 6\)
\(\epsilon_t = -0.17\)
Então,
\(y_t = 3,83\)
Que tal nós simularmos no R um série com padrão autorregressivo de \(p = 1\)? Na simulação vamos usar a função arima.sim do pacote stats, do R-base.
# Simulação de um processo AR(1)
set.seed(123) # Valores aleatórios, mas fixos.
phi <- 0.8 # Esamos determinando o coeficiente autorregressivo.
n <- 100 # Tamanho da série
y <- arima.sim(model = list(ar = phi), n = n)
# Plotar a série simulada
plot(y, type = "l", col = "blue", lwd = 2,
main = paste("Série AR(1) Simulada (φ =", phi, ")"),
xlab = "Tempo", ylab = expression(y[t]))
grid()
Esclarecendo alguns pontos do código acima. Usamos a função arima.sim para simular um modelo autorregressivo, com lag = 1. Como argumento da função arima.sim a função model pede que um conjunto de argumentos (uma lista, na verdade, daí o uso da função list). E que argumentos são esses em nosso exemplo? O termo de defasagem e o número de observações. A função cria, então o conjunto de dados, gerados aleatoriamente, com um erro normalmente distribuído, ou seja,
\(\epsilon_t \sim \mathcal{N}(0,1)\)
Ou seja, o erro segue um padrão de ruído branco.
A rigor, portanto, cada ponto da linha azul foi gerado levando em conta a seguinte equação:
\(y_t = 0,8y_{t-1} + \epsilon_{t}\)
É importante esclarecer que nós não geramos um modelo Arima de séries temporais. Geramos, artificilmente, um série temporal: valores ao longo do tempo. Esses valores são unidos pela linha azul, o que dá origem ao gráfico acima.
Quando queremos estimar os parâmetros \(\phi\) do modelo AR que representa a série temporal, os pacotes utilizam o método da máxima verossimilhança: o software encontra os parâmetros que maximizam a verossimilhança, ou seja, que tornam os dados observados mais prováveis.
Vamos visualizar as consequências de escolhas diferentes para \(\phi\). Para isso, utilizamos o seguinte código:
set.seed(123)
n <- 100
# Três séries com diferentes níveis de autorregressão
y1 <- arima.sim(model = list(ar = 0.2), n = n)
y2 <- arima.sim(model = list(ar = 0.8), n = n)
y3 <- arima.sim(model = list(ar = -0.8), n = n)
#png("graficoAR.png", width = 1000, height = 800) # Para exportar o gráfico em PNG.
par(mfrow = c(3,1))
plot(y1, type = "l", col = "darkgreen", lwd = 2,
main = "AR(1) com φ = 0.2 (fraca dependência)",
xlab = "Tempo", ylab = expression(y[t]))
grid()
plot(y2, type = "l", col = "blue", lwd = 2,
main = "AR(1) com φ = 0.8 (forte dependência positiva)",
xlab = "Tempo", ylab = expression(y[t]))
grid()
plot(y3, type = "l", col = "red", lwd = 2,
main = "AR(1) com φ = -0.8 (forte dependência negativa)",
xlab = "Tempo", ylab = expression(y[t]))
grid()
Analisando esses resultados:
- \(\phi = 0,2\) : a influência de valores passados é pequena. A série é quase um ruído branco.
- \(\phi = 0,8\) : a série fica menos errática, se move com mais lentidão. Há memória persistente na sua composição.
- \(\phi = -0,8\) : valores positivos e negativos se alterna.
A intuição por trás do modelo AR: As séries econômicas geralmente são do tipo AR, ou seja, os resultados ou números de hoje são reflexo dos de ontem. É como se em Economia, a Lei de Brandford é bem presente.
Componente de média móvel, MA(q)
O termo média móvel tem dois significados distintos em estátística das séries temporais ou econometria.
O primeiro significado é o mais comum: média móvel é um processo de suavização de uma série temporal. Por exemplo, durante a pandemia do Covid-19, era sempre informada a média móvel de casos registrados. Média móvel aqui é a média de \(k\) observações da série. Então, cada \(k\) observações eram substituídas por por sua média, recriando uma série suavizada e com menos valores. O outro sentido é o de uma estimativa dos erros usada em modelos do tipo ARIMA.
Vamos ilustrar a média móvel como a suavização de uma série usando dados sobre a Covid-19 no Brasil. No site https://covid.saude.gov.br/, clique em Painel Interativo. Em seguida, clique em Baixar Dados de Casos novos por semana. Infelizmente o arquivo vem em formato proprietário. Eu abri o arquivo, exclui a coluna de datas, nomeie a coluna de valores como “casos” e o salvei com o nome “covid.csv”.
Vamos carregar esses dados:
covid <-read.csv("covid.csv", header = T)Vamos acrescentar uma coluna de datas semanais ao nosso data frame.
covid$semana <- seq(
as.Date("2025-01-01"),
by = "week",
length.out = nrow(covid)
)Vamos acrescentar a média móvel de 3 semanas ao nosso data frame.
covid_mm <- covid %>%
mutate(media_movel = rollmean(casos, k = 3, fill = NA, align = "right"))Explicando o código acima: Usamos um dos mais antigos e importantes pacotes do R que trabalha com séries temporais, o pacote zoo. A função é a rollmean que produz médias móveis. Escolhemos fazer médias de 3 semanas. Quando não houver valor, a informação fica com NA. E de fato, como precisamos de 3 observações para se calcular a média móvel, as duas primeiras ficam com NA. Calculamos uma média e passamos para a seguinte, avançando uma semana (daí o termo móvel). O alinhamento à direita significa que estamos representando os dados da última semana sendo considerada.
Vamos fazer o gráfico da série original com sua média móvel.
ggplot(covid_mm, aes(x = semana)) +
geom_line(aes(y = casos), color = "darkgreen", linewidth = 1.2) +
geom_line(aes(y = media_movel), color = "blue", linewidth = 1.2) +
labs(
title = "Série com Média Móvel",
x = "Data",
y = "Valor"
) +
theme_minimal()Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).
Notem o aviso de que dois valores foram desconsiderados.
O outro tipo de erro é o estocástico e será tratado no contexto de um modelo de séries temporais ARMA e ARIMA.
Modelos ARMA, prelúdio para os modelos ARIMA
Modelos ARMA são usados para representar séries temporais consideradas estacionárias (mais sobre isso depois). São dois os componentes de um modelo ARMA:
Um componente autorregressivo (AR) — que, como vimos, usa os valores passados da própria série.
Um componente de média móvel (MA) — que usa os erros (ou choques) passados do modelo.
A estrutura do modelo é a seguinte:
\(y_t = \mu + \phi_1y_{t-1} + ... + \phi_py_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + ... + \theta_q\epsilon_{t-q}\)
Onde:
\(y_t\) é o valor da série no momento t.
\(\mu\) é o valor médio da série.
\(\phi_1\) é o parâmetro estimado relacionado à observação da série no tempo t-1.
\(\epsilon_t\) o erro da série no momento t.
\(\theta_1\) o parâmetro associado ao erro defasado de um período.
O componente MA(1)
O fator de média móvel se refere à importância de estimativas dos erros para o valor de uma série. O valor atual é resultado de uma combinação do erro atual e do erro passado:
\(y_t = \mu + \epsilon_t + \theta_1\epsilon_{t-1}\)
Portanto, o termo média móvel quer dizer aqui média de erros e não média de observações, como no caso anterior.
No caso de um ARMA(1,1) estamos considerando que a série pode ser representada pelo seguinte modelo:
\(y_t = \mu + \phi_1y_{t-1} + \epsilon_t + \theta_1\epsilon_{t-1}\)
Intuitivamente, estamos afirmando que um bom modelo a representar uma série temporal é o leva em conta o passado recente da série e um componente de correção do erros de previsão cometidos no passado também recente.
Se \(\theta_1 > 0\) erro de ontem “reforça” o erro de hoje. Se \(\theta < 0\), os erros de ontem amenizam a ocorrência do erro de hoje. É como aprender com seus próprios erros. Vamos estudar esse assunto com ajuda de simulações feitas no R, mais uma vez.
Intuição por trás de um modelo com coeficiente de MA negativo (\(\theta<0\)): aqui está a crença de que os agentes econômicos aprendem com seus erros. Ou seja, suas expectativas estão sendo adaptadas conforme erros de previsão aconteçam. É um caso de expectativas adaptativas.
Simulando modelos de média móvel no R.
set.seed(123)
y1 <- arima.sim(model = list(ma = 0.8), n = 100)
y2 <- arima.sim(model = list(ma = -0.8), n = 100)
#png("graficoMA.png", width = 1000, height = 800) # Para exportar o gráfico em PNG.
par(mfrow = c(2,1))
plot(y1, type = "l", col = "blue", main = "MA(1) com θ = 0.8")
plot(y2, type = "l", col = "red", main = "MA(1) com θ = -0.8")
Os dois gráficos mostram o seguinte:
\(\theta = 0,8\) : o erro passado reforça o erro atual.
\(\theta = -0,8\) : o erro passado corrige o erro atual.
Para estimar os parâmetros \(\theta\) com base em um conjunto de dados, usamos o método da máxima verossimilhança mais uma vez.
Condições para aplicação do modelo ARMA:
A série precisa ser estacionária (mais sobre isso depois).
O erro é do tipo erro branco.
Procedimentos estatísticos são usados para definir as ordens dos componentes AR(p) e MA(q).
Simulando um modelo ARMA
Exemplo no R. Acompanhe o seguinte código, usando, mais uma vez, a função arima.sim. Agora temos dois argumentos, dando mais sentido à lista no argumento da função.
set.seed(123)
# parâmetros
phi <- 0.6 # AR(1)
theta <- 0.8 # MA(1)
n <- 100 # tamanho da série
# simulação
y <- arima.sim(model = list(ar = phi, ma = theta), n = n)
# gráfico
plot(y, type = "l", col = "purple", lwd = 2,
main = paste("Série ARMA(1,1): φ =", phi, ", θ =", theta),
xlab = "Tempo", ylab = expression(y[t]))
grid()
Vamos comparar os três modelos que estudamos até aqui: AR, MA e ARMA.
set.seed(123)
# AR(1)
y_ar <- arima.sim(model = list(ar = 0.6), n = 100)
# MA(1)
y_ma <- arima.sim(model = list(ma = 0.8), n = 100)
# ARMA(1,1)
y_arma <- arima.sim(model = list(ar = 0.6, ma = 0.8), n = 100)
#png("graficoARMA.png", width = 1000, height = 800) # Para exportar o gráfico em PNG.
par(mfrow = c(3,1))
plot(y_ar, type = "l", col = "blue", main = "AR(1) φ=0.6")
plot(y_ma, type = "l", col = "red", main = "MA(1) θ=0.8")
plot(y_arma, type = "l", col = "purple", main = "ARMA(1,1) φ=0.6 θ=0.8")
Resultados: no primeiro gráfico vê-se uma forte persistência. No segundo gráfico note a elevada flutuação de valores. E no terceiro gráfico temos a combinação dos dois valores.
Estimando um modelo ARMA(1,1)
Digamos que tenhamos escolhidos os parâmetros \(p = 1, q = 1\) para nosso modelo ARMA, ou seja, acreditamos (por razões estatísticas) que nossa série é razoavelmente estacionária e que informações recentes sobre os valores da variável e sobre erros (choques) a explicam bem.
O que significa, então, estimar um modelo ARMA?
O modelo a ser estimado é o seguinte:
\(y_t = \mu + \phi_1y_{t-1} + \epsilon_t + \theta_1\epsilon_{t-1}\)
Precisamos estimar os parâmetros: \(\mu\), \(\phi_1\) e \(\theta_1\).
Uma dificuldade prática com essa estimação é que, ao contrário do que acontece com \(y_t\), não conhecemos previamente os valores do erro. Então não podemos simplesmente aplicar um método de regressão linear padrão.
Então a ideia é escolher os parâmetros que tornam o modelo o mais próximo possível do conjunto de dados originais. Ou seja, parâmetros que têm a maior chance de representar de forma satisfatória os dados disponíveis. Esse é sentido de um modelo de máxima verossimilhança.
Trecho técnico
No estimação do modelo pelo método da máxima verossimilhança (ML), supomos que os erros são normalmente distribuídos.
\(\epsilon_t \sim \mathcal{N}(0,1)\)
Portanto, a função densidade de probabilidade (FDP) dos erros é uma típica função normal. Uma FDP estima a probabilidade de ocorrência de eventos para o caso de variáveis contínuas.
\(f(\epsilon_t) = \frac{1}{\sqrt{2\pi\sigma^2}} exp^{(-\frac{\epsilon^2_i}{2\sigma^2})}\)
A fórmula acima expressa a probabilidade de ocorrência de um dado erro (\(\epsilon_t\)). Levando em conta todos os erros previstos, a probabilidade total será dada pelo produto dessas probabilidades individuais, considerando que os erros são independentemente distribuídos.
\(L(μ,ϕ,θ,σ^2)= \prod_{t=1}^{T}(2\pi\sigma^2)^{-1/2}exp(-\frac{\epsilon^2_t}{2\sigma^2})\)
Cada \(\epsilon_t\) é resultado da diferença entre o valor verdadeiro da variável e o valor estimado pelo modelo (por isso, os parâmetros estão implícitos na função \(L(.)\) acima.
Maximizar a verossimilhança consiste em determinar os valores dos parâmetros (\(\mu\), \(\phi\), \(\theta\) e \(\sigma^2\)) que maximizam a função \(L(.)\).
Intuitivamente, estamos dizendo que a função \(L(.)\) expressa a probabilidade de ocorrência dos valores observados, dado os conjunto de parâmetros, ou seja:
\(L(\mu, \phi, \theta, \sigma^2) = P(y_1, ...y_T| \mu, \phi, \theta, \sigma^2)\)
Se o modelo está bem ajustados, os \(\epsilon_t\) serão pequenos. (Observe que na fórmula do \(L(.)\) o erro é um expoente com sinal negativo. Portanto, minimizar o erro é o mesmo que maximizar a função de máxima verossimilhança.
Abrindo \(\epsilon_t\)
\(\epsilon_t = \hat{y}_t - y_t\)
onde:
\(y_t = \mu + \phi_1y_{t-1} + \epsilon_t + \theta_1\epsilon_{t-1}\)
e
\(\hat{y}_t = \mu + \phi_1y_{t-1} + \theta_1\epsilon_{t-1}\)
Método interativo
A escolha dos parâmetros do modelo, ou seja, aqueles que maximizam a função \(L(.)\) e minimizam os \(\epsilon_t\) se dá por interação.
Tecnicamente falando, o algoritmo que, por exemplo, o R usará para definir os parâmetros do modelo é um método numérico interativo de otimização*
O processo incia com o algoritmo lançando um conjunto de valores inciais para os parâmeros, digamos:
\(\mu^{(0)}\), \(\phi^{(0)}\), \(\theta^{(0)}\), \(\sigma^{2(0)}\)
Os resíduos dessa primeira rodada são calculados:
\(\epsilon^{(0)}_t = y_t - \mu^{(0)} - \phi^{(0)}y_{t-1} - \theta^{(0)}\epsilon^{(0)}_{t-1}\)
A verossimilhança é calculada com esses resíduos (na verdade, por facilidade, o log da verossimilhança. Os parâmetros são levemente ajustados e a verossimilhança é calculada outra vez.
Em cada interação é calculado o gradiente da função no ponto. Esse gradiente aponta na direção para que os novos parâmetros sejam escolhidos de forma a aumentar o valor de \(L(.)\). Isso é feito até que não haja novos aumentos significativos, ou seja, se chegou à cominação ótima de parâmetros que otimiza \(L(.)\) e o gradiente, nesse ponto máximo, é aproximadamente igual a zero. Intuitivamente, é como se o algorítimo usasse o gradiente com bússola, ela sinaliza para qual direção a função deve crescer mais.
Componente de integração (d) ou diferenciação.
O que são séries não estacionárias? (próxima aula!).