Dados

A Figura abaixo mostra as vendas mensais (em quilolitros) de vinho tinto por vinicultores australianos de janeiro de 1980 a outubro de 1991. 

wine = read.table(file = "http://leg.ufpr.br/~lucambio/CE017/20212S/wine.tsm")
head(wine)
##     V1
## 1  464
## 2  675
## 3  703
## 4  887
## 5 1139
## 6 1077
wine = ts(wine, start = c(1980,1), frequency = 12)
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0), pch=19)
plot(wine, main="Vendas mensais de vinho tinto por vinicultores australianos", xlab="", ylab="Quilolitros", type="b")
grid()

Observações

Verifica-se no gráfico que as vendas têm uma tendência ascendente e um padrão sazonal com um pico em julho e um vale em janeiro. Vamos visualizar isso melhor.

Uma primeira trnasformação sugerida é aplicarmos o logaritmo, isto pode provocar uma diminuição da variabilidade dos dados e, como consequência, podemos melhor visualizar o possível padrão zasonal. 

par(mfrow=c(1,1), mar=c(1,3,4,2)+.5, mgp=c(1.6,.6,0), pch=19)
lwine = log(wine)
plot(lwine, main="Vendas mensais de vinho tinto por vinicultores australianos \n escala logarítmica", xlab="", ylab="log(Quilolitros)", type="b")
grid()

adj1 = lm(lwine ~ time(lwine))
par(mfrow=c(1,1), mar=c(1,3,4,2)+.5, mgp=c(1.6,.6,0), pch=19)
plot(lwine, main="Vendas mensais de vinho tinto por vinicultores australianos \n escala logarítmica", xlab="", ylab="log(Quilolitros)", type="b")
lines(ts(fitted(adj1), start = c(1980,1), frequency = 12), col="red", lwd=2)
grid()

Dados transformados

lwine1 = lwine - fitted(adj1)
par(mfrow=c(1,1), mar=c(1,3,4,2)+.5, mgp=c(1.6,.6,0), pch=19)
plot(lwine1, main="Vendas mensais de vinho tinto por vinicultores australianos \n escala logarítmica transformada", xlab="", ylab="log(Quilolitros) - tendência", type="b")
abline(h=0, col="red", lwd=2)
grid()

library(astsa)
monthplot(lwine1, xlab="Meses", ylab="Quilolitros", type="l")
grid()

Na figura acima percebemos que o comportamento resultante de cada ano, em cada mês, é mais devido à variabilidade devido à produção do que algum comportamento de tendência, ou seja, algum crescimento ou decrescimento. Por outro lado, percebe-se que entre janeiro e agosto existe crescimento das vendas, as quais diminuem fortemente em setembro e se estabilizam até dezembro. 

acf2(lwine1)

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12]
## ACF  0.47  0.20  0.05 -0.16 -0.40 -0.61 -0.43 -0.18 0.01  0.14  0.37  0.77
## PACF 0.47 -0.02 -0.05 -0.21 -0.31 -0.43 -0.05  0.08 0.11 -0.02  0.10  0.60
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.38  0.16  0.04 -0.16 -0.41 -0.59 -0.43 -0.22 -0.05  0.07  0.31  0.67
## PACF -0.29 -0.07  0.01 -0.03 -0.19  0.01 -0.02 -0.09 -0.07 -0.02 -0.02  0.11
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF   0.35  0.13  0.03 -0.13 -0.38 -0.53 -0.39 -0.20 -0.04  0.06  0.28  0.62
## PACF -0.08 -0.12 -0.05  0.06 -0.05  0.02  0.01 -0.05  0.00  0.00 -0.01  0.07
##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF   0.33  0.15  0.05 -0.11 -0.33 -0.46 -0.32 -0.18 -0.04  0.07  0.26  0.54
## PACF -0.07  0.01 -0.04 -0.07  0.05  0.08  0.02 -0.09 -0.01  0.04 -0.06 -0.05
lag1.plot(lwine1, 12)

Percebemos, tanto no gráfico das funções de autocorrelação e autocorrelação parcial quanto no gráfico das 12 defasagens acima. que existe uma fraca dependência de ordem 1 e uma forte dependência de ordem 12. Propomos então um modelo autoregressivo de ordem 12.

Modelos

mod01 = arima(lwine1, order = c(12,0,0))
mod01
## 
## Call:
## arima(x = lwine1, order = c(12, 0, 0))
## 
## Coefficients:
##          ar1     ar2     ar3      ar4      ar5      ar6      ar7     ar8
##       0.1151  0.0250  0.0685  -0.0371  -0.0255  -0.1384  -0.0117  0.0319
## s.e.  0.0583  0.0591  0.0588   0.0591   0.0591   0.0558   0.0591  0.0590
##          ar9     ar10     ar11    ar12  intercept
##       0.0423  -0.0120  -0.0266  0.7080    -0.0099
## s.e.  0.0584   0.0603   0.0602  0.0572     0.0360
## 
## sigma^2 estimated as 0.01734:  log likelihood = 80.63,  aic = -133.26

Caso seja aceito este modelo como adequado podemos, entre outras coisas, proceder a predizer valores de vendas futuros. Primeiro fazer ista previsão nos dados transformados, posteriormente apresentamos esta previsão nos dados originais. 

# horizonte de provisão
m = 12
preditos = predict(mod01, n.ahead = m)
par(mfrow = c(1,1), mar=c(2,3,4,1), mgp=c(1.6,.6,0), pch=19)
ts.plot(lwine1, preditos$pred, col=1:2, xlim=c(1986,1993), lwd=2, ylab="log(Quilolitros) - tendência", xlab="", main = "Vendas mensais de vinho tinto por vinicultores australianos \n escala logarítmica transformada", type="b")
U = preditos$pred+preditos$se; L = preditos$pred-preditos$se
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()

Procedemos agora a aplicar a transformação inversa nos dados, ou seja, acrescentrar a tendência estimada aos dados transformados e ainda predizemos os valores de tendência para os valores preditos serem transformados à escala origina. 

n = length(lwine)
novas.vendas = ts(c(wine,rep(0,m)), start = c(1980,1), frequency = 12)
lwine2 = coef(adj1)[1]+coef(adj1)[2]*time(novas.vendas)
novos.preditos = exp(lwine2[(n+1):(n+m)]+preditos$pred)
par(mfrow = c(1,1), mar=c(4,3,1,1), mgp=c(1.6,.6,0), pch=19)
ts.plot(wine, novos.preditos, col=1:2, xlim=c(1980,1993), lwd=2, ylab="Vendas de xampu", xlab="Tempo", type="b")
lines(novos.preditos, type="p", col=2)
grid()