Script

library(readr)
library(mFilter)
library(strucchange)
library(seasonal)
library(ggplot2)
library(urca)
library(forecast)
library(astsa)
library(BETS)

data.set <- read.csv(file.choose(), sep = ";", dec = ".", header = T)
st.rabies <- ts(data.set[3], start = c(2006, 1), end = c(2019, 12), frequency = 12)

decomp <- decompose(st.rabies)
plot(decomp)

beta <- time(st.rabies)
reg <- lm(st.rabies ~ beta)
summary(reg)
plot(st.rabies, xlab <- "Time", ylab = "N. cases", font.lab = 8)
abline(reg, col = "red")

mosum.rabies <- efp(diff(st.rabies)~lag(diff(st.rabies)), type = "Rec-CUSUM")
sctest(mosum.rabies)
bound.rabies <- boundary(mosum.rabies, alpha = 0.05)

ma4 <- filter(x = st.rabies, filter = rep(x = 1/4, times = 4), sides = 2)
ma3 <- filter(x = st.rabies, filter = rep(x = 1/3, times = 3), sides = 2)
plot(st_rabies, xlab = "Time", ylab = "Cases", font.lab = 8)
lines(ma4, col = "blue", lwd = 2)
lines(ma3, col = "red", lwd = 2)

monthplot(st.rabies, ylab = "N. cases", xlab = "Time", font.lab = 8, 
          col.base = "red", lty.base = 1, lwd.base = 2)
abline(h = c(0, 10, 20, 30, 40), col = "darkgray", lty = 2)
legend("topright", bty = "n", cex = 0.8, legend = c("Avarage"), 
       col = c("red"), lwd = 2)

seasonal <- seas(st.rabies)
summary(seasonal)
qs(seasonal)

spec.orig <- data.frame(series(seasonal, "sp0"))
ggplot(aes(x = 0:60, y = X10.Log.Spectrum_AdjOri.), data = spec.orig, 
       colour = "black") + geom_line() + 
       geom_vline(colour = "red", xintercept = c(10, 20, 30, 40, 50), linetype = 5) + 
       geom_vline(colour = "blue", xintercept = c(42, 52), linetype = 3) + 
                  ylab(" ") + xlab(" ") + theme_bw() + 
                  ggtitle("Spectral plot of the first-differenced original series") + 
                  theme(plot.title = element_text(lineheight = 2, face = "bold", size = 8))

zivot.t <- ur.za(st.rabies, model = c("trend"), lag = 12)
summary(zivot.t)
zivot.b <- ur.za(st.rabies, model = c("both"), lag = 12)
summary(zivot.b)
acf2(st.rabies)

ndiffs(st.rabies)
dif.st.rabies <- diff(st.rabies, lag = 1, differences = 1)

zivot.t.dif <- ur.za(dif.st.rabies, model = c("trend"), lag = 12)
summary(zivot.t.dif)
zivot.b.dif <- ur.za(dif.st.rabies, model = c("both"), lag = 12)
summary(zivot.b.dif)
acf2(dif.st.rabies)

arima.1 <- Arima(dif_st_rabies, order = c(6,1,8), seasonal = c(0,0,0))
acf2(arima.1$residuals)
t_test(arima.1)
arima.2 <- Arima(st_rabies, order = c(4,1,4), seasonal = c(0,0,0))
diag <- tsdiag(arima.2, gof.lag = 20)
Box.test(arima.2$residuals, lag = 24, type = "Ljung-Box", fitdf = 1)
shapiro.test(arima.2$residuals)

training <- window(st.rabies, start = c(2006, 1), end = c(2017, 12), frequency = 12)
test <- window(st.rabies, start = c(2018, 1), end = c(2019, 12), frequency = 12)
arima.3 <- Arima(training, order = c(4,1,4), seasonal = c(0,0,0))
fct.arima.3 <- forecast(fct.arima.3, h = 24, level = 95)
accuracy <- accuracy(fct.arima.3, teste)
accuracy

model <- Arima(st.rabies, order = c(4,1,4), seasonal = c(0,0,0), 
               include.mean = FALSE, include.drift = TRUE,
               include.constant = FALSE, lambda = "auto", biasadj = TRUE,
               method = c("CSS-ML"), transform.pars = FALSE, model = NULL,
               fixed = NULL, init = NULL, SSinit = c("Gardner1980"),
               optim.method = "BFGS", optim.control = list(), kappa = 1e6,
               hessian = FALSE)
forecasting <- forecast(model, h = 36, level = 95)
forecasting