Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
562 views
Tweet
Share
Share
Kernel: R (R-Project)
# 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)
Loading required package: fit.models Attaching package: ‘tseries’ The following object is masked from ‘package:imputeTS’: na.remove Loading required package: timeDate Loading required package: timeSeries Loading required package: fBasics Loading required package: parallel 
# 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

Análise da Temperatura

options(repr.plot.width=16, repr.plot.height=8)

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

jpeg("g1.jpg", width = 1000, height = 700)
plot(T_INMET)
dev.off()

plot(T_INMET)
Image in a Jupyter notebook
T_ICEA <- ts(T$T_ICEA,start=c(1961,1),f=12)

jpeg("g2.jpg", width = 1000, height = 700)
plot(T_ICEA)
dev.off()

plot(T_ICEA)
Image in a Jupyter notebook

Observe que as séries possuem falhas. Inicialmente as preencheremos com uma regressão robusta entre INMET e ICEA.

options(repr.plot.width=10, repr.plot.height=10)

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

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

jpeg("g3.jpg", width = 1000, height = 700)
plot(bin)
dev.off()

plot(bin)
Call: lmRob(formula = T_INMET ~ T_ICEA, data = T) Residuals: Min 1Q Median 3Q Max -3.8137 -0.7333 -0.1029 0.6328 5.7032 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.39915 0.77383 5.685 2.24e-08 *** T_ICEA 0.80588 0.02864 28.142 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9937 on 495 degrees of freedom Multiple R-Squared: 0.4702 Test for Bias: statistic p-value M-estimate 9.596 8.246e-03 LS-estimate 45.666 1.213e-10 194 observations deleted due to missingness
Image in a Jupyter notebook
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)
Image in a Jupyter notebook
# As falhas que sobraram foram preenchidas pelo estimador de Kalman.
T_INMET <- na.kalman(T_INMET, model = "auto.arima")

options(repr.plot.width=16, repr.plot.height=8)

jpeg("g4.jpg", width = 1000, height = 700)
plot(T_INMET)
dev.off()

plot(T_INMET)
Image in a Jupyter notebook

Análise da Umidade Relativa

UR_INMET <- ts(T$UR_INMET,start=c(1961,1),f=12)
plot(UR_INMET)
Image in a Jupyter notebook
UR_ICEA <- ts(T$UR_ICEA,start=c(1961,1),f=12)
plot(UR_ICEA)
Image in a Jupyter notebook

Observe que as séries possuem falhas. Inicialmente as preencheremos com uma regressão entre INMET e ICEA.

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

bin <- hexbin(UR_INMET, UR_ICEA, xbins=30)
plot(bin)
Call: lmRob(formula = UR_INMET ~ UR_ICEA, data = T) Residuals: Min 1Q Median 3Q Max -14.8431 -2.7160 -0.1321 2.7989 16.2860 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 14.02245 1.56939 8.935 <2e-16 *** UR_ICEA 0.86527 0.02302 37.593 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.127 on 432 degrees of freedom Multiple R-Squared: 0.5769 Test for Bias: statistic p-value M-estimate 7.301 2.598e-02 LS-estimate 26.278 1.967e-06 257 observations deleted due to missingness
Image in a Jupyter notebook
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)
Image in a Jupyter notebook
# As falhas que sobraram foram preenchidas pelo estimador de Kalman.
UR_INMET <- na.kalman(UR_INMET, model = "auto.arima")
plot(UR_INMET)
Image in a Jupyter notebook

Cálculo do índice térmico

install.packages('weathermetrics')
Installing package into ‘/home/user/R/x86_64-pc-linux-gnu-library/3.4’ (as ‘lib’ is unspecified) 
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)

jpeg("g9.jpg", width = 1000, height = 700)
plot(it)
dev.off()

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

jpeg("g10.jpg", width = 1000, height = 700)
plot(decompose(it))
dev.off()

plot(decompose(it))
Image in a Jupyter notebook

ARIMA

Etapas: a) Verificar se existe a necessidade de uma transformação na série original, com objetivo de estabilizar a variância; b) Tornar a série estacionária por meio de diferenças, de modo que o processo dZt seja reduzido a um ARMA(p,q) c) Identificar o processo ARMA(p,q) resultante. d) Verificação da estacionariedade e da invertibilidade.

FAC : correlação simples entre Zt e Zt – k em função da defasagem k. FACP: correlação entre Zt e Zt – k em função da defasagem k, filtrado o efeito de todas as outras defasagens sobre Zt e Zt – k.

FACP -> AR

D -> I

FAC -> MA

1- Número de AR (auto-regressivo) termos (p): termos AR são apenas defasagens da variável dependente. Por exemplo, se o símbolo p representa 5, os preditores de x (t) irá ser X (t-1) … .x (T-5). 2- Número de MA (média móvel) termos (q): termos MA estão defasados erros de previsão na equação de predição. Por exemplo, se q é 5, os preditores para x (t) será E (t-1) … .e (t-5) onde e (i) é a diferença entre a média móvel ao valor imediato e real. 3- Número de Diferenças (d): Estes são o número de diferenças não sazonal, ou seja, neste caso, tomamos a primeira diferença de ordem. Assim, ou nós podemos passar essa variável e colocar d = 0, ou passar a variável original e coloca -d = 1. Ambos irão gerar mesmos resultados.

Notação: arima(p, d, q), sendo p relacionado a autocorrelação parcial, d a diferença entre os valores, q associado a autocorrelação).

Testes da estacionariedade da série

urkpssTest(it, type = c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)
Title: KPSS Unit Root Test Test Results: NA Description: Thu Oct 18 18:23:32 2018 by user:
Image in a Jupyter notebook
adf.test(it, alternative="stationary")
Warning message in adf.test(it, alternative = "stationary"): “p-value smaller than printed p-value”
Augmented Dickey-Fuller Test data: it Dickey-Fuller = -14.159, Lag order = 8, p-value = 0.01 alternative hypothesis: stationary 
jpeg("g11.jpg", width = 1000, height = 700)
acf(it)
dev.off()
acf(it)

jpeg("g12.jpg", width = 1000, height = 700)
pacf(it)
dev.off()
pacf(it)
Image in a Jupyter notebookImage in a Jupyter notebook

Se rejeitarmos a hipótese da estacionariedade, desta forma, partimos de d=1.

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

# Calculando as correlações para uma diferença de ordem 1
jpeg("g12.jpg", width = 1000, height = 700)
acf(diff(it, differences=1))
dev.off()
acf(diff(it, differences=1))

jpeg("g13.jpg", width = 1000, height = 700)
pacf(diff(it, differences=1))
dev.off()
pacf(diff(it, differences=1))
Image in a Jupyter notebookImage in a Jupyter notebook

Como o ACF cai depois do primeiro lag, podemos partir de p = 1-3. Para o PACF o q poderia ser igual a 2, temos, então, um ARIMA(p, d, q) = ARIMA (1, 1, 2) com os mesmos parâmetros para sazonalidade.

Ajustando os modelos

ajuste_0 <- Arima(it, c(1, 1, 2)) # Sem sazonalidade
ajuste_0
Series: it ARIMA(1,1,2) Coefficients: ar1 ma1 ma2 -0.2744 0.3753 0.0765 s.e. 0.3108 0.3057 0.0477 sigma^2 estimated as 4.389: log likelihood=-1487.86 AIC=2983.72 AICc=2983.78 BIC=3001.87
ajuste_1 <- Arima(it, c(1, 1, 2), seasonal=list(order = c(1, 1, 2), period = 12)) # Com sazonalidade
ajuste_1
Series: it ARIMA(1,1,2)(1,1,2)[12] Coefficients: ar1 ma1 ma2 sar1 sma1 sma2 0.5930 -1.3374 0.3538 -0.9781 -0.0084 -0.9405 s.e. 0.1312 0.1503 0.1430 0.0489 0.0696 0.0716 sigma^2 estimated as 1.513: log likelihood=-1119.62 AIC=2253.23 AICc=2253.4 BIC=2284.87

Utilizando via auto.arima

# Automatizando
ajuste_auto <- auto.arima(it, stepwise=FALSE, approximation=FALSE, trace=TRUE)
ajuste_auto
ARIMA(0,0,0)(0,1,0)[12] : 2686.114 ARIMA(0,0,0)(0,1,0)[12] with drift : 2688.049 ARIMA(0,0,0)(0,1,1)[12] : 2324.46 ARIMA(0,0,0)(0,1,1)[12] with drift : Inf ARIMA(0,0,0)(0,1,2)[12] : 2326.375 ARIMA(0,0,0)(0,1,2)[12] with drift : Inf ARIMA(0,0,0)(1,1,0)[12] : 2477.808 ARIMA(0,0,0)(1,1,0)[12] with drift : 2479.523 ARIMA(0,0,0)(1,1,1)[12] : 2326.363 ARIMA(0,0,0)(1,1,1)[12] with drift : Inf ARIMA(0,0,0)(1,1,2)[12] : Inf ARIMA(0,0,0)(1,1,2)[12] with drift : Inf ARIMA(0,0,0)(2,1,0)[12] : 2394.978 ARIMA(0,0,0)(2,1,0)[12] with drift : 2396.522 ARIMA(0,0,0)(2,1,1)[12] : 2326.994 ARIMA(0,0,0)(2,1,1)[12] with drift : Inf ARIMA(0,0,0)(2,1,2)[12] : 2330.21 ARIMA(0,0,0)(2,1,2)[12] with drift : Inf ARIMA(0,0,1)(0,1,0)[12] : 2640.773 ARIMA(0,0,1)(0,1,0)[12] with drift : 2642.737 ARIMA(0,0,1)(0,1,1)[12] : 2276.591 ARIMA(0,0,1)(0,1,1)[12] with drift : Inf ARIMA(0,0,1)(0,1,2)[12] : 2278.55 ARIMA(0,0,1)(0,1,2)[12] with drift : Inf ARIMA(0,0,1)(1,1,0)[12] : 2439.248 ARIMA(0,0,1)(1,1,0)[12] with drift : 2441.059 ARIMA(0,0,1)(1,1,1)[12] : 2278.543 ARIMA(0,0,1)(1,1,1)[12] with drift : Inf ARIMA(0,0,1)(1,1,2)[12] : Inf ARIMA(0,0,1)(1,1,2)[12] with drift : Inf ARIMA(0,0,1)(2,1,0)[12] : 2366.268 ARIMA(0,0,1)(2,1,0)[12] with drift : 2367.953 ARIMA(0,0,1)(2,1,1)[12] : Inf ARIMA(0,0,1)(2,1,1)[12] with drift : Inf ARIMA(0,0,1)(2,1,2)[12] : 2281.766 ARIMA(0,0,1)(2,1,2)[12] with drift : Inf ARIMA(0,0,2)(0,1,0)[12] : 2640.657 ARIMA(0,0,2)(0,1,0)[12] with drift : 2642.63 ARIMA(0,0,2)(0,1,1)[12] : 2270.698 ARIMA(0,0,2)(0,1,1)[12] with drift : Inf ARIMA(0,0,2)(0,1,2)[12] : 2272.719 ARIMA(0,0,2)(0,1,2)[12] with drift : Inf ARIMA(0,0,2)(1,1,0)[12] : 2435.692 ARIMA(0,0,2)(1,1,0)[12] with drift : 2437.532 ARIMA(0,0,2)(1,1,1)[12] : 2272.718 ARIMA(0,0,2)(1,1,1)[12] with drift : Inf ARIMA(0,0,2)(1,1,2)[12] : Inf ARIMA(0,0,2)(1,1,2)[12] with drift : Inf ARIMA(0,0,2)(2,1,0)[12] : 2360.922 ARIMA(0,0,2)(2,1,0)[12] with drift : 2362.653 ARIMA(0,0,2)(2,1,1)[12] : Inf ARIMA(0,0,2)(2,1,1)[12] with drift : Inf ARIMA(0,0,3)(0,1,0)[12] : 2626.19 ARIMA(0,0,3)(0,1,0)[12] with drift : 2628.18 ARIMA(0,0,3)(0,1,1)[12] : Inf ARIMA(0,0,3)(0,1,1)[12] with drift : Inf ARIMA(0,0,3)(0,1,2)[12] : Inf ARIMA(0,0,3)(0,1,2)[12] with drift : Inf ARIMA(0,0,3)(1,1,0)[12] : 2422.329 ARIMA(0,0,3)(1,1,0)[12] with drift : 2424.215 ARIMA(0,0,3)(1,1,1)[12] : Inf ARIMA(0,0,3)(1,1,1)[12] with drift : Inf ARIMA(0,0,3)(2,1,0)[12] : 2351.283 ARIMA(0,0,3)(2,1,0)[12] with drift : 2353.079 ARIMA(0,0,4)(0,1,0)[12] : 2627.224 ARIMA(0,0,4)(0,1,0)[12] with drift : 2629.224 ARIMA(0,0,4)(0,1,1)[12] : Inf ARIMA(0,0,4)(0,1,1)[12] with drift : Inf ARIMA(0,0,4)(1,1,0)[12] : 2423.745 ARIMA(0,0,4)(1,1,0)[12] with drift : 2425.646 ARIMA(0,0,5)(0,1,0)[12] : 2627.556 ARIMA(0,0,5)(0,1,0)[12] with drift : 2629.565 ARIMA(1,0,0)(0,1,0)[12] : 2634.957 ARIMA(1,0,0)(0,1,0)[12] with drift : 2636.929 ARIMA(1,0,0)(0,1,1)[12] : 2265.29 ARIMA(1,0,0)(0,1,1)[12] with drift : Inf ARIMA(1,0,0)(0,1,2)[12] : Inf ARIMA(1,0,0)(0,1,2)[12] with drift : Inf ARIMA(1,0,0)(1,1,0)[12] : 2431.531 ARIMA(1,0,0)(1,1,0)[12] with drift : 2433.377 ARIMA(1,0,0)(1,1,1)[12] : Inf ARIMA(1,0,0)(1,1,1)[12] with drift : Inf ARIMA(1,0,0)(1,1,2)[12] : Inf ARIMA(1,0,0)(1,1,2)[12] with drift : Inf ARIMA(1,0,0)(2,1,0)[12] : 2359.552 ARIMA(1,0,0)(2,1,0)[12] with drift : 2361.286 ARIMA(1,0,0)(2,1,1)[12] : Inf ARIMA(1,0,0)(2,1,1)[12] with drift : Inf ARIMA(1,0,0)(2,1,2)[12] : Inf ARIMA(1,0,0)(2,1,2)[12] with drift : Inf ARIMA(1,0,1)(0,1,0)[12] : 2629.704 ARIMA(1,0,1)(0,1,0)[12] with drift : 2631.698 ARIMA(1,0,1)(0,1,1)[12] : Inf ARIMA(1,0,1)(0,1,1)[12] with drift : Inf ARIMA(1,0,1)(0,1,2)[12] : Inf ARIMA(1,0,1)(0,1,2)[12] with drift : Inf ARIMA(1,0,1)(1,1,0)[12] : 2421.485 ARIMA(1,0,1)(1,1,0)[12] with drift : 2423.41 ARIMA(1,0,1)(1,1,1)[12] : Inf ARIMA(1,0,1)(1,1,1)[12] with drift : Inf ARIMA(1,0,1)(1,1,2)[12] : Inf ARIMA(1,0,1)(1,1,2)[12] with drift : Inf ARIMA(1,0,1)(2,1,0)[12] : 2347.271 ARIMA(1,0,1)(2,1,0)[12] with drift : 2349.139 ARIMA(1,0,1)(2,1,1)[12] : Inf ARIMA(1,0,1)(2,1,1)[12] with drift : Inf ARIMA(1,0,2)(0,1,0)[12] : 2630.53 ARIMA(1,0,2)(0,1,0)[12] with drift : 2632.532 ARIMA(1,0,2)(0,1,1)[12] : Inf ARIMA(1,0,2)(0,1,1)[12] with drift : Inf ARIMA(1,0,2)(0,1,2)[12] : Inf ARIMA(1,0,2)(0,1,2)[12] with drift : Inf ARIMA(1,0,2)(1,1,0)[12] : 2423.202 ARIMA(1,0,2)(1,1,0)[12] with drift : 2425.136 ARIMA(1,0,2)(1,1,1)[12] : Inf ARIMA(1,0,2)(1,1,1)[12] with drift : Inf ARIMA(1,0,2)(2,1,0)[12] : 2349.278 ARIMA(1,0,2)(2,1,0)[12] with drift : 2351.154 ARIMA(1,0,3)(0,1,0)[12] : 2626.242 ARIMA(1,0,3)(0,1,0)[12] with drift : 2628.245 ARIMA(1,0,3)(0,1,1)[12] : Inf ARIMA(1,0,3)(0,1,1)[12] with drift : Inf ARIMA(1,0,3)(1,1,0)[12] : 2422.703 ARIMA(1,0,3)(1,1,0)[12] with drift : 2424.623 ARIMA(1,0,4)(0,1,0)[12] : 2627.554 ARIMA(1,0,4)(0,1,0)[12] with drift : 2629.567 ARIMA(2,0,0)(0,1,0)[12] : 2634.301 ARIMA(2,0,0)(0,1,0)[12] with drift : 2636.284 ARIMA(2,0,0)(0,1,1)[12] : Inf ARIMA(2,0,0)(0,1,1)[12] with drift : Inf ARIMA(2,0,0)(0,1,2)[12] : Inf ARIMA(2,0,0)(0,1,2)[12] with drift : Inf ARIMA(2,0,0)(1,1,0)[12] : 2427.388 ARIMA(2,0,0)(1,1,0)[12] with drift : 2429.271 ARIMA(2,0,0)(1,1,1)[12] : Inf ARIMA(2,0,0)(1,1,1)[12] with drift : Inf ARIMA(2,0,0)(1,1,2)[12] : Inf ARIMA(2,0,0)(1,1,2)[12] with drift : Inf ARIMA(2,0,0)(2,1,0)[12] : 2353.586 ARIMA(2,0,0)(2,1,0)[12] with drift : 2355.382 ARIMA(2,0,0)(2,1,1)[12] : Inf ARIMA(2,0,0)(2,1,1)[12] with drift : Inf ARIMA(2,0,1)(0,1,0)[12] : 2630.941 ARIMA(2,0,1)(0,1,0)[12] with drift : 2632.943 ARIMA(2,0,1)(0,1,1)[12] : Inf ARIMA(2,0,1)(0,1,1)[12] with drift : Inf ARIMA(2,0,1)(0,1,2)[12] : Inf ARIMA(2,0,1)(0,1,2)[12] with drift : Inf ARIMA(2,0,1)(1,1,0)[12] : 2423.271 ARIMA(2,0,1)(1,1,0)[12] with drift : 2425.205 ARIMA(2,0,1)(1,1,1)[12] : Inf ARIMA(2,0,1)(1,1,1)[12] with drift : Inf ARIMA(2,0,1)(2,1,0)[12] : Inf ARIMA(2,0,1)(2,1,0)[12] with drift : Inf ARIMA(2,0,2)(0,1,0)[12] : Inf ARIMA(2,0,2)(0,1,0)[12] with drift : Inf ARIMA(2,0,2)(0,1,1)[12] : Inf ARIMA(2,0,2)(0,1,1)[12] with drift : Inf ARIMA(2,0,2)(1,1,0)[12] : 2422.676 ARIMA(2,0,2)(1,1,0)[12] with drift : 2424.618 ARIMA(2,0,3)(0,1,0)[12] : 2627.22 ARIMA(2,0,3)(0,1,0)[12] with drift : 2629.234 ARIMA(3,0,0)(0,1,0)[12] : 2624.389 ARIMA(3,0,0)(0,1,0)[12] with drift : 2626.391 ARIMA(3,0,0)(0,1,1)[12] : Inf ARIMA(3,0,0)(0,1,1)[12] with drift : Inf ARIMA(3,0,0)(0,1,2)[12] : Inf ARIMA(3,0,0)(0,1,2)[12] with drift : Inf ARIMA(3,0,0)(1,1,0)[12] : 2419.92 ARIMA(3,0,0)(1,1,0)[12] with drift : 2421.846 ARIMA(3,0,0)(1,1,1)[12] : Inf ARIMA(3,0,0)(1,1,1)[12] with drift : Inf ARIMA(3,0,0)(2,1,0)[12] : 2347.963 ARIMA(3,0,0)(2,1,0)[12] with drift : 2349.818 ARIMA(3,0,1)(0,1,0)[12] : 2625.213 ARIMA(3,0,1)(0,1,0)[12] with drift : 2627.217 ARIMA(3,0,1)(0,1,1)[12] : Inf ARIMA(3,0,1)(0,1,1)[12] with drift : Inf ARIMA(3,0,1)(1,1,0)[12] : 2420.888 ARIMA(3,0,1)(1,1,0)[12] with drift : 2422.81 ARIMA(3,0,2)(0,1,0)[12] : 2626.973 ARIMA(3,0,2)(0,1,0)[12] with drift : 2628.986 ARIMA(4,0,0)(0,1,0)[12] : 2625.894 ARIMA(4,0,0)(0,1,0)[12] with drift : 2627.9 ARIMA(4,0,0)(0,1,1)[12] : Inf ARIMA(4,0,0)(0,1,1)[12] with drift : Inf ARIMA(4,0,0)(1,1,0)[12] : 2421.564 ARIMA(4,0,0)(1,1,0)[12] with drift : 2423.489 ARIMA(4,0,1)(0,1,0)[12] : 2627.148 ARIMA(4,0,1)(0,1,0)[12] with drift : 2629.159 ARIMA(5,0,0)(0,1,0)[12] : 2626.384 ARIMA(5,0,0)(0,1,0)[12] with drift : 2628.399 Best model: ARIMA(1,0,0)(0,1,1)[12]
Series: it ARIMA(1,0,0)(0,1,1)[12] Coefficients: ar1 sma1 0.2988 -0.8859 s.e. 0.0375 0.0245 sigma^2 estimated as 1.592: log likelihood=-1129.63 AIC=2265.25 AICc=2265.29 BIC=2278.82

No caso, obteve um AIC maior, mas com SE menores. Pelo princípio da parcimônia, e sabendo que existe sazonalidade, utilizaremos, de início, o modelo (0,0,2)(0,1,1)12.

acf(ajuste_auto$residuals)
pacf(ajuste_auto$residuals)
Image in a Jupyter notebookImage in a Jupyter notebook
jpeg("g14.jpg", width = 1000, height = 700)
qqnorm(ajuste_auto$residuals)
qqline(ajuste_auto$residuals)
dev.off()

qqnorm(ajuste_auto$residuals)
qqline(ajuste_auto$residuals)

lillie.test(ajuste_auto$residuals)
Lilliefors (Kolmogorov-Smirnov) normality test data: ajuste_auto$residuals D = 0.04776, p-value = 0.0007485
Image in a Jupyter notebook

Não adere a Normal, mas também não apresenta comportamento anômalo. É o caso de tentar uma transformação de Box-Cox.

l <- BoxCox.lambda(it)
l

ajuste_auto_bc <- auto.arima(it, stepwise=FALSE, approximation=FALSE, trace=TRUE, lambda=l)
ajuste_auto_bc
ARIMA(0,0,0)(0,1,0)[12] : -5184.219 ARIMA(0,0,0)(0,1,0)[12] with drift : -5182.264 ARIMA(0,0,0)(0,1,1)[12] : -5533.35 ARIMA(0,0,0)(0,1,1)[12] with drift : Inf ARIMA(0,0,0)(0,1,2)[12] : -5531.396 ARIMA(0,0,0)(0,1,2)[12] with drift : Inf ARIMA(0,0,0)(1,1,0)[12] : -5382.485 ARIMA(0,0,0)(1,1,0)[12] with drift : -5380.732 ARIMA(0,0,0)(1,1,1)[12] : -5531.404 ARIMA(0,0,0)(1,1,1)[12] with drift : Inf ARIMA(0,0,0)(1,1,2)[12] : Inf ARIMA(0,0,0)(1,1,2)[12] with drift : Inf ARIMA(0,0,0)(2,1,0)[12] : -5456.414 ARIMA(0,0,0)(2,1,0)[12] with drift : -5454.814 ARIMA(0,0,0)(2,1,1)[12] : -5531.065 ARIMA(0,0,0)(2,1,1)[12] with drift : Inf ARIMA(0,0,0)(2,1,2)[12] : Inf ARIMA(0,0,0)(2,1,2)[12] with drift : Inf ARIMA(0,0,1)(0,1,0)[12] : -5218.002 ARIMA(0,0,1)(0,1,0)[12] with drift : -5216.026 ARIMA(0,0,1)(0,1,1)[12] : -5567.573 ARIMA(0,0,1)(0,1,1)[12] with drift : Inf ARIMA(0,0,1)(0,1,2)[12] : -5566.239 ARIMA(0,0,1)(0,1,2)[12] with drift : Inf ARIMA(0,0,1)(1,1,0)[12] : -5408.399 ARIMA(0,0,1)(1,1,0)[12] with drift : -5406.568 ARIMA(0,0,1)(1,1,1)[12] : -5566.321 ARIMA(0,0,1)(1,1,1)[12] with drift : Inf ARIMA(0,0,1)(1,1,2)[12] : -5563.524 ARIMA(0,0,1)(1,1,2)[12] with drift : Inf ARIMA(0,0,1)(2,1,0)[12] : -5474.445 ARIMA(0,0,1)(2,1,0)[12] with drift : -5472.733 ARIMA(0,0,1)(2,1,1)[12] : Inf ARIMA(0,0,1)(2,1,1)[12] with drift : Inf ARIMA(0,0,1)(2,1,2)[12] : Inf ARIMA(0,0,1)(2,1,2)[12] with drift : Inf ARIMA(0,0,2)(0,1,0)[12] : -5216.614 ARIMA(0,0,2)(0,1,0)[12] with drift : -5214.63 ARIMA(0,0,2)(0,1,1)[12] : -5571.615 ARIMA(0,0,2)(0,1,1)[12] with drift : Inf ARIMA(0,0,2)(0,1,2)[12] : -5570 ARIMA(0,0,2)(0,1,2)[12] with drift : Inf ARIMA(0,0,2)(1,1,0)[12] : -5410.267 ARIMA(0,0,2)(1,1,0)[12] with drift : -5408.412 ARIMA(0,0,2)(1,1,1)[12] : Inf ARIMA(0,0,2)(1,1,1)[12] with drift : Inf ARIMA(0,0,2)(1,1,2)[12] : Inf ARIMA(0,0,2)(1,1,2)[12] with drift : Inf ARIMA(0,0,2)(2,1,0)[12] : -5478.05 ARIMA(0,0,2)(2,1,0)[12] with drift : -5476.297 ARIMA(0,0,2)(2,1,1)[12] : Inf ARIMA(0,0,2)(2,1,1)[12] with drift : Inf ARIMA(0,0,3)(0,1,0)[12] : -5226.15 ARIMA(0,0,3)(0,1,0)[12] with drift : -5224.152 ARIMA(0,0,3)(0,1,1)[12] : -5579.704 ARIMA(0,0,3)(0,1,1)[12] with drift : Inf ARIMA(0,0,3)(0,1,2)[12] : Inf ARIMA(0,0,3)(0,1,2)[12] with drift : Inf ARIMA(0,0,3)(1,1,0)[12] : -5418.451 ARIMA(0,0,3)(1,1,0)[12] with drift : -5416.558 ARIMA(0,0,3)(1,1,1)[12] : Inf ARIMA(0,0,3)(1,1,1)[12] with drift : Inf ARIMA(0,0,3)(2,1,0)[12] : -5483.142 ARIMA(0,0,3)(2,1,0)[12] with drift : -5481.336 ARIMA(0,0,4)(0,1,0)[12] : -5224.486 ARIMA(0,0,4)(0,1,0)[12] with drift : -5222.48 ARIMA(0,0,4)(0,1,1)[12] : Inf ARIMA(0,0,4)(0,1,1)[12] with drift : Inf ARIMA(0,0,4)(1,1,0)[12] : -5416.737 ARIMA(0,0,4)(1,1,0)[12] with drift : -5414.831 ARIMA(0,0,5)(0,1,0)[12] : -5223.889 ARIMA(0,0,5)(0,1,0)[12] with drift : -5221.875 ARIMA(1,0,0)(0,1,0)[12] : -5220.45 ARIMA(1,0,0)(0,1,0)[12] with drift : -5218.47 ARIMA(1,0,0)(0,1,1)[12] : -5574.769 ARIMA(1,0,0)(0,1,1)[12] with drift : Inf ARIMA(1,0,0)(0,1,2)[12] : Inf ARIMA(1,0,0)(0,1,2)[12] with drift : Inf ARIMA(1,0,0)(1,1,0)[12] : -5412.812 ARIMA(1,0,0)(1,1,0)[12] with drift : -5410.959 ARIMA(1,0,0)(1,1,1)[12] : Inf ARIMA(1,0,0)(1,1,1)[12] with drift : Inf ARIMA(1,0,0)(1,1,2)[12] : -5570.727 ARIMA(1,0,0)(1,1,2)[12] with drift : Inf ARIMA(1,0,0)(2,1,0)[12] : -5478.239 ARIMA(1,0,0)(2,1,0)[12] with drift : -5476.494 ARIMA(1,0,0)(2,1,1)[12] : Inf ARIMA(1,0,0)(2,1,1)[12] with drift : Inf ARIMA(1,0,0)(2,1,2)[12] : Inf ARIMA(1,0,0)(2,1,2)[12] with drift : Inf ARIMA(1,0,1)(0,1,0)[12] : -5222.121 ARIMA(1,0,1)(0,1,0)[12] with drift : -5220.123 ARIMA(1,0,1)(0,1,1)[12] : Inf ARIMA(1,0,1)(0,1,1)[12] with drift : Inf ARIMA(1,0,1)(0,1,2)[12] : Inf ARIMA(1,0,1)(0,1,2)[12] with drift : Inf ARIMA(1,0,1)(1,1,0)[12] : -5420.563 ARIMA(1,0,1)(1,1,0)[12] with drift : -5418.635 ARIMA(1,0,1)(1,1,1)[12] : Inf ARIMA(1,0,1)(1,1,1)[12] with drift : Inf ARIMA(1,0,1)(1,1,2)[12] : Inf ARIMA(1,0,1)(1,1,2)[12] with drift : Inf ARIMA(1,0,1)(2,1,0)[12] : -5487.814 ARIMA(1,0,1)(2,1,0)[12] with drift : -5485.94 ARIMA(1,0,1)(2,1,1)[12] : Inf ARIMA(1,0,1)(2,1,1)[12] with drift : Inf ARIMA(1,0,2)(0,1,0)[12] : -5222.038 ARIMA(1,0,2)(0,1,0)[12] with drift : -5220.031 ARIMA(1,0,2)(0,1,1)[12] : Inf ARIMA(1,0,2)(0,1,1)[12] with drift : Inf ARIMA(1,0,2)(0,1,2)[12] : Inf ARIMA(1,0,2)(0,1,2)[12] with drift : Inf ARIMA(1,0,2)(1,1,0)[12] : -5418.903 ARIMA(1,0,2)(1,1,0)[12] with drift : -5416.966 ARIMA(1,0,2)(1,1,1)[12] : Inf ARIMA(1,0,2)(1,1,1)[12] with drift : Inf ARIMA(1,0,2)(2,1,0)[12] : -5485.803 ARIMA(1,0,2)(2,1,0)[12] with drift : -5483.921 ARIMA(1,0,3)(0,1,0)[12] : -5225.087 ARIMA(1,0,3)(0,1,0)[12] with drift : -5223.078 ARIMA(1,0,3)(0,1,1)[12] : Inf ARIMA(1,0,3)(0,1,1)[12] with drift : Inf ARIMA(1,0,3)(1,1,0)[12] : -5418.182 ARIMA(1,0,3)(1,1,0)[12] with drift : -5416.243 ARIMA(1,0,4)(0,1,0)[12] : -5223.891 ARIMA(1,0,4)(0,1,0)[12] with drift : -5221.873 ARIMA(2,0,0)(0,1,0)[12] : -5219.32 ARIMA(2,0,0)(0,1,0)[12] with drift : -5217.331 ARIMA(2,0,0)(0,1,1)[12] : Inf ARIMA(2,0,0)(0,1,1)[12] with drift : Inf ARIMA(2,0,0)(0,1,2)[12] : Inf ARIMA(2,0,0)(0,1,2)[12] with drift : Inf ARIMA(2,0,0)(1,1,0)[12] : -5415.145 ARIMA(2,0,0)(1,1,0)[12] with drift : -5413.26 ARIMA(2,0,0)(1,1,1)[12] : Inf ARIMA(2,0,0)(1,1,1)[12] with drift : Inf ARIMA(2,0,0)(1,1,2)[12] : Inf ARIMA(2,0,0)(1,1,2)[12] with drift : Inf ARIMA(2,0,0)(2,1,0)[12] : -5482.265 ARIMA(2,0,0)(2,1,0)[12] with drift : -5480.467 ARIMA(2,0,0)(2,1,1)[12] : Inf ARIMA(2,0,0)(2,1,1)[12] with drift : Inf ARIMA(2,0,1)(0,1,0)[12] : -5221.426 ARIMA(2,0,1)(0,1,0)[12] with drift : -5219.42 ARIMA(2,0,1)(0,1,1)[12] : Inf ARIMA(2,0,1)(0,1,1)[12] with drift : Inf ARIMA(2,0,1)(0,1,2)[12] : Inf ARIMA(2,0,1)(0,1,2)[12] with drift : Inf ARIMA(2,0,1)(1,1,0)[12] : -5418.844 ARIMA(2,0,1)(1,1,0)[12] with drift : -5416.907 ARIMA(2,0,1)(1,1,1)[12] : Inf ARIMA(2,0,1)(1,1,1)[12] with drift : Inf ARIMA(2,0,1)(2,1,0)[12] : Inf ARIMA(2,0,1)(2,1,0)[12] with drift : Inf ARIMA(2,0,2)(0,1,0)[12] : Inf ARIMA(2,0,2)(0,1,0)[12] with drift : Inf ARIMA(2,0,2)(0,1,1)[12] : Inf ARIMA(2,0,2)(0,1,1)[12] with drift : Inf ARIMA(2,0,2)(1,1,0)[12] : -5418.63 ARIMA(2,0,2)(1,1,0)[12] with drift : -5416.687 ARIMA(2,0,3)(0,1,0)[12] : -5224.318 ARIMA(2,0,3)(0,1,0)[12] with drift : -5222.298 ARIMA(3,0,0)(0,1,0)[12] : -5226.97 ARIMA(3,0,0)(0,1,0)[12] with drift : -5224.964 ARIMA(3,0,0)(0,1,1)[12] : Inf ARIMA(3,0,0)(0,1,1)[12] with drift : Inf ARIMA(3,0,0)(0,1,2)[12] : Inf ARIMA(3,0,0)(0,1,2)[12] with drift : Inf ARIMA(3,0,0)(1,1,0)[12] : -5420.527 ARIMA(3,0,0)(1,1,0)[12] with drift : -5418.604 ARIMA(3,0,0)(1,1,1)[12] : Inf ARIMA(3,0,0)(1,1,1)[12] with drift : Inf ARIMA(3,0,0)(2,1,0)[12] : -5485.637 ARIMA(3,0,0)(2,1,0)[12] with drift : -5483.787 ARIMA(3,0,1)(0,1,0)[12] : -5226.009 ARIMA(3,0,1)(0,1,0)[12] with drift : -5223.999 ARIMA(3,0,1)(0,1,1)[12] : Inf ARIMA(3,0,1)(0,1,1)[12] with drift : Inf ARIMA(3,0,1)(1,1,0)[12] : -5418.923 ARIMA(3,0,1)(1,1,0)[12] with drift : -5417.005 ARIMA(3,0,2)(0,1,0)[12] : -5224.479 ARIMA(3,0,2)(0,1,0)[12] with drift : -5222.461 ARIMA(4,0,0)(0,1,0)[12] : -5225.24 ARIMA(4,0,0)(0,1,0)[12] with drift : -5223.23 ARIMA(4,0,0)(0,1,1)[12] : Inf ARIMA(4,0,0)(0,1,1)[12] with drift : Inf ARIMA(4,0,0)(1,1,0)[12] : -5418.561 ARIMA(4,0,0)(1,1,0)[12] with drift : -5416.635 ARIMA(4,0,1)(0,1,0)[12] : -5224.07 ARIMA(4,0,1)(0,1,0)[12] with drift : -5222.054 ARIMA(5,0,0)(0,1,0)[12] : -5225.74 ARIMA(5,0,0)(0,1,0)[12] with drift : -5223.72 Best model: ARIMA(0,0,3)(0,1,1)[12]
Series: it ARIMA(0,0,3)(0,1,1)[12] Box Cox transformation: lambda= -0.7557019 Coefficients: ma1 ma2 ma3 sma1 0.2236 0.1118 0.1139 -0.8862 s.e. 0.0388 0.0381 0.0353 0.0249 sigma^2 estimated as 1.526e-05: log likelihood=2794.9 AIC=-5579.79 AICc=-5579.7 BIC=-5557.19
qqnorm(ajuste_auto_bc$residuals)
qqline(ajuste_auto_bc$residuals)
lillie.test(ajuste_auto_bc$residuals)
hist(ajuste_auto_bc$residuals)
Lilliefors (Kolmogorov-Smirnov) normality test data: ajuste_auto_bc$residuals D = 0.048052, p-value = 0.0006699
Image in a Jupyter notebookImage in a Jupyter notebook

Melhorou as estatísticas, mas não a questão da aderência a Normal. Tem outliers...

# 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(0, 0, 3), seasonal=list(order=c(0, 1, 1), period=12)), types = c("AO","LS","TC"), maxit.iloop=50)
Warning message in locate.outliers.oloop(y = y, fit = fit, types = types, cval = cval, : “stopped when 'maxit.oloop = 4' was reached”
outliers
Call: structure(list(method = NULL), .Names = "method") Coefficients: ma1 ma2 ma3 sma1 AO507 0.2195 0.1154 0.1241 -0.8833 -0.0170 s.e. 0.0388 0.0378 0.0355 0.0252 0.0036 sigma^2 estimated as 1.468e-05: log likelihood = 2805.82, aic = -5599.64 Outliers: type ind time coefhat tstat 1 AO 507 2003:03 -0.01704 -4.709
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
Call: structure(list(method = NULL), .Names = "method") Coefficients: ar1 ma1 ma2 sma1 TC88 0.9910 -0.7919 -0.1484 -0.9594 -0.0157 s.e. 0.0118 0.0388 0.0396 0.0317 0.0036 sigma^2 estimated as 1.922e-05: log likelihood = 2709.07, aic = -5406.15 Outliers: type ind time coefhat tstat 1 TC 88 1968:04 -0.01565 -4.379
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
Call: structure(list(method = NULL), .Names = "method") Coefficients: ar1 ma1 ma2 sma1 TC89 0.9907 -0.7903 -0.1505 -0.9571 -0.0093 s.e. 0.0119 0.0388 0.0395 0.0306 0.0023 sigma^2 estimated as 8.025e-06: log likelihood = 3005.94, aic = -5999.88 Outliers: type ind time coefhat tstat 1 TC 89 1968:05 -0.00934 -4.052

Reajustando

# Reajustando
ajuste_auto_bc_outliers <- auto.arima(it, stepwise=FALSE, approximation=FALSE, trace=TRUE, lambda=l)
ajuste_auto_bc_outliers
ARIMA(0,0,0)(0,1,0)[12] : -5595.614 ARIMA(0,0,0)(0,1,0)[12] with drift : -5593.661 ARIMA(0,0,0)(0,1,1)[12] : -5935.103 ARIMA(0,0,0)(0,1,1)[12] with drift : Inf ARIMA(0,0,0)(0,1,2)[12] : -5933.397 ARIMA(0,0,0)(0,1,2)[12] with drift : Inf ARIMA(0,0,0)(1,1,0)[12] : -5784.9 ARIMA(0,0,0)(1,1,0)[12] with drift : -5783.151 ARIMA(0,0,0)(1,1,1)[12] : -5933.43 ARIMA(0,0,0)(1,1,1)[12] with drift : Inf ARIMA(0,0,0)(1,1,2)[12] : -5931.091 ARIMA(0,0,0)(1,1,2)[12] with drift : Inf ARIMA(0,0,0)(2,1,0)[12] : -5857.596 ARIMA(0,0,0)(2,1,0)[12] with drift : -5856.001 ARIMA(0,0,0)(2,1,1)[12] : -5932.721 ARIMA(0,0,0)(2,1,1)[12] with drift : Inf ARIMA(0,0,0)(2,1,2)[12] : Inf ARIMA(0,0,0)(2,1,2)[12] with drift : Inf ARIMA(0,0,1)(0,1,0)[12] : -5624.147 ARIMA(0,0,1)(0,1,0)[12] with drift : -5622.173 ARIMA(0,0,1)(0,1,1)[12] : -5965.234 ARIMA(0,0,1)(0,1,1)[12] with drift : Inf ARIMA(0,0,1)(0,1,2)[12] : -5964.415 ARIMA(0,0,1)(0,1,2)[12] with drift : Inf ARIMA(0,0,1)(1,1,0)[12] : -5806.584 ARIMA(0,0,1)(1,1,0)[12] with drift : -5804.76 ARIMA(0,0,1)(1,1,1)[12] : -5964.537 ARIMA(0,0,1)(1,1,1)[12] with drift : Inf ARIMA(0,0,1)(1,1,2)[12] : Inf ARIMA(0,0,1)(1,1,2)[12] with drift : Inf ARIMA(0,0,1)(2,1,0)[12] : -5872.149 ARIMA(0,0,1)(2,1,0)[12] with drift : -5870.45 ARIMA(0,0,1)(2,1,1)[12] : Inf ARIMA(0,0,1)(2,1,1)[12] with drift : Inf ARIMA(0,0,1)(2,1,2)[12] : Inf ARIMA(0,0,1)(2,1,2)[12] with drift : Inf ARIMA(0,0,2)(0,1,0)[12] : -5622.49 ARIMA(0,0,2)(0,1,0)[12] with drift : -5620.509 ARIMA(0,0,2)(0,1,1)[12] : -5969.018 ARIMA(0,0,2)(0,1,1)[12] with drift : Inf ARIMA(0,0,2)(0,1,2)[12] : -5967.759 ARIMA(0,0,2)(0,1,2)[12] with drift : Inf ARIMA(0,0,2)(1,1,0)[12] : -5808.343 ARIMA(0,0,2)(1,1,0)[12] with drift : -5806.496 ARIMA(0,0,2)(1,1,1)[12] : Inf ARIMA(0,0,2)(1,1,1)[12] with drift : Inf ARIMA(0,0,2)(1,1,2)[12] : -5965.001 ARIMA(0,0,2)(1,1,2)[12] with drift : Inf ARIMA(0,0,2)(2,1,0)[12] : -5875.818 ARIMA(0,0,2)(2,1,0)[12] with drift : -5874.077 ARIMA(0,0,2)(2,1,1)[12] : Inf ARIMA(0,0,2)(2,1,1)[12] with drift : Inf ARIMA(0,0,3)(0,1,0)[12] : -5634.082 ARIMA(0,0,3)(0,1,0)[12] with drift : -5632.085 ARIMA(0,0,3)(0,1,1)[12] : -5978.213 ARIMA(0,0,3)(0,1,1)[12] with drift : Inf ARIMA(0,0,3)(0,1,2)[12] : Inf ARIMA(0,0,3)(0,1,2)[12] with drift : Inf ARIMA(0,0,3)(1,1,0)[12] : -5816.817 ARIMA(0,0,3)(1,1,0)[12] with drift : -5814.928 ARIMA(0,0,3)(1,1,1)[12] : Inf ARIMA(0,0,3)(1,1,1)[12] with drift : Inf ARIMA(0,0,3)(2,1,0)[12] : -5880.987 ARIMA(0,0,3)(2,1,0)[12] with drift : -5879.188 ARIMA(0,0,4)(0,1,0)[12] : -5633.091 ARIMA(0,0,4)(0,1,0)[12] with drift : -5631.085 ARIMA(0,0,4)(0,1,1)[12] : -5977.304 ARIMA(0,0,4)(0,1,1)[12] with drift : Inf ARIMA(0,0,4)(1,1,0)[12] : -5815.636 ARIMA(0,0,4)(1,1,0)[12] with drift : -5813.731 ARIMA(0,0,5)(0,1,0)[12] : -5631.492 ARIMA(0,0,5)(0,1,0)[12] with drift : -5629.478 ARIMA(1,0,0)(0,1,0)[12] : -5626.03 ARIMA(1,0,0)(0,1,0)[12] with drift : -5624.052 ARIMA(1,0,0)(0,1,1)[12] : -5971.74 ARIMA(1,0,0)(0,1,1)[12] with drift : Inf ARIMA(1,0,0)(0,1,2)[12] : Inf ARIMA(1,0,0)(0,1,2)[12] with drift : Inf ARIMA(1,0,0)(1,1,0)[12] : -5810.389 ARIMA(1,0,0)(1,1,0)[12] with drift : -5808.544 ARIMA(1,0,0)(1,1,1)[12] : Inf ARIMA(1,0,0)(1,1,1)[12] with drift : Inf ARIMA(1,0,0)(1,1,2)[12] : Inf ARIMA(1,0,0)(1,1,2)[12] with drift : Inf ARIMA(1,0,0)(2,1,0)[12] : -5875.374 ARIMA(1,0,0)(2,1,0)[12] with drift : -5873.644 ARIMA(1,0,0)(2,1,1)[12] : Inf ARIMA(1,0,0)(2,1,1)[12] with drift : Inf ARIMA(1,0,0)(2,1,2)[12] : Inf ARIMA(1,0,0)(2,1,2)[12] with drift : Inf ARIMA(1,0,1)(0,1,0)[12] : -5629.227 ARIMA(1,0,1)(0,1,0)[12] with drift : -5627.228 ARIMA(1,0,1)(0,1,1)[12] : Inf ARIMA(1,0,1)(0,1,1)[12] with drift : Inf ARIMA(1,0,1)(0,1,2)[12] : Inf ARIMA(1,0,1)(0,1,2)[12] with drift : Inf ARIMA(1,0,1)(1,1,0)[12] : -5819.324 ARIMA(1,0,1)(1,1,0)[12] with drift : -5817.396 ARIMA(1,0,1)(1,1,1)[12] : Inf ARIMA(1,0,1)(1,1,1)[12] with drift : Inf ARIMA(1,0,1)(1,1,2)[12] : Inf ARIMA(1,0,1)(1,1,2)[12] with drift : Inf ARIMA(1,0,1)(2,1,0)[12] : -5885.866 ARIMA(1,0,1)(2,1,0)[12] with drift : -5883.994 ARIMA(1,0,1)(2,1,1)[12] : Inf ARIMA(1,0,1)(2,1,1)[12] with drift : Inf ARIMA(1,0,2)(0,1,0)[12] : -5628.545 ARIMA(1,0,2)(0,1,0)[12] with drift : -5626.539 ARIMA(1,0,2)(0,1,1)[12] : Inf ARIMA(1,0,2)(0,1,1)[12] with drift : Inf ARIMA(1,0,2)(0,1,2)[12] : Inf ARIMA(1,0,2)(0,1,2)[12] with drift : Inf ARIMA(1,0,2)(1,1,0)[12] : -5817.44 ARIMA(1,0,2)(1,1,0)[12] with drift : -5815.504 ARIMA(1,0,2)(1,1,1)[12] : Inf ARIMA(1,0,2)(1,1,1)[12] with drift : Inf ARIMA(1,0,2)(2,1,0)[12] : -5883.839 ARIMA(1,0,2)(2,1,0)[12] with drift : -5881.961 ARIMA(1,0,3)(0,1,0)[12] : -5633.48 ARIMA(1,0,3)(0,1,0)[12] with drift : -5631.472 ARIMA(1,0,3)(0,1,1)[12] : Inf ARIMA(1,0,3)(0,1,1)[12] with drift : Inf ARIMA(1,0,3)(1,1,0)[12] : -5816.716 ARIMA(1,0,3)(1,1,0)[12] with drift : -5814.777 ARIMA(1,0,4)(0,1,0)[12] : -5632.07 ARIMA(1,0,4)(0,1,0)[12] with drift : -5630.052 ARIMA(2,0,0)(0,1,0)[12] : -5624.965 ARIMA(2,0,0)(0,1,0)[12] with drift : -5622.978 ARIMA(2,0,0)(0,1,1)[12] : -5976.243 ARIMA(2,0,0)(0,1,1)[12] with drift : Inf ARIMA(2,0,0)(0,1,2)[12] : Inf ARIMA(2,0,0)(0,1,2)[12] with drift : Inf ARIMA(2,0,0)(1,1,0)[12] : -5813.072 ARIMA(2,0,0)(1,1,0)[12] with drift : -5811.194 ARIMA(2,0,0)(1,1,1)[12] : Inf ARIMA(2,0,0)(1,1,1)[12] with drift : Inf ARIMA(2,0,0)(1,1,2)[12] : -5972.219 ARIMA(2,0,0)(1,1,2)[12] with drift : Inf ARIMA(2,0,0)(2,1,0)[12] : -5879.922 ARIMA(2,0,0)(2,1,0)[12] with drift : -5878.135 ARIMA(2,0,0)(2,1,1)[12] : Inf ARIMA(2,0,0)(2,1,1)[12] with drift : Inf ARIMA(2,0,1)(0,1,0)[12] : -5628.116 ARIMA(2,0,1)(0,1,0)[12] with drift : -5626.11 ARIMA(2,0,1)(0,1,1)[12] : Inf ARIMA(2,0,1)(0,1,1)[12] with drift : Inf ARIMA(2,0,1)(0,1,2)[12] : Inf ARIMA(2,0,1)(0,1,2)[12] with drift : Inf ARIMA(2,0,1)(1,1,0)[12] : -5817.417 ARIMA(2,0,1)(1,1,0)[12] with drift : -5815.482 ARIMA(2,0,1)(1,1,1)[12] : Inf ARIMA(2,0,1)(1,1,1)[12] with drift : Inf ARIMA(2,0,1)(2,1,0)[12] : Inf ARIMA(2,0,1)(2,1,0)[12] with drift : Inf ARIMA(2,0,2)(0,1,0)[12] : Inf ARIMA(2,0,2)(0,1,0)[12] with drift : Inf ARIMA(2,0,2)(0,1,1)[12] : Inf ARIMA(2,0,2)(0,1,1)[12] with drift : Inf ARIMA(2,0,2)(1,1,0)[12] : -5816.598 ARIMA(2,0,2)(1,1,0)[12] with drift : -5814.659 ARIMA(2,0,3)(0,1,0)[12] : -5632.052 ARIMA(2,0,3)(0,1,0)[12] with drift : -5630.034 ARIMA(3,0,0)(0,1,0)[12] : -5635.615 ARIMA(3,0,0)(0,1,0)[12] with drift : -5633.61 ARIMA(3,0,0)(0,1,1)[12] : Inf ARIMA(3,0,0)(0,1,1)[12] with drift : Inf ARIMA(3,0,0)(0,1,2)[12] : Inf ARIMA(3,0,0)(0,1,2)[12] with drift : Inf ARIMA(3,0,0)(1,1,0)[12] : -5819.236 ARIMA(3,0,0)(1,1,0)[12] with drift : -5817.316 ARIMA(3,0,0)(1,1,1)[12] : Inf ARIMA(3,0,0)(1,1,1)[12] with drift : Inf ARIMA(3,0,0)(2,1,0)[12] : -5883.556 ARIMA(3,0,0)(2,1,0)[12] with drift : -5881.712 ARIMA(3,0,1)(0,1,0)[12] : -5633.852 ARIMA(3,0,1)(0,1,0)[12] with drift : -5631.841 ARIMA(3,0,1)(0,1,1)[12] : Inf ARIMA(3,0,1)(0,1,1)[12] with drift : Inf ARIMA(3,0,1)(1,1,0)[12] : -5817.198 ARIMA(3,0,1)(1,1,0)[12] with drift : -5815.272 ARIMA(3,0,2)(0,1,0)[12] : -5632.358 ARIMA(3,0,2)(0,1,0)[12] with drift : -5630.341 ARIMA(4,0,0)(0,1,0)[12] : -5633.71 ARIMA(4,0,0)(0,1,0)[12] with drift : -5631.7 ARIMA(4,0,0)(0,1,1)[12] : Inf ARIMA(4,0,0)(0,1,1)[12] with drift : Inf ARIMA(4,0,0)(1,1,0)[12] : -5817.205 ARIMA(4,0,0)(1,1,0)[12] with drift : -5815.274 ARIMA(4,0,1)(0,1,0)[12] : -5632.009 ARIMA(4,0,1)(0,1,0)[12] with drift : -5629.992 ARIMA(5,0,0)(0,1,0)[12] : -5632.617 ARIMA(5,0,0)(0,1,0)[12] with drift : -5630.598 Best model: ARIMA(0,0,3)(0,1,1)[12]
Series: it ARIMA(0,0,3)(0,1,1)[12] Box Cox transformation: lambda= -0.8397411 Coefficients: ma1 ma2 ma3 sma1 0.2070 0.1108 0.1212 -0.8823 s.e. 0.0387 0.0377 0.0356 0.0253 sigma^2 estimated as 8.5e-06: log likelihood=2994.15 AIC=-5978.3 AICc=-5978.21 BIC=-5955.7
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)
Lilliefors (Kolmogorov-Smirnov) normality test data: ajuste_auto_bc_outliers$residuals D = 0.045034, p-value = 0.002032
Image in a Jupyter notebookImage in a Jupyter notebook

Testando o modelo

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

Previsões para os próximos 3 meses

previsoes_auto_bc_outliers <- forecast(ajuste_auto_bc_outliers, h=3)
previsoes_auto_bc_outliers
accuracy(previsoes_auto_bc_outliers)
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 Aug 2018 26.33039 24.87485 27.95070 24.16290 28.88467 Sep 2018 28.98351 27.21961 30.96980 26.36381 32.12565 Oct 2018 30.73892 28.76825 32.97308 27.81667 34.28051
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)

jpeg("g15.jpg", width = 1000, height = 700)
plot(bin)
dev.off()
plot(bin)

options(repr.plot.width=16, repr.plot.height=8)
jpeg("g16.jpg", width = 1000, height = 700)
ts.plot(ajuste_auto_bc_outliers$fitted,ajuste_auto_bc_outliers$x, lty=c(2,1), gpars = list(col = c("blue", "black")))
dev.off()
ts.plot(ajuste_auto_bc_outliers$fitted,ajuste_auto_bc_outliers$x, lty=c(2,1), gpars = list(col = c("blue", "black")))
Call: lmRob(formula = ajuste_auto_bc_outliers$x ~ ajuste_auto_bc_outliers$fitted, data = T) Residuals: Min 1Q Median 3Q Max -3.915993 -0.759128 0.005429 0.725391 4.303154 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.59538 0.57686 1.032 0.302 ajuste_auto_bc_outliers$fitted 0.98361 0.02091 47.032 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.133 on 689 degrees of freedom Multiple R-Squared: 0.5979 Test for Bias: statistic p-value M-estimate 0.1828 9.127e-01 LS-estimate 39.6429 2.464e-09
Image in a Jupyter notebookImage in a Jupyter notebook
jpeg("g17.jpg", width = 1000, height = 700)
plotarimapred(it, ajuste_auto_bc_outliers, xlim=c(2010, 2019))
dev.off()
plotarimapred(it, ajuste_auto_bc_outliers, xlim=c(2010, 2019))
Image in a Jupyter notebook

Bootstrap

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
STATIONARY BOOTSTRAP FOR TIME SERIES Average Block Length of 20 Call: tsboot(tseries = it, statistic = AIC_Boot, R = 1000, l = 20, sim = "geom") Bootstrap Statistics : original bias std. error t1* -5971.776 673.5559 48.29323
boot.ci(boot_ajuste_auto_bc_outliers)
Warning message in boot.ci(boot_ajuste_auto_bc_outliers): “bootstrap variances needed for studentized intervals”Warning message in boot.ci(boot_ajuste_auto_bc_outliers): “BCa intervals not defined for time series bootstraps”
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 1000 bootstrap replicates CALL : boot.ci(boot.out = boot_ajuste_auto_bc_outliers) Intervals : Level Normal Basic Percentile 95% (-6740, -6551 ) (-6731, -6548 ) (-5395, -5213 ) Calculations and Intervals on Original Scale