beer
Fernando Lucambio Pérez
20/03/2022
Dados
Produção mensal de cerveja na Austrália em milhões de litros: Janeiro de 1956 - Agosto de 1995. Queremos estimar a produção um ano à frente.
beer = ts(read.table('http://leg.ufpr.br/~lucambio/STemporais/beer.txt', header = TRUE), start = c(1956,1), frequency = 12)
#
min(beer); max(beer); length(beer)## [1] 64.8## [1] 217.8## [1] 476Apresenção dos dados
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0))
plot(beer, main = "Produção mensal de cerveja na Austrália", ylab = "Milhões de litros", xlab = "", ylim = c(60, 220))
grid()
Observamos estabilidade até meados da década de 60, logo depois começa o crescimento até mais ou menos metade da década eos anos 70, estabiliza-se a produção até o ano 1990 e depois começa um leve decrescimento. Claro que nos referimos à tendência, dentro de cada ano acontece sazonalidade na produção.
Observemos com detalhes esses intervalo de tempo.
par(mfrow=c(2,2), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0))
plot(beer, main = "Produção mensal entre 1956 a 1965", ylab = "Milhões de litros", xlab = "", xlim = c(1956,1965), ylim = c(60, 220))
grid()
plot(beer, main = "Produção mensal entre 1965 a 1975", ylab = "Milhões de litros", xlab = "", xlim = c(1965,1975), ylim = c(60, 220))
grid()
#
plot(beer, main = "Produção mensal ente 1975 a 1990", ylab = "Milhões de litros", xlab = "", xlim = c(1975,1990), ylim = c(60, 220))
grid()
plot(beer, main = "Produção mensal entre 1990 a 1995", ylab = "Milhões de litros", xlab = "", xlim = c(1990,1995), ylim = c(60, 220))
grid()
Nossa última observaão não se mostrou realista, ou seja, os períodos 1975 a 1990 e 1990 a 1995 são, a princípio semelhantes. Pensamos então utilizar uma regressão segmentada para estimar os instantes de mudanças de comportamento e o próprio comportamento de tendência.
library(segmented)
tempo = time(beer)
ajuste = lm(beer ~ tempo, na.action=NULL)
ajuste1 = segmented(ajuste, seg.Z = ~ tempo, psi = c(1965,1975))
summary(ajuste1)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = ajuste, seg.Z = ~tempo, psi = c(1965, 1975))
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.tempo 1963.395 1.500
## psi2.tempo 1976.833 0.639
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3530.0584 1835.9200 -1.923 0.0551 .
## tempo 1.8473 0.9369 1.972 0.0492 *
## U1.tempo 3.3349 1.0115 3.297 NA
## U2.tempo -6.0945 0.4470 -13.633 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.92 on 470 degrees of freedom
## Multiple R-Squared: 0.6888, Adjusted R-squared: 0.6855
##
## Boot restarting based on 10 samples. Last fit:
## Convergence attained in 5 iterations (rel. change -1.2754e-09)par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0))
plot.segmented(ajuste1, lwd=2, main = "Produção mensal de cerveja na Austrália", ylab = "Milhões de litros", xlab = "", ylim = c(60, 220))
points(tempo, beer, pch=19, col="red", type="l")
grid()
Eliminemos o comportanmento de tendência com o intuito de transformar os dados numa série temporal estacionária (fracamente).
beer.sim = beer - fitted(ajuste1)
min(beer.sim); max(beer.sim)## [1] -50.19011## [1] 58.81409par(mfrow=c(1,1), mar=c(1,3,4,2)+.5, mgp=c(1.6,.6,0))
plot(beer.sim, main = "Produção mensal de cerveja na Austrália \n dados transformados", ylab = "Milhões de litros (transformados)", xlab = "", ylim = c(-60, 60))
grid()
library(astsa)
acf2(beer.sim)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## ACF 0.46 0.16 -0.01 -0.28 -0.38 -0.53 -0.44 -0.21 0.00 0.16 0.53 0.77
## PACF 0.46 -0.06 -0.07 -0.30 -0.17 -0.39 -0.16 -0.10 -0.03 -0.11 0.40 0.57
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF 0.42 0.22 -0.05 -0.29 -0.35 -0.56 -0.41 -0.21 -0.07 0.16 0.50 0.70
## PACF 0.02 0.14 0.00 -0.06 0.03 -0.08 0.04 -0.10 -0.14 -0.20 -0.08 0.17
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF 0.43 0.21 -0.09 -0.26 -0.37 -0.61 -0.39 -0.20 -0.11 0.20 0.46 0.65
## PACF -0.03 0.02 -0.03 -0.02 -0.01 -0.18 0.04 0.04 -0.11 0.03 -0.04 0.09
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF 0.47 0.14 -0.11 -0.22 -0.40 -0.60 -0.34 -0.26 -0.09 0.23 0.39 0.65
## PACF 0.06 -0.13 -0.01 0.06 0.03 -0.08 0.12 -0.08 -0.11 0.06 -0.07 0.02Modelos
Pensamos em modelos autoregressivos de ordem até 12.
modelo01 = ar.ols(beer.sim, order.max = 12)
modelo01##
## Call:
## ar.ols(x = beer.sim, order.max = 12)
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.0386 -0.0290 0.0284 -0.1483 0.1285 -0.1352 -0.0949 0.0066
## 9 10 11 12
## 0.0198 -0.1265 0.2402 0.6102
##
## Intercept: -0.04467 (0.4779)
##
## Order selected 12 sigma^2 estimated as 105.9modelo02 = ar.yw(beer.sim, order.max = 12)
modelo02##
## Call:
## ar.yw.default(x = beer.sim, order.max = 12)
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.0681 -0.0356 0.0292 -0.1485 0.1235 -0.1521 -0.0892 0.0034
## 9 10 11 12
## 0.0270 -0.1225 0.2320 0.5680
##
## Order selected 12 sigma^2 estimated as 115.9modelo03 = ar.mle(as.numeric(beer.sim), order.max = 12)
modelo03##
## Call:
## ar.mle(x = as.numeric(beer.sim), order.max = 12)
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.0412 -0.0267 0.0280 -0.1481 0.1256 -0.1379 -0.0950 0.0062
## 9 10 11 12
## 0.0199 -0.1222 0.2374 0.5989
##
## Order selected 12 sigma^2 estimated as 104Vamos então procurar as predições um ano à frente.
# horizonte de provisão
m = 12
preditos = predict(modelo01, n.ahead = m)
par(mfrow = c(1,1),mar=c(4,3,1,1),mgp=c(1.6,.6,0), pch=19)
ts.plot(beer.sim, preditos$pred, col=1:2, xlim=c(1992,1996), lwd=2, ylab="Produção mensal de cerveja", xlab="Tempo", 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()
Transformando os valores utilizados para modelar beer.sim e os valores preditos na escala original.
slope(ajuste1)## $tempo
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 1.8473 0.93685 1.9718 0.0063792 3.68830
## slope2 5.1823 0.38147 13.5850 4.4327000 5.93190
## slope3 -0.9122 0.23306 -3.9141 -1.3702000 -0.45424Acima observamos os diferentes coeficientes estimados pela regressão segmentada. Fizemos a supossição de dois pontos de corte: 1965 e 1975 (chutes iniciais). A regresão segmentada estimou como: 1963.395 e 1976.833.
Esses coeficientes permitirão a predição da reta de regressão segmentada para novos valores; faremos isso considerando dois cenários: estabilidade e crescimento de 10% das vendas.
valores.trans = beer.sim + fitted(ajuste1)
n = length(valores.trans)
# cenário de estabilidade
novos.preditos = list(pred = rep(fitted(ajuste1)[n],m)+preditos$pred, se = preditos$se)
par(mfrow = c(1,1),mar=c(4,3,1,1),mgp=c(1.6,.6,0), pch=19)
ts.plot(beer, novos.preditos$pred, col=1:2, xlim=c(1992,1997), lwd=2, ylab="Produção mensal de cerveja", xlab = "", type="b", ylim = c(100, 220))
mtext(side=3,"Cenário de estabilidade")
U = novos.preditos$pred+novos.preditos$se; L = novos.preditos$pred-novos.preditos$se
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
grid()
# cenário de crescimento
# Considerando crescimento de de 10% ao ano
beta = fitted(ajuste1)[n]/time(fitted(ajuste1))[n]
novos.preditos = list(pred = (1.1*beta*time(preditos$pred))+preditos$pred, se = preditos$se)
par(mfrow = c(1,1),mar=c(4,3,1,1),mgp=c(1.6,.6,0), pch=19)
ts.plot(beer, novos.preditos$pred, col=1:2, xlim=c(1992,1997), lwd=2, ylab="Produção mensal de cerveja", xlab = "", type="b", ylim = c(60, 220))
mtext(side=3,"Cenário de crescimento de 10")
U = novos.preditos$pred+novos.preditos$se; L = novos.preditos$pred-novos.preditos$se
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
grid()