SARIMA para temperaturas-3.r

# https://cran.r-project.org/web/views/TimeSeries.html
# https://www.datacamp.com/tracks/time-series-with-r
# https://www.datascience.com/blog/introduction-to-forecasting-with-arima-in-r-learn-data-science-tutorials
# https://www.otexts.org/fpp/8/7

# install.packages('robust')
#install.packages('TSPred')

library(TSPred)
library(hexbin)
library(robust)
library(imputeTS)
library(tseries)
library(forecast)
library(fUnitRoots)
library(portes)
library(nortest)
library(tsoutliers)
library(boot)

# T <- read.xls("Dados das estacoes com 2018.xlsx", sheet = 1, header = TRUE, stringsAsFactors = FALSE)
T <- read.csv("Dados das estacoes com 2018.csv")

head(T)

sapply(T, class) # Verificando os tipos

T_INMET <- ts(T$T_INMET,start=c(1961,1),f=12)
plot(T_INMET)

T_ICEA <- ts(T$T_ICEA,start=c(1961,1),f=12)
plot(T_ICEA)

mod1 = lmRob(T_INMET ~ T_ICEA, data = T)
summary(mod1)

bin <- hexbin(T_INMET, T_ICEA, xbins=30)
plot(bin)

aux <- T[is.na(T$T_INMET)&!is.na(T$T_ICEA),]
head(aux)

aux <- predict(mod1, newdata=T[is.na(T$T_INMET)&!is.na(T$T_ICEA),]) # Estima INMET com base no ICEA

head(aux)

# T[is.na(T$T_INMET)&!is.na(T$T_ICEA),] <- aux

T[is.na(T$T_INMET)&!is.na(T$T_ICEA),'T_INMET'] <- aux

T_INMET <- ts(T$T_INMET,start=c(1961,1),f=12)
plot(T_INMET)

# As falhas que sobraram foram preenchidas pelo estimador de Kalman.
T_INMET <- na.kalman(T_INMET, model = "auto.arima")
plot(T_INMET)

UR_INMET <- ts(T$UR_INMET,start=c(1961,1),f=12)
plot(UR_INMET)

UR_ICEA <- ts(T$UR_ICEA,start=c(1961,1),f=12)
plot(UR_ICEA)

mod2 = lmRob(UR_INMET ~ UR_ICEA, data = T)
summary(mod2)

bin <- hexbin(UR_INMET, UR_ICEA, xbins=30)
plot(bin)

aux <- T[is.na(T$UR_INMET)&!is.na(T$UR_ICEA),]
head(aux)

aux <- predict(mod2, newdata=T[is.na(T$UR_INMET)&!is.na(T$UR_ICEA),]) # Estima INMET com base no ICEA

head(aux)

T[is.na(T$UR_INMET)&!is.na(T$UR_ICEA),'UR_INMET'] <- aux

UR_INMET <- ts(T$UR_INMET,start=c(1961,1),f=12)
plot(UR_INMET)

# As falhas que sobraram foram preenchidas pelo estimador de Kalman.
UR_INMET <- na.kalman(UR_INMET, model = "auto.arima")
plot(UR_INMET)

install.packages('weathermetrics')

library(weathermetrics)

it <- heat.index(t=T_INMET, rh=UR_INMET, temperature.metric="celsius", output.metric="celsius", round=2)
it <- ts(it,start=c(1961,1),f=12)
plot(it)

# Decomposição da série temporal (aditiva, possivelmente baseada na transformáda rápida de Fourier, e apenas para fins de análise preliminar)
plot(decompose(it))

urkpssTest(it, type = c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)

adf.test(it, alternative="stationary")

acf(it)
pacf(it)

ndiffs(it) # Valor de d para fazer a série estacionária.

# Calculando as correlações para uma diferença de ordem 1
acf(diff(it, differences=1))
pacf(diff(it, differences=1))

ajuste_0 <- Arima(it, c(1, 1, 2)) # Sem sazonalidade
ajuste_0

ajuste_1 <- Arima(it, c(1, 1, 2), seasonal=list(order = c(1, 1, 2), period = 12)) # Com sazonalidade
ajuste_1

# Automatizando
ajuste_auto <- auto.arima(it, stepwise=FALSE, approximation=FALSE, trace=TRUE)
ajuste_auto

acf(ajuste_auto$residuals)
pacf(ajuste_auto$residuals)

qqnorm(ajuste_auto$residuals)
qqline(ajuste_auto$residuals)
lillie.test(ajuste_auto$residuals)

l <- BoxCox.lambda(it)
l

ajuste_auto_bc <- auto.arima(it, stepwise=FALSE, approximation=FALSE, trace=TRUE, lambda=l)
ajuste_auto_bc

qqnorm(ajuste_auto_bc$residuals)
qqline(ajuste_auto_bc$residuals)
lillie.test(ajuste_auto_bc$residuals)
hist(ajuste_auto_bc$residuals)

# A função não possui o parâmetro lambda
# outliers <- tso(BoxCox(it,l),tsmethod = "auto.arima", args.tsmethod = list(lambda=l), types = c("AO","LS","TC"))

outliers <- tso(BoxCox(it,l),tsmethod = "arima", args.tsmethod = list(order=c(1, 0, 1), seasonal=list(order=c(0, 1, 1), period=12)), types = c("AO","LS","TC"), maxit.iloop=50)

outliers

it[outliers$outliers$ind] <- NA

# Substituindo os outliers por estimativas...
it <- na.kalman(it, model = "auto.arima")

l <- BoxCox.lambda(it)
l

outliers2 <- tso(BoxCox(it,l),tsmethod = "arima", args.tsmethod = list(order=c(1, 0, 2), seasonal=list(order=c(0, 1, 1), period=12)), types = c("AO","LS","TC"), maxit.iloop=50)
outliers2

it[outliers2$outliers$ind] <- NA
it <- na.kalman(it, model = "auto.arima")
l <- BoxCox.lambda(it)
l

outliers3 <- tso(BoxCox(it,l),tsmethod = "arima", args.tsmethod = list(order=c(1, 0, 2), seasonal=list(order=c(0, 1, 1), period=12)), types = c("AO","LS","TC"), maxit.iloop=50)
outliers3

# Reajustando
ajuste_auto_bc_outliers <- auto.arima(it, stepwise=FALSE, approximation=FALSE, trace=TRUE, lambda=l)
ajuste_auto_bc_outliers

qqnorm(ajuste_auto_bc_outliers$residuals)
qqline(ajuste_auto_bc_outliers$residuals)
lillie.test(ajuste_auto_bc_outliers$residuals)
hist(ajuste_auto_bc_outliers$residuals)

LjungBox(ajuste_auto_bc_outliers, lag=c(5, 10, 20, 25, 39, 35, 40, 45, 50), season=12) # Testa se existe autocorrelação entre os dados considerando a sazonalidade

previsoes_auto_bc_outliers <- forecast(ajuste_auto_bc_outliers, h=3)
previsoes_auto_bc_outliers
accuracy(previsoes_auto_bc_outliers)

mod2 = lmRob(ajuste_auto_bc_outliers$x ~ ajuste_auto_bc_outliers$fitted, data = T)
summary(mod2)

options(repr.plot.width=8, repr.plot.height=8)
bin <- hexbin(ajuste_auto_bc_outliers$fitted, ajuste_auto_bc_outliers$x, xbins=30)
plot(bin)

options(repr.plot.width=16, repr.plot.height=8)
ts.plot(ajuste_auto_bc_outliers$fitted,ajuste_auto_bc_outliers$x, lty=c(2,1), gpars = list(col = c("blue", "black")))

plotarimapred(it, ajuste_auto_bc_outliers, xlim=c(2010, 2021))

AIC_Boot <- function(ts) {
    aux <- Arima(ts, lambda=l, order=c(1, 0, 0), seasonal=list(order=c(0, 1, 1), period=12))
    return(aux$aic)
}

aux <- Arima(it, lambda=l, order=c(1, 0, 0), seasonal=list(order=c(0, 1, 1), period=12))

AIC_Boot(it)

# Bootstrap estacionário com bloco de 20
boot_ajuste_auto_bc_outliers <- tsboot(it, AIC_Boot, R = 1000, l = 20, sim = "geom")

boot_ajuste_auto_bc_outliers

boot.ci(boot_ajuste_auto_bc_outliers)