Dessazonalizando o PIB brasileiro pelo R – Parte II

No post anterior, abordamos rapidamente a importância de aplicar o ajustamento sazonal em séries temporais. Agora focaremos apenas nos passos para replicar a metodologia do IBGE para a dessazonalização da série trimestral do PIB brasileiro. Vamos lá?

  1. Baixando e extraindo o X-13-ARIMA-SEATS no site do U.S. Census Bureau
    • A versão do X-13 que precisamos é a “x13ashtmlall.zip” (no momento que escrevo a versão mais atual é “_V1.1_B19”), disponível neste link;
    • Depois de feito o download, extraia os arquivos em qualquer pasta do seu computador (OBS: a pasta será utilizada no código);
    • Para facilitar, manteremos no código do R o nome original da pasta do X-13: “x13ashtml”;
  2. Baixando a série histórica do PIB Total das Contas Nacionais Trimestrais, disponível neste link
    • No site do IBGE, precisamos das “Tabelas Completas”, com link disponível no canto esquerdo da página;
    • Aconselhamos salvar e extrair a planilha na mesma pasta do X-13, pois ela também servirá de output para os dados dessazonalizados pelo código R.
    • Os dados que nos interessam nesta planilha estão na aba “Série Encadeada”, coluna R (referente ao PIB Total, em número índice).
  3. Verificando os parâmetros de dessazonalização do IBGE, disponível neste link
    • O IBGE divulga os parâmetros de dessazonalização do PIB na “Publicação completa” das Contas Nacionais Trimestrais” (em formato PDF, também no canto esquerdo da página). Dentro do documento, procure pelo quadro de “modelos adotados no ajuste sazonal”, disponível no anexo de notas metodológicas.
    • Importante destacar que os parâmetros de ajustamento sazonal não são fixos, de modo que sempre que houver nova divulgação as Contas Nacionais Trimestrais, o IBGE estimará novamente o modelo ARIMA e gerará diferentes outliers (efeitos de intervenção).
      Pronto, já temos todos os elementos para executar a rotina de dessazonalização no R. Mas antes cabem mais algumas ressalvas sobre as linhas do código: (i) atenção para a AMOSTRA que se pretende dessazonalizar – no exemplo abaixo, “ ,2015.1” significa ajustar do início da amostra até o 1º trimestre de 2015 e (ii) atenção para as pastas onde foram salvos o programa X-13 e os dados do IBGE – no exemplo, os caminhos estão definidos como “C:\Users\Tho\Google Drive\Blog\201508\x13ashtml”

O código precisará de três bibliotecas: seasonal, XLConnect e dplyr. Se algum erro surgir na execução do código, tente reinstalar as bibliotecas utilizando o install.packages(“XLConnect”, “dplyr”, “seasonal”)

library(XLConnect)
library(dplyr)

# Especificar o caminho onde foi salvo o programa X-13
# Para rodar da primeira vez, necessita apontar o caminho do x13
Sys.setenv(X13_PATH = "C:\\Users\\Tho\\Google Drive\\Blog\\201508\\x13ashtml")
library(seasonal)

# Mudar a AMOSTRA de acordo com a base de dados
AMOSTRA <- ",2015.1"

#especificar o caminho onde foi salva a planilha "Tab_Compl_CNT_1T15.xls" com a série do PIB. Neste exemplo, utilizamos a mesma pasta do X-13
arq <- 'C:\\Users\\Tho\\Google Drive\\Blog\\201508\\x13ashtml\\Tab_Compl_CNT_1T15.xls'
Original = readWorksheetFromFile(arq, sheet = "Série Encadeada",
                                 startRow = 4, startCol = 18, endCol = 18, useCachedValues = TRUE)

PIBT <- ts(Original[,1], start = 1996, frequency = 4)

plot(PIBT)

plot of chunk unnamed-chunk-2

#abaixo são definidas a decomposição do modelo (mult para Multiplicativa ou add para Aditiva), o modelo ARIMA "(0 1 1)(0 1 1), por exemplo" e os efeitos de intervenção ("TC1996.3", "AO1998.1", ...). Se a decomposição for Aditiva (x11.mode = "add"), colocar transform.function = NULL. Se x11.mode = "mult", então transform.function = "log"

PIBTSA <- seas(x = PIBT,
               series.modelspan = AMOSTRA,
               arima.model = "(0 1 1)(0 1 1)",
               regression.variables = c("LS2008.4", "Easter[1]"),
               forecast.maxlead = 6,
               forecast.lognormal = "yes",
               regression.aictest = NULL,
               outlier = NULL,
               transform.function = "log",
               automdl = NULL,
               x11.mode = "mult")

O resultados finais estão destacados no gráfico e nas tabelas abaixo.

#testes, gráficos e resultados
plot(PIBTSA)

plot of chunk unnamed-chunk-3

summary(PIBTSA)
##
## Call:
## seas(x = PIBT, transform.function = "log", regression.aictest = NULL,
##     outlier = NULL, automdl = NULL, series.span = AMOSTRA, regression.variables = c("TC1996.3",
##         "AO1998.1", "LS2008.4", "TC2009.1", "Easter[15]"), arima.model = "(0 2 2)(0 1 1)",
##     x11.mode = "mult")
##
## Coefficients:
##                    Estimate Std. Error z value Pr(&amp;amp;amp;amp;gt;|z|)
## TC1996.3           0.036471   0.006370   5.726 1.03e-08 ***
## AO1998.1          -0.025803   0.005112  -5.047 4.49e-07 ***
## LS2008.4          -0.053941   0.008096  -6.663 2.69e-11 ***
## TC2009.1          -0.028206   0.006147  -4.589 4.46e-06 ***
## Easter[15]        -0.006106   0.001515  -4.029 5.60e-05 ***
## MA-Nonseasonal-01  0.635040   0.113684   5.586 2.32e-08 ***
## MA-Nonseasonal-02  0.237159   0.113764   2.085   0.0371 *
## MA-Seasonal-04     0.738180   0.092061   8.018 1.07e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## X11 adj.  ARIMA: (0 2 2)(0 1 1)  Obs.: 77  Transform: log
## AICc: 257.1, BIC: 274.6  QS (no seasonality in final):    0
## Box-Ljung (no autocorr.): 27.65   Shapiro (normality): 0.9656 *
final(PIBTSA)
##          Qtr1     Qtr2     Qtr3     Qtr4
## 1996 100.0209 100.3977 104.7880 103.3985
## 1997 103.8248 104.8589 106.6537 107.1998
## 1998 104.3627 106.8397 107.1789 105.7283
## 1999 105.5060 105.9662 106.7138 108.0361
## 2000 109.5622 110.6865 111.6985 112.9109
## 2001 113.2017 113.0359 112.1191 112.2626
## 2002 114.3034 115.0038 116.8060 118.2141
## 2003 116.6877 116.6858 117.6694 119.1089
## 2004 121.2424 123.6196 125.3305 126.5121
## 2005 126.5316 128.4487 128.1028 129.3265
## 2006 130.9385 132.2814 134.0533 135.5494
## 2007 138.0144 140.4844 141.9269 144.2904
## 2008 146.7550 149.0367 151.7199 145.5084
## 2009 142.3680 146.2933 149.7956 153.1179
## 2010 155.9768 158.3498 160.1333 162.0052
## 2011 163.3604 166.4313 165.5810 166.1209
## 2012 166.5143 167.2685 169.3465 169.9862
## 2013 171.2721 173.3229 173.4634 173.5229
## 2014 174.7563 172.3464 172.5598 173.0335
## 2015 172.7567
resid(PIBTSA)
##               Qtr1          Qtr2          Qtr3          Qtr4
## 1996 -7.631667e-04 -2.201111e-03  3.141728e-03  2.574970e-03
## 1997  3.174729e-03  6.212903e-03  9.204532e-03 -7.465839e-03
## 1998 -9.728653e-03 -6.811626e-03 -5.964617e-03 -1.765721e-02
## 1999  4.723447e-05  4.385868e-03 -3.390666e-03  1.649382e-02
## 2000  9.088939e-03  4.091440e-03 -5.358513e-03  1.031379e-02
## 2001 -5.678862e-03 -6.629682e-03 -1.881175e-02  4.415687e-03
## 2002  1.588737e-02 -2.421547e-03  9.973459e-03  4.824630e-03
## 2003 -2.221758e-02  1.505260e-03 -1.293940e-04  6.493953e-03
## 2004  1.705166e-02  9.154115e-03 -1.592008e-03 -1.634295e-03
## 2005 -3.946405e-03  4.692298e-03 -1.941391e-02  5.489389e-03
## 2006  1.143963e-02 -7.772136e-03  4.263996e-03  3.790999e-04
## 2007  1.489767e-02 -3.828167e-03 -2.912046e-03  7.137433e-03
## 2008  5.731716e-03 -7.618234e-03  6.515553e-03 -2.574858e-03
## 2009 -6.548068e-03 -1.033460e-02  6.591033e-03  6.130106e-03
## 2010  3.502380e-03 -5.098606e-03 -3.128757e-03 -3.658016e-03
## 2011 -4.866644e-03  3.833883e-03 -2.047175e-02 -5.243749e-03
## 2012 -5.751823e-03 -7.500150e-03  9.345213e-03 -8.699906e-03
## 2013  4.248774e-03  8.307355e-04 -6.875811e-03 -4.441129e-03
## 2014  4.045059e-03 -2.492166e-02  5.948884e-03  3.909982e-04
## 2015 -4.024602e-03
coef(PIBTSA)
##          TC1996.3          AO1998.1          LS2008.4          TC2009.1
##       0.036471448      -0.025802666      -0.053941263      -0.028206039
##        Easter[15] MA-Nonseasonal-01 MA-Nonseasonal-02    MA-Seasonal-04
##      -0.006105615       0.635040284       0.237159253       0.738180493
fivebestmdl(PIBTSA)
## Error in seas(x = PIBT, transform.function = "log", regression.aictest = NULL, : object 'AMOSTRA' not found
spc(PIBTSA)
## series{
##   title = "PIBT"
##   file = "C:\Users\Tho\AppData\Local\Temp\RtmpgZIDc2/x13/data.dta"
##   format = "datevalue"
##   period = 4
##   span = (,2015.1)
## }
##
## transform{
##   function = log
##   print = aictransform
## }
##
## regression{
##   variables = (TC1996.3 AO1998.1 LS2008.4 TC2009.1 Easter[15])
## }
##
## arima{
##   model = (0 2 2)(0 1 1)
## }
##
## x11{
##   mode = mult
##   save = (d10 d11 d12 d13 d16 e18)
## }
##
## estimate{
##   save = (model estimates lkstats residuals)
## }
##
## spectrum{
##   print = qs
## }

Construindo o output em Excel, sendo que a primeira linha define o vetor de datas

Data &amp;amp;amp;amp;lt;- seq(as.POSIXlt("1996-3-1"), as.POSIXlt("2015-3-1"), by = "quarter")
PIBTSA <- as.numeric(series(PIBTSA, "d11"))
PIBTSA <- as.data.frame(PIBTSA)
PIBT <- as.numeric(PIBT)
PIBT <- as.data.frame(PIBT)

Matriz <- cbind(Data, PIBT, PIBTSA)

# Tire o "#" de comentário na linha abaixo para gravar a série dessazonalizada na nova aba "Dessaz"
# writeWorksheetToFile("Tab_Compl_CNT_1T15.xls", Matriz, sheet = "Dessaz", startCol = 1, clearSheets = TRUE)
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s