Índice de Atividade Econômica do Banco Central do Brasil

library(fpp2)
library(tidyverse)
library(TSstudio)
library(tseries)
library(BETS) 
library(dplyr)

Procurando a identificação

df.res <- BETSsearch("Central Bank Economic Activity Index",
                     view = FALSE)

Verificando o resultado e selecionando o índice

df.res %>%
  kable() %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"))

code	description	unit	periodicity	start	last_value	source
24363	Central Bank Economic Activity Index	Index	M	01/01/2003	feb/2018	BCB-Depec
24364	Central Bank Economic Activity Index (IBC-Br) - seasonally adjusted	Index	M	01/01/2003	feb/2018	BCB-Depec

id <- 24364 # série ajustada sazonalmente

Define as datas de início e término da série

inicio = '2003-01-01' 

fim = as.character(Sys.Date() - 1) # ontem

Obtém e visualiza os dados

df.BETS <- BETSget(code = id, 
                   data.frame = TRUE, 
                   from = inicio, 
                   to = fim)
df <- df.BETS

head(df)

##         date  value
## 1 2003-01-01 100.39
## 2 2003-02-01 102.06
## 3 2003-03-01 102.03
## 4 2003-04-01 100.92
## 5 2003-05-01  99.76
## 6 2003-06-01 100.31

tail(df)

##           date  value
## 206 2020-02-01 140.01
## 207 2020-03-01 131.76
## 208 2020-04-01 119.41
## 209 2020-05-01 121.63
## 210 2020-06-01 128.10
## 211 2020-07-01 130.85

ibc <- df$value

ts.ibc <- ts(ibc, frequency =  12, start = c(2003,01), end = c(2020,06))
ts.ibc

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2003 100.39 102.06 102.03 100.92  99.76 100.31  99.06  99.72 101.77 101.85
## 2004 103.51 104.52 105.99 106.43 106.31 107.25 107.82 108.20 109.41 108.71
## 2005 109.54 110.21 110.45 110.98 110.95 110.60 110.63 111.26 110.91 111.09
## 2006 113.37 113.63 113.54 114.31 115.63 114.70 116.94 116.60 116.83 118.09
## 2007 119.65 120.18 119.42 121.24 122.33 123.38 124.26 124.26 125.35 126.29
## 2008 127.23 126.75 126.46 127.80 129.41 131.13 130.93 130.99 130.92 128.12
## 2009 120.90 122.19 123.01 123.64 124.88 125.92 126.23 128.09 129.08 129.80
## 2010 133.61 135.18 136.57 137.17 136.46 136.19 136.94 137.71 139.26 138.67
## 2011 140.22 140.77 141.08 140.88 141.30 142.14 142.25 141.95 141.79 141.64
## 2012 138.97 140.38 140.06 141.02 142.61 143.55 144.05 144.29 143.70 144.88
## 2013 144.96 143.82 144.94 146.29 146.81 146.33 147.18 147.14 148.05 148.01
## 2014 148.29 147.99 147.80 146.95 146.25 143.55 144.84 145.61 146.19 145.41
## 2015 145.28 144.33 144.07 142.46 141.19 139.93 138.98 139.12 137.84 138.74
## 2016 135.65 135.59 134.67 134.72 134.16 134.67 134.43 133.83 133.65 133.20
## 2017 133.73 136.01 135.84 135.86 134.99 136.10 136.44 135.62 135.44 135.56
## 2018 137.35 137.36 137.35 138.45 133.05 137.63 137.95 139.02 137.92 137.73
## 2019 139.29 137.76 137.54 137.32 138.19 138.92 137.82 138.42 139.05 139.52
## 2020 139.30 140.01 131.76 119.41 121.63 128.10                            
##         Nov    Dec
## 2003 101.97 101.64
## 2004 109.62 108.79
## 2005 111.81 112.36
## 2006 118.86 119.91
## 2007 126.18 126.90
## 2008 124.96 120.88
## 2009 130.03 131.26
## 2010 140.11 138.17
## 2011 142.28 140.43
## 2012 144.61 144.11
## 2013 148.37 148.65
## 2014 145.64 146.10
## 2015 136.91 136.21
## 2016 133.07 132.09
## 2017 135.89 137.47
## 2018 138.32 138.63
## 2019 139.53 139.17
## 2020

autoplot(ts.ibc, lwd = 1.2, col = "steelblue") +
  labs(y = 'IBC-BR', 
       x = 'Data',
       title = 'Índice de Atividade do Banco Central do Brasil',
       caption = 'Fonte: BACEN')

O impacto da Covid-19 em números

ibc.2020 <- window(ts.ibc, 2020)

covid <- data.frame(Data = time(ibc.2020), 
                    Valor = as.matrix(ibc.2020))

arit.ret <- function(precos) {
####  Calcula o retorno aritmético de uma série de preços ####    
    c(NA, precos[2:length(precos)]/ precos[1:length(precos)-1] - 1)
}


covid$delta <- arit.ret(as.numeric(covid$Valor)) * 100

covid

##       Data  Valor      delta
## 1 2020.000 139.30         NA
## 2 2020.083 140.01  0.5096913
## 3 2020.167 131.76 -5.8924363
## 4 2020.250 119.41 -9.3731026
## 5 2020.333 121.63  1.8591408
## 6 2020.417 128.10  5.3194113

autoplot(ibc.2020, lwd = 1.2, col = "steelblue") +
  labs(y = 'IBC-BR', 
       x = 'Data',
       title = 'Índice de Atividade do Banco Central do Brasil em 2020',
       caption = 'Fonte: BACEN')

Usando o pacote `TSstudio`

ts_seasonal(ts.ibc, 
            type = 'box', 
            title = 'Sazonalidade - Índice de Atividade do BCB')

Testes

acf <- ggAcf(ts.ibc)
acf

pacf <- ggPacf(ts.ibc)
pacf

A série é claramente não estacionária. O Pacf sugere que a séria ficaria estacionária no lag = 3.

Utlizaremos os testes de estacionariedade com o pacote tseries.

Teste ADF - no nível:

$H_0$: A série não é estacionária

$H_1$: A série é estacionária

adf.test(ts.ibc)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts.ibc
## Dickey-Fuller = -0.5433, Lag order = 5, p-value = 0.9791
## alternative hypothesis: stationary

Rejeita-se Ho. p-valor: 0.756
Note que o ADF está sugerindo um lag = 5

Teste KPSS

Note a diferença da hipótese básica do teste

$H_0$: A série é estacionária
$H_1$: A série não é estacionária

kpss.test(ts.ibc)

## 
##  KPSS Test for Level Stationarity
## 
## data:  ts.ibc
## KPSS Level = 2.9662, Truncation lag parameter = 4, p-value = 0.01

Rejeita-se a hipótese nula (p-value smaller than printed p-value).

Teste Philipps - Perron

$H_0$: A série não é estacionária
$H_1$: A série é estacionária

pp.test(ts.ibc)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  ts.ibc
## Dickey-Fuller Z(alpha) = -1.6139, Truncation lag parameter = 4, p-value
## = 0.9774
## alternative hypothesis: stationary

Não se pode rejeitar a hipótese nula a menos de $\alpha$ = 21%: p-valor = 0.2136

Efetuando a diferença

diff <- diff(ts.ibc)

lim <- 1.96/sqrt(length(ibc))

autoplot(diff, col ="blue", lwd = 1.1) +
  geom_hline(yintercept = 0,
             colour = 'darkred', 
             linetype = 'solid', 
             size = 0.75) +
  xlab('Data') +
  ylab('Diff')

Efetuando o correlograma

acf.diff <- ggAcf(diff, lag.max = 24)
acf.diff

Uma diferenciação de lag = 2, mas observa-se a sazonalidade nos instantes 6 e 12 e 18 e 24. O teste sugere 4 lags (5-1).

O Acf é útil para medir o Modelo MA(q).

pacf.diff <- ggPacf(diff, lag.max = 24)
pacf.diff

Conte os primeiros spikes do PACF para o lag do AR(p).

Testes com a série diferenciada

Teste ADF - na diferença

$H_0$: A série não é estacionária
$H_1$: A série é estacionária

adf.test(diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff
## Dickey-Fuller = -5.8767, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Teste KPSS na diferença

$H_0$: A série é estacionária
$H_1$: A série não é estacionária

kpss.test(diff)

## 
##  KPSS Test for Level Stationarity
## 
## data:  diff
## KPSS Level = 0.58481, Truncation lag parameter = 4, p-value = 0.02402

Teste Philipps - Perron na diferença

$H_0$: A série não é estacionária
$H_1$: A série é estacionária

pp.test(diff)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff
## Dickey-Fuller Z(alpha) = -145.61, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary

AR(p)

ar <- arima(ts.ibc, order = c(2,1,0))
ar

## 
## Call:
## arima(x = ts.ibc, order = c(2, 1, 0))
## 
## Coefficients:
##          ar1      ar2
##       0.2421  -0.1646
## s.e.  0.0710   0.0711
## 
## sigma^2 estimated as 2.343:  log likelihood = -385.59,  aic = 777.18

Efetuando as previsões com o modelo AR(p)

fit <- forecast(ar, h = 12) #  IC level = 95

autoplot(ts.ibc, series = "Série_Observada") +
    autolayer(fit$fitted, series = "Série Predita_AR") +
    autolayer(fit, series = "Projeção h = 12") +
    xlab('Data')

MA(q)

ma <- arima(ts.ibc, order = c(0,1,2))
ma

## 
## Call:
## arima(x = ts.ibc, order = c(0, 1, 2))
## 
## Coefficients:
##          ma1      ma2
##       0.2228  -0.0512
## s.e.  0.0772   0.0773
## 
## sigma^2 estimated as 2.372:  log likelihood = -386.86,  aic = 779.71

Efetuando as previsões com o modelo MA(q)

fit <- forecast(ma, h = 12) # SE QUISER FIXAR IC level = 95)

autoplot(ts.ibc, series = "Série_Observada") +
    autolayer(fit$fitted, series = "Série Predita_MA") +
    autolayer(fit, series = "Projeção h = 12") +
    xlab('Data')

#$ ARIMA(p,d,q)

arima <- arima(ts.ibc, order = c(2,1,2))
arima

## 
## Call:
## arima(x = ts.ibc, order = c(2, 1, 2))
## 
## Coefficients:
##          ar1      ar2     ma1     ma2
##       0.0112  -0.8059  0.2272  0.7382
## s.e.  0.3419   0.2230  0.3075  0.3748
## 
## sigma^2 estimated as 2.275:  log likelihood = -382.75,  aic = 775.49

Efetuando as previsões com o modelo ARIMA(p,d,q)

fit <- forecast(arima, h = 12) # IC 95% padrão

autoplot(ts.ibc, series = "Série_Observada") +
  autolayer(fit$fitted, series = "Série Predita_ARIMA") +
  autolayer(fit, series = "Projeção h = 12") +
  xlab('Data')

sarima <-  arima(ts.ibc, 
                 order= c(2,1,0),
                 seasonal = list(order = c(0,1,1), 
                 period = 12))    

sarima

## 
## Call:
## arima(x = ts.ibc, order = c(2, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##          ar1      ar2     sma1
##       0.1953  -0.1432  -0.7691
## s.e.  0.0738   0.0730   0.0814
## 
## sigma^2 estimated as 2.579:  log likelihood = -378.25,  aic = 764.51

Efetuando as previsões com SARIMA

fit <- forecast(sarima, h = 12) 

autoplot(ts.ibc, series = "Serie_Observada") +
  autolayer(fit$fitted, series = "Série Predita_ SARIMA") +
  autolayer(fit, series = "Projeção h = 12", showgap = FALSE) +
  xlab('Data')

Usando o Auto.Arima

auto.arima <- auto.arima(ts.ibc, seasonal = TRUE,
                         stepwise = FALSE, 
                         approximation = FALSE)
auto.arima

## Series: ts.ibc 
## ARIMA(0,2,4) 
## 
## Coefficients:
##           ma1      ma2      ma3     ma4
##       -0.7566  -0.3146  -0.1080  0.2055
## s.e.   0.0681   0.0855   0.0946  0.0722
## 
## sigma^2 estimated as 2.34:  log likelihood=-383.14
## AIC=776.28   AICc=776.58   BIC=792.97

print('Modelo -> SARIMA(3,1,0)(0,1,1)[12]')

## [1] "Modelo -> SARIMA(3,1,0)(0,1,1)[12]"

fit <- forecast(auto.arima, h = 12) 

autoplot(ts.ibc, series = "Série_Observada") +
  autolayer(fit$fitted, series = "Série Predita_AUTO.ARIMA") +
  autolayer(fit, series = "Projeção h = 12") +
  xlab('Data')

round(accuracy(ar),  4)

##                  ME   RMSE    MAE    MPE   MAPE   MASE   ACF1
## Training set 0.1229 1.5271 0.9583 0.1015 0.7473 1.0108 -0.017

round(accuracy(ma),  4)

##                  ME   RMSE    MAE    MPE   MAPE   MASE   ACF1
## Training set 0.1157 1.5365 0.9543 0.0955 0.7444 1.0066 0.0027

round(accuracy(arima), 4)

##                  ME   RMSE   MAE    MPE   MAPE   MASE  ACF1
## Training set 0.1187 1.5048 0.925 0.0984 0.7235 0.9756 7e-04

round(accuracy(sarima), 4)

##                  ME   RMSE    MAE     MPE   MAPE   MASE    ACF1
## Training set -0.123 1.5556 0.9669 -0.0916 0.7419 1.0199 -0.0129

round(accuracy(auto.arima), 4)

##                   ME   RMSE    MAE     MPE   MAPE   MASE    ACF1
## Training set -0.1099 1.5077 0.9185 -0.0832 0.7145 0.1946 -0.0071

Melhor modelo encontrado (pela função `auto.arima`): SARIMA(3,1,0)(0,1,1)[12]

Checando a significância estatística dos parâmetros

t_test(model = auto.arima)

##         Coeffs Std.Errors         t Crit.Values Rej.H0
## ma1 -0.7565757 0.06810458 11.109028    1.971661   TRUE
## ma2 -0.3145982 0.08554356  3.677637    1.971661   TRUE
## ma3 -0.1080443 0.09458919  1.142248    1.971661  FALSE
## ma4  0.2054761 0.07221053  2.845514    1.971661   TRUE

Os parâmetros são significativos.

Diagnóstico

Teste de ausência de autocorrelação serial

Teste Ljung-Box

$H_0$: Os resíduos são independentes e identicamente distribuídos (inexiste aurocorrelação), ou seja, não há dependência serial para todas as defasagens.
$H_1$: não $H_0$

O modelo SARIMA(3,1,0)(0,1,1)[12]

IMPORTANTE:

\[fitdf = p + q + P + Q = 3 + 0 + 0 + 1 = 4\]

Box.test(x = auto.arima$residuals, 
         lag = 24, 
         type = "Ljung-Box", 
         fitdf = 4)

## 
##  Box-Ljung test
## 
## data:  auto.arima$residuals
## X-squared = 15.607, df = 20, p-value = 0.7407

Não podemos rejeitar $H_0$. A série aparenta ser iid.

Resíduos

checkresiduals(auto.arima, col ="blue")

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,2,4)
## Q* = 15.607, df = 20, p-value = 0.7407
## 
## Model df: 4.   Total lags used: 24

Checando a Heterocedasticidade

Teste proposto por Engle (1982) - GARCH

$H_0$: Os resíduos ao quadrado são uma sequência de ruídos brancos, ou seja, os resíduos são homocedásticos.
$H_1$: A série é heterocedástica

library(FinTS)

ArchTest(auto.arima$residuals, lags = 36)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  auto.arima$residuals
## Chi-squared = 133.21, df = 36, p-value = 4.442e-13

Não rejeitamos a hipótese nula.

NORMALIDADE

$H_0$: $\sim N (\mu, \sigma^2)$
$H_1$: não $H_0$

library(normtest)

jb.norm.test(auto.arima$residuals, nrepl = 2000)

## 
##  Jarque-Bera test for normality
## 
## data:  auto.arima$residuals
## JB = 2084.1, p-value < 2.2e-16

shapiro.test(auto.arima$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  auto.arima$residuals
## W = 0.79599, p-value = 7.544e-16

hist(auto.arima$residuals,
     main = 'Histograma dos resíduos',
     xlab = 'resíduos',
     ylab = 'Frequências')

round(summary(auto.arima$residuals), digits = 4)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -10.2346  -0.6319   0.0336  -0.1099   0.6295   6.0551

Calculando a moda

moda <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

result <- moda(as.numeric(auto.arima$residuals))
print(result)

## [1] 0.04489577

mean(auto.arima$residuals, trim = 0.05)

## [1] -0.02832715

Estamos violando esse diagóstico, porém se efetuarmos o arrendondamento, \[média \approx mediana \approx moda \approx zero\]

fit

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jul 2020       130.2386 128.2783 132.1990 127.2405 133.2368
## Aug 2020       128.8342 125.7061 131.9623 124.0502 133.6182
## Sep 2020       127.9480 124.0665 131.8295 122.0118 133.8842
## Oct 2020       127.7077 123.3654 132.0500 121.0667 134.3487
## Nov 2020       127.4674 122.6874 132.2474 120.1571 134.7778
## Dec 2020       127.2272 122.0262 132.4281 119.2730 135.1813
## Jan 2021       126.9869 121.3776 132.5962 118.4082 135.5656
## Feb 2021       126.7466 120.7385 132.7548 117.5579 135.9353
## Mar 2021       126.5064 120.1067 132.9061 116.7189 136.2939
## Apr 2021       126.2661 119.4805 133.0517 115.8885 136.6437
## May 2021       126.0258 118.8588 133.1929 115.0648 136.9869
## Jun 2021       125.7856 118.2404 133.3307 114.2463 137.3248

Índice de Atividade Econômica do Banco Central do Brasil

Alvaro Villarinho e Rogério Guerra da Silveira

26 de setembro de 2020

Procurando a identificação

Verificando o resultado e selecionando o índice

Define as datas de início e término da série

Obtém e visualiza os dados

Usando o pacote `TSstudio`

Testes

Teste ADF - no nível:

Teste KPSS

Teste Philipps - Perron

Efetuando a diferença

Efetuando o correlograma

Testes com a série diferenciada

Teste ADF - na diferença

Teste KPSS na diferença

Teste Philipps - Perron na diferença

AR(p)

Efetuando as previsões com o modelo AR(p)

MA(q)

Efetuando as previsões com o modelo MA(q)

Efetuando as previsões com o modelo ARIMA(p,d,q)

Efetuando as previsões com SARIMA

Usando o Auto.Arima

Melhor modelo encontrado (pela função `auto.arima`): SARIMA(3,1,0)(0,1,1)[12]

Checando a significância estatística dos parâmetros

Diagnóstico

Teste de ausência de autocorrelação serial

O modelo SARIMA(3,1,0)(0,1,1)[12]

Checando a Heterocedasticidade

NORMALIDADE

Calculando a moda

Índice de Atividade Econômica do Banco Central do Brasil

Alvaro Villarinho e Rogério Guerra da Silveira

26 de setembro de 2020

Procurando a identificação

Verificando o resultado e selecionando o índice

Define as datas de início e término da série

Obtém e visualiza os dados

Usando o pacote TSstudio

Testes

Teste ADF - no nível:

Teste KPSS

Teste Philipps - Perron

Efetuando a diferença

Efetuando o correlograma

Testes com a série diferenciada

Teste ADF - na diferença

Teste KPSS na diferença

Teste Philipps - Perron na diferença

AR(p)

Efetuando as previsões com o modelo AR(p)

MA(q)

Efetuando as previsões com o modelo MA(q)

Efetuando as previsões com o modelo ARIMA(p,d,q)

Efetuando as previsões com SARIMA

Usando o Auto.Arima

Melhor modelo encontrado (pela função auto.arima): SARIMA(3,1,0)(0,1,1)[12]

Checando a significância estatística dos parâmetros

Diagnóstico

Teste de ausência de autocorrelação serial

O modelo SARIMA(3,1,0)(0,1,1)[12]

Checando a Heterocedasticidade

NORMALIDADE

Calculando a moda

Usando o pacote `TSstudio`

Melhor modelo encontrado (pela função `auto.arima`): SARIMA(3,1,0)(0,1,1)[12]