Análise da série temporal do ICMS: um modelo em R

Em diversos sites nacionais e internacionais, podem ser encontradas rotinas para a construção de modelos, bem como explicações sobre testes de verificação relativos à estacionariedade, normalidade, diagnóstico e previsão de modelos ARIMA (p, d, q) para séries temporais. O termo “p” refere-se à ordem autoregressiva do modelo, “d” é relativo à diferenciação, e “q” refere-se às médias móveis.

Neste exercício, usaremos como exemplo a arrecadação do Imposto sobre Circulação de Mercadorias e Serviços (ICMS) no período de janeiro de 2000 a abril de 2024. Apresentaremos um script detalhado para estimação e previsão de uma série temporal, seguindo as melhores práticas estatísticas. Este script é especialmente útil para aqueles interessados em aplicar algoritmos de R para análise de séries temporais univariadas. No Excel, o conjunto de dados, denominado icms_data, deve ser organizado em apenas duas colunas: date e valor. Realizada esta operação, o script realiza as demais tarefas.

# Carregar bibliotecas necessárias
library(tidyverse)
library(ggfortify)
library(xts)
library(tseries)
library(timeSeries)
library(changepoint)
library(readxl)
library(TSstudio)
library(forecast)
library(seasonal)
library(ggplot2)
library(urca)

# Carregar o conjunto de dados "icmsn" do arquivo Excel
icmsn_data <- read_excel("/Users/robertocalazans/icmsn.xlsx")

# Selecionar a coluna relevante para a série temporal
icmsn <- ts(icmsn_data$icmsn, start = c(2000, 1), frequency = 12)
dicmsn <- diff(icmsn)

# Decomposição STL da série temporal
stl_decomp <- stl(icmsn, s.window = "periodic", robust = FALSE)

# Detecção de mudanças na média e variância usando changepoint
cpt_meanvar <- changepoint::cpt.meanvar(icmsn)

# Criar um data frame para visualização 
icmsn_df <- data.frame(date = as.Date(time(icmsn)), value = as.numeric(icmsn))

# Visualização dos pontos de mudança da variância
p <- autoplot(cpt_meanvar) +
  ggtitle("Detecção de mudanças na variância") +
  theme_minimal() +
  scale_x_continuous(breaks = seq(2000, end(icmsn)[1], by = 1))
print(p)

# Plotar a decomposição STL separadamente com ggplot2
autoplot(stl_decomp) +
  ggtitle("Decomposição STL do ICMS nominal do RS") +
  theme_minimal()

# Teste de normalidade
qqnorm(icmsn, main = "QQ Plot do ICMS RS", xlab = "Quantis teóricos N(0,1)", pch = 20, ylab = "ICMS RS")
qqline(icmsn, lty = 2, col = "red")

# Gráficos da série em nível e diferenciada
Grafico.icmsn_df <- data.frame(date = as.Date(time(icmsn[-1])), 
                               nivel = as.numeric(icmsn[-1]), 
                               diff = as.numeric(dicmsn))

ggplot(Grafico.icmsn_df, aes(x = date)) +
  geom_line(aes(y = nivel, color = "Nível"), size = 1) +
  geom_line(aes(y = diff, color = "Diferenciado"), linewidth = 1, linetype = "dashed") +
  labs(title = "ICMS nominal do RS, nível e diff", y = "R$ milhões", x = "Meses") +
  scale_x_date(date_labels = "%Y-%m", date_breaks = "1 year") +
  scale_color_manual(values = c("blue", "red")) +
  theme_minimal()

# Dessazonalizando os dados
icmsn_sazo <- seasonal::seas(icmsn)
plot(icmsn, main = "ICMS nominal do RS com e sem sazonalidade", ylab = "ICMS", xlab = "Tempo")
lines(seasonal::final(icmsn_sazo), col = 2)
legend("topleft", legend = c("com sazonalidade", "sem sazonalidade"), lty = c(1, 1), col = c(1, 2))

# Teste de autocorrelação Ljung-Box Test
Box.test(icmsn, lag = 24, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  icmsn
## X-squared = 5233.7, df = 24, p-value < 2.2e-16

# p-valor < 0,05 sugere a rejeição da hipótese nula, indicando que os resíduos não são independentes e, 
# portanto, há autocorrelação.

# Dickey Fuller Test
tseries::adf.test(diff(icmsn, 1))

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(icmsn, 1)
## Dickey-Fuller = -10.547, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

summary(urca::ur.df(icmsn, type = "trend", selectlags = "BIC", lags = 1))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -759.91 -102.58    3.79   79.76 1048.92 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 28.76185   24.15185   1.191  0.23469    
## z.lag.1     -0.29161    0.04930  -5.915 9.48e-09 ***
## tt           3.70408    0.62923   5.887 1.10e-08 ***
## z.diff.lag  -0.17117    0.05858  -2.922  0.00376 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 201.4 on 286 degrees of freedom
## Multiple R-squared:  0.1975, Adjusted R-squared:  0.1891 
## F-statistic: 23.47 on 3 and 286 DF,  p-value: 1.315e-13
## 
## 
## Value of test-statistic is: -5.9145 12.4553 17.6412 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

# [tau3]  >  5% e é significativo, rejeita-se Ho, a série é estacionária.
# F > phi2 (5%), rejeita-se drift = 0, há drift na regressão.
# F > phi3  (5%), rejeita-se trend = 0, há tendência determinística.


# Phillips Perron Test
tseries::pp.test(icmsn)

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

# KPSS Test
kpps_icmsn <- urca::ur.kpss(icmsn)
summary(kpps_icmsn)

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 4.8188 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

tseries::kpss.test(icmsn)

## 
##  KPSS Test for Level Stationarity
## 
## data:  icmsn
## KPSS Level = 4.8188, Truncation lag parameter = 5, p-value = 0.01

# p-value < 0,05 e é significativo, logo rejeita-se Ho e a série apresenta raiz unitária.

# Remoção de tendência e sazonalidade
icmsn_des <- stl_decomp$time.series[, "remainder"]


# Aplicar os testes na série ajustada

# Remoção de tendência e sazonalidade
# Após a remoção da tendência e sazonalidade, todos os testes (ADF, PP e KPSS) 
# indicam que a série dessazonalizada é estacionária.

tseries::adf.test(icmsn_des)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  icmsn_des
## Dickey-Fuller = -11.216, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

tseries::pp.test(icmsn_des)

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

tseries::kpss.test(icmsn_des)

## 
##  KPSS Test for Level Stationarity
## 
## data:  icmsn_des
## KPSS Level = 0.0087624, Truncation lag parameter = 5, p-value = 0.1

# Ajuste do modelo ARIMA com drift
fit <- forecast::auto.arima(icmsn, seasonal = TRUE, stepwise = FALSE, approximation = FALSE)
summary(fit)

## Series: icmsn 
## ARIMA(2,1,1)(2,0,0)[12] with drift 
## 
## Coefficients:
##          ar1     ar2      ma1    sar1    sar2    drift
##       0.4441  0.2444  -0.9799  0.3635  0.1763  12.3613
## s.e.  0.0597  0.0609   0.0164  0.0603  0.0618   1.6212
## 
## sigma^2 = 32255:  log likelihood = -1922.35
## AIC=3858.7   AICc=3859.09   BIC=3884.41
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -2.475312 177.4305 112.4798 -1.656375 6.379175 0.5146687
##                     ACF1
## Training set -0.03021488

# RESULTADO:--->>>>>> ARIMA(2,1,1)(2,0,0)[12] with drift

# Verificação da estacionaridade dos resíduos do modelo ARIMA
checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,0,0)[12] with drift
## Q* = 39.234, df = 19, p-value = 0.004119
## 
## Model df: 5.   Total lags used: 24