Código.
Faz a leitura de pacotes.
library(xts)
## Carregando pacotes exigidos: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
##
## first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(forecast)
library(tsbox)
## Warning: package 'tsbox' was built under R version 4.1.3
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.1.3
##
## Attaching package: 'TSstudio'
## The following object is masked from 'package:tsbox':
##
## ts_plot
Seta a pasta e e importa dados.
# Seta a pasta
setwd('C:/Users/Usuario/Documents/covid_clima_R')
# Importa dados.
covid_barret <- read.table("dados_covid_Barretos.csv", header = TRUE, quote = "\"", sep = ";", dec = ",")
Efetua análises descritivas nos dados para identificar possíveis inconsistências.
# Identifica a existência de NAs na coluna numérica.
sum(is.na(covid_barret$casosNovos))
## [1] 0
## Também vemos dados faltantes da coluna de tipo character.
sum(covid_barret$Data_Medicao == "")
## [1] 0
# Visualiza summary para checar se existem valores suspeitos. De acordo com outras fontes de dados de COVID, os dados não aparentam suspeitos.
summary(covid_barret)
## Data_Medicao casosNovos
## Length:148 Min. : 0.00
## Class :character 1st Qu.: 7.50
## Mode :character Median : 30.00
## Mean : 33.26
## 3rd Qu.: 48.25
## Max. :193.00
Trata do DF covid_barret para obtermos objeto que possibilite modelagem temporal.
# Transforma a coluna que trata de datas para formato de data.
covid_barret[,1] <- as.Date(covid_barret[,1], '%d/%m/%Y')
# Temos melhor noção da extensão temporal da série.
summary(covid_barret$Data_Medicao)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "2020-05-25" "2020-06-30" "2020-08-06" "2020-08-06" "2020-09-12" "2020-10-19"
# Transforma o DF para formato de xts.
covid_barret_xts <- xts(covid_barret$casosNovos, covid_barret$Data_Medicao)
# Renomeia coluna de covid_barret_ts.
colnames(covid_barret_xts) <- "casos_novos"
Visualiza a série temporal.
covid_barret_xts %>% plot.ts(main="Novos casos diários de COVID-19 ")
Neste ponto, podemos selecionar o modelo ARIMA mais adequado à série temporal de duas formas: uma mais manual e outro mais programática/algorítmica.
Primeira descrevendo a manual, é investigada a estacionaridade da série, isso formalmente por meio de testes (como Dickey-Fuller aumentado; Ljung-Box e KPSS). Caso não constatada, a série deve ser transformada (usualmente por diferenças ou log) até que seja obtida estacionaridade. Caso contatada estacionaridade, manter a série temporal que a apresenta.
Em sequência, com base em plots de ACF e PACF, são selecionadas ordens dos termos AR e MA. Os modelos selecionados são estimados. O melhor modelo é selecionado critério de informação (como AIC, BIC, AICc).
A validade deste modelo selecionado é conferida por teste que investiga propriedades de ruído branco, usualmente focalizando autocorrelação nos resíduos. Caso a última seja existente, deve-se novamente olhar ACF e PACF para selecionar lags e reestimar modelo até que os resíduos dele não apresentem autocorrelação. Após isso ser satisfeito, é feita previsão adotando este último modelo (HYNDMAN; ATHANASOPOULOS, 2018; DEHESH; MARDANI-FARD; DEHESH, 2020).
Todavia, este processo iterativo pode tornar-se bastante longo devido à grande quantidade de possibilidades existentes de combinações de AR, MA e de integração. Ainda, a seleção de lags de AR e MA com base nos mencionados plots usualmente não segue uma regra formal consolidada, tornando a escolha do modelo não tão simples e precisa.
Nesse contexto, é sugerido usar a função auto.arima no R. Tal função é um algoritmo que combina testes de raiz unitária, minimização do AICc e MLE para obter um modelo ARIMA (HYNDMAN; ATHANASOPOULOS, 2018). Essa função adota procedimentos similares ao descrito para selecionar o melhor modelo. Todavia, tal função não cobre a investigação de autocorrelação nos resíduos, assim, esta testagem deve ser efetuada a parte.
Empregamos a função auto.arima e testamos o modelo selecionado.
# O modelo ARIMA é estimado. Os parâmetros são selecionados para buscar em todos os modelos (stepwise = FALSE) e também são evitadas aproximações que poderiam impossibilitar a escolha de modelo de menor AICc (approximation = FALSE). A função está programada inclusive para considerar modelos ARIMA sazonais.
aut_arima <- auto.arima(covid_barret_xts, stepwise = FALSE, approximation = FALSE)
# Vimos summary do modelo selecionado. Ele é um ARIMA(4,1,1).
summary(aut_arima)
## Series: covid_barret_xts
## ARIMA(4,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1
## 0.0809 -0.1880 -0.1527 -0.2573 -0.8804
## s.e. 0.0922 0.0865 0.0919 0.0962 0.0575
##
## sigma^2 estimated as 900: log likelihood=-707.29
## AIC=1426.57 AICc=1427.17 BIC=1444.51
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 2.880538 29.38614 20.4653 -Inf Inf 0.7749611 -0.04931025
# Teste de Ljung-Box. É um teste de portmanteau (Hyndman; Athanasopoulos, 2018). Sua hipótese nula é que se autocorrelações em conjunto são iguais a 0. Como seu valor p > 0.1, a hipótese nula não é rejeitada.
Box.test(residuals(aut_arima), 20)
##
## Box-Pierce test
##
## data: residuals(aut_arima)
## X-squared = 14.426, df = 20, p-value = 0.8083
# Visualiza raízes unitárias. Um indicativo da validade do modelo é elas estarem situadas dentro do círculo unitário.
autoplot(aut_arima)
O teste de Ljung-Box indica a aceitabilidade do modelo (LEVY et al, 2021). Reforçando tal aceitabilidade, constatamos que as raízes AR e MA encontram-se dentro do círculo unitário, indicando que o processo ARIMA é (dada covariância) estacionário e invertível (POURGHASEMI et al., 2020).
Assim, podemos visualizar previsões (neste caso, selecionamos um mês) usando o modelo ARIMA(4,1,1).
autoplot(forecast(aut_arima), h = 30)
