Introdução

Os cursos de econometria de séries de tempo, usualmente, começam pelo ensino de modelos lineares univariados para séries estacionárias. Estes modelos são da família ARMA e tentam representar uma série de tempo \(y_{t}\) em função de suas defasagens \(y_{t-1}, y_{t-2}, \dots, y_{t-n}\) e de choques aleatórios (inovações) \(\epsilon_{t}, \epsilon_{t-1}, \epsilon_{t-2}, \dots, y_{t-n}\). Contudo, pode ser mais interessante relacionar duas séries de tempo \(y_{t}\) e \(x_{t}\) diferentes via um modelo linear. Em alguns casos isto pode ser equivalente a um VAR ou VARMA, mas o modelo linear tem a vantagem de ser mais simples de implementar e de interpretar. O tipo de modelo linear que estamos interessados é da forma

\[ y_{t} = \beta_{0} + \beta_{1}x_{t} + w_{t} \]

onde \(y_{t}\) é a série que queremos “explicar” em função da série \(x_{t}\). É evidente que podemos estender este modelo para incluir defasagens das variáveis \(x_{t}\) e \(y_{t}\), além de incluir outras séries, dummies, efeitos sazonais e tendências temporais.

Quando se usa dados em forma de séries de tempo numa regressão linear, é bastante comum que se enfrente algum nível de autocorrelação nos resíduos. Uma das hipóteses do modelo “clássico” de regressão linear é de que as observações são i.i.d., isto é, que os dados são independentes e identicamente distribuídos. Isto obviamente não se aplica no contexto de séries de tempo (os dados não são independentes), então é preciso algum cuidado no uso de modelos de regressão linear. Neste sentido, o diagnósito dos resíduos é o mais importante passo para verificar a qualidade do modelo. Idealmente, os resíduos do modelo devem se comportar como ruído branco (i.e., não devem apresentar autocorrelação).

Outro problema típico deste tipo de análise, chamado de “regressão espúria”, acontece quando se faz a regressão de séries não-estacionárias. O último exemplo deste post mostra exatamente isto. Quaisquer duas séries com tendência serão fortemente linearmente relacionadas. Isto leva a uma regressão com \(R^2\) altíssimo e estatísticas-t muito significativas e a vários modelos sem sentido. Exemplos disto podem ser vistos neste site (em inglês).

Pacotes

# Bases de dados
library("AER")
library("astsa")
# Funções para séries de tempo
library("forecast")
library("tseries")
library("dyn")
# Funções para facilitar a manipulação de dados com datas
library("xts")
library("lubridate")
# Visualização
library("ggplot2")

Exemplo: PIB e Consumo

Como primeiro exemplo vamos analisar a variação do PIB. A base de dados USMacroG é um conjunto de 12 séries macroeconômicas dos EUA, disponibilizadas pelo livro de econometria do Greene. Como as séries do PIB e do consumo estão em nível precisamos fazer alguma transformação para convertê-los em taxas percentuais. O importante do código abaixo é notar o uso da função ts.instersect que serve para emparelhar as séries e transformá-las em colunas de um data.frame. O modelo que queremos estimar relaciona as variações do PIB com as variações do consumo da seguinte forma.

\[ \Delta PIB_{t} = \beta_{0} + \beta_{1} \Delta C_{t} + u_{t} \]

onde \(\Delta\) representa a variação percentual, isto é, \(\Delta x_{t} = \frac{x_{t} - x_{t-1}}{x_{t-1}}\). Para estimar o modelo uso a função lm. A saída abaixo é essencialmente idêntica à de uma regressão linear de dados em cross-section

\[ \hat{PIB}_{t} = \underset{(0.00007)}{0.00264} + \underset{(0.06316)}{0.68428}C_{t} \]

Dependent variable:
PIB (variação)
Consumo (variação) 0.684***
(0.063)
constante 0.003***
(0.001)
Observations 203
R2 0.369
Adjusted R2 0.366
Residual Std. Error 0.008 (df = 201)
F Statistic 117.364*** (df = 1; 201)
Note: p<0.1; p<0.05; p<0.01
# Carrega os dados
data("USMacroG")
d <- USMacroG
# Variação percentual do PIB
pib <- d[, "gdp"] / lag(d[, "gdp"], -1) - 1
# Variação percentual do Consumo
cons <- d[, "consumption"] / lag(d[, "consumption"], -1) - 1
dados <- ts.intersect(pib, cons, dframe = TRUE)
# Regressão do PIB contra o Consumo (contas nacionais)
summary(fit <- lm(pib ~ cons, data = dados))

Como há apenas duas séries podemos visualizar a sua relação num gráfico de dispersão ou mesmo num gráfico de linha.

# Scatterplot
plot(pib ~ cons, data = dados)
# Linha de regressão em vermelho
abline(coef(fit), col = "red")

# Gráfico de linha
plot.ts(pib)
lines(cons, col = "red")
legend("topright", legend = c("PIB", "Consumo"), col = 1:2, lty = 1)

Ainda que práticas, as funções base do R para visualização não são muito elegantes e há vários pacotes que oferecem opções mais interessantes. O popular pacote ggplot2, em particular, foi atualmente integrado ao pacote forecast nas funções autoplot e autolayer.

ggplot(dados, aes(x = cons, y = pib)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

autoplot(pib) +
  autolayer(cons)

É sempre importante verificar a presença de autocorrelação nos resíduos do modelo - este problema é especialmente presente em contextos de séries de tempo.

u <- resid(fit)
par(mfrow = c(2, 1))
forecast::ggAcf(u)

forecast::ggPacf(u)

Neste caso parece não haver grande problema já que a maior parte das autocorrelações estão dentro do intervalo de confiança (com exceção do segundo lag). Há testes mais formais, como o Breusch-Godfrey ou Ljung-Box, para verificar a presença de autocorrelação nos resíduos de um modelo. Em ambos os casos, queremos testar se os resíduos do modelo se comportam como ruído branco e, idealmente, queremos não-rejeitar a hipótese nula (isto é, busca-se um p-valor “grande”). O teste BG está implementado no pacote lmtest que é carregado automaticamente junto com o pacote AER. Ele essencialmente faz uma regressão linear entre os resíduos do modelo com defasagens da variável explicada (\(y\)) e dos próprios resíduos e verifica se os coeficientes estimados são significativos (diferentes de zero). O código abaixo mostra como fazer o teste para a segunda ordem.

\[ \hat{u}_{t} = \alpha_{0} + \alpha_{1}C_{t} + \alpha_{2}\hat{u}_{t-1} + \alpha_{3}\hat{u}_{t-2} + \epsilon_{t} \]

bgtest(fit, order = 2)
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  fit
## LM test = 4.8858, df = 2, p-value = 0.08691

A tabela abaixo reporta o valor do teste para diversas ordens. Como se vê, não há evidências de que exista correlação serial nos resíduos do modelo. Caso houvesse autocorrelação, uma saída seria modelar os resíduos como um processo ARMA e estimar o modelo conjuntamente por máxima verossimilhança. Isto é, seria um modelo de regressão linear com erros ARMA. Discuto mais sobre este tipo de modelo e como estimá-lo neste outro post.

Ordem da defasagem Estatística LM p-valor
1 0.1074866 0.7430239
2 4.8857641 0.0869100
3 4.8859138 0.1803443
4 7.3978074 0.1163009
5 8.3037142 0.1402726
6 8.3038739 0.2166753
7 8.3473089 0.3029668
8 8.4821106 0.3878461
9 8.5054023 0.4841218
10 8.8502628 0.5463708
11 8.8803034 0.6329399
12 11.6353655 0.4753901

Também seria importante verificar se as séries das variações trimestrais do PIB e do consumo são estacionárias. Caso elas não fossem, os resultados da regressão provavelmente seriam espúrios. Contudo, em geral, quando tiramos a diferença percentual de uma série ela se torna estacionária, então deixo esta questão de lado. Discuto mais sobre não-estacionaridade e regressão espúria neste outro post.

Decomposição de tendência e sazonalidade

Pode-se propor um modelo linear que separa efeitos de tendência e sazonalidade da seguinte forma:

\[ y_{t} = T_{t} + S_{t} + u_{t} \]

onde \(T_{t}\) é uma componente de tendência (em geral, um polinômio ou média móvel) e \(S_{t}\) é um componente sazonal. Talvez a forma mais simples deste modelo seja assumir uma tendência linear e incluir dummies sazonais.

\[ y_{t} = \beta_{0} + \beta_{1}t + \sum_{k = 2}^{k = 12}\alpha_{k}d_{k} + u_{t} \]

Para facilitar a notação a equação acima assume dados mensais, onde \(d_{k}\) é uma variável dummy igual a \(1\) quando a observação for referente ao mês \(k\). Isto é, \(d_{2} = 1\) quando os dados forem do mês de fevereiro, logo, \(\alpha_{2}\) captura o efeito marginal do mês de fevereiro. Não se pode colocar uma dummy para cada mês pois haveria um problema de singularidade na matriz de regressores, mas é claro que quando o dado for do mês de janeiro, todas as dummies terão valor igual a zero.

Exemplo: Tendência linear

\[ y_{t} = \beta_{0} + \beta_{1}t + \epsilon_{t} \]

gdp <- d[, "gdp"]
time_trend <- seq_along(gdp)
fit <- tslm(gdp ~ time_trend)
Dependent variable:
PIB
t 35.236***
(0.449)
constante 950.919***
(53.065)
Observations 204
R2 0.968
Adjusted R2 0.968
Residual Std. Error 377.571 (df = 202)
F Statistic 6,161.458*** (df = 1; 202)
Note: p<0.1; p<0.05; p<0.01

O gráfico abaixo mostra o ajuste com a linha de tendência em vermelho. Fica claro que os dados nas extremidades da série fogem bastante da tendência. Isto pode tanto ser indício de um mau ajuste, como de que o crescimento recente foi “acima da tendência”.

autoplot(gdp) +
  autolayer(fitted(fit)) +
  scale_color_discrete(name = "", label = "Tendência Linear") +
  labs(title = "PIB (tendência linear)", x = "", y = "")

Parte do motivo de se estimar uma tendência para a série é para poder removê-la da série original. Isto é, queremos fazer

\[ y_{t}^{\text{detrend}} = y_{t} - \hat{y_{t}} = y_{t} - \hat{\beta_{0}} + \hat{\beta_{1}}t \]

onde \(y_{t}^{\text{detrend}}\) é a série sem tendência. Isto é particularmente útil quando estamos lidando com séries “tendência-estacionárias”, i.e., séries que se tornam estacionárias quando subtraímos a sua tendência. No caso da série do PIB há algum debate sobre como torná-la estacionária. Em geral, é mais comum tratá-la como “diferença-estacionária” e tirar a primeira diferença da série (implicitamente assumindo que ela possui uma raiz unitária). Contudo, há alguma evidência de que a série possa ser “tendência-estacionária” quando se considera uma tendência com quebra.

De qualquer forma, para este modelo simples podemos facilmente computar a série sem tendência. Visualmente, a série não parece estacionária, mas seria necessário fazer algum teste formal para concluir isto com mais certeza.

detrend <- gdp - fitted(fit)
autoplot(detrend)

Exemplo: Tendências polinomiais (ordens mais altas)

Uma forma de estender o modelo acima é considerar ordens polinomiais mais elevadas. De forma geral, o modelo teria a seguinte forma:

\[ y_{t} = \sum_{i = 0}^{k}\beta_{k}t^{k} + \epsilon_{t} \]

onde \(k\) é a ordem que estamos considerando. Para uma tendência cúbica, por exemplo, teríamos:

\[ y_{t} = \beta_{0} + \beta_{1}t + \beta_{2}t^{2} + \beta_{3}t^{3}+ \epsilon_{t} \]

O código abaixo estima justamente este modelo e os resultados são reportados na tabela.

trend <- seq_along(time_trend)
fit_poly <- tslm(gdp ~ poly(trend, 3, raw = TRUE))
Dependent variable:
PIB
t 22.411***
(1.566)
t2 -0.011
(0.018)
t3 0.0004***
(0.0001)
constante 1,563.528***
(37.166)
Observations 204
R2 0.996
Adjusted R2 0.996
Residual Std. Error 130.279 (df = 200)
F Statistic 17,749.800*** (df = 3; 200)
Note: p<0.1; p<0.05; p<0.01
autoplot(gdp) +
  autolayer(fitted(fit_poly)) +
  scale_color_discrete(name = "", label = "Tendência Cúbica")

Como se vê no gráfico, a tendência cúbica oferece uma aproximação melhor da série original, contudo, há motivos para ter cuidado com este tipo de modelo. À medida que se aumenta a ordem do polinômio, o ajuste da tendência vai progressivamente melhorando. Isto é problemático por pelo menos três motivos: (1) enquanto talvez exista algum argumento que justifique uma tendência cúbica nos dados do PIB, é difícil imaginar algum para explicar por que o PIB exibe uma tendência polinomial de ordem 10, 17, 21, etc; (2) ordens altas geralmente implicam num sobreajuste dos dados, isto é, o modelo vai se tornar muito bom em explicar o conjunto atual de dados, mas será péssimo para fazer previsões; (3) cada parâmetro do modelo é estimado com algum erro, então quanto mais parâmetros forem incluídos no modelo, maior será este erro.

O gráfico abaixo mostra o ajuste com polinômios de ordens diferentes. Note como as tendências com ordens de polinômios mais elevados apresentam um ajuste melhor aos dados. Isto é especialmente verdadeiro no começo e no final da série.

Exemplo: tendência e sazonalidade (aditiva) na produção industrial

O código abaixo estima um modelo como da equação acima para o índice de produção industrial mensal do FED. Os resultados estão apresentados na tabela. A regressão encontra efeitos significativos e positivos para os meses de junho, setembro e outubro.

data(prodn)
fit <- tslm(prodn ~ trend + season)
Dependent variable:
Production
Trend 0.291***
(0.003)
Feb 1.806
(1.610)
Mar 2.128
(1.610)
Apr 1.950
(1.610)
May 1.936
(1.610)
Jun 3.649**
(1.610)
Jul -1.429
(1.610)
Aug 1.776
(1.611)
Sep 4.105**
(1.611)
Oct 4.207***
(1.611)
Nov 2.294
(1.611)
Dec -0.320
(1.611)
Intercept 28.727***
(1.266)
Observations 372
R2 0.962
Adjusted R2 0.961
Residual Std. Error 6.340 (df = 359)
F Statistic 755.242*** (df = 12; 359)
Note: p<0.1; p<0.05; p<0.01

Ao visualizar o ajuste dos dados vê-se também a limitação desta abordagem. Implicitamente se assume que o efeito sazonal é sempre o mesmo.

autoplot(prodn) +
  autolayer(fitted(fit))

Exemplo: tendência e sazonalidade (multiplicativa)

Com os mesmos dados do exemplo anterior queremos estimar um modelo da forma

\[ y_{t} = T_{t}\cdot S_{t}\cdot u_{t} \]

O código abaixo estima este modelo usando a função decompose do pacote base do R e apresenta os resultados visualmente. O primeiro quadrante do gráfico mostra a série original. Os quadrantes seguintes mostram a parte sazonal, a tendência (polinomial) e o resíduo, i.e., a série dessazonalidada.

fit <- decompose(prodn, type = "multiplicative") 
autoplot(fit)

forecast::Acf(na.omit(fit)$random)

autoplot(stl(prodn, s.window = "periodic"))

Regressão com variáveis defasadas

Uma especificação bastante comum é descrever o comportamento de uma variável em função de seus valores passados. Isto pode refletir, por exemplo, algum componente inercial nos dados em que o valor da inflação do mês de outubro está correlacionada ao valor da inflação no mês de setembro. Também é possível estender este modelo incluindo os valores contemporâneos e defasados de outras variáveis. Pode-se especular, por exemplo, que o desempenho das vendas de um setor seja função da renda no mesmo período (mês, trimestre, etc.) e em períodos anteriores. Ou seja, queremos “explicar” a série \(y\) em função de seus valores passados (\(y_{t-1}, y_{t-2}, \dots\)) e dos valores contemporâneos e passados de outras séries (\(x_{t-1}, z_{t-1}, \dots\)).

Exemplo simples

Um modelo para explicar \(y\) em função de seu valor defasado em um período, do valor contemporâneo de \(x\) e do valor defasado de \(x\) em um período.

\[ y_{t} = \beta_{0} + \beta_{1}y_{t - 1} + \beta_{2}x_{t} + \beta_{3}x_{t - 1} + w_{t} \]

Note que o modelo acima não seria muito útil para gerar previsões de \(y_{t + 1}\) pois ele exigiria conhecimento de \(x_{t + 1}\). Então, seria necessário primeiro prever o valor de \(x_{t + 1}\) para computar uma estimativa para \(y_{t + 1}\).

\[ \mathbb{E}(y_{t + 1} | \mathbb{I}_{t}) = \beta_{0} + \beta_{1}y_{t} + \beta_{2}x_{t + 1} + \beta_{3}x_{t} \]

Exemplo: Índice de Produção Industrial

Para o exemplo abaixo uso o pacote ribge para carregar a série do Índice de Produção Industrial (IPI). Este pacote é hospedado num Github (ao invés do CRAN) e então é preciso usar devtools::install_github() para instalar o pacote. Depois de instalado, ele funciona como qualquer outro pacote.

library(ribge)
# Baixa os dados
ipi <- bcb_serie_carrega(21859)
# Converte a série para ts
prod <- ts(ipi$valor, start = c(2002, 01), frequency = 12)
# Gráfico da série
autoplot(prod) + labs(title = "Índice de Produção Industiral - geral (2012 = 100)")

Pode-se visualizar a relação linear entre valores correntes e defasados do IPI usando a função lag.plot. Na imagem abaixo, cada quadrado mostra um gráfico de dispersão dos valores do IPI em \(t\) contra os valores do IPI em \(t-k\). Alguns lags parecem não exibir muita relação como o lag 6. Já o primeiro e último lag parecem apresentar uma relação linear mais acentuada.

gglagplot(prod, 12, do.lines = FALSE, colour = FALSE)

Como exemplo, podemos propor o modelo abaixo para o IPI. A escolha dos lags aqui foi um tanto arbitrária e há métodos mais apropriados para escolhê-los.

\[ IPI_{t} = \beta_{0} + \beta_{1}IPI_{t - 1} + \beta_{2}IPI_{t - 2} + \beta_{3}IPI_{t - 4} + \beta_{2}IPI_{t - 12} \]

Para estimar este modelo no R há um pequeno inconveniente da função lag que, na verdade, funciona como um operador foreward. Além disso, agora a função base lm (e por conseguinte, também a função tslm) se prova um tanto inconveniente, pois ela não funciona bem com variáveis defasadas. Para usar a função lm seria necessário primeiro “emparelhar” as diferentes defasagens da série, isto é, seria necessário criar um data.frame (ou ts) em que cada coluna mostra os valores das defasagens escolhidas. Por motivo de completude, deixo um código de exemplo que faz isto. Na prática, vale mais a pena escolher alguma outra função como dynlm::dynlm ou dyn::dyn$lm. O código abaixo usa o forecast::tslm, mas nos exemplos seguintes uso o dyn::dyn$lm.

# Exemplo usando forecast::tslm (tb funcionaria com stats::lm)
df <- ts.intersect(prod, lag(prod, -1), lag(prod, -2), lag(prod, -4),
                  lag(prod, -12), dframe = TRUE)
fit <- tslm(prod ~ ., data = df)

Pode-se contornar o problema da função lag definindo uma nova função, L, que funciona da maneira desejada. O código abaixo estima a regressão usando dyn::dyn$lm. A sintaxe dentro da função é praticamente idêntica à que vimos acima com as funções lm e tslm.

# Define uma função L
L <- function(x, k) {lag(x, -k)}
fit <- dyn$lm(prod ~ L(prod, 1) + L(prod, 2) + L(prod, 4) + L(prod, 12))

Os resultados da regressão acima estão resumidos na tabela abaixo.

Dependent variable:
ipi
ipi [1] 0.409***
(0.062)
ipi [4] 0.168**
(0.067)
ipi [8] -0.123***
(0.046)
ipi [12] 0.452***
(0.044)
constante 9.044**
(3.718)
Observations 199
R2 0.786
Adjusted R2 0.782
Residual Std. Error 4.457 (df = 194)
F Statistic 178.551*** (df = 4; 194)
Note: p<0.1; p<0.05; p<0.01

Podemos combinar a informação de outros indicadores industriais para adicionar informação potencialmente relevante para nossa regressão. Neste exemplo, uso outros indicadores industriais para encontrar aqueles que “ajudam a explicar” o indicador geral. O pacote ribge funciona melhor para importação de séries individuais, mas é fácil contornar isto com um loop. Note que no código abaixo não tive muito cuidado com eficiência pois estamos importando apenas 9 séries. O loop importa cada série individualmente, converte ela em um objeto xts e grava ela numa lista. No final a lista l tem cada uma das nove séries como elementos. Por fim uso o do.call para aplicar a função merge em todos os elementos da lista l. Isto é, junto todas as séries num único xts. O resto do código renomeia as variáveis manualmente e depois seleciona alguma das variáveis para a análise. Os códigos utilizados são diretamente copiados do sistema de séries temporais do BCB.

codigos = c(21859, 21861:21868)

l = list()
# Loop para baixar cada série
for (i in 1:length(codigos)) {
  x = bcb_serie_carrega(codigos[i])
  datas = as.Date(x$data, format = "%d/%m/%Y")
  index = as.POSIXct(datas)
  serie = xts(x[, 2], order.by = index)
  l[[as.character(codigos[i])]] = serie
}
indus_ts = do.call(merge, l)
nomes_var =
  c("geral", "extrativa_mineral", "transformacao", "capital", "intermediarios",
    "consumo", "consumo_duraveis", "semiduraveis_e_nao_duraveis",
    "insumos_da_construcao_civil")
names(indus_ts) = nomes_var
df <- indus_ts[, c("transformacao", "extrativa_mineral",
                "capital", "intermediarios", "consumo",
                "consumo_duraveis", "semiduraveis_e_nao_duraveis")]
mes_inicio = month(index(df))[1]
ano_inicio = year(index(df))[1]
series <- ts(as.matrix(df), start = c(ano_inicio, mes_inicio), frequency = 12)
series <- na.omit(series)

Feito o trabalho de importação dos dados podemos propor um modelo simples que toma o valor defasado das variáveis. Novamente a escolha das defasagens e dos regressores foi completamente arbitrária. Discuto formas melhores de escolha neste outro post. O modelo estimado usa defasagens de outras séries para modelar o comportamento da série do índice de produção da indústria de tranformação

\[ Transf_{t} = \beta_{0} + \beta_{1}Durav_{t - 1} + \beta_{2}Durav_{t - 12} + \beta_{3}Capital_{t - 1} + \beta_{4}Capital_{t - 6} + \beta_{5}Transf_{t - 1} + \alpha_{1}t + \sum_{k = 2}^{12}\alpha_{k}d_{k} \]

fit <- dyn$lm(transformacao ~ L(consumo_duraveis, 1) + L(consumo_duraveis, 12) +
                              L(capital, 1) + L(capital, 6) +
                              L(intermediarios, 12) +
                              L(transformacao, 1) +
                              time(transformacao) + as.factor(cycle(transformacao)),
              data = series)
autoplot(series[, 1]) +
  autolayer(fitted(fit))

prod_recente <- window(prod, start = c(2015, 1))
summary(fit <- tslm(prod_recente ~ trend + season))
## 
## Call:
## tslm(formula = prod_recente ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3932 -2.1167  0.5397  2.0039  6.9737 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 82.23575    1.57774  52.122  < 2e-16 ***
## trend       -0.04943    0.02721  -1.816 0.076466 .  
## season2     -2.61057    2.01333  -1.297 0.201832    
## season3      5.23886    2.01388   2.601 0.012768 *  
## season4      3.32829    2.01480   1.652 0.106008    
## season5      8.13772    2.01608   4.036 0.000225 ***
## season6      8.32715    2.01774   4.127 0.000170 ***
## season7     12.57658    2.01975   6.227 1.87e-07 ***
## season8     15.42443    2.13543   7.223 6.98e-09 ***
## season9     11.24886    2.13595   5.266 4.46e-06 ***
## season10    13.72329    2.13682   6.422 9.80e-08 ***
## season11     7.77272    2.13803   3.635 0.000751 ***
## season12    -2.40285    2.13959  -1.123 0.267798    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.183 on 42 degrees of freedom
## Multiple R-squared:  0.8137, Adjusted R-squared:  0.7605 
## F-statistic: 15.29 on 12 and 42 DF,  p-value: 1.168e-11
plot.ts(prod_recente)
lines(fitted(fit), col = 2)

Dummies

Na análise de sazonalidade já lidamos com dummies, mas vale a pena discutí-las num contexto mais amplo. Pode-se usar variáveis dummy para estimar o efeito de certos eventos extraordinários (greves, atentados, crise externa, etc.) ou mesmo o efeito de eventos recorrentes e previsíveis (feriados, dias úteis, alguma estação do ano, etc.).

Exemplo: Efeito de greves e preços administrados no IPCA

\[ \text{IPCA}_{t} = \beta_{0} + \beta_{1}\text{Greve_cam}_{t} + \beta_{2}\text{Greve_2013}_{t} + \beta_{3}\text{Preços_adm}_{t} + \text{Dummies_sazonais} + \epsilon_{t} \]

ipca <- bcb_serie_carrega(433)
ipca$data <- as.Date(ipca$data, format = "%d/%m/%Y")
index <- as.POSIXct(ipca$data)
ipca <-  xts(x = ipca$valor, order.by = index)
ipca <- ipca["2003-05/"]
# dummy greve caminhoneiro
greve_caminhao <- xts(rep(0, length(ipca)), order.by = index(ipca))
greve_caminhao["2018-05"] <- 1

# greves 2013
greve_2013 <- xts(rep(0, length(ipca)), order.by = index(ipca))
greve_2013["2018-07"] <- 1

# precos_adm
precos_adm <- xts(rep(0, length(ipca)), order.by = index(ipca))
precos_adm["2015-01/2015-04"] <- 1
x = cbind(ipca, greve_caminhao, greve_2013, precos_adm)
x = ts(coredata(x), start = c(year(index(x))[1], month(index(x))[1]), frequency = 12)

ipca = x[, 1]
greve_caminhao = x[, 2]
greve_2013 = x[, 3]
precos_adm = x[, 4]

tslm(ipca ~ greve_caminhao + greve_2013 + precos_adm + season)
## 
## Call:
## tslm(formula = ipca ~ greve_caminhao + greve_2013 + precos_adm + 
##     season)
## 
## Coefficients:
##    (Intercept)  greve_caminhao      greve_2013      precos_adm  
##       0.568698        1.040000       -0.366250        0.510833  
##        season2         season3         season4         season5  
##      -0.081250       -0.089375       -0.157500       -0.348698  
##        season6         season7         season8         season9  
##      -0.268698       -0.292448       -0.188698       -0.096823  
##       season10        season11        season12  
##      -0.112448       -0.001198        0.083802
fit <- dyn$lm(ipca ~ L(ipca, 1) + greve_caminhao + greve_2013 + precos_adm +
                     as.factor(cycle(ipca)))

Na tabela abaixo omito os interceptos das dummies sazonais

Dependent variable:
IPCA
ipca [1] 0.508***
(0.061)
Greve (2018) 1.039***
(0.205)
Greve (2013) -0.382*
(0.204)
Preços adm 0.212*
(0.108)
constante 0.108
(0.070)
Constant 0.256***
(0.062)
Observations 194
R2 0.522
Adjusted R2 0.482
Residual Std. Error 0.198 (df = 178)
F Statistic 12.982*** (df = 15; 178)
Note: p<0.1; p<0.05; p<0.01

Previsão (alguns problemas)

Modelos de regressão linear podem também ser usados para fazer previsões. A princípio os valores futuros de uma série podem ser computados usando os valores estimados dos \(\beta_{i}\). Num modelo linear de \(y_{t}\) em função de \(k\) séries de tempo, como o abaixo, podemos calcular o valor de \(\hat{y_{t}}\)

\[ \hat{y_{t}} = \hat{\beta_{0}} + \hat{\beta_{1}}x_{1, t} + \hat{\beta_{2}}x_{2, t} + \dots + \hat{\beta_{k}}x_{k, t} \]

contudo, o valor de \(y_{t + 1}\) exige o conhecimento do valor de todas as \(k\) séries em \(t+1\). O mesmo vale quando usamos variáveis defasadas como regressores. Há diversas formas de resolver este problema. Uma solução simples seria tomar o valor médio dos regressores para formar as previsões. Outra solução seria propor modelos ARMA simples para cada uma das séries \(x_{k, t}\) e usar estas previsões dentro do modelo linear para \(y_{t}\). Alternativamente, pode-se também montar cenários futuros hipotéticos definindo valores pré-estabelecidos para as variáveis regressoras.

Há também alguns casos mais simples, em que o valor futuro das variáveis regressoras é conhecido. Por exemplo, pode-se modelar a demanda por energia elétrica numa região como função de uma dummy de estação (verão x não-verão) e de dia útil (dia útil x dia não-útil). É bastante simples fazer previsões futuras neste caso pois sabe-se de antemão todos os valores futuros das séries que se está usando para “explicar” a demanda por energia elétrica. Já se o mesmo modelo considerasse também o PIB desta região como variável na regressão, seria mais complexo prever os valores futuros da demanda por energia, pois seria necessário saber também os valores futuros do PIB.

Exemplo: tendência e sazonalidade

É relativamente simples prever os valores futuros de modelos de tendência e sazonalidade determinísticas. Como o adjetivo “determinístico” sugere sabe-se de antemão todos os valores que esta série vai exibir. O exemplo abaixo estima um modelo simples para a demanda por passagens aéreas (voos internacionais).

Vale notar que não se costuma fazer previsões de longo prazo com este tipo de modelo, pois a hipótese de que a sazonalidade/tendência continua exatamente igual ao longo do tempo vai se tornando cada vez mais frágil. A curto prazo, contudo, pode ser razoável supor que este modelo linear simples ofereça uma boa aproximação da realidade.

fit <- tslm(AirPassengers ~ trend + season)
autoplot(forecast(fit, h = 24), include = 24)

Exemplo: previsão de cenário

No caso da regressão acima do IPCA, pode-se estimar o impacto de uma nova greve dos caminhoneiros. Claro que é preciso cuidado

dummies = cbind(greve_caminhao, greve_2013, precos_adm)
fit <- Arima(ipca, order = c(1, 0, 0), xreg = coredata(dummies))
greve_caminhao_2020 = c(rep(0, 9), 1, 0, 0)
novo_xreg = cbind(greve_caminhao = greve_caminhao_2020, greve_2013 = rep(0, 12), precos_adm = rep(0, 12))
autoplot(forecast(fit, xreg = novo_xreg), include = 20)

Exemplo: variáveis defasadas

Pode ser um tanto difícil fazer previsões com modelos que usam informação de outras séries. Num modelo simples como \(y_{t} = \beta_{0} + \beta_{1}x_{t - 1}\) para prever valores futuros de \(y_{t}\) é preciso fazer previsõs para a série \(x_{t}\), pois, \(y_{t + 2}\) é função linear de \(x_{t + 1}\). Há muitas maneiras de abordar este problema e eu provavelmente vou discutir mais sobre as alternativas num post futuro. O exemplo abaixo mostra como usar informação disponível de outras séries

df_ajustada <-
  ts.intersect(transf = series[, "transformacao"],
               lag(series[, "consumo_duraveis"], -1),
               lag(series[, "consumo_duraveis"], -12),
               lag(series[, "capital"], -1),
               lag(series[, "capital"], -6),
               lag(series[, "intermediarios"], -12),
               lag(series[, "transformacao"], -1),
               dframe = TRUE
               )
fit <- tslm(transf ~ ., data = df_ajustada)
sub <- df_ajustada[(length(df_ajustada[, "transf"]) - 12):length(df_ajustada[, "transf"]), ]
autoplot(forecast(fit, sub), include = 36)

Pode-se, como de costume, separar os dados em train e test para avaliar a qualidade das previsões dentro da amostra.