Vamos inicialmente ler o dados de óbitos por COVID no estado do Amazonas:

# leitura dos dados doencas respiratorias estado AM
dados = read.csv("../Dados/obitosCOVID_AM.csv", sep=",", h=T)
names(dados)
##  [1] "date"                                "state"                              
##  [3] "epidemiological_week_2020"           "deaths_total_2019"                  
##  [5] "deaths_total_2020"                   "deaths_covid19"                     
##  [7] "new_deaths_total_2020"               "deaths_indeterminate_2019"          
##  [9] "deaths_indeterminate_2020"           "deaths_others_2019"                 
## [11] "deaths_others_2020"                  "deaths_pneumonia_2019"              
## [13] "deaths_pneumonia_2020"               "deaths_respiratory_failure_2019"    
## [15] "deaths_respiratory_failure_2020"     "deaths_sars_2019"                   
## [17] "deaths_sars_2020"                    "deaths_septicemia_2019"             
## [19] "deaths_septicemia_2020"              "new_deaths_covid19"                 
## [21] "new_deaths_indeterminate_2020"       "new_deaths_others_2020"             
## [23] "new_deaths_pneumonia_2020"           "new_deaths_respiratory_failure_2020"
## [25] "new_deaths_sars_2020"                "new_deaths_septicemia_2020"
dados2=dados[order(dados$date),]

Selecionamos agora somente os dias que não têm dados faltantes. Em tese seria melhor inputar os dados faltantes utilizando algum método adequado, mas como nosso objetivo nesta aula é apenas entender os passos na modelagem de ST, não discutiremos isto aqui.

#serie Covid
covid19=dados2$new_deaths_covid19[!is.na(dados2$new_deaths_covid19)]

Plotando a série, notamos que a série é não estacionária e possui heterocedasticidade.

# Descritivas
plot.ts(dados2$new_deaths_covid19)

plot.ts(covid19)

#ACF, PACF
ppacf(covid19, new.window = FALSE)

Vamos tomar o log na tentativa de controlar a heterocedasticidade. Em seguida, vamos testar se existe raiz unitária nesta série do log:

# tomando o log

lcovid=log(covid19)
ppacf(lcovid, new.window = FALSE)

ru.lcovid=ur.df(lcovid, type = "none", lags = 2)
print.adftest(ru.lcovid)
## *** Model: diff(X[t]) = a0 + d.x[t-1] + e[t], e[t] ~ WN ***
## 
## Testing H01: d = 0, there is 1 unit root (tau1).
## 
##      Statistics tau_5pct Result
## tau1   -0.96873    -1.95  Ac H0
## 
##  Signif. codes: '*' Rej H0 (sig=10%)    '**' Rej H0 (sig=5%)    '***' Rej H0 (sig=1%) 
## 
## For more details, see ENDERS, W. (2015), Applied Econometric Time Series.

Como identificamos R.U. na série, tomemos a diferença e se o processo resultante for estacionário, procedemos para a identificação do processo:

dlcovid=diff(lcovid)
ppacf(dlcovid,lags = 1:60, new.window = FALSE)

Inicialmente, escolhemos o modelo MA(1) para a série de log-diferenças do número de óbito por COVID no Amazonas, com base no ACF e PACF que vimos acima. \[ \Delta log\left(covi19_t\right)=\theta \epsilon_{t-1} \]

Ajustaremos também alguns outros modelos razoáveis para esta série, e compararemos seu ajuste para escolher o melhor modelo.

fit1 = arima(dlcovid, order=c(0,0,1))
fit2 = arima(lcovid, order=c(0,1,1))
fit3 = arima(lcovid, order=c(0,1,3))
fit4 = arima(lcovid, order=c(1,1,1))
fit5 = arima(dlcovid, order=c(0,0,2))
predfit2=predict(fit2,n.ahead = 10)
predfit5=predict(fit5,n.ahead = 10)

Alternativamente, poderiamo utilizar o pacote forecast e a função auto.arima

auto.arima(lcovid)
## Series: lcovid 
## ARIMA(2,1,3) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2      ma3
##       1.2150  -0.5384  -2.0361  1.4643  -0.3010
## s.e.  0.1667   0.1476   0.1779  0.3177   0.1643
## 
## sigma^2 estimated as 0.2444:  log likelihood=-193.58
## AIC=399.16   AICc=399.47   BIC=420.81

Poderiamos também separar a base em duas para avaliar a adequação da previsão:

lcovid_mod=lcovid[1:250]
lcovid_val=lcovid[251:274]

fit2.1 = arima(lcovid_mod, order=c(0,1,1))

predfit2.1=predict(fit2.1,n.ahead = 24)

auto.arima(lcovid_mod)
## Series: lcovid_mod 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1     ma2     ma3
##       0.7024  -1.5210  0.4818  0.2051
## s.e.  0.0825   0.0913  0.1383  0.0747
## 
## sigma^2 estimated as 0.2374:  log likelihood=-173.36
## AIC=356.72   AICc=356.97   BIC=374.31
autofit=arima(lcovid_mod, order=c(1,1,3))
predautofit=predict(autofit,n.ahead = 24)

Gráficos das previsões:

plot.ts(cbind(lcovid_val,predfit2.1$pred,predautofit$pred), 
               plot.type = "single",
               col=c("black", "red", "blue"), lty=c(1,2,2))
abline(h=mean(lcovid_val))