A série temporal plotada na abaixo consiste em 57 overshorts diários consecutivos de um tanque subterrâneo de gasolina em um posto de gasolina no Colorado. Se \(y_t\) é a quantidade medida de combustível no tanque no final do dia \(t\) e \(a_t\) é a quantidade medida vendida menos a quantidade entregue durante o curso do dia \(t\), então o overshort no final do dia \(t\) é definido como \(x_t = y_t - y_{t−1} +a_t\).
Devido ao erro na medição da quantidade atual de combustível no tanque, a quantidade vendida e a quantidade entregue ao posto, vemos \(Y_t\), \(a_t\) e \(X_t\) como valores observados de algum conjunto de variáveis aleatórias \(Y_t\), \(A_t\) e \(X_t\) para \(t = 1,\cdots,57\). Na ausência de qualquer erro de medição e qualquer vazamento no tanque, cada \(x_t\) seria zero.
oshorts = read.table(file = "http://leg.ufpr.br/~lucambio/CE017/20212S/oshorts.tsm")
head(oshorts)
## V1
## 1 78
## 2 -58
## 3 53
## 4 -65
## 5 13
## 6 -6
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0), pch=19)
X = ts(oshorts$V1)
plot(X, ylab="Galões", xlab="Dias", col="violet")
points(X)
abline(h=0, lty=1, lwd=2)
grid()
Como sempre devemos ter uma ideia de qual modelo propôr. Vejamos as funções autocorrelação e autocorrelação parcial o que nos informam, assim como outras correlações.
library(astsa)
lag1.plot(X, 4)
acf2(X)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## ACF -0.5 0.12 -0.21 0.08 0.02 0.12 -0.22 0.25
## PACF -0.5 -0.18 -0.32 -0.26 -0.15 0.04 -0.19 0.12
Parece plausível um modelo \(AR(1)\).
modelo01 = arima(X, order=c(1,0,0))
modelo01
##
## Call:
## arima(x = X, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## -0.5129 -4.4509
## s.e. 0.1141 4.4196
##
## sigma^2 estimated as 2518: log likelihood = -304.22, aic = 614.44
modelo02 = arima(X, order=c(2,0,0))
modelo02
##
## Call:
## arima(x = X, order = c(2, 0, 0))
##
## Coefficients:
## ar1 ar2 intercept
## -0.6026 -0.1779 -4.6034
## s.e. 0.1308 0.1322 3.7097
##
## sigma^2 estimated as 2438: log likelihood = -303.33, aic = 614.66
Aceitamos como modelo adequado \(AR(1)\) (menor AIC). Podemos então proceder a predizer overshorts futuros.
# horizonte de provisão
m = 3
preditos = predict(modelo01, n.ahead = m)
preditos
## $pred
## Time Series:
## Start = 58
## End = 60
## Frequency = 1
## [1] 1.986064 -7.752596 -2.757274
##
## $se
## Time Series:
## Start = 58
## End = 60
## Frequency = 1
## [1] 50.17688 56.39276 57.91744
par(mfrow = c(1,1), mar=c(4,3,1,1), mgp=c(1.6,.6,0), pch=19)
ts.plot(X, preditos$pred, col=1:2, lwd=2, ylab="Overshorts", xlab="Dias")
U = preditos$pred+preditos$se; U; L = preditos$pred-preditos$se; L
## Time Series:
## Start = 58
## End = 60
## Frequency = 1
## [1] 52.16295 48.64016 55.16017
## Time Series:
## Start = 58
## End = 60
## Frequency = 1
## [1] -48.19082 -64.14536 -60.67472
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
lines(preditos$pred, type="p", col=2)
grid()