Leitura dos dados

Este conjunto de dados descreve o número mensal de vendas de xampu em um período de 3 anos. As unidades são uma contagem de vendas e há 36 observações. O conjunto de dados original é creditado a Makridakis, Wheelwright e Hyndman (1998). 

shampoo = read.table(file = "http://leg.ufpr.br/~lucambio/CE017/20212S/shampoo.txt", sep=",", header = TRUE)
head(shampoo)
##   Month Sales
## 1  1-01 266.0
## 2  1-02 145.9
## 3  1-03  83.1
## 4  1-04 119.3
## 5  1-05 180.3
## 6  1-06 168.5
vendas = ts(shampoo$Sales, start = c(2001,1), frequency = 12)
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0), pch=19)
plot(vendas, ylab="Número mensal de vendas de xampu", xlab="Tempo")
grid()

Observação

O conjunto de dados mostra uma tendência crescente e possivelmente algum componente sazonal. Pensamos em tendência não linear, possivelmente uma função exponencial.

Pensamos sempre como refrência alguma curva matemática conhecida: logaritmo, inversa, exponencial, linear, polinômio, etc. Visualizemos algumas destas curvas a continuação. 

par(mfrow=c(2,2), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0), pch=19)
x = seq(0.01,5,by=0.01)
plot(x,log(x), type="l", col="red",xlab="x",ylab = "logaritmo")
grid()
x = seq(0.01,5,by=0.01)
plot(x,exp(x), type="l", col="red",xlab="x",ylab = "exponencial")
grid()
x = seq(0.01,5,by=0.01)
plot(x,1/x, type="l", col="red",xlab="x",ylab = "inversa")
grid()
x = seq(0.01,5,by=0.01)
plot(x,x+x^2, type="l", col="red",xlab="x",ylab = "polinômio")
grid()

Eliminando a tendência

Transformando os dados originais podemos obter como resposta uma série sem tendência e mesmo um série estacionária. Faremos a continuação a eleiminação da tendência visualizada e verificamos a estaiocionariedade da série obtida. 

trend = lm(vendas ~ time(vendas)+I(time(vendas)^2))
(summary(trend))
## 
## Call:
## lm(formula = vendas ~ time(vendas) + I(time(vendas)^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -103.816  -44.693   -7.286   33.564  145.622 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.571e+08  6.515e+07   3.947 0.000391 ***
## time(vendas)      -2.570e+05  6.507e+04  -3.949 0.000389 ***
## I(time(vendas)^2)  6.420e+01  1.625e+01   3.951 0.000386 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.27 on 33 degrees of freedom
## Multiple R-squared:  0.8269, Adjusted R-squared:  0.8164 
## F-statistic:  78.8 on 2 and 33 DF,  p-value: 2.713e-13
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0), pch=19)
plot(vendas, ylab="Número mensal de vendas de xampu \n tendência estimada", xlab="Tempo")
lines(ts(fitted(trend),start = c(2001,1), frequency = 12), col="red", lwd=2)
grid()

Dados transformados

vendas1 = vendas - fitted(trend)
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0), pch=19)
plot(vendas1, ylab="Número mensal de vendas de xampu", xlab="Tempo", col="violet", type="b")
abline(h=0, lwd=2)
grid()

library(astsa)
acf2(vendas1)

##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
## ACF  -0.45  0.16 -0.19  0.03  0.03 -0.11  0.08 -0.23
## PACF -0.45 -0.05 -0.17 -0.16 -0.03 -0.17 -0.08 -0.30

Modelos

Devemos ter alguma ideia de qual modelo propôr. Uma forma de fazer isso é verificando as possíveis influências das observações atrás na resposta no momento atual. 

lag1.plot(vendas1, 6)

Percebemos que este modelo apresenta-se como autoregressivo de ordem máximo 2. 

mod1 =  arima(vendas1, order=c(1,0,0))
mod1
## 
## Call:
## arima(x = vendas1, order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.4653    -1.0452
## s.e.   0.1492     6.3574
## 
## sigma^2 estimated as 3065:  log likelihood = -195.7,  aic = 397.41
mod2 =  arima(vendas1, order=c(2,0,0))
mod2
## 
## Call:
## arima(x = vendas1, order = c(2, 0, 0))
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.4939  -0.0605    -1.1361
## s.e.   0.1687   0.1670     5.9993
## 
## sigma^2 estimated as 3053:  log likelihood = -195.64,  aic = 399.28

Podemos investir em mais parâmetros. 

mod3 =  arima(vendas1, order=c(3,0,0))
mod3
## 
## Call:
## arima(x = vendas1, order = c(3, 0, 0))
## 
## Coefficients:
##           ar1      ar2      ar3  intercept
##       -0.4893  -0.1313  -0.1475    -0.8507
## s.e.   0.1674   0.1865   0.1791     5.2756
## 
## sigma^2 estimated as 2993:  log likelihood = -195.3,  aic = 400.61

Caso seja aceito o modelo no objeto mod1 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 = 3
preditos = predict(mod1, n.ahead = m)
par(mfrow = c(1,1), mar=c(4,3,1,1), mgp=c(1.6,.6,0), pch=19)
ts.plot(vendas1, preditos$pred, col=1:2, xlim=c(2001,2004.2), lwd=2, ylab="Vendas de xampú", xlab="Tempo")
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(vendas)
novas.vendas = ts(c(shampoo$Sales,rep(0,m)), start = c(2001,1), frequency = 12)
vendas2 = coef(trend)[1]+coef(trend)[2]*time(novas.vendas)+coef(trend)[3]*I(time(novas.vendas)^2)
novos.preditos = vendas2[(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(vendas, novos.preditos, col=1:2, xlim=c(2001,2004.2), lwd=2, ylab="Vendas de xampu", xlab="Tempo")
U = novos.preditos+preditos$se; L = novos.preditos-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(novos.preditos, type="p", col=2)
grid()