Cleaning

rm(list=ls())

packages

Instalando alguns pacotes padrões

install.packages("readxl", repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'readxl' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'readxl'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\Thomas\AppData\Local\R\win-library\4.2\00LOCK\readxl\libs\x64\readxl.dll
## to C:\Users\Thomas\AppData\Local\R\win-library\4.2\readxl\libs\x64\readxl.dll:
## Permission denied
## Warning: restored 'readxl'
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages("readr",repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'readr' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages("dplyr", repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'dplyr' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'dplyr'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\Thomas\AppData\Local\R\win-library\4.2\00LOCK\dplyr\libs\x64\dplyr.dll
## to C:\Users\Thomas\AppData\Local\R\win-library\4.2\dplyr\libs\x64\dplyr.dll:
## Permission denied
## Warning: restored 'dplyr'
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('tidyverse', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'tidyverse' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('BETS', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'BETS' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('urca', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'urca' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'urca'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\Thomas\AppData\Local\R\win-library\4.2\00LOCK\urca\libs\x64\urca.dll
## to C:\Users\Thomas\AppData\Local\R\win-library\4.2\urca\libs\x64\urca.dll:
## Permission denied
## Warning: restored 'urca'
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('TSA', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'TSA' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('lmtest', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'lmtest' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'lmtest'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\Thomas\AppData\Local\R\win-library\4.2\00LOCK\lmtest\libs\x64\lmtest.dll
## to C:\Users\Thomas\AppData\Local\R\win-library\4.2\lmtest\libs\x64\lmtest.dll:
## Permission denied
## Warning: restored 'lmtest'
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('aTSA', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'aTSA' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
install.packages('forecast', repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/Thomas/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)
## package 'forecast' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'forecast'
## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## problem copying C:\Users\Thomas\AppData\Local\R\win-
## library\4.2\00LOCK\forecast\libs\x64\forecast.dll to C:
## \Users\Thomas\AppData\Local\R\win-library\4.2\forecast\libs\x64\forecast.dll:
## Permission denied
## Warning: restored 'forecast'
## 
## The downloaded binary packages are in
##  C:\Users\Thomas\AppData\Local\Temp\RtmpacPDMy\downloaded_packages
library(readxl)
## Warning: package 'readxl' was built under R version 4.2.2
library(readr)
## Warning: package 'readr' was built under R version 4.2.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.2
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ stringr 1.4.1
## ✔ tidyr   1.2.1     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(BETS)
## Warning: package 'BETS' was built under R version 4.2.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## 
## Attaching package: 'BETS'
## 
## The following object is masked from 'package:stats':
## 
##     predict
library(urca)
## Warning: package 'urca' was built under R version 4.2.2
library(TSA)
## Warning: package 'TSA' was built under R version 4.2.2
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## 
## The following object is masked from 'package:readr':
## 
##     spec
## 
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## 
## The following object is masked from 'package:utils':
## 
##     tar
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.2.2
## Carregando pacotes exigidos: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(aTSA)
## 
## Attaching package: 'aTSA'
## 
## The following object is masked from 'package:graphics':
## 
##     identify
library(forecast)
## Warning: package 'forecast' was built under R version 4.2.2
## 
## Attaching package: 'forecast'
## 
## The following object is masked from 'package:aTSA':
## 
##     forecast

importing data

importando a planilha para o R

df_projeto <- read_excel("C:/Users/Thomas/Desktop/df projeto.xlsx")
## New names:
## • `` -> `...1`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`

Filtering data

Vou enxutar a base deixando só os dados que me interessam por ora.

aproveito também para tirar a primeira linha que pode atrapalhar.

dfpo1 <- select(df_projeto,3)
dfpo1 <- dfpo1 %>% slice(-1)

Time series

Preciso transformar em serie temporal para poder rodar os códigos que quero

Obs: Já faço a filtragem para a data determinada pelo exercício.

Ploto a série para ter uma primeira ideia do que estou lidando.

tsdf <- ts(dfpo1, start=c(2012, 03), end=c(2019, 12), frequency=12)
plot.ts(tsdf,col=2,main="População ocupada",ylab="Ocupação",xlab="Mês")

Inicios de análise

Podemos ver que há tendência de crescimento.

Não consigo afirmar se há heterocidacidade.

Me parece haver sazonalidade, pois há picos de tempos em tempos.

Utilizarei o “monthplot” para averiguar a sazonalidade.

monthplot(tsdf)

### Vejo que março é um mês onde a ocupação tem queda relevante, portanto há sazonalidade. PArece tambem que do primeiro para o segundo trimestre do ano a ocupação é menor, isso dialoga com o comercio no Brasil que é mais forte no final do ano devido aos feriados.

Vou decompor as plotagens para averiguar os outros fatores.

plot(decompose(tsdf))

A serie aparenta ter tendencia de crescimento e drift, pois ela oscila de tempos em tempos.

A sazonalidade é bem característica e determinada.

O erro está longe de ser um Ruído branco, teremos de tratar isso.

adf.test(tsdf)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 1.041   0.918
## [2,]   1 0.126   0.679
## [3,]   2 0.222   0.706
## [4,]   3 0.556   0.801
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -0.976   0.704
## [2,]   1 -1.792   0.410
## [3,]   2 -1.651   0.464
## [4,]   3 -1.274   0.600
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -1.22   0.900
## [2,]   1 -2.13   0.518
## [3,]   2 -1.98   0.580
## [4,]   3 -1.56   0.756
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Tratamento de estacionariedade

Como estudado na matéria não podemos rodar um SARIMA em cima desses dados pois aparentam não serem estacionários.

Vamos tratar isso futuramente.

acf(tsdf, lag.max = 24)

As autocorrelações caem de forma lenta, portanto há mais um indicativo de que não é estacionário. Uma vez que um dos pontos da estacionariedade é que “t” não deve explicar “Y”.

pacf(tsdf, lag.max = 24)

Vou fazer um teste ADF para averiguar se há raiz unitária, pois temos muitas evidência que há.

Em resumo aqui estamos testando se a série é ou não estacionária, se em Yt = a +pYt-1 +Et, p = ou > 1 então não é estacionário.

adf.trend <- ur.df(tsdf, type = c('trend'), lags = 12)
adf.trend@teststat
##                tau3     phi2     phi3
## statistic -3.271578 3.730011 5.360429
adf.trend@cval
##       1pct  5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2  6.50  4.88  4.16
## phi3  8.73  6.49  5.47

Utilizou se lag = 12, pois vimos que a sazonalidade é anual, com queda em março.

Podemos ver que rejeitamos a hipótese nula a 10%. Ou seja, Podemos dizer que não há raiz unitária a 10% de confiança. Entretanto vendo as autocorrelações sabemos que há raiz unitária. O que deve explicar esse resultado é que a tendencia varia, e isto deve confundir o teste, criando a impressão que é estacionário. No futuro teremos que olhar as diferenças para anular a tendencia.

Para maior compreensão vamos olhar agora com drift.

adf.drift <- ur.df(tsdf, type = c('drift'), lags = 12)
adf.drift@teststat
##                tau2     phi1
## statistic -3.029467 4.821761
adf.drift@cval
##       1pct  5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1  6.70  4.71  3.86

Neste cenário rejeitamos a hipótese nula a 5%.

Vou analisar as diferenciações da serie para anular o efeito da tendência. Vou plotar até 4 diferenciações e ver como s ecomportam

Importante lembrar que queremos sempre trabalhar com o mais parcimonioso, então vamos tentar utilizar menos diferenciações possiveis..

ts.plot(diff(tsdf))

ts.plot(diff(diff(tsdf)))

ts.plot(diff(diff(diff(tsdf))))

ts.plot(diff(diff(diff(diff(tsdf)))))

A primeria diferenciação já resolveu o problema da tendencia, enportanto vamos utilizar ela.

Vou aplicar log nessas séries para ver se corrige a heterocidacidade.

ts.plot(diff(log(tsdf)))

Me parece ter melhorado um pouco em relação a anterior. Vou portanto utilizar a diferença do log.

Vamos agora olhar as autocorrelações para ver se corrigimos o problema.

acf(diff(log(tsdf)), lag.max = 36)

A primeira diferença parece ter corrigido o problema da autocorrelação, mas não corrigiu a autocorrelação sazonal. Já a terceira diferença corrige os dois. Vou fazer o teste dos dois e ver se resolvem o problema da raiz unitária.

adf.trend.tsdf <- ur.df(diff(log(tsdf)), type = "trend" , lag = 12,  selectlags = "AIC")
adf.trend.tsdf@teststat
##                tau3     phi2     phi3
## statistic -5.604845 10.55467 15.80068
adf.trend.tsdf@cval
##       1pct  5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2  6.50  4.88  4.16
## phi3  8.73  6.49  5.47

O resultado é o esperado, corrigimos o problema da raiz unitária, podemos rejeitar a hipótese nula com significÃncia de 99%

Modelando a série

Vou agora analisar o ACF e PACF para ter uma noção de qual o melhor modelo SARIMA usar.

acf(diff(log(tsdf)), lag.max = 36)

pacf(diff(log(tsdf)), lag.max = 36)

Chego a conclusão que o modelo a ser usado é SARIMA (1,1,0)(1,0,0)[12]

Vou utilizar o comando auto.arima para ver qual conclusão ele chega

forecast::auto.arima(log(tsdf))
## Series: log(tsdf) 
## ARIMA(2,1,0)(2,0,0)[12] 
## 
## Coefficients:
##          ar1      ar2    sar1    sar2
##       0.6492  -0.1723  0.2816  0.1986
## s.e.  0.1064   0.1060  0.1048  0.1077
## 
## sigma^2 = 0.01202:  log likelihood = 74.2
## AIC=-138.4   AICc=-137.71   BIC=-125.74

Ele descreve que o modelo (2,1,0)(2,0,0)[12]

vou analisar sobre os dois modelos

fit1.sarima <- Arima(log(tsdf), order = c(1,1,0), seasonal = c(1,0,0))
fit1.sarima
## Series: log(tsdf) 
## ARIMA(1,1,0)(1,0,0)[12] 
## 
## Coefficients:
##          ar1    sar1
##       0.5291  0.4005
## s.e.  0.0872  0.0963
## 
## sigma^2 = 0.01245:  log likelihood = 71.77
## AIC=-137.54   AICc=-137.27   BIC=-129.94
BETS::t_test(fit1.sarima, alpha = 0.05)
##         Coeffs Std.Errors        t Crit.Values Rej.H0
## ar1  0.5291498 0.08721198 6.067398    1.986377   TRUE
## sar1 0.4005210 0.09632992 4.157805    1.986377   TRUE
fit2.sarima <- Arima(log(tsdf), order = c(2,1,0), seasonal = c(2,0,0))
fit2.sarima
## Series: log(tsdf) 
## ARIMA(2,1,0)(2,0,0)[12] 
## 
## Coefficients:
##          ar1      ar2    sar1    sar2
##       0.6492  -0.1723  0.2816  0.1986
## s.e.  0.1064   0.1060  0.1048  0.1077
## 
## sigma^2 = 0.01202:  log likelihood = 74.2
## AIC=-138.4   AICc=-137.71   BIC=-125.74
BETS::t_test(fit1.sarima, alpha = 0.05)
##         Coeffs Std.Errors        t Crit.Values Rej.H0
## ar1  0.5291498 0.08721198 6.067398    1.986377   TRUE
## sar1 0.4005210 0.09632992 4.157805    1.986377   TRUE

A soma do ar1 com ar2 é muito próximo do ar1 do modelo que eu estimei,isso vale para a parte sazonal, portanto vou utilizar o que eu estimei pois é mais parcimonioso.

Resíduos

Vamos analisar os reíduos para ver se são similares a um ruido branco, caso isso não seja verdade, podemos estar excluindo fatores importantes do modelo.

plot(fit1.sarima$residuals)

plot(fit2.sarima$residuals)

Eles são similares a um ruído branco.

Projeção

plot(forecast(object = fit1.sarima, h=12, level=0.95))

tsdf2 <- ts(dfpo1, start=c(2012, 03), end=c(2022, 08), frequency=12)
plot(tsdf2)

Fiz o forecast da série e comparei com a realidade, e a projeção erra no curto prazo, isso ocorre por conta da pandemia, que é um fator exógeno, o domelo SARIMA não considera fatores exógenos. Entretanto se considerarmos que após o fim da pandemia (ao menos considerando que a pior fase da pandemia durou até fim de 2021)o nivel de ocupação voltou ao que era o modelo acertou, um detalhe relevante é que para uma confiança alta de 95% temos que projetar um intervalo d econfiança muito amplo, o que não é positivo, mas isso se deve a série ser relativamente pequena. Outra questão, é que fora o choque exógeno da oandemia outros choque ocorreram, e o modelo tambem não considera. portanto, temos vantagens em trabalhar com o SARIMA, e aparentemente no longo prazo ele possui boa precisão, mas para periodos mais curtos, ou em séries que sofrem efeito de outras variáveis que não da própria série ele não é eficiente.