MODELAGEM E PREVISÃO DE SÉRIE ECONÔMICA UTILIZANDO MODELAGEM BOX & JENKINS

Macroeconometria - Mestrado Profissional em Economia - Insper

Professora Dra. Juliana Inhasz

Lucas Marin Avelleda

Parte I: Pacotes Utilizados

Foram utilizados os pacotes fredr, zoo, xts, lubridate, ggplot2, tidyr, reshape2,dplyr,kableExtra,knitr,FitAR,seastests,dplyr,aTSA,urca e forecast.

Parte II: Download e tratamento dos dados

Para a obtenção dos dados, utilizei o pacote fredr, o qual possibilita o downolad de séries disponíveis na base de dados do FRED (Federal Reserve Bank of St. Louis) diretamente através do R. Para tanto, é necessário um cadastro no site e a obtenção de uma API, a qual precisa ser validada antes da obtenção dos dados. Por se tratar de uma informação confidencial, essa parte do código foi ocultada deste script. No entanto, o processo de obtenção da API é simples e imediato.

A série escolhida para análise, modelagem e previsão nesse trabalho é a “Industrial Production: Total Index”, cujo código no FRED é “IPB50001N”. Trata-se de um indicador que mensura o produto real de todas as indústrias localizadas nos Estados Unidos (excluindo os territórios). Mais informações sobre a série e sua mensuração podem ser encontradas no link https://fred.stlouisfed.org/series/IPB50001N. Cabe notar que essa série não possui ajuste sazonal (em muitos casos de séries no FRED há tal ajuste).

O processo de obtenção dos dados é simples, conforme as linhas de código abaixo. Basta informar o código da série à função fredr e a data de início desejada.

indpro<-fredr(series_id = "IPB50001N",
                observation_start = as.Date("1990-01-01"))
str(indpro)
## Classes 'tbl_df', 'tbl' and 'data.frame':    362 obs. of  3 variables:
##  $ date     : Date, format: "1990-01-01" "1990-02-01" ...
##  $ series_id: chr  "IPB50001N" "IPB50001N" "IPB50001N" "IPB50001N" ...
##  $ value    : num  63.2 64.1 64.5 63.2 63.8 ...

Para a análise e tratamento de séries temporais no R, é desejável transformar a base de dados em um formato xts para obter valiosas ferramentas de análises que doravante serão implementadas. Para transformar a série, originalmente em tbl_df, tbl, data.frame para xts, utilizou-se o código abaixo, ordenando os dados pelas datas, as quais precisam ser declaradas como “as.Date”.

indpro_xts<-xts(indpro[,3],order.by=as.Date(indpro$date),frequency = 12)
names(indpro_xts)<-c("INDPRO")

Parte III: Análise Preliminar dos Dados

Com a série agora em formato adequado, é possível extrair, através de funções dos pacotes zoo/xts, informações como as apresentadas na Tabela 1, abaixo.

sumario<-function(x){
  c(Início=as.character(start(x)),Fim=as.character(end(x)),
    Observações=length(indpro_xts),Frequência="Mensal")
}

tabela1<-as.data.frame(sapply(indpro_xts,sumario))

knitr::kable(tabela1,caption="Tabela 1: Informações da Base de Dados",digits=3,align="c", format.args = list(big.mark=",",scientific=FALSE)) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width =TRUE,position = "center")
Tabela 1: Informações da Base de Dados
INDPRO
Início 1990-01-01
Fim 2020-02-01
Observações 362
Frequência Mensal

Também é importante nessa análise inicial observar o gráfico da série. Para tanto, foi utilizado o pacote ggplot2 com o código abaixo. Nota-se que a série possui tendência, aparentemente estocástica (quiçá um passeio aleatório com drift). Através da análise visual, ademais, a não convergência a uma média constante no tempo provê indídicos de não estacionariedade.

gg1<-ggplot(indpro,aes(x=date,y=value))+
  geom_line(size=1.0)+
  theme_minimal()+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
      axis.line = element_line(color="black",size=1),
      panel.border=element_rect(colour="black",fill=NA,size=0.1),
      legend.background = element_blank(),
      legend.box.background = element_rect(colour="black",size=1),
      strip.text = element_text(size = 16, color = "black"),
      axis.title = element_text(color = "black", hjust = 0, face = "italic"),
      axis.text = element_text(color = "black"), 
      plot.title = element_text(face = "bold", size = (10)),
      plot.subtitle = element_text(face = "italic", size = (8)))+
  labs(y="Industrial Production Index",x="Tempo",title="Gráfico 1: Industrial Production Index",subtitle = "Dados obtidos em https://fred.stlouisfed.org/series/IPB50001N ")
gg1

Para extrair informações sobre a sazonalidade, selecionei um período de 4 anos para plotar em um gráfico, apresentado abaixo. Percebem-se no gráfico picos de crescimento nos meses de verão e menor atividade industrial nos meses de inverno, sugerindo sazonalidade. Cabe ainda mencionar que no período selecionado (2016 a 2020) há uma tendência aparentemente determinística na série.

indpro4<-tail(indpro,4*12)

gg2<-ggplot(indpro4,aes(x=date,y=value))+
  geom_line(size=1.0)+
  theme_minimal()+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
      axis.line = element_line(color="black",size=1),
      panel.border=element_rect(colour="black",fill=NA,size=0.1),
      legend.background = element_blank(),
      legend.box.background = element_rect(colour="black",size=1),
      strip.text = element_text(size = 16, color = "black"),
      axis.title = element_text(color = "black", hjust = 0, face = "italic"),
      axis.text = element_text(color = "black"), 
      plot.title = element_text(face = "bold", size = (10)),
      plot.subtitle = element_text(face = "italic", size = (8)))+
  labs(y="Industrial Production Index",x="Tempo",title="Gráfico 2: Industrial Production Index",subtitle = "Dados obtidos em https://fred.stlouisfed.org/series/IPB50001N ")
gg2

No entanto, para obter maior rigor científico acerca da existência da sazonalidade, realizaram-se dois testes de sazonalidade através do código abaixo. O primeiro é o Teste-WO de sazonalidade em séries temporais, cuja hipótese nula é a inexistência de sazonalidade. Como o p-valor registrado é praticamente nulo, rejeita-se a hipótese nula e o teste acaba por informar que identificou sazonalidade na série. Foi realizado também o Teste-QS que ao retornar “TRUE”, confirma a existência de sazonalidade na série.

summary(wo(indpro_xts,freq=12))
## Test used:  WO 
##  
## Test statistic:  1 
## P-value:  0 0 0 
##  
## The WO - test identifies seasonality
isSeasonal(indpro_xts,test="qs",freq=12)
## [1] TRUE

Ademais, de forma a sumarizar o exposto até aqui, realiza-se abaixo a decomposição da série em seus componentes aleatório, sasonal e de tendência, os quais tornam-se deveras evidentes.

indpro_ts<-ts(indpro_xts,frequency=12)
indpro_decomp<-decompose(indpro_ts)
plot(indpro_decomp)

Parte IV: Processo Gerador dos Dados

Para melhor entender o processo gerador dos dados, recorre-se abaixo aos gráficos de Autocorrelação e Autocorrelação Parcial. Uma vez que a ACF decai lentamente e a PACF trunca no lag=2 e há ocorrência de autocorrelações significativas a cada 12 lags (o que já era esperado), um candidato processo gerador da série original é um modelo SAR(p)s, SAR(1)12

plot(acf(indpro_xts))

plot(pacf(indpro_xts))

Entretanto, do exposto inicialmente, sabe-se que aparentemente a série não é estacionária. Testaremos essa hipótese a seguir com o código abaixo, do pacote urca, que reproduz o Augmented Dickey-Fuller Test de raiz unitária para uma série temporal.O teste abaixo é para a série no nível (sem diferenciação).

Sob a primeira hipótese nula do teste (H0: Existência de raiz unitária), conclui-se que não é possível rejeitá-la uma vez que o valor da estatística tau3 encontrado (-2.75) é maior que o valor crítico tauc3 (-3.42). Dessa forma, uma vez que o processo gerador possui raiz dentro do círculo unitário, a série não pode ser considerada estacionária.

ur0<-ur.df(indpro_xts, type="trend", lags=0, selectlags="AIC")
ur0@teststat
##                tau3     phi2     phi3
## statistic -2.748348 3.297445 3.922288
ur0@cval
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36
plot(ur0)

Procedemos, então, com a primeira diferenciação da série e seu respectivo teste de raiz unitária.

Agora, sob a primeira hipótese nula do teste (H0: Existência de raiz unitária), conclui-se que é possível rejeitá-la uma vez que o valor da estatística tau3 encontrado (-31.34) é menor que o valor crítico tauc3 (-3.98). Dessa forma, uma vez que o processo gerador possui raiz fora do círculo unitário, o processo pode ser considerado estacionário.

indpro_1diff_xts<-xts(diff(indpro_xts)[-1])
ur1<-ur.df(indpro_1diff_xts, type="trend", lags=0, selectlags="AIC")
ur1@teststat
##                tau3     phi2     phi3
## statistic -31.34378 327.4809 491.2214
ur1@cval
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36
plot(ur1)

Agora, considerando que a série é estacionária na primeira diferença, podemos novamente buscar seu processo gerador através dos gráficos de ACF e PACF. Percebe-se que tanto a ACF quanto a PACF truncam e apresentam correlação expressiva a cada 12 lags (mantendo a característica sazonal). Portanto, um candidato a processo gerador dessa série é um processo SARIMA(P,Q,D)(p,q,d)s da forma SARIMA(2,1,2)(2,1,2)12.

plot(acf(indpro_1diff_xts))

plot(pacf(indpro_1diff_xts))

Podemos, então, estimar o processo gerador:

genproc1<-arima(indpro_xts, order = c(2, 1, 2),
      seasonal = list(order = c(2, 1, 2), period = 12))
summary(genproc1)
## 
## Call:
## arima(x = indpro_xts, order = c(2, 1, 2), seasonal = list(order = c(2, 1, 2), 
##     period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sar1     sar2     sma1     sma2
##       1.2444  -0.3848  -1.2706  0.5562  -0.5386  -0.1066  -0.0803  -0.3082
## s.e.  0.3435   0.3074   0.3107  0.2281   0.3035   0.0952   0.3051   0.2243
## 
## sigma^2 estimated as 0.5097:  log likelihood = -381.09,  aic = 780.19
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.00492328 0.7010817 0.5215845 0.01061082 0.5697955 0.4003625
##                      ACF1
## Training set -0.007789471

E, em seguida, realizar uma análise dos resíduos. Sua distribuição não aparenta ser gaussiana, o que também é confirmado pelo Teste-Shapiro.

plot(genproc1$residuals)

hist(genproc1$residuals,breaks=100,density = 20,prob=TRUE)

shapiro.test(genproc1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  genproc1$residuals
## W = 0.96537, p-value = 1.459e-07

No entanto, não há indícios de autocorrelação nem direta nem indireta na série dos resíduos, conforme constatam os gráficos de ACF e PACF abaixo.

acf(genproc1$residuals)

pacf(genproc1$residuals)

A tabela abaixo apresenta um sumário estatístico dos resíduos do modelo. Importante notar a média próxima a zero e o fato dela aparentar ser constante, conforme visto mais acima no gráfico da série dos resíduos. Embora não possua distribuição normal, ele se aproxima de um ruido branco.

sumario<-function(x){
  c(mínimo=min(x),média=mean(x),máximo=max(x),mediana=median(x),desvpad=sd(x))
}

tabela1<-sapply(genproc1$residuals,sumario)
knitr::kable(tabela1,caption="Sumário Estatístico dos Resíduos",digits=3,align ="c", format.args = list(big.mark=",",scientific=FALSE)) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE,position = "center")
Sumário Estatístico dos Resíduos
mínimo -4.088
média 0.005
máximo 3.063
mediana -0.004
desvpad 0.702

Parte V: Previsão Utilizou-se o pacote Forecast para gerar a previsão de 12 meses a frente com o processo gerador da série encontrado:

genproc1_forecast<-forecast(genproc1,h=12)
plot(genproc1_forecast,include=12*6)

Por fim, é feita uma comparação gráfica entre os valores originais da série, em preto, e os valores gerados pelo modelo (fitted), em vermelho. Percebe-se que visualmente o modelo escolhido parece se adequar bem à série.

dates<-indpro$date
fitted_prc<-fitted(genproc1)

mari<-as.data.frame(dates)
mari$fitted<-fitted_prc
mari$actual<-indpro$value

gg3<-ggplot(mari,aes(x=dates))+
  geom_line(aes(y=actual), colour="black")+
  geom_line(aes(y=fitted), colour="red")  +
  theme_minimal()+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
      axis.line = element_line(color="black",size=1),
      panel.border=element_rect(colour="black",fill=NA,size=0.1),
      legend.background = element_blank(),
      legend.box.background = element_rect(colour="black",size=1),
      strip.text = element_text(size = 16, color = "black"),
      axis.title = element_text(color = "black", hjust = 0, face = "italic"),
      axis.text = element_text(color = "black"), 
      plot.title = element_text(face = "bold", size = (10)),
      plot.subtitle = element_text(face = "italic", size = (8)))+
  labs(y="Industrial Production Index",x="Tempo",title="Gráfico 3: Industrial Production Index, actual x fitted",subtitle = "Dados obtidos em https://fred.stlouisfed.org/series/IPB50001N ")
gg3