Séries Temporais - PIB Paraíba

ICMS - Paraíba

Vamos estar utilizando uma série temporal, encontrada no site do IPEA DATA, sobre a arrecadação de ICMS na Paraíba. O link onde pode ser encontrado a série temporal é: Série Temporal

Comentário: Imposto sobre Operações Relativas à Circulação de Mercadorias e sobre Prestações de Serviços de Transporte Interestadual e Intermunicipal e de Comunicação (ICMS). É um dos Impostos sobre a Produção e a Circulação, sendo um tributo indireto e regressivo. De competência dos Estados, tem como fato gerador as operações citadas, ainda que as operações e as prestações se iniciem no exterior. Incide ainda sobre a entrada de mercadoria importada do exterior. Números referente a um dado estado. Nota: a partir de março de 2014, dados do “Demonstrativo da Receita Corrente Líquida” do “Relatório Resumido de Execução Orçamentária” (RREO), obtidos no “Sistema de Informações Contábeis e Fiscais do Setor Público Brasileiro” (Siconfi). Mais informações: Função Social dos Tributos.

Carregando a série Temporal

setwd("C:\\Users\\Pessoal\\Desktop\\ESTATÍSTICA\\UFPB\\8º PERÍODO\\SÉRIES TEMPORAIS")

dados <- read.csv2("dados2.csv", header = F, dec = ",")

Fazendo algumas manipulações no banco de dados

names(dados) <- c("Data", "ICMS")

Verificando o tipo de variável

glimpse(dados)

Rows: 362
Columns: 2
$ Data <chr> "1993.01", "1993.02", "1993.03", "1993.04", "1993.05", "1993.06",…
$ ICMS <dbl> 187.90, 214.49, 263.32, 334.06, 426.48, 632.42, 872.97, 1181.24, …

Como o real entrou em vigor a partir de 1994, vamos omitir os dois primeiros anos da série e trabalhar apenas de 1995 em diante.

dados <- dados[-c(1:24),]

Visualizando as primeiras 10 linhas

dados %>%
  slice_head(n=10) %>%
  gt::gt()

Data	ICMS
1995.01	30896
1995.02	28483
1995.03	29352
1995.04	26762
1995.05	26876
1995.06	27771
1995.07	25503
1995.08	28720
1995.09	32856
1995.10	29159

Deflacionando a série via ipc (Indíce de preço ao consumidor)

times <- seq(as.Date('1995-01-01'), as.Date('2023-02-01'),by='1 month')
icms <- deflateBR::deflate(dados$ICMS, times,"02/2023" ,index = "ipc")

Downloading necessary data from IPEA's API
...

Transformando em série Temporal

icms_PB <- ts(icms, start = c(1995,1), frequency = 12)

icms_PB

           Jan       Feb       Mar       Apr       May       Jun       Jul
1995 203899.03 184967.58 186932.90 165896.71 161913.77 163690.32 143994.59
1996 185537.78 183152.05 179851.07 175097.05 177124.56 177596.72 166912.29
1997 212482.48 194474.63 165775.90 182484.63 172570.52 187964.77 180973.32
1998 219391.19 186447.38 182247.56 177461.56 191591.68 189045.38 187507.88
1999 229657.40 209932.18 207016.17 217105.33 200809.60 190914.93 191899.87
2000 229003.28 221787.81 215433.95 214762.53 220815.74 235071.27 235134.12
2001 322768.15 260657.07 270032.00 263447.62 266731.99 278307.12 269748.77
2002 273319.67 295517.12 242564.82 261138.72 203546.35 239113.14 226018.12
2003 319154.07 282557.94 240745.49 224400.46 213777.74 260119.45 265406.67
2004 257786.47 266556.24 235536.36 245477.64 258746.17 232353.29 232830.33
2005 301129.31 267135.70 270002.99 282620.00 281029.65 288502.47 270858.68
2006 342472.04 317493.92 287776.07 307182.33 305108.38 305229.93 308402.04
2007 354235.44 333083.81 303469.41 310119.74 326391.73 334968.69 340400.54
2008 396287.64 367856.88 348449.00 353390.22 356867.15 362192.77 351707.46
2009 393448.07 368972.29 355871.49 336826.67 348890.50 361317.01 394604.59
2010 459608.71 413507.63 399448.69 430980.60 408065.42 436925.05 446163.89
2011 508425.37 476718.96 441466.61 445730.28 447724.26 459146.20 468139.49
2012 518326.94 481169.05 479639.38 472424.15 482979.87 493300.50 498118.49
2013 620927.43 545395.01 471983.41 534918.43 531177.08 556251.30 599703.75
2014 637707.79 559601.48 583446.19 568134.68 614247.37 583915.26 554951.32
2015 685436.84 570419.10 523430.60 581592.26 577583.98 546551.40 549466.69
2016 601946.20 552610.19 486134.68 543434.91 524943.94 541629.94 541336.35
2017 633362.68 545049.55 528011.25 558457.05 537703.02 536646.71 537118.98
2018 626620.46 574339.26 541742.34 550492.88 566026.35 514114.45 577379.00
2019 662558.26 596537.23 567038.44 581187.19 609405.66 582043.03 511409.55
2020 695670.16 603384.25 563265.07 516278.45 429554.74 533006.83 553502.79
2021 750374.75 690605.80 641276.72 604020.25 598533.21 651385.32 706911.44
2022 839304.72 641169.95 611303.06 658104.53 663508.97 671881.48 661400.64
2023 708411.01 634287.65                                                  
           Aug       Sep       Oct       Nov       Dec
1995 158008.62 179431.42 158174.69 160708.25 191885.76
1996 173176.30 186253.72 188010.80 182002.22 213369.37
1997 172529.50 178686.69 192282.83 178056.27 200813.09
1998 185814.58 198815.37 208879.26 222130.43 217804.04
1999 221512.71 196899.80 198050.33 204373.96 239634.07
2000 224338.17 248547.21 257320.21 259965.52 264868.18
2001 268252.10 254971.22 260656.71 269563.52 295611.79
2002 257812.71 306419.47 249605.66 274246.37 257230.09
2003 229347.18 234123.41 267342.97 231119.04 266857.41
2004 237336.92 313664.73 254111.29 305243.02 291320.24
2005 281878.63 299454.58 288547.62 312647.16 313031.06
2006 329171.28 338984.82 318558.10 333561.18 356834.46
2007 340218.10 336093.60 357328.26 349206.66 388251.66
2008 379262.10 369229.40 384644.31 392046.39 377439.14
2009 376728.28 394285.73 398528.38 446306.14 417180.25
2010 432473.88 446414.98 449441.04 478620.56 468045.35
2011 468041.68 466213.98 443575.35 452581.42 463047.98
2012 489703.88 491242.90 516556.06 538149.15 565930.91
2013 542380.47 538870.23 553659.59 573670.25 574673.76
2014 564608.65 588289.51 617436.98 647719.43 621369.65
2015 526635.96  43242.99  44527.95 550050.73 555644.86
2016 539636.24 525643.27 561526.28 561817.28 571130.69
2017 536484.96 567071.11 582729.66 597634.79 609277.82
2018 567301.47 640898.18 623451.55 621135.64 666382.98
2019 518965.83 585538.51 593474.88 643324.66 667361.62
2020 628134.29 655216.78 662137.98 698101.90 755371.78
2021 755547.84 713269.30 737088.34 708873.42 770463.45
2022 580198.90 623991.16 613980.04 644417.83 667701.18
2023                                                  

Estatística Descritiva da Série

which.max(icms_PB)

[1] 325

Tivemos que o maior valor de arrecadação de ICMS na Paraíba foi em 2022-01-01, e que o valor que foi arrecadado foi de R$ 8.3930472^{5}.

which.min(icms_PB)

[1] 249

Tivemos que o menor valor de arrecadação de ICMS na Paraíba foi em 2015-09-01, e que o valor que foi arrecadado foi de R$ 4.3242994^{4}.

skimr::skim(icms_PB)

Data summary

Name	icms_PB
Number of rows	338
Number of columns	1
_______________________
Column type frequency:
ts	1
________________________
Group variables	None

Variable type: ts

skim_variable	n_missing	complete_rate	start	end	frequency	deltat	mean	sd	min	max	median	line_graph
x	0	1	1995	2023	12	0.08	401908.2	172366.8	43242.99	839304.7	372978.8	⣀⣀⣀⡠⠒⠉⠉⠉

Visualização Gráfica da Série

forecast::autoplot(icms_PB)+
  ggplot2::labs(title = "Arrecadação de ICMS na Paraíba", subtitle = "De Janeiro de 1955 à Fevereiro de 2023") +
  ggplot2::xlab("Tempo") + ggplot2::ylab("Arrecadação em R$")

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Características da série

Aqui, vamos testar algumas características que estão presentes em algumas séries temporais, tais como: estacionariedade, tendência, sazonalidade.

Testando a estacionariedade

As hipóteses são:

\[H_0: A \hspace{0.1cm} série \hspace{0.1cm} é \hspace{0.1cm} não \hspace{0.1cm} estacionária\] \[H_1: A \hspace{0.1cm} série \hspace{0.1cm} é \hspace{0.1cm} estacionária\]

tseries::adf.test(icms_PB)

Warning in tseries::adf.test(icms_PB): p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  icms_PB
Dickey-Fuller = -4.8425, Lag order = 6, p-value = 0.01
alternative hypothesis: stationary

Resposta: Rejeitamos \(H_0\), ou seja, com 95% de confiança a série, segundo o teste de Dickey-Fuller é não estacionária. O resultado do p-valor foi de 0.01.

Sendo assim, precisamos tomar diferença para tornar a série estacionária, vamos saber quantas diferenças vai ser preciso tomar.

forecast::ndiffs(icms_PB)

[1] 1

Segundo o comando ndiffs precisamos tomar uma diferença para tornar a série estacionária.

Verificando se é preciso tomar diferença na componente sazonal

forecast::nsdiffs(icms_PB)

[1] 0

Conforme o resultado do teste não vai ser preciso tomar diferença na componente sazonal

Testando a tendência

• Como pelo gráfico, é observado que possivelmente tem se uma tendência crescente, esse será o argumento a ser testado.

• As hipóteses são as seguintes:

\[H_0: A \hspace{0.1cm} série \hspace{0.1cm} não \hspace{0.1cm} tem \hspace{0.1cm} tendência \hspace{0.1cm} crescente\] \[H_1: A \hspace{0.1cm} série \hspace{0.1cm} tem \hspace{0.1cm} tendência \hspace{0.1cm} crescente\]

randtests::cox.stuart.test(icms_PB, alternative = "right.sided")

    Cox Stuart test

data:  icms_PB
statistic = 167, n = 169, p-value < 2.2e-16
alternative hypothesis: increasing trend

Resposta: Rejeitamos \(H_0\), ou seja, com 95% de confiança a série tem tendência crescente. O resultado do p-valor foi menor que \(2.2 \times 10^{-16}\).

Testando a sazonalidade

Visualização gráfica

forecast::ggseasonplot(icms_PB) + 
  ggplot2::labs(title = "Gráfico de Sazonalidade da Arrecadação do ICSM na Paraíba", subtitle = "De Janeiro de 1955 à Fevereiro de 2023") + ggplot2::ylab("Arrecadação em R$") + ggplot2::xlab("Mês")

Comentário: De acordo com a visualização gráfica, podemos observar que a série tem a presença da componente sazonal, vamos utilizar um teste para comprovar.

• As hipóteses do teste são:

\[H_0: A \hspace{0.1cm} série \hspace{0.1cm} não \hspace{0.1cm} possui \hspace{0.1cm} sazonalidade\] \[H_1: A \hspace{0.1cm} série \hspace{0.1cm} possui \hspace{0.1cm} sazonalidade\]

seastests::combined_test(icms_PB)

Test used:  WO 
 
Test statistic:  1 
P-value:  0.000196446 0.0003388343 1.676437e-14

seastests::isSeasonal(icms_PB)

[1] TRUE

Resposta: Rejeitamos \(H_0\), ao nível de confiança de 95% rejeitamos a hipótese nula, ou seja, com 95% de confiança a série possui sazonalidade.

Então a série de Arrecadação de ICMS da Paraíba de 1955 à 2023 possui tendência crescente, é não estacionária (sendo necessário tomar uma diferença) e possui sazonalidade.

Gráficos de Autocorrelação

Gráfico de Autocorrelação Amostral

• O gráfico de correlação amostral é útil na identificação da ordem do processo se ele for um MA, e se for um AR, é esperado um decaimento exponencial.

forecast::Acf(icms_PB, main = "Gráfico de correlação amostral", xlab = "Defasagem", ylab = "Autocorrelações")

É observado um decaimento exponencial, sugerindo assim que o processo seja do tipo AR

Gráfico de Correlação Parcial

• Acima, verificamos que a série possivelmente seja do tipo AR, com o gráfico de Correlação Parcial, é sugerido a ordem do processo, então vamos verificar.

forecast::Pacf(icms_PB, main = "Gráfico de correlação Parcial", xlab = "Defasagem", ylab = "Autocorrelações")

forecast::ggtsdisplay(icms_PB)

Verificação da série

O objetivo é fazer a verificação da série para ajustar um processo que mais se adeque a série, temos já como evidência de que existe a parte autoregressora da série, através do gráfico de correlação parcial, pode ser comprovado.

• Utilizaremos a função auto.arima do pacote forecast.

forecast::auto.arima(icms_PB)

Series: icms_PB 
ARIMA(1,1,2)(0,0,2)[12] with drift 

Coefficients:
         ar1      ma1      ma2    sma1    sma2      drift
      0.3173  -0.5500  -0.3806  0.1349  0.1075  1519.8133
s.e.  0.1094   0.1107   0.0882  0.0557  0.0503   322.2786

sigma^2 = 1.985e+09:  log likelihood = -4083.4
AIC=8180.8   AICc=8181.14   BIC=8207.54

Ajustando um processo

fit <- astsa::sarima(icms_PB,p=1,d=1,q=2,P=0,D=0,Q=2,S=12)

initial  value 10.832313 
iter   2 value 10.740252
iter   3 value 10.717259
iter   4 value 10.712812
iter   5 value 10.709726
iter   6 value 10.708581
iter   7 value 10.705633
iter   8 value 10.703939
iter   9 value 10.700861
iter  10 value 10.697729
iter  11 value 10.697040
iter  12 value 10.696675
iter  13 value 10.696481
iter  14 value 10.696461
iter  15 value 10.696455
iter  16 value 10.696454
iter  17 value 10.696453
iter  18 value 10.696453
iter  19 value 10.696452
iter  20 value 10.696452
iter  20 value 10.696452
iter  20 value 10.696452
final  value 10.696452 
converged
initial  value 10.698024 
iter   2 value 10.698019
iter   3 value 10.698003
iter   4 value 10.697992
iter   5 value 10.697987
iter   6 value 10.697985
iter   7 value 10.697981
iter   8 value 10.697981
iter   8 value 10.697981
iter   8 value 10.697981
final  value 10.697981 
converged

fit$ttable

          Estimate       SE t.value p.value
ar1         0.3173   0.1094  2.9003  0.0040
ma1        -0.5500   0.1107 -4.9668  0.0000
ma2        -0.3806   0.0882 -4.3136  0.0000
sma1        0.1349   0.0557  2.4220  0.0160
sma2        0.1075   0.0503  2.1381  0.0332
constant 1519.8133 322.2786  4.7158  0.0000

Resposta: De acordo com o teste é observado que todos os coeficientes são significativos para o modelo e foram iguais aos encontrados na função auto.arima, sendo assim, manteremos este ajuste.

De acordo com a visualização gráfica do ajuste, observamos que os resíduos são não correlacionados, pelo teste de Ljung-box.

Análise Diagnóstico dos Resíduos

• Precisamos padronizar os resíduos, para que os valores discrepantes não interfiram muito por conta da escala

residuos <- fit$fit$residuals
residuos_padronizados <- c((residuos - mean(residuos))/sd(residuos))

• Realizando o teste de maneira mais genérica, através da função checkresiduals, da biblioteca forecast.

forecast::checkresiduals(residuos_padronizados)

    Ljung-Box test

data:  Residuals
Q* = 3.7984, df = 10, p-value = 0.956

Model df: 0.   Total lags used: 10

Resposta: De acordo com a saída, do teste de Ljung-Box, e com a análise gráfica residual é observado que os resíduos são não correlacionados, pois a hipótese nula de que os resíduos são não correlacionados não foi rejeitada.

Verificando ouliers

ggplot2::ggplot(icms_PB, ggplot2::aes(x=as.vector(1:338), y = residuos_padronizados))+
  ggplot2::geom_point(size = 4, color = "#e7ad52") +
  ggplot2::geom_abline(intercept = -2.5, slope = 0,color = "#001132")+
  ggplot2::geom_abline(intercept = 2.5, slope = 0,color = "#001132") + 
  ggplot2::ylim(-11.6,6) + 
  ggplot2::xlab("Observações") + 
  ggplot2::ylab("Resíduos Padronizados")+
  ggplot2::labs(title = "Gráfico dos resíduos padronizados")

Traçando uma reta em 2.5

Tivemos a quantidade de 5 outliers.

Que foram:

banco <- data.frame(times,icms)
banco[which(residuos_padronizados>2.5 | residuos_padronizados< -2.5),]

         times      icms
249 2015-09-01  43242.99
250 2015-10-01  44527.95
251 2015-11-01 550050.73
305 2020-05-01 429554.74
326 2022-02-01 641169.95

Testando a normalidade

• Para o teste de normalidade dos resíduos vamos utilizar três testes

Teste de Bera-Jarque

• As hipóteses são:

\[H_0: Os \hspace{0.1cm} resíduos \hspace{0.1cm} possuem \hspace{0.1cm} distribuição \hspace{0.1cm} normal\]

\[H_1: Os \hspace{0.1cm} resíduos \hspace{0.1cm} não \hspace{0.1cm} possuem \hspace{0.1cm} distribuição \hspace{0.1cm} normal\]

tseries::jarque.bera.test(residuos_padronizados)

    Jarque Bera Test

data:  residuos_padronizados
X-squared = 41033, df = 2, p-value < 2.2e-16

Resposta: Rejeitamos \(H_0\), ou seja, com 95% de confiança os resíduos não estão distribuidos normalmente.

Teste de Lillie-Fors

• As hipóteses são:

\[H_0: Os \hspace{0.1cm} resíduos \hspace{0.1cm} possuem \hspace{0.1cm} distribuição \hspace{0.1cm} normal\]

\[H_1: Os \hspace{0.1cm} resíduos \hspace{0.1cm} não \hspace{0.1cm} possuem \hspace{0.1cm} distribuição \hspace{0.1cm} normal\]

nortest::lillie.test(fit$fit$residuals)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  fit$fit$residuals
D = 0.13894, p-value < 2.2e-16

Resposta: Rejeitamos \(H_0\), ou seja, com 95% de confiança os resíduos não estão distribuidos normalmente.

Teste de Shapiro-Wilk

• As hipóteses são:

\[H_0: Os \hspace{0.1cm} resíduos \hspace{0.1cm} possuem \hspace{0.1cm} distribuição \hspace{0.1cm} normal\]

\[H_1: Os \hspace{0.1cm} resíduos \hspace{0.1cm} não \hspace{0.1cm} possuem \hspace{0.1cm} distribuição \hspace{0.1cm} normal\]

shapiro.test(residuos_padronizados)

    Shapiro-Wilk normality test

data:  residuos_padronizados
W = 0.70652, p-value < 2.2e-16

Resposta: Rejeitamos \(H_0\), ou seja, com 95% de confiança os resíduos não estão distribuidos normalmente.

Sintese: Os resíduos não estão distribuidos normalmente.

Alisamento Exponencial

• Como a série tem a presença de tendência e de sazonalidade, vamos utilizar o alisamento exponencial de HoltWinters, que é utilizado na presença dessas componentes na série temporal.

fit2 <- HoltWinters(icms_PB, seasonal = "multiplicative")
fit2

Holt-Winters exponential smoothing with trend and multiplicative seasonal component.

Call:
HoltWinters(x = icms_PB, seasonal = "multiplicative")

Smoothing parameters:
 alpha: 0.3159753
 beta : 0
 gamma: 0.10589

Coefficients:
            [,1]
a   6.342858e+05
b   1.063932e+03
s1  9.554144e-01
s2  9.764387e-01
s3  9.684284e-01
s4  9.906426e-01
s5  9.774125e-01
s6  9.726577e-01
s7  9.600750e-01
s8  9.433302e-01
s9  1.035607e+00
s10 1.074824e+00
s11 1.122869e+00
s12 1.005939e+00

Fazendo previsões para os 6 próximos meses

prev <- predict(fit2, n.ahead = 6, prediction.interval = T)
prev

              fit      upr      lwr
Mar 2023 607022.3 645292.6 568751.9
Apr 2023 621418.9 671188.8 571649.0
May 2023 617351.4 675780.5 558922.3
Jun 2023 632566.4 699787.6 565345.2
Jul 2023 625158.3 698515.3 551801.4
Aug 2023 623152.0 702544.2 543759.8

Visualizando graficamente

plot(fit2,prev,main = "Arrecadação de ICMS na Paraíba", xlab = "Tempo", ylab = "Arrecadaçaõ em R$")

Decomposição da Série

• Vamos aplicar o método de decomposição stl para a série, pois ele é sensível a outliers.

decomp <- stl(icms_PB, s.window = "periodic")
forecast::autoplot(decomp) +
    ggplot2::labs(title = "Decomposição via método STL da série ICMS - PB")+
    ggplot2::xlab("Ano")