if (!require(pacman)) install.packages("pacman")
pacman::p_load(gtrendsR, ggplot2, ggthemes, patchwork, RColorBrewer, magrittr, dplyr, gridExtra, ggfortify, fpp2, snpar, dygraphs, tseries, forecast, xts, seasonalview, seastests, astsa, dygraphs, fpp, lmtest, FitAR, deflateBR, sidrar, tidyverse, writexl, rbcb)

primeiro passo- ler o banco e olhar a serie

# Periodo JAN/95 - AGO/2019

# Preço constante - AGO/2019

 rm(list=ls())

 dados = read.table("icms.txt", header=F, stringsAsFactors = FALSE)

 igpdi=dados[,1]

 icmspe=dados[,2]

 str(igpdi)
##  int [1:296] 1087850 1100390 1120350 1146140 1150710 1180900 1207330 1222890 1209670 1212410 ...
 str(icmspe)
##  int [1:296] 98137000 104690000 91726000 98261000 87281000 92092000 93807000 97990000 97117000 94190000 ...
 tam = length(igpdi)

 tam_icmspe = length(icmspe)

 

 # Aqui dividimos todas as observações da série de preços

 # por seu ultimo valor. Dessa forma, a série deflacionada

 # fica a preços constantes da ultima observaçãoo.

 igpdi_d = igpdi/igpdi[tam]

 igpdi_d
##   [1] 0.1501736 0.1519047 0.1546601 0.1582203 0.1588512 0.1630188 0.1666674
##   [8] 0.1688154 0.1669904 0.1673686 0.1695898 0.1700550 0.1731058 0.1744256
##  [15] 0.1748038 0.1760214 0.1789838 0.1811719 0.1831515 0.1831584 0.1833930
##  [22] 0.1837961 0.1843152 0.1859331 0.1888666 0.1896617 0.1918705 0.1929983
##  [29] 0.1935809 0.1949309 0.1951007 0.1950152 0.1961651 0.1968360 0.1984704
##  [36] 0.1998426 0.2016000 0.2016400 0.2021107 0.2018388 0.2022985 0.2028603
##  [43] 0.2020969 0.2017463 0.2017007 0.2016345 0.2012673 0.2032468 0.2055798
##  [50] 0.2147005 0.2189413 0.2190062 0.2182511 0.2204750 0.2239842 0.2272407
##  [57] 0.2305759 0.2349298 0.2408852 0.2438545 0.2463490 0.2468267 0.2472795
##  [64] 0.2475956 0.2492604 0.2515699 0.2572533 0.2619372 0.2637359 0.2647216
##  [71] 0.2657473 0.2677683 0.2690797 0.2699908 0.2721595 0.2752283 0.2764390
##  [78] 0.2804644 0.2849964 0.2875710 0.2886699 0.2928444 0.2950821 0.2956081
##  [85] 0.2961575 0.2966986 0.2970341 0.2991089 0.3024179 0.3076664 0.3139737
##  [92] 0.3213965 0.3298932 0.3437931 0.3638623 0.3736801 0.3818055 0.3878878
##  [99] 0.3943153 0.3959373 0.3932937 0.3905507 0.3897722 0.3921962 0.3963045
## [106] 0.3980384 0.3999448 0.4023523 0.4055702 0.4099642 0.4137881 0.4185341
## [113] 0.4246523 0.4301189 0.4350016 0.4407043 0.4428357 0.4451881 0.4488546
## [120] 0.4511806 0.4526743 0.4545048 0.4589899 0.4613091 0.4601371 0.4580691
## [127] 0.4562207 0.4526356 0.4520379 0.4549024 0.4564071 0.4567053 0.4600004
## [134] 0.4597354 0.4576661 0.4577710 0.4594883 0.4625474 0.4633342 0.4652310
## [141] 0.4663436 0.4701040 0.4727835 0.4740287 0.4760524 0.4771596 0.4782018
## [148] 0.4788520 0.4796016 0.4808537 0.4826483 0.4893670 0.4950793 0.4987721
## [155] 0.5040068 0.5114406 0.5164848 0.5184437 0.5220826 0.5279081 0.5378074
## [162] 0.5479800 0.5541259 0.5520055 0.5540168 0.5600632 0.5604429 0.5579622
## [169] 0.5580436 0.5573437 0.5526722 0.5529166 0.5538857 0.5521380 0.5485861
## [176] 0.5490899 0.5504428 0.5502178 0.5506071 0.5499859 0.5555325 0.5616080
## [183] 0.5651599 0.5692212 0.5781528 0.5801269 0.5813872 0.5877843 0.5942241
## [190] 0.6003382 0.6098248 0.6121343 0.6181213 0.6240338 0.6278412 0.6309527
## [197] 0.6309955 0.6301672 0.6298470 0.6337081 0.6384762 0.6410163 0.6437524
## [204] 0.6427239 0.6446469 0.6451011 0.6486930 0.6552820 0.6612677 0.6658122
## [211] 0.6759033 0.6846389 0.6906646 0.6884904 0.6902160 0.6947632 0.6968988
## [218] 0.6982820 0.7004121 0.7000145 0.7022619 0.7076222 0.7086093 0.7118879
## [225] 0.7215538 0.7260762 0.7280862 0.7331028 0.7360584 0.7422787 0.7532962
## [232] 0.7566935 0.7532520 0.7484784 0.7443591 0.7448284 0.7449651 0.7493881
## [239] 0.7579290 0.7608404 0.7659288 0.7700260 0.7793648 0.7865018 0.7896714
## [246] 0.7950607 0.7996922 0.8029017 0.8143306 0.8286487 0.8385494 0.8422587
## [253] 0.8551633 0.8619137 0.8656327 0.8687871 0.8786201 0.8929769 0.8895092
## [260] 0.8933703 0.8936561 0.8948336 0.8953140 0.9027547 0.9066573 0.9072274
## [267] 0.9038011 0.8925697 0.8879962 0.8794636 0.8768669 0.8789597 0.8843987
## [274] 0.8852615 0.8923612 0.8989764 0.9041586 0.9055488 0.9106427 0.9190704
## [281] 0.9341519 0.9479580 0.9521684 0.9586124 0.9757577 0.9783254 0.9671533
## [288] 0.9627979 0.9634564 0.9755175 0.9859855 0.9948923 0.9988708 1.0051726
## [295] 1.0050925 1.0000000
 # Agora dividimos os valores da arrecadados do ICMS

 # pelo IGP-DI a fim de obter uma série livre de

 # dinâmina inflacionaria

 icmspe_d = icmspe/igpdi_d

 

 

 # Aqui formatamos os dados como uma série temporal

 icmspe_d.ts = ts(icmspe_d,start=c(1995,1),frequency=12)

summary(icmspe_d.ts)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 5.495e+08 7.036e+08 9.017e+08 9.840e+08 1.284e+09 1.581e+09
plot(icmspe_d.ts) #simples

autoplot(icmspe_d.ts)+xlab("Ano")+ylab("ICMs")# autoplot

segundo passo - Analisando Sazonalidade

isSeasonal(icmspe_d.ts) #fala se tem sazonalidade ou não
## [1] TRUE
ggseasonplot(icmspe_d.ts)+xlab("Meses")

resultado <- combined_test(icmspe_d.ts)
# Verificar tendencia da série
p_valor <- resultado$Pval
tendencia <- ifelse(p_valor < 0.05, "tem tendencia", "não tem tendencia")
mensagem <- paste("A série é", tendencia, "com um valor-p do teste cox stuart de", p_valor)
print(mensagem)
## [1] "A série é tem tendencia com um valor-p do teste cox stuart de 2.22044604925031e-16"
## [2] "A série é tem tendencia com um valor-p do teste cox stuart de 0"                   
## [3] "A série é tem tendencia com um valor-p do teste cox stuart de 0"

terceiro passo - Analisando Tendência

randtests::cox.stuart.test(icmspe_d.ts, alternative = "right.sided") 
## 
##  Cox Stuart test
## 
## data:  icmspe_d.ts
## statistic = 148, n = 148, p-value < 2.2e-16
## alternative hypothesis: increasing trend
resultado <- randtests::cox.stuart.test(icmspe_d.ts,alternative = "right.sided")
# Verificar tendencia da série
p_valor <- resultado$p.value
tendencia <- ifelse(p_valor < 0.05, "tem tendencia", "não tem tendencia")
mensagem <- paste("A série é", tendencia, "com um valor-p do teste cox stuart de", p_valor)
print(mensagem)
## [1] "A série é tem tendencia com um valor-p do teste cox stuart de 0"

quarto passo - testa a estacionalidade

adf.test(icmspe_d.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  icmspe_d.ts
## Dickey-Fuller = -2.4431, Lag order = 6, p-value = 0.3895
## alternative hypothesis: stationary
# Realizar o teste ADF
resultado <- adf.test(icmspe_d.ts)

# Verificar estacionariedade da série
p_valor <- resultado$p.value
estacionaria <- ifelse(p_valor < 0.05, "estacionária", "não estacionária")

# Mensagem sobre estacionariedade da série
mensagem <- paste("A série é", estacionaria, "com um valor-p do teste ADF de", p_valor)
print(mensagem)
## [1] "A série é não estacionária com um valor-p do teste ADF de 0.389487029595882"

quarto passo - 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(icmspe_d.ts, main = "Gráfico de correlação amostral", xlab = "Defasagem", ylab = "Autocorrelações")

quinto passo - 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(icmspe_d.ts, main = "Gráfico de correlação Parcial", xlab = "Defasagem", ylab = "Autocorrelações")

forecast::ggtsdisplay(icmspe_d.ts)

sexto passo - 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.

forecast::auto.arima(icmspe_d.ts)
## Series: icmspe_d.ts 
## ARIMA(2,0,0)(1,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2     sar1     sma1      drift
##       0.4533  0.4443  -0.1612  -0.6590  2868388.3
## s.e.  0.0536  0.0531   0.0841   0.0722   831809.9
## 
## sigma^2 = 3.223e+15:  log likelihood = -5476.29
## AIC=10964.57   AICc=10964.88   BIC=10986.47
fit <- astsa::sarima(icmspe_d.ts,p=2,d=0,q=0,P=1,D=1,Q=1,S=12)
## initial  value 18.304201 
## iter   2 value 18.013900
## iter   3 value 17.887100
## iter   4 value 17.884408
## iter   5 value 17.868205
## iter   6 value 17.861322
## iter   7 value 17.859677
## iter   8 value 17.859468
## iter   9 value 17.859458
## iter  10 value 17.859437
## iter  11 value 17.859421
## iter  12 value 17.859417
## iter  12 value 17.859417
## final  value 17.859417 
## converged
## initial  value 17.864546 
## iter   2 value 17.864533
## iter   3 value 17.863917
## iter   4 value 17.863883
## iter   5 value 17.863795
## iter   6 value 17.863775
## iter   7 value 17.863766
## iter   8 value 17.863764
## iter   9 value 17.863763
## iter   9 value 17.863763
## iter   9 value 17.863763
## final  value 17.863763 
## converged

setimo passo - analise de residuos

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

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 19.402, df = 10, p-value = 0.03545
## 
## Model df: 0.   Total lags used: 10
resultados = forecast::checkresiduals(residuos_padronizados)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 19.402, df = 10, p-value = 0.03545
## 
## Model df: 0.   Total lags used: 10
p_valor <- resultado$p.value
resi <- ifelse(p_valor < 0.05, " não correlacionados", " são correlacionados")

# Mensagem sobre a série
mensagem <- paste("A série é", resi, "com um valor-p do teste de residuos é de", p_valor)
print(mensagem)
## [1] "A série é  são correlacionados com um valor-p do teste de residuos é de 0.389487029595882"

oitavo passo - teste de normalidade

tseries::jarque.bera.test(residuos_padronizados)
## 
##  Jarque Bera Test
## 
## data:  residuos_padronizados
## X-squared = 225.81, df = 2, p-value < 2.2e-16
resultado=jarque.bera.test(residuos_padronizados)

p_valor <- resultado$p.value
normalidade <- ifelse(p_valor < 0.05, " tem distribuição normal", " não tem distribuição normal")

# Mensagem sobre a série
mensagem <- paste("A série", normalidade, "com um valor-p do teste de residuos é de", p_valor)
print(mensagem)
## [1] "A série  tem distribuição normal com um valor-p do teste de residuos é de 0"

nono passo -normal

nortest::lillie.test(fit$fit$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fit$fit$residuals
## D = 0.075546, p-value = 0.0003297
resultados = nortest::lillie.test(fit$fit$residuals)
residuos <- ifelse(p_valor < 0.05, " tem distribuição normal", " não tem distribuição normal")

# Mensagem sobre a série
mensagem <- paste("os residos", residuos, "com um valor-p do teste de residuos é de", p_valor)
print(mensagem)
## [1] "os residos  tem distribuição normal com um valor-p do teste de residuos é de 0"
shapiro.test(residuos_padronizados)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_padronizados
## W = 0.95463, p-value = 5.984e-08
resultados = shapiro.test(residuos_padronizados)
residuos <- ifelse(p_valor < 0.05, " tem distribuição normal", " não tem distribuição normal")

# Mensagem sobre a série
mensagem <- paste("os residos ", residuos, "com um valor-p do teste de residuos é de", p_valor)
print(mensagem)
## [1] "os residos   tem distribuição normal com um valor-p do teste de residuos é de 0"

decimo passo - alisamento exponencial

fit2 <- HoltWinters(icmspe_d.ts, seasonal = "multiplicative")
fit2
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = icmspe_d.ts, seasonal = "multiplicative")
## 
## Smoothing parameters:
##  alpha: 0.3196109
##  beta : 0.02094918
##  gamma: 0.3025145
## 
## Coefficients:
##             [,1]
## a   1.469974e+09
## b   3.001412e+06
## s1  1.008624e+00
## s2  9.851449e-01
## s3  1.038383e+00
## s4  1.049252e+00
## s5  1.083190e+00
## s6  9.144453e-01
## s7  9.110835e-01
## s8  9.867417e-01
## s9  9.380593e-01
## s10 9.510049e-01
## s11 9.692830e-01
## s12 9.601856e-01
prev <- predict(fit2, n.ahead = 6, prediction.interval = T)
prev
##                 fit        upr        lwr
## Sep 2019 1485678633 1543722639 1427634628
## Oct 2019 1454050905 1521937987 1386163822
## Nov 2019 1535745476 1614782884 1456708067
## Dec 2019 1554970269 1643031320 1466909218
## Jan 2020 1608516496 1706585431 1510447562
## Feb 2020 1360678357 1455162225 1266194489
plot(fit2,prev,main = "Arrecadação de ICMS na Paraíba", xlab = "Tempo", ylab = "Arrecadaçaõ em R$")

decimo primeiro passo -Decomposição da Série

medias moveis

#ma(auscafe,5)
#mean(auscafe[1:5])
autoplot(icmspe_d.ts,series = "Data")+autolayer(ma(icmspe_d.ts,5),series = "5-MA")+
  xlab("Ano")+ylab("icms")+ggtitle("icms")+
  scale_colour_manual(values=c("Data"="grey50","5-MA"="red"),breaks = c("Data","5-MA"))
## Warning: Removed 4 rows containing missing values (`geom_line()`).

Decomposição Clássica

#aditiva
decomp_auscafe.A=decompose(icmspe_d.ts,type = "additive")
autoplot(decomp_auscafe.A)

#multiplicativa
decomp_auscafe.M=decompose(icmspe_d.ts,type = "multiplicative")

autoplot(decomp_auscafe.M)

stl

library(seasonal)
stl_1<-seas(icmspe_d.ts)
#summary(stl_1)

library(dplyr)
icmspe_d.ts%>%
stl(s.window = "periodic")%>%
  plot