library(tidyverse)
library(forecast)LDE-Séries Temporais
Análise de Séries Temporais com ARIMA no R
Primeiro, vamos acionar os pacotes que serão utilizados nesse 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.
library(sidrar)
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 %>%
select(Mes = `Mês (Código)`,
Valor = Valor) %>%
mutate(Mes = as.numeric(Mes),
Valor = as.numeric(Valor)) %>%
arrange(Mes)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),
frequency = 12)
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) +
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 por 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 = c + \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.
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:
\(c = 0,1\)
\(\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,93\)
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 # 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 o termo model pede que uma 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 = \phi_1y_{t-1} + \epsilon_{t}\)
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)
Enquanto no AR usamos valores passados da própria variável para explicá-la, no componente de média móvel vamos usar informações sobre os “erros” passados para deduzir o erro hoje.
Então, a equação de uma previsão com médias móveis seria do tipo:
\(y_t = \mu + \epsilon_t + \theta_{t-1}\)
Onde:
\(\mu\) é a média de valores observados da séries.
\(\theta\) é o parâmetro que precisa a influência do erro de ontem no resultado de hoje.
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.
Vamos simular um modelo 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.
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.
Modelos ARMA (prelúdio para o ARIMA)
O modelo ARMA(p,q) sugere que os valores atuais de uma séries dependem de dois componentes ou memórias:
AR(p): valores atuais dependem de valores passados da própria série.
MA(q): erros passados.
Uma expressão genérica para esse tipo de modelo seria como a seguinte:
\(y_t = \phi_{t-1}y_{t-1} + ... + \phi_{t-p}y_{t-p} + \theta_{t-1} \epsilon_{t-1} + ... + \theta_{t-q} \epsilon_{t-q}\)
Notem na equação acima as referências aos parâmetros p e q.
Vamos usar uma analogia com o tempo para entender a importância dos componentes AR e MA.
O componente AR sugere que o tempo feito ontem influencia o tempo hoje. Já o componente MA no dias choque climáticos, ou mudanças bruscas no tempo, também afetam o tempo hoje.
O modelo ARMA combina estes dois efeitos.
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.
Componente de integração (d) ou diferenciação.
O que são séries não estacionárias? (próxima aula!).