AT2_EconometriaSARIMA
Guilherme Ziegler Trintinalia
2022-10-12
Importação da base de dados
# importando a base de dados
rm(list=ls());
library(readxl)
dados<-read_xlsx("..\\Dados\\AT2_Grupo_9.xlsx",sheet = "Dados", col_names = T); #Dados em excel
dados## # A tibble: 163 × 2
## Data PEAm
## <chr> <dbl>
## 1 2009.01 20167.
## 2 2009.02 20206.
## 3 2009.03 20127.
## 4 2009.04 20174.
## 5 2009.05 20356.
## 6 2009.06 20641.
## 7 2009.07 20675.
## 8 2009.08 20767.
## 9 2009.09 20852.
## 10 2009.10 20614.
## # … with 153 more rowsTransformação e manipulação dos dados
# Transformar os dados
peam <-as.vector(dados$PEAm)
logpeam <-as.vector(log(peam))
logpeam## [1] 9.911800 9.913729 9.909797 9.912149 9.921152 9.935037 9.936659
## [8] 9.941101 9.945223 9.933711 9.950589 9.953904 9.961057 9.965266
## [15] 9.973821 9.976353 9.973425 9.984111 9.995189 9.991123 9.993666
## [22] 9.981427 9.969557 9.973474 9.986844 9.996863 9.994969 9.993523
## [29] 9.997002 10.001645 10.006296 10.003517 10.009002 9.998345 9.990747
## [36] 9.992533 10.001112 10.001339 10.006533 9.998074 9.997563 10.001163
## [43] 10.012212 10.012121 10.014336 10.005772 10.004153 10.008705 10.012070
## [50] 10.008190 10.011383 10.022312 10.029842 10.039129 10.043993 10.039995
## [57] 10.040292 10.028459 10.026595 10.028332 10.038267 10.034614 10.035368
## [64] 10.042946 10.042262 10.053163 10.056401 10.055275 10.056963 10.044824
## [71] 10.045503 10.048078 10.052403 10.056463 10.053651 10.064110 10.064055
## [78] 10.065322 10.072914 10.078615 10.075152 10.068841 10.067733 10.061037
## [85] 10.067666 10.064481 10.067204 10.067076 10.074637 10.081071 10.080542
## [92] 10.077825 10.080478 10.084450 10.079687 10.084115 10.090325 10.090151
## [99] 10.094848 10.090135 10.096171 10.099278 10.100504 10.101721 10.100352
## [106] 10.098007 10.089687 10.097848 10.103371 10.104736 10.108930 10.105341
## [113] 10.107674 10.113746 10.115358 10.113961 10.114591 10.105670 10.103957
## [120] 10.110930 10.117288 10.118469 10.127959 10.122556 10.115806 10.121529
## [127] 10.131690 10.139504 10.139414 10.135627 10.131855 10.125934 10.124158
## [134] 10.125032 10.129286 10.129785 10.133012 10.133414 10.133645 10.134720
## [141] 10.129503 10.128239 10.123666 10.121489 10.118285 10.117113 10.117013
## [148] 10.121120 10.117901 10.126953 10.124028 10.129602 10.136097 10.123746
## [155] 10.124351 10.120464 10.121684 10.126389 10.130352 10.131367 10.137479
## [162] 10.136822 10.134792# declarar séries temporais
graphics.off()
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo# Plotando PEA nas regiões metropolitanas
logpeam <-ts(logpeam, start=c(2009,1), frequency=12)
peam <-ts(peam, start=c(2009,1), frequency=12)
par(mfrow=c(1,2))
plot(peam, col=6, main = "População Economicamente Ativa nas regiões metropolitanas")
plot(logpeam, col=4, main= "PEAm em Ln")
# Juntar os dados em um ?nico data-frame
dados<-as.data.frame(cbind(logpeam, peam))Análise da estacionariedade pela FAC e FACP
#Elaborar FAC e FACP
library(forecast)
tsdisplay(peam, lag.max = 24)
# Gráficos dividindo a tela
par(mfrow=c(1,2))
acf(logpeam, ci.type = "ma" , main = " Função de Autocorrelação",
ylab = "Autocorrelação", xlab = "Defasagens", lag.max = 12)
pacf(logpeam, main = "Função de Autocorrelação Parcial",
ylab = "Autocorrela??o", xlab = "Defasagens", lag.max = 12)
# Diferenciação das séries (d = 1)
# Percebe-se que as séries são não estacionárias, logo necessita-se diferencia-las
par(mfrow=c(1,2))
Acf(diff(logpeam), main='Correlograma da Diferen?a desemprego no Brasil',
ylab='rho', ylim=c(-0.5, 0.6), lag.max = 12)
# p do arima = 1
Pacf(diff(logpeam), main='Correlograma da Diferen?a do desemprego do Brasil',
ylab='rho', ylim=c(-0.5, 0.6), lag.max = 12)
# q do arima = 1
Acf(diff(logpeam), main='Correlograma da Diferen?a desemprego no Brasil',
ylab='rho', ylim=c(-0.5, 0.6), lag.max = 48)
# P = 1, mem?ria longa da FAQP sendo projetada na FAQ nos picos de 12.
Pacf(diff(logpeam), main='Correlograma da Diferen?a do desemprego do Brasil',
ylab='rho', ylim=c(-0.5, 0.6), lag.max = 48)
# pela Análise das séries diferenciadas, pode-se atribuir MODELO SARIMA, tal que:
# Q = 1
# s = 12
# (p, d, q) X (P, D, Q)s = (1, 1, 1)X(1,0,1)12Pacote auto.arima
# A seguir realiza-se a função auto.arima para estabelecer os melhores parâmetros e comparar com aqueles analisados pelo cientista
# modelo retornado pela função auto.arima abaixo
arima010200s12 <-auto.arima(logpeam, d=1, D=0, max.p=4,
max.q = 4, max.P = 2, max.Q = 2,test = c("adf"),
ic=c("bic"), allowdrift = TRUE, method = 'ML')
arima010200s12## Series: logpeam
## ARIMA(0,1,0)(2,0,0)[12]
##
## Coefficients:
## sar1 sar2
## 0.3154 0.2336
## s.e. 0.0838 0.0913
##
## sigma^2 = 2.783e-05: log likelihood = 620.8
## AIC=-1235.6 AICc=-1235.45 BIC=-1226.33accuracy(arima010200s12)## ME RMSE MAE MPE MAPE
## Training set 0.0006358402 0.005226679 0.004104881 0.006350898 0.04085532
## MASE ACF1
## Training set 0.2371783 -0.01312832# Modelo sugerido
arima111101s12 <-Arima(logpeam, order=c(1,1,1), seasonal=c(1,0,1),
include.drift=TRUE, method = 'ML')
arima111101s12## Series: logpeam
## ARIMA(1,1,1)(1,0,1)[12] with drift
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produzidos## ar1 ma1 sar1 sma1 drift
## -0.0052 -0.0335 0.9437 -0.7723 9e-04
## s.e. 0.0016 NaN 0.0002 NaN 9e-04
##
## sigma^2 = 2.605e-05: log likelihood = 626.67
## AIC=-1241.35 AICc=-1240.81 BIC=-1222.82accuracy(arima111101s12)## ME RMSE MAE MPE MAPE
## Training set -2.317061e-05 0.005008705 0.003927331 -0.0001957287 0.03908663
## MASE ACF1
## Training set 0.2269195 0.0003510027#Análise de Estabilidade (RU)
# Modelo sugerido tem MAPA ligeriramente melhor e menor do que o modelo do autoarima, mas envia mensagem de erro, que deve ser inspecionada
autoplot(arima010200s12)
autoplot(arima111101s12)
#ambos são estáveis com raízes unitárias (RU) dentro do círculo unitário.#Análise de autocorrelação
# teste de Ljung Box
Box.test(arima010200s12$residuals, lag = 1, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: arima010200s12$residuals
## X-squared = 0.028614, df = 1, p-value = 0.8657Box.test(arima010200s12$residuals, lag = 3, type='Ljung-Box')##
## Box-Ljung test
##
## data: arima010200s12$residuals
## X-squared = 1.3084, df = 3, p-value = 0.7271Box.test(arima010200s12$residuals, lag = 6, type='Ljung-Box')##
## Box-Ljung test
##
## data: arima010200s12$residuals
## X-squared = 4.2757, df = 6, p-value = 0.6394Box.test(arima010200s12$residuals, lag = 12, type='Ljung-Box')##
## Box-Ljung test
##
## data: arima010200s12$residuals
## X-squared = 10.42, df = 12, p-value = 0.5792# Não há autocorreção para o modelo da função auto.arima para as defasagens 1, 3, 6 e 12, para nivel de confiança de 95%, hipótese nula não é rejeitada, matendo-se hipótese nula de ausência de autocorrelação.
Box.test(arima111101s12$residuals, lag = 1, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: arima111101s12$residuals
## X-squared = 2.0454e-05, df = 1, p-value = 0.9964Box.test(arima111101s12$residuals, lag = 3, type='Ljung-Box')##
## Box-Ljung test
##
## data: arima111101s12$residuals
## X-squared = 1.5455, df = 3, p-value = 0.6718Box.test(arima111101s12$residuals, lag = 6, type='Ljung-Box')##
## Box-Ljung test
##
## data: arima111101s12$residuals
## X-squared = 5.0087, df = 6, p-value = 0.5427Box.test(arima111101s12$residuals, lag = 12, type='Ljung-Box')##
## Box-Ljung test
##
## data: arima111101s12$residuals
## X-squared = 11.729, df = 12, p-value = 0.4677# o mesmo é verificado para o modelo sugerido pela análise do pesquisador.# Teste BG para autocorrelação
checkresiduals(arima010200s12, test="LB", lag = 12)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)(2,0,0)[12]
## Q* = 10.42, df = 10, p-value = 0.4045
##
## Model df: 2. Total lags used: 12library(lmtest)## Carregando pacotes exigidos: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericres_arima010200s12<-lm(residuals(arima010200s12) ~ 1)
bgtest(res_arima010200s12)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: res_arima010200s12
## LM test = 0.028114, df = 1, p-value = 0.8668checkresiduals(arima111101s12, test="LB", lag = 12)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,0,1)[12] with drift
## Q* = 11.729, df = 8, p-value = 0.1637
##
## Model df: 4. Total lags used: 12library(lmtest)
res_arima111101s12<-lm(residuals(arima111101s12) ~ 1)
bgtest(res_arima111101s12)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: res_arima111101s12
## LM test = 2.0165e-05, df = 1, p-value = 0.9964# Novamente, p-valor maior que 0,05 para nível de confiança de 95%, NÃO rejeitando H0 de ausência de autocorrelação. STeste de heterocedasticia
# Teste ARCH-LM com defasagens 1, 3 e 6.
library(FinTS)##
## Attaching package: 'FinTS'## The following object is masked from 'package:forecast':
##
## AcfArchTest(arima010200s12$residuals, lags=1)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: arima010200s12$residuals
## Chi-squared = 0.48325, df = 1, p-value = 0.487ArchTest(arima010200s12$residuals, lags=3)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: arima010200s12$residuals
## Chi-squared = 0.88776, df = 3, p-value = 0.8284ArchTest(arima010200s12$residuals, lags=6)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: arima010200s12$residuals
## Chi-squared = 6.194, df = 6, p-value = 0.4018# Modeloauto.arima passou
ArchTest(arima111101s12$residuals, lags=1)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: arima111101s12$residuals
## Chi-squared = 0.93836, df = 1, p-value = 0.3327ArchTest(arima111101s12$residuals, lags=3)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: arima111101s12$residuals
## Chi-squared = 1.1968, df = 3, p-value = 0.7538ArchTest(arima111101s12$residuals, lags=6)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: arima111101s12$residuals
## Chi-squared = 5.3558, df = 6, p-value = 0.499# Modelo sugerido também passa para os lags 1, 3 e 6, todos os p-valores são maiores que 0,05, ficando com H0. Normalidade dos resíduos
# Teste de normalidade, Ho = normalidade p-valor > 5%, H1 = Não há normalidade, p-valor < 5%.
jarque.bera.test(arima010200s12$residuals)##
## Jarque Bera Test
##
## data: arima010200s12$residuals
## X-squared = 2.3497, df = 2, p-value = 0.3089shapiro.test(arima010200s12$residuals)##
## Shapiro-Wilk normality test
##
## data: arima010200s12$residuals
## W = 0.99345, p-value = 0.6761jarque.bera.test(arima111101s12$residuals)##
## Jarque Bera Test
##
## data: arima111101s12$residuals
## X-squared = 0.9833, df = 2, p-value = 0.6116shapiro.test(arima111101s12$residuals)##
## Shapiro-Wilk normality test
##
## data: arima111101s12$residuals
## W = 0.99315, p-value = 0.6382Projeção do modelo autoarima
arima010200s12_f<-arima(logpeam, order = c(0,1,0), seasonal = c(2,0,0))
update(arima010200s12_f, method ='ML')##
## Call:
## arima(x = logpeam, order = c(0, 1, 0), seasonal = c(2, 0, 0), method = "ML")
##
## Coefficients:
## sar1 sar2
## 0.3154 0.2336
## s.e. 0.0838 0.0913
##
## sigma^2 estimated as 2.688e-05: log likelihood = 620.8, aic = -1235.6arima010200s12_f_cv<-forecast(arima010200s12_f, h = 5)
# Prevendo valores 5 passos (meses) a frente.
arima010200s12_f_cv## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Aug 2022 10.13680 10.13016 10.14345 10.12664 10.14696
## Sep 2022 10.13763 10.12824 10.14703 10.12326 10.15200
## Oct 2022 10.13344 10.12193 10.14495 10.11584 10.15104
## Nov 2022 10.13256 10.11927 10.14585 10.11224 10.15289
## Dec 2022 10.13083 10.11597 10.14569 10.10811 10.15355autoplot(arima010200s12_f_cv, type = "l", fcol = 3, flty=1, flwd=1)
Projeção do modelo sugerido
arima111101s12_f<-arima(logpeam, order = c(1,1,1), seasonal = c(1,0,1))
update(arima111101s12_f, method ='ML')##
## Call:
## arima(x = logpeam, order = c(1, 1, 1), seasonal = c(1, 0, 1), method = "ML")
##
## Coefficients:
## ar1 ma1 sar1 sma1
## 0.5364 -0.6073 0.9658 -0.7916
## s.e. 0.3196 0.2960 0.0298 0.0918
##
## sigma^2 estimated as 2.419e-05: log likelihood = 626.37, aic = -1242.74arima111101s12_f_cv<-forecast(arima111101s12_f, h = 5) # Prevendo valores 5 passos (meses) a frente.
arima111101s12_f_cv## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Aug 2022 10.13695 10.13065 10.14326 10.12731 10.14659
## Sep 2022 10.13784 10.12924 10.14645 10.12468 10.15100
## Oct 2022 10.13255 10.12228 10.14283 10.11684 10.14827
## Nov 2022 10.13046 10.11881 10.14211 10.11264 10.14828
## Dec 2022 10.13041 10.11756 10.14326 10.11076 10.15006autoplot(arima111101s12_f_cv, type = "l", fcol = 3, flty=1, flwd=1)
# Gráficos comparativos
par(mfrow=c(1,2))
plot(arima111101s12_f_cv, type = "l", fcol = 5, flty=1, flwd=1)
plot(arima010200s12_f_cv, type = "l", fcol = 3, flty=1, flwd=1)
# Comparação entre os modelos
accuracy(arima111101s12)## ME RMSE MAE MPE MAPE
## Training set -2.317061e-05 0.005008705 0.003927331 -0.0001957287 0.03908663
## MASE ACF1
## Training set 0.2269195 0.0003510027accuracy(arima010200s12)## ME RMSE MAE MPE MAPE
## Training set 0.0006358402 0.005226679 0.004104881 0.006350898 0.04085532
## MASE ACF1
## Training set 0.2371783 -0.01312832Conclui-se que o modeo sugerido superior a função auto.arima pelos critérios de erro para a previsão com exceção de ME, logo a escolha do SARIMA(1, 1, 1)X(1,0,1)12 é superior para a previsão.