rm(list=ls()) # Remoção de arquivos no espaço de trabalho.
library(tidyverse) # Data science.
library(forecast) # Previsão de séries temporais no R.
library(sidrar) # Raspagem de dados do SIDRA-IBGE.
library(zoo) # Tratamento de séries temporais no R.
library(quantmod) # Dados financeiros em R.
library(tseries) # Pacote para tratar séries temporais.
library(lmtest) # Realiza testes com os modelos estimados.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. Para abrir chunks: ctrl+Alt+I.
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.9 0.5 1.0 -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.6 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.1 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.7 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.7 -19.8 7.8 10.0 9.2 2.9 2.7 1.5 0.4 -0.5
2021 0.7 -1.5 -2.1 -1.8 1.1 -0.4 -1.3 -0.1 -0.3 -0.5 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.2 0.8 -0.4 0.6 -0.4 -0.5 0.4 -0.4 0.4 1.1 1.5
2024 -1.2 0.2 0.5 -0.4 -1.2 4.4 -1.4 -0.3 1.0 0.0 -0.7 -0.4
2025 0.2 0.0 1.8 -0.7 -0.5 0.1 -0.1 0.7 -0.4 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, artificialmente, 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: o componente de erro
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. E assim por diante. Até o ponto em que um novo valor de parâmetros não produza resultados significativos.
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? Tecnicamente, dizemos que uma série é estacionária quando atende aos seguintes pré-requisito:
A média da série é constante.
A variância da série é contante.
A autocovariância da série depende apenas do intervalo de tempo considerado e não de uma “memória” entre valores passados e presente.
Trecho técnico: autocovariância
A covariância entre duas variáveis é uma medida de quão relacionadas estão. Por exemplo, se uma variável \(X\) cresce no momento em que outra variável, \(Y\), podemos dizer que essas duas séries apresentam covariância positiva. A fórmula da covariância é a seguinte:
\(Cov(X, Y) = E[(X - E[Y])(Y - E[Y])]\)
Onde E(.) significa esperança matemática ou uma média de valores. Se os valores de X são maiores que sua média no momento em que os valores de Y são, também, maiores que sua média; os valores de X são menores que sua média exatamente quando os valores de Y são, também, menores que sua média, o resultado da multiplicação será positivo. Isso evidencia que os dados variam na mesma direção.
No caso da autocovariância, comparamos não os dados de duas séries, mas dados de uma mesma série em momentos diferentes no tempo. Considere uma série \(y_t\). Vamos então apresentar a fórmula da covariância para um intervalo (lag) de tempo igual a \(k\):
\(\gamma_k = Cov(y_t, y_{t-k}) = E[(y_t - \mu)(y_{t-k} - \mu)]\)
Em que \(\mu\) representa o valor médio da série.
Para entender melhor, digamos que a série em questão é o PIB per capita (PIBpc), de 2001 até 2024. Se \(k=3\), então estamos comparando os seguintes valores: PIBpc de 2024 com o de 2021; PIBpc de 2023 com o de 2020 e assim por diante. Ou seja, criamos duas séries. uma de valores “atuais” e outra com valores defasado. Mas a média dessas “duas” séries é uma só, como dissemos, \(\mu\).
Então, para que uma série seja considerada estacionária, é suposto que a autocovariância depende mais das defasagens, por exemplo o valor, o valor do PIBpc de 2024 está relacionado (varia no mesmo sentido) com o PIBpc de 2021, assim como o de 2023 está relacionado (mesma covariância) com o de 2020 e assim por diante. Ou seja, a relação está determinada pelo tamanho da defasagem e a autocovariância é constante.
O que seria um indicativo de não estacionariedade? Se a autocovariância não se mantêm constante (a cada \(k\) períodos de defasagem). Ou seja, a autocovariância se desvincula do lag e passa a ser determinada pelo tempo, t, mesmo.
Vamos usar uma analogia com temperaturas para entender melhor esse ponto. Considere que temos uma série de temperaturas diárias e queremos calcular a autocovariância para um lag de 90 dias ou três meses.
Se entre os dias 21 e 22 de setembro de 2020 a temperatura subiu 3 graus e o mesmo aconteceu entre os dias 21 e 22 de junho dos mesmo ano, podemos considerar que essa série é estacionária. Agora se a temperatura aumentou 1 grau por dia em um mês de verão e caiu 1 grau por dia em um mês de inverno (o lag de três meses os separando), dizemos que a série é não estacionária.
White noise & Random walk
Para visualizar a questão, vamos comparar o gráfico de uma série estacionária (ruído branco) com uma série não estacionária (passeio aleatório).
Ruído branco:
set.seed(123)
x <- rnorm(100) # 100 valores aleatórios N(0,1)
# Visualização
plot(x, type="l", main="Série Estacionária: Ruído Branco", col="blue", lwd = 2)
Passeio aletatório:
# Série não estacionária: passeio aleatório
set.seed(123)
y <- cumsum(rnorm(100)) # soma cumulativa
# Visualização
plot(y, type="l", main="Série Não Estacionária: Passeio Aleatório", col="red", lwd = 2)
Questão: para você, qual gráfico representa melhor uma série com preços de uma ação?
Diferenciação
Considere a seguinte série de PIBpc a preços do ano anterior:
### Pegar a série.
pibpc <- get_sidra(api = "/t/6784/n1/all/v/9809/p/all")All others arguments are desconsidered when 'api' is informed### Limpar a série.
pib_pc <- pibpc %>% # Vamos trabalhar com pipes, para canalizar o que queremos do df salvo.
select(Ano = `Ano (Código)`, # Usamos essa função para selecionar as colunas de nosso interesse.
Valor = Valor) %>%
mutate(Ano = as.numeric(Ano), # Estamos transformando as colunas selecionadas.
Valor = as.numeric(Valor)) %>% # Queremos que elas represente números.
arrange(Ano) # Essa função é usada para ordenar os dados.Criando uma série temporal.
pib_pc_ts <- ts(pib_pc$Valor,
start = c(1996), # Queremos que a série comece nessa data.
frequency = 1) # Frequência, no caso, 12 meses em um ano.Gráfico:
plot(pib_pc_ts,
main = "PIB per capita Brasil (1996-2022)",
xlab = "Ano",
ylab = "PIB per capita",
col = "blue",
lwd = 2)
grid(col = "darkgreen", lty = "dotted") # grade leve
Claramente, a série não é estacionária, apresenta uma forte tendência de crescimento.
Vamos diferenciar a série para tentar torná-la mais estacionária. Em seguida, vamos ao seu gráfico.
dpib_pc_ts <- diff(pib_pc_ts)plot(dpib_pc_ts,
main = "PIB per capita Brasil - série diferenciada",
xlab = "Ano",
ylab = "PIB per capita",
col = "blue",
lwd = 2)
grid(col = "darkgreen", lty = "dotted") # grade leve
Vamos diferenciar duas vezes a série.
dpib_pc_ts2 <- diff(pib_pc_ts, differences = 2)plot(dpib_pc_ts2,
main = "PIB per capita Brasil - série diferenciada",
xlab = "Ano",
ylab = "PIB per capita",
col = "blue",
lwd = 2)
grid(col = "darkgreen", lty = "dotted") # grade leve
A série acima está mais estacionária.
Série de preços de ações
Baixar dados da ação PETR4 (Petrobras - preferencial nominativa é a ação preferencial da Petrobras, com prioridade em dividendos, mas sem direito a voto. O pacote utilizado será o quantmod.
getSymbols("PETR4.SA", src = "yahoo", from = "2020-01-01", to = Sys.Date())[1] "PETR4.SA"class(PETR4.SA)[1] "xts" "zoo"No código acima, a fonte dos dados é o site “Yahoo” e Sys.Date() representa a data de hoje. No resultado, xts significa eXtensible Time Series, ou seja, um objeto de tipo série temporal que é extensível. É uma estrutura de dados especial usada no R para lidar com informações que variam ao longo do tempo — como preços diários de ações, taxas de câmbio, PIB mensal, temperaturas, etc.
Vamos transforma o objeto de tipo xts em data frame.
petra_df <- data.frame(data = index(PETR4.SA), fechamento = as.numeric(Cl(PETR4.SA)))
class(petra_df)[1] "data.frame"No código acima, Cl(.) é uma função do pacote quantmod para pegar o preço de fechamento das ações. Vamos ao gráfico dessa ação.
ggplot(petra_df, aes(x = data, y = fechamento)) +
geom_line(color = "blue", size = 1.1) +
labs(title = "Preço de Fechamento - PETR4 (Petrobras)",
x = "Data",
y = "Preço (R$)") +
theme_minimal()Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Vamos calcular uma média móvel de 15 dias para inserir no gráfico. Vamos usar a função SMA do pacote TTR. Esse pacote já vem junto com o quantmod e serverllll .
mm_15 <- SMA(Cl(PETR4.SA), n = 15)Agora vamos criar um data frame com média móvel de 15 dias.
petra_dfm <- data.frame(
data = index(PETR4.SA),
fechamento = as.numeric(Cl(PETR4.SA)),
mm15 = as.numeric(mm_15)
)O gráfico mais completo fica, então na seguinte forma.
ggplot(petra_dfm, aes(x = data)) +
geom_line(aes(y = fechamento), color = "blue", size = 1.2) +
geom_line(aes(y = mm15), color = "red", size = 1) +
labs(title = "PETR4 - Preço de Fechamento e Média Móvel de 15 dias",
x = "Data",
y = "Preço (R$)") +
theme_minimal()Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).
Diferenciando a série
Vamos diferenciar a série original com os preços de fechamento, segundo a fórmula abaixo:
\(variacaopct= \frac{P_t - P_{t-1}}{P_{t-1}}\)
Para isso, criamos um data frame só com os preços de fechamento. e, em seguida, usamos a função Delt() do pacote quantmod.
fechamento <- Cl(PETR4.SA)Em seguida aplicamos ao fechamento a função Delt(), do pacote quantmod.
variacao <- Delt(fechamento) * 100Agora estruturamos um data frame com o preço de fechamento e sua variação percentual.
petra_dfpct <- data.frame(
data = index(PETR4.SA),
fechamento = as.numeric(fechamento),
variacao_pct = as.numeric(variacao)
)E, finalmente, o seu gráfico.
ggplot(petra_dfpct, aes(x = data, y = variacao_pct)) +
geom_line(color = "darkgreen", size = 1) +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Variação Percentual Diária - PETR4",
x = "Data",
y = "Variação (%)") +
theme_minimal()Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).
Compare os dois gráficos: os dados originais que esboçam um Random walk e os dados da série diferenciada, mais estacionária.
ARIMA
Já vimos com certo grau de detalhes os componentes de um modelo ARIMA: autorregressivo, média móvel e integração.
Chegou a hora de entender e desenvolver em sua inteireza o modelo ARIMA.
O modelo foi proposto no início dos anos 70 por dois estatísticos britânicos, George P. E. Box (1919-2013) e Gwilym M. (1932 – 1982). Vamos, mais uma vez, usar um exemplo prático para explicar como se produz previsões a partir de um modelo ARIMA.
Vamos usar a de taxas de desemprego do Brasil, conforme publicada mensalmente na pesquisa PNAD Contínua (PNAD-c) do IBGE. Primeiro raspamos os dados que nos interessa. Vocês verão que esta é um série que informa o dado de cada mês como uma média de valores trimestrais.
dn <- get_sidra(api="/t/6381/n1/all/v/4099/p/all/d/v4099%201")All others arguments are desconsidered when 'api' is informedVamos separar do data frame as colunas que realmente nos interessam.
dn_df <- dn %>% # Vamos trabalhar com pipes, para canalizar o que queremos do df salvo.
select(trim_m = `Trimestre Móvel (Código)`, Valor = Valor) %>%
mutate(trim_m = as.numeric(trim_m), # Estamos transformando as colunas selecionadas.
Valor = as.numeric(Valor)) %>% # Queremos que elas represente números.
arrange(trim_m) # Essa função é usada para ordenar os dados.Vamos transformar o data frame dn_df em uma série temporal.
dn_ts <- ts(dn_df$Valor,
start = c(2012, 03), frequency = 12)
dn_ts Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2012 8.0 7.8 7.7 7.6 7.5 7.3 7.1 6.9 6.8 6.9
2013 7.2 7.8 8.0 7.9 7.6 7.5 7.4 7.2 7.0 6.8 6.5 6.2
2014 6.5 6.8 7.2 7.2 7.0 6.9 7.0 7.0 6.8 6.7 6.6 6.6
2015 6.9 7.5 8.0 8.1 8.2 8.4 8.6 8.8 9.0 9.1 9.1 9.1
2016 9.6 10.3 11.0 11.3 11.3 11.4 11.7 11.9 11.9 11.9 12.0 12.1
2017 12.7 13.3 13.9 13.7 13.4 13.1 12.9 12.7 12.5 12.3 12.1 11.9
2018 12.3 12.7 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.7 11.7
2019 12.2 12.5 12.8 12.6 12.4 12.1 11.9 11.9 11.9 11.7 11.3 11.1
2020 11.3 11.7 12.4 12.7 13.1 13.6 14.1 14.8 14.9 14.6 14.3 14.2
2021 14.4 14.6 14.9 14.8 14.7 14.2 13.7 13.2 12.7 12.1 11.6 11.1
2022 11.2 11.2 11.1 10.5 9.8 9.3 9.1 8.9 8.7 8.3 8.1 7.9
2023 8.4 8.6 8.8 8.5 8.3 8.0 7.9 7.8 7.7 7.6 7.5 7.4
2024 7.6 7.8 7.9 7.5 7.1 6.9 6.8 6.6 6.4 6.2 6.1 6.2
2025 6.5 6.8 7.0 6.6 6.2 5.8 5.6 5.6 5.6 O nosso objetivo é realizar previsões estatisticamente significativas com base nos dados disponíveis. Vamos seguir um roteiro que de maneira alguma é o único ou, quem sabe, o melhor. Mas é um que tenho adotado em meus trabalhos de previsão).
O modelo ARIMA é montado em duas etapas. Na primeira, escolhemos os argumentos de seus três componentes: AR(p), I(d) e MA(q). Os argumentos p, d e q são números inteiros, por exemplo, ARIMA (1, 1, 1). Na segunda etapa, vamos estimar os parâmetros do modelo. Tais parâmetros são números reais.
Um bom e tradicional início é o de visualizar o gráfico da série.
plot(dn_ts, main = "Taxa de Desemprego - Brasil (IBGE)", ylab = "%", xlab = "Trimestre Móvel")
O gráfico já revela que a série não tem nada de estacionária. Mas vamos confirmar essa não estacionariedade usando um teste estatístico bem padronizado na montagem de modelos ARIMA, o teste Dickey-Fuller Aumentado (ADF). Esse teste investiga a existência de raiz unitária na série. Se a série apresenta raiz unitária, ela é não estacionária.
(Ver trecho técnico depois).
A \(H_0\) do teste afirma que a série é não-estacionária. Logo, se rejeitamos a \(H_0\) a série pode ser considerada estacionária. Para realizar o teste, usamos o pacote tseries.
adf.test(dn_ts) # Teste de Dickey-Fuller Aumentado
Augmented Dickey-Fuller Test
data: dn_ts
Dickey-Fuller = -1.0005, Lag order = 5, p-value = 0.9354
alternative hypothesis: stationaryNo resultado acima, note que o valor da estatística p calculada (o p-value) é maior que 0,05. Isso nos leva a acreditar que a série é não estacionária, pois não podemos rejeitar a \(H_0\).
Para tentar corrigir essa não-estacionariedade, vamos diferenciar a série em um período.
dn_diff <- diff(dn_ts)
plot(dn_diff, main="1ª diferença da taxa de desemprego")
O gráfico acima já revela uma série (diferenciada) bem mais estacionária que a série original. Vamos aplicar, outra vez, o ADF.
adf.test(dn_diff)Warning in adf.test(dn_diff): p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: dn_diff
Dickey-Fuller = -4.672, Lag order = 5, p-value = 0.01
alternative hypothesis: stationaryO p-value agora está bem menor que 0,05. Logo, a serie original diferenciada em um período fica estacionária.
O que concluímos com essa análise? Que um bom argumento para o I do modelo ARIMA é \(d = 1\). Falta definir, agora, os argumentos p e q.
Vamos usar dois testes estatístico bem famosos para definir os argumentos do AR e do MA, os ACF e PACF.
ACF é uma função de autocorrelação, cujo gráfico será empregado para definir o parâmetro q do MA. A função apresenta um gráfico da correlação entre os valores da variável (\(y_t\)) com seus valores defasados (\(y_{t-k}\)), onde \(k\) é a ordem de defasagem. Cada ponto no gráfico é o valor da correlação entre a variável e uma ordem de defasagem, \(k = 1, 2, 3 ...\).
PACF, por sua vez, é uma função de autocorrelação parcial. Mede a correlação parcial entre \(y_t\) e seus valores defasados (\(y_{t - k}\)), mas ignorando os valores intermediários. Exemplo, no cálculo da correlação entre \(y_t\) e \(y_{t-4}\), as observações \(y_{t-2}\) e \(y_{t-3}\) são ignoradas. Vamos usar o gráfico da função de correlação parcial para definir o termo do componente AR, o p.
Gráficos da ACF e PACF
Primeiro, o gráfico da função de autocorrelação.
acf(dn_diff)
Vamos analisar o gráfico acima em detalhes. O eixo x (Lag) mostra o número de defasagens (\(k\) ) — ou seja, em quanto tempo defasamos a série. O eixo vertical é a correlação da série com cada uma de suas defasagens.
As linhas azuis traçejadas são os limites de significância (o mais comum é o de 95%). As barras que ultrapassam essas linhas azuis denotam correlação estatisticamente significante.
As três primeiras defasagens evidenciam uma correlação forte. Ademais, observe que as barras caem de forma abrupta e não suave. A forte correlação observada nos primeiros lags indica que a série guarda muita memória recente, os valores recentes têm forte influência sobre o valor atual da série.
Agora o gráfico da autocorrelação parcial.
pacf(dn_ts)
No gráfico acima, apenas a primeira defasagem apresenta forte autocorrelação parcial. A autocorrelação da primeira defasagem corta a linha azul de forma brusca. Isso indicar um componente AR de ordem \(p=1\).
É comum usar-se uma tabela com a de baixo para definir os parâmetros p e q.
| Comportamento da ACF | Comportamento da PACF | Modelo sugerido | Interpretação geral |
|---|---|---|---|
| Decai exponencialmente ou de forma senoidal | Corta bruscamente após lag p | AR(p) | A série tem forte dependência autorregressiva — os valores passados explicam os atuais. |
| Corta bruscamente após lag q | Decai exponencialmente ou de forma senoidal | MA(q) | A série é explicada por choques anteriores (médias móveis). |
| Decai lentamente (não estacionária) | Decai lentamente (não estacionária) | Necessita diferenciação (d>0) | A série apresenta tendência — precisa diferenciar antes de ajustar ARIMA. |
| Decai exponencialmente | Decai exponencialmente | ARMA(p,q) | A série tem componentes tanto AR quanto MA — ambos influenciam. |
| ACF e PACF cortam no mesmo lag | Ambos decaem rapidamente | ARMA(p,q) possível | Ambos os componentes podem estar presentes, mas precisa testar. |
| Padrão sazonal (picos regulares a cada s lags) | Padrão sazonal nos mesmos pontos | SARIMA | Há sazonalidade periódica (ex.: dados trimestrais ou mensais). |
No caso dos gráficos que encontramos, uma análise seria com a da tabela a seguir.
| Gráfico | Padrão observado | Interpretação comum |
|---|---|---|
| ACF | corta após lag 1 (ou cai rapidamente) | componente MA(1) |
| PACF | corta após lag 1 (ou é significativa só no lag 1) | componente AR(1) |
Portanto, vamos trabalhar com o modelo ARIMA (1,1,1).
Estimando o modelo ARIMA (1,1,1)
Podemos não só estimar, fácil e rapidamente o modelo ARIMA(1,1,1) como compará-lo com outras versões: ARIMA(0,1,1) e ARIMA(1,1,0).
mod_111 <- Arima(dn_ts, order = c(1,1,1))
summary(mod_111)Series: dn_ts
ARIMA(1,1,1)
Coefficients:
ar1 ma1
0.6028 0.1985
s.e. 0.0747 0.0737
sigma^2 = 0.04583: log likelihood = 20.49
AIC=-34.97 AICc=-34.82 BIC=-25.71
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.004593877 0.2121033 0.1568367 -0.04458902 1.645907 0.1053485
ACF1
Training set 0.04614341No quadro acima nós temos os parâmetros do modelo estimado (\(ar1= 0,6028\) e \(ma1 = 0,19\)), valores do teste AIC (Akaike Information Criterion: quanto menor, melhor) e medidas dos erros médios (diferença média entre os valores observados e os valores estimado com base nos parâmetros encontrados).
Esses valores são úteis quando comparamos o modelo estimado com alternativas, ARIMA(0,1,1) e ARIMA(1,1,0). Por exemplo, o melhor modelo é o que apresenta menores estatísticas de AIC e menores erros médios.
Não vamos fazer essas comparações agora. Ao invés, vamos produzir um gráfico com os valores estimados e os valores observados. Quanto mais o gráfico estimado (o “ajustado”) aderir ao gráfico com valores observados, melhor.
Gráfico valores estimados x valores observados
Primeiro criamos uma série temporal com os valores estimados pelo modelo.
# Valores ajustados
ajustados <- fitted(mod_111)Agora o gráfico. Notem que acrescentamos com a função lines o traçado dos valores ajustados pelo modelo.
# Criação do gráfico
plot(dn_ts, col = "darkgreen", lwd = 2,
main = "Taxa de Desemprego: Observado vs Estimado (ARIMA(1,1,1))",
ylab = "Taxa de desemprego (%)",
xlab = "Tempo")
lines(ajustados, col = "blue", lwd = 2)
legend("bottom",
legend = c("Observado", "Estimado (ARIMA(1,1,1))"),
col = c("darkgreen", "blue"),
lwd = 2, lty = c(1,1))
O modelo ficou tão bem ajustado que as duas linhas praticamente se confundem.
Valores previstos
Vamos fazer previsões para os 6 meses (trimestres móveis) seguintes. Para isso usamos a função forecast, do pacote de mesmo nome. O gráfico será feito no ggplot. Como temos uma série temporal e não um data frame, devemos usar a função autoplot, assim como explicamos antes. No código abaixo, a função autolayer acrescenta uma camada de dados em séries temporais no gráfico ggplot.
# Previsão de 12 períodos à frente
previsao <- forecast(mod_111, h = 6)
# Gráfico incluindo previsão
autoplot(previsao) +
autolayer(fitted(mod_111), series = "Ajustado") +
autolayer(dn_ts, series = "Observado") +
ggtitle("Previsão ARIMA(1,1,1) da taxa de desemprego") +
ylab("Taxa de desemprego (%)") +
xlab("Tempo")
Os dados previstos juntamente com seus intervalos de confiança.
previsao Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Oct 2025 5.595800 5.321442 5.870157 5.176205 6.015394
Nov 2025 5.593267 5.028001 6.158534 4.728767 6.457767
Dec 2025 5.591741 4.747728 6.435753 4.300935 6.882546
Jan 2026 5.590821 4.490034 6.691607 3.907312 7.274329
Feb 2026 5.590266 4.255536 6.924996 3.548972 7.631559
Mar 2026 5.589931 4.042043 7.137819 3.222641 7.957222Podemos, ainda, organizar os dados acima em um data frame mais apresentável. Primeiro, para evitar problemas de datas no data frame a ser criado, vamos converter as datas que estão em série temporais em datas de data frame. Para isso, vamos usar a função as.yermon do pacote zoo.
datas_formatadas <- as.yearmon(time(previsao$mean))E agora, montamos a tabela.
tabela_prev <- data.frame(
Data = datas_formatadas,
Previsao = as.numeric(previsao$mean),
Limite_Inferior_95 = as.numeric(previsao$lower[,2]),
Limite_Superior_95 = as.numeric(previsao$upper[,2])
)
print(tabela_prev) Data Previsao Limite_Inferior_95 Limite_Superior_95
1 out 2025 5.595800 5.176205 6.015394
2 nov 2025 5.593267 4.728767 6.457767
3 dez 2025 5.591741 4.300935 6.882546
4 jan 2026 5.590821 3.907312 7.274329
5 fev 2026 5.590266 3.548972 7.631559
6 mar 2026 5.589931 3.222641 7.957222E, finalmente, exportamos os resultados em dados do tipo CSV.
write.csv(tabela_prev, "previsoes_desemprego_br.csv", row.names = FALSE)Na função acima escolhermos row.names=FALSE para impedir que o programa insira uma coluna só com os índices da tabela. Vamos acompanhar as próximas divulgações das taxas de desemprego para ver como as previsões acertaram (ou não!).
Pacotes como o X-13ARIMA-SEATS, do U.S Census Bureau, junto com o pacote seasonal, podem escolher os valores de p, d e q e já estimar os parâmetros do modelo. Esses pacotes são especialmente indicados para as séries com forte componente sazonal.