Os dados a seguir correspondem à produçãoo total mensal de electricidade na Alemanha entre 1955 e 1979, em milhões de quilowatts-hora.
electric = ts(c(6410,5861,6471,5779,5815,5649,5844,6108,6352,6975,7124,7393,
7374,7209,7114,6688,6465,6406,6524,6649,6751,7604,7706,7776,
8187,7261,7542,6993,7364,6521,7069,7187,7542,8318,8316,8584,
8737,7623,8205,7551,7226,6979,7367,7329,7720,8449,8370,8655,
8660,7922,8154,7964,7501,7691,7896,8433,8888,9809,9832,10173,
10079,9637,9971,8902,9061,8387,8966,9083,9462,10158,10270,10916,
11141,9783,10560,9420,9574,9235,9352,9474,9951,11130,11365,11478,
11658,10673,11685,10294,10582,9784,10288,10437,10819,12203,12410,12975,
13620,12169,12973,11454,11307,10208,11024,11189,11588,13263,12970,13508,
14413,13374,13694,12991,11789,11895,12909,12714,13618,14820,14667,15077,
15157,14141,15158,13602,13488,12510,12772,12798,13419,14833,15764,16041,
16111,14170,15699,13959,13412,12907,13094,13098,14347,15761,16404,16529,
16521,14557,15288,14757,13795,13833,13336,13524,14736,16720,17396,17776,
18069,16657,17339,15360,15675,13967,15371,15779,16400,18745,19677,19878,
20238,18581,19939,17956,16965,16382,17020,16887,18133,20497,20917,22533,
22490,20284,21433,20559,18375,17731,18030,17387,18937,21896,22052,23124,
23907,21992,24282,20835,19492,19405,19297,18995,20504,22659,23932,24330,
25468,23425,23896,21849,21457,20054,19519,20451,22029,24581,25451,26595,
27907,25190,26571,24532,23406,21240,21416,22098,23132,27008,27974,28526,
28995,26005,27825,24563,25163,22552,23397,22985,24668,28974,28473,28054,
28656,26084,27344,25755,22487,21826,20867,20802,22455,27429,28378,29719,
30462,29363,30630,26602,25574,24509,24401,24642,26682,28241,30195,32351,
32423,28242,30005,27903,25843,24277,23349,23929,26475,29055,30948,32870,
33266,30695,30594,28529,26841,25652,24849,25419,28088,31290,33259,34944,
38139,32933,34351,30217,29565,26236,27078,27633,28220,31829,33604,32945),
start = c(1955,1), frequency = 12)
Fonte: A First Course on Time Series Analysis Chair of Statistics, University of Wurzburg, September 18, 2006
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0))
plot(electric, main = "Produção total mensal de electricidade na Alemanha", ylab = "milhões de Kw/h",
xlab = "", type = "l")
grid()
Observemos a média e a função de autocorrelação.
library(astsa)
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0))
electric1 = stl(electric, s.window = "periodic")
plot(electric1, main = "Produção total mensal de electricidade na Alemanha")
grid()
acf1(electric, col = 2:7, lwd=4, gg=TRUE)
## [1] 0.98 0.95 0.91 0.88 0.86 0.85 0.85 0.85 0.87 0.88 0.89 0.88 0.86 0.84 0.80
## [16] 0.77 0.75 0.74 0.74 0.75 0.76 0.77 0.78 0.78 0.76 0.73 0.70 0.67 0.65 0.64
## [31] 0.64 0.65 0.66 0.67 0.67 0.67 0.65 0.63 0.60 0.57 0.55 0.54 0.54 0.54 0.55
## [46] 0.56 0.56 0.56
Podemos utilizar a tendência estimada acima para remover o crescimento desta série ou podemos diferenciar a série. Apresentamos o resultado de aplicar as duas possibilidades a continuação.
diff1 = diff(electric)
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
plot(diff1, main = expression(paste("Primeira diferença: ", (1-B)*X[t], " = ", X[t]-X[t-1])),
ylab = "", xlab = "", type = "l")
abline(h=mean(diff1), col = "red", lwd = 4)
grid()
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
electric2 = electric - electric1$time.series[,2]
plot(electric2, main = "Produção total mensal de electricidade na Alemanha \n eliminando a tendência",
ylab = "", xlab = "", type = "l")
abline(h=mean(electric2), col = "red", lwd = 4)
grid()
Observamos que estas séries com média zero aparentam forte sazonalidade anual, ainda aparenta crescimento da variabilidade, isso faz pensar em modelos multiplicativos.
log.electric = log(electric)
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
plot(log.electric, main = "Produção total mensal de electricidade na Alemanha \n escala logaritmica",
ylab = "log(milhoes de Kw/h)", xlab = "", type = "l")
grid()
diff2 = diff(log.electric)
plot(diff2, main = expression(paste("Primeira diferença: ", (1-B)*Y[t], " = ", Y[t]-Y[t-1], ", ",Y[t], "= ",log(X[t]))),
ylab = "", xlab = "", type = "l")
grid()
par(mfrow=c(1,1), mar=c(3,3,1,2)+.5, mgp=c(1.6,.6,0))
electric2 = stl(log.electric, s.window = "periodic")
plot(electric2, main = "Produção total mensal de electricidade na Alemanha \n escala logaritmica")
grid()
acf1(log.electric, col = 2:7, lwd=4, gg=TRUE)
## [1] 0.98 0.96 0.94 0.91 0.90 0.88 0.88 0.88 0.88 0.89 0.89 0.88 0.87 0.85 0.83
## [16] 0.80 0.79 0.77 0.77 0.77 0.77 0.78 0.78 0.77 0.76 0.74 0.72 0.70 0.68 0.67
## [31] 0.66 0.66 0.66 0.67 0.67 0.66 0.65 0.63 0.61 0.59 0.57 0.56 0.55 0.55 0.55
## [46] 0.55 0.55 0.55
Vamos utilizar agora o estimador não paramétrico da tendência para construirmos uma série com média que não dependa do tempo, isso para satisfazer a primeira exigência da definição de série estacionária.
log.electric1 = log.electric - electric2$time.series[,2]
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
plot(log.electric1, main = "Produção total mensal de electricidade na Alemanha \n escala logaritmica sem tendência")
abline(h=mean(log.electric1), col = "red", lwd = 4)
grid()
acf1(log.electric1, col = 2:7, lwd=4, gg=TRUE)
## [1] 0.71 0.44 -0.02 -0.45 -0.74 -0.84 -0.72 -0.44 -0.01 0.42 0.69 0.90
## [13] 0.68 0.42 -0.01 -0.43 -0.68 -0.80 -0.69 -0.43 -0.02 0.39 0.65 0.86
## [25] 0.66 0.42 0.00 -0.40 -0.66 -0.76 -0.66 -0.41 -0.03 0.37 0.61 0.82
## [37] 0.63 0.40 0.01 -0.36 -0.62 -0.72 -0.63 -0.40 -0.04 0.33 0.57 0.78
Vejamos as primeiras 8 autocorrelações na nova série estabilizada utllizando como tendência um modelo não paramétrico.
lag1.plot(log.electric1,9)
Percebemos a concordância do gráfico de autocorrelações acima com o gráfico da função de autocorrelação (ACF).
Um problema quanto à estimação não paramétrica da tendência e na estimação não paramétrica em geral é que não dispomos de um modelo para reconstruirmos os dados. Em algumas situações vamos utilizar somente estimação não paramétrica mas, na maioria das situações preferimos estimar paramétricamente a tendência. Faremos isso a continuação utilizando o logaritmo da série original.
tempo = time(log.electric)
modelo1 = lm(log.electric ~ tempo)
modelo2 = lm(log.electric ~ tempo + tempo^2)
modelo3 = lm(log.electric ~ tempo + tempo^2 + tempo^3)
AIC(modelo1, modelo2, modelo3)
## df AIC
## modelo1 3 -478.0799
## modelo2 3 -478.0799
## modelo3 3 -478.0799
BIC(modelo1, modelo2, modelo3)
## df BIC
## modelo1 3 -466.9686
## modelo2 3 -466.9686
## modelo3 3 -466.9686
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
plot(log.electric, main = "Produção total mensal de electricidade na Alemanha \n escala logaritmica",
ylab = "", xlab = "")
lines(ts(fitted(modelo1),start = c(1955,1), frequency = 12), col = "red", lwd = 4)
lines(ts(fitted(modelo2),start = c(1955,1), frequency = 12), col = "blue", lwd = 4)
lines(ts(fitted(modelo3),start = c(1955,1), frequency = 12), col = "green", lwd = 4)
grid()
Estes resultados indicam, a princípio, que o modelo linear simples seria suficiente para identifcar a tendência de crescimento linear. Observemos o comportamento da série corrigida por este modelo.
log.electric2 = log.electric - ts(fitted(modelo1),start = c(1955,1), frequency = 12)
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
plot(log.electric2, main = "Produção total mensal de electricidade na Alemanha \n escala logaritmica sem tendência")
abline(h=mean(log.electric2), col = "red", lwd = 4)
grid()
Percebemos que este séria ainda não é estável, queremos dizer que percebemos crescimento da produção até o ano 1963, depois estabilidade, decrescimento até 1968, aumento e estabilidade até 1974 e depois novamento decrescimento da produção.
Utilizaremos regressão segmentada para estimarmos este comportamento médio.
library(segmented)
seg.electric = segmented(modelo1, seg.Z = ~ tempo, psi = c(1963,1968,1975))
summary(seg.electric)
##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = modelo1, seg.Z = ~tempo, psi = c(1963, 1968,
## 1975))
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.tempo 1965.167 0.290
## psi2.tempo 1965.395 0.435
## psi3.tempo 1972.935 0.670
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.533e+02 5.694e+00 -26.917 <2e-16 ***
## tempo 8.285e-02 2.905e-03 28.519 <2e-16 ***
## U1.tempo -3.026e-01 7.990e-01 -0.379 NA
## U2.tempo 2.978e-01 7.990e-01 0.373 NA
## U3.tempo -4.256e-02 6.796e-03 -6.262 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09417 on 292 degrees of freedom
## Multiple R-Squared: 0.9655, Adjusted R-squared: 0.9647
##
## Boot restarting based on 6 samples. Last fit:
## Convergence attained in 7 iterations (rel. change -5.7829e-08)
slope(seg.electric)
## $tempo
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.082847 0.0029049 28.51900 0.077130 0.088564
## slope2 -0.219790 0.7990200 -0.27508 -1.792400 1.352800
## slope3 0.078055 0.0045095 17.30900 0.069180 0.086931
## slope4 0.035498 0.0050848 6.98130 0.025491 0.045506
intercept(seg.electric)
## $tempo
## Est.
## intercept1 -153.260
## intercept2 441.480
## intercept3 -143.910
## intercept4 -59.948
plot.segmented(seg.electric, lwd=4, xlab = "", col = "red")
lines(log.electric)
grid()
Percebemos uma boa representação do comportamento da tendência descrito. Esperamos então que os dados sem o comportamento da média dependendo do tempo sejam estáveis.
log.electric3 = log.electric - ts(fitted(seg.electric),start = c(1955,1), frequency = 12)
par(mfrow=c(1,1), mar=c(3,3,2,2)+.5, mgp=c(1.6,.6,0))
plot(log.electric3, main = "Produção total mensal de electricidade na Alemanha \n escala logaritmica sem tendência", xlab = "", ylab = "")
abline(h=mean(log.electric3), col = "red", lwd = 4)
grid()
acf1(log.electric3, col = 2:7, lwd=4, gg=TRUE)
## [1] 0.73 0.47 0.04 -0.37 -0.65 -0.76 -0.66 -0.41 -0.02 0.38 0.63 0.83
## [13] 0.63 0.38 -0.01 -0.41 -0.65 -0.76 -0.66 -0.42 -0.05 0.34 0.59 0.79
## [25] 0.60 0.38 0.00 -0.37 -0.61 -0.72 -0.63 -0.40 -0.05 0.33 0.55 0.76
## [37] 0.58 0.38 0.02 -0.34 -0.57 -0.68 -0.60 -0.39 -0.05 0.30 0.53 0.73
lag1.plot(log.electric3, 9)
library(forecast)
auto.arima(log.electric3)
## Series: log.electric3
## ARIMA(2,0,2)(1,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1
## 1.3224 -0.4753 -0.9089 0.3715 -0.0600 -0.7988
## s.e. 0.2773 0.2356 0.2736 0.1228 0.0701 0.0439
##
## sigma^2 = 0.0008322: log likelihood = 608.52
## AIC=-1203.04 AICc=-1202.64 BIC=-1177.4
Uma forma de verificarmos a significância dos coeficientes sugeridos pelo modelo de escolha automática é ajustar novamente o mesmo modelo com o comando sarima.
modelo.sarima = sarima(log.electric3, p = 2, d = 0, q = 2, P = 1, D = 1, Q = 1, S = 12,
no.constant = TRUE)
## initial value -3.027103
## iter 2 value -3.280457
## iter 3 value -3.449984
## iter 4 value -3.507181
## iter 5 value -3.521695
## iter 6 value -3.529273
## iter 7 value -3.534104
## iter 8 value -3.534657
## iter 9 value -3.534984
## iter 10 value -3.535226
## iter 11 value -3.535282
## iter 12 value -3.535376
## iter 13 value -3.535729
## iter 14 value -3.535985
## iter 15 value -3.536248
## iter 16 value -3.536503
## iter 17 value -3.536551
## iter 18 value -3.536625
## iter 19 value -3.536731
## iter 20 value -3.536982
## iter 21 value -3.537056
## iter 22 value -3.537134
## iter 23 value -3.537139
## iter 24 value -3.537140
## iter 24 value -3.537140
## iter 24 value -3.537140
## final value -3.537140
## converged
## initial value -3.523791
## iter 2 value -3.523928
## iter 3 value -3.524229
## iter 4 value -3.524264
## iter 5 value -3.524292
## iter 6 value -3.524323
## iter 7 value -3.524337
## iter 8 value -3.524349
## iter 9 value -3.524361
## iter 10 value -3.524412
## iter 11 value -3.524616
## iter 12 value -3.524732
## iter 13 value -3.524905
## iter 14 value -3.525492
## iter 15 value -3.526227
## iter 16 value -3.526861
## iter 17 value -3.526945
## iter 18 value -3.527004
## iter 19 value -3.527204
## iter 20 value -3.527500
## iter 21 value -3.527947
## iter 22 value -3.528563
## iter 23 value -3.529494
## iter 24 value -3.530166
## iter 25 value -3.530427
## iter 26 value -3.530757
## iter 27 value -3.531222
## iter 28 value -3.531707
## iter 29 value -3.531789
## iter 30 value -3.531801
## iter 31 value -3.531816
## iter 32 value -3.531823
## iter 33 value -3.531829
## iter 34 value -3.531832
## iter 35 value -3.531833
## iter 36 value -3.531835
## iter 37 value -3.531838
## iter 38 value -3.531843
## iter 39 value -3.531847
## iter 40 value -3.531848
## iter 40 value -3.531848
## final value -3.531848
## converged
modelo.sarima$ttable
## Estimate SE t.value p.value
## ar1 1.3224 0.2773 4.7690 0.0000
## ar2 -0.4753 0.2356 -2.0174 0.0446
## ma1 -0.9089 0.2736 -3.3216 0.0010
## ma2 0.3715 0.1228 3.0245 0.0027
## sar1 -0.0600 0.0701 -0.8558 0.3929
## sma1 -0.7988 0.0439 -18.2094 0.0000
modelo.sarima$AIC
## [1] -4.177208
Observamos que o coefciente autoregressivo sazonal não pe significativo então, uma maneira de desconsiderar o ajuste deste coeficiênte no modelo seria da seguinte forma.
modelo.sarima1 = sarima(log.electric3, p = 2, d = 0, q = 2, P = 1, D = 1, Q = 1, S = 12,
fixed = list(NA,NA,NA,NA,0,NA), no.constant = TRUE)
## initial value -3.027103
## iter 2 value -3.246809
## iter 3 value -3.474216
## iter 4 value -3.505585
## iter 5 value -3.516406
## iter 6 value -3.526796
## iter 7 value -3.529548
## iter 8 value -3.531613
## iter 9 value -3.532142
## iter 10 value -3.532336
## iter 11 value -3.532486
## iter 12 value -3.532889
## iter 13 value -3.533082
## iter 14 value -3.533172
## iter 15 value -3.533437
## iter 16 value -3.533689
## iter 17 value -3.533834
## iter 18 value -3.533964
## iter 19 value -3.533986
## iter 20 value -3.533997
## iter 21 value -3.533998
## iter 22 value -3.533998
## iter 22 value -3.533998
## iter 22 value -3.533998
## final value -3.533998
## converged
## initial value -3.522415
## iter 2 value -3.522582
## iter 3 value -3.522666
## iter 4 value -3.522693
## iter 5 value -3.522731
## iter 6 value -3.522750
## iter 7 value -3.522777
## iter 8 value -3.522828
## iter 9 value -3.523312
## iter 10 value -3.523696
## iter 11 value -3.524036
## iter 12 value -3.524160
## iter 13 value -3.524891
## iter 14 value -3.525485
## iter 15 value -3.525935
## iter 16 value -3.525936
## iter 17 value -3.526053
## iter 18 value -3.526115
## iter 19 value -3.526516
## iter 20 value -3.527316
## iter 21 value -3.528444
## iter 22 value -3.529964
## iter 23 value -3.530183
## iter 24 value -3.530246
## iter 25 value -3.530289
## iter 26 value -3.530302
## iter 27 value -3.530316
## iter 28 value -3.530323
## iter 29 value -3.530379
## iter 30 value -3.530459
## iter 31 value -3.530543
## iter 32 value -3.530587
## iter 33 value -3.530604
## iter 34 value -3.530604
## iter 34 value -3.530604
## final value -3.530604
## converged
modelo.sarima1$ttable
## Estimate SE t.value p.value
## ar1 1.3285 0.2719 4.8869 0.0000
## ar2 -0.4828 0.2297 -2.1020 0.0364
## ma1 -0.9146 0.2671 -3.4239 0.0007
## ma2 0.3784 0.1174 3.2238 0.0014
## sma1 -0.8160 0.0376 -21.6797 0.0000
modelo.sarima1$AIC
## [1] -4.181665
Percebemos uma pequena melhoria na qualidade do ajuste e ainda observamos que o valor do AIC o demonstra. Selecionado o modelo procedemos a obtermos predições.
par(pch=19)
prev.electric = sarima.for(log.electric3, n.ahead = 12, p = 2, d = 0, q = 2, P = 1, D = 1, Q = 1, S = 12,
fixed = list(NA,NA,NA,NA,0,NA), no.constant = TRUE)
grid()
Utilizamos diversas transformações para obtermos uma série fracamente estacionária à qual ajustamos um modelo SARIMA, do qual mostramos acima as previsões para o ano de 1980.
Faremos a suposição de que os dados, em 1980, devem comportar-se como uma continuação do modelo de regressão segmentada e assim mostrarmos as previsões na escala origial.
prev.electric.org = ts(c(electric,exp(prev.electric$pred)), start = c(1955,1), frequency = 12)
novo.tempo = time(prev.electric.org)
tendencia = c(intercept(seg.electric)$tempo[1]
+ slope(seg.electric)$tempo[1]*novo.tempo[novo.tempo < seg.electric$psi[1,2]],
intercept(seg.electric)$tempo[2]
+ slope(seg.electric)$tempo[2]*novo.tempo[seg.electric$psi[1,2] < novo.tempo
& novo.tempo < seg.electric$psi[2,2]],
intercept(seg.electric)$tempo[3]
+ slope(seg.electric)$tempo[3]*novo.tempo[seg.electric$psi[2,2] < novo.tempo
& novo.tempo < seg.electric$psi[3,2]],
intercept(seg.electric)$tempo[4]
+ slope(seg.electric)$tempo[4]*novo.tempo[seg.electric$psi[3,2] < novo.tempo])
tendencia.prev = ts(tendencia, start = c(1955,1), frequency = 12)
electric.prev = ts(c(electric,exp(prev.electric$pred+tendencia.prev[novo.tempo>=1980])), start = c(1955,1), frequency = 12)
par(mfrow=c(1,1), mar=c(3,3,3,2)+.5, mgp=c(1.6,.6,0))
plot(electric.prev,
main = "Produção total mensal de electricidade na Alemanha \n valores previstos para 1980",
ylab = "milhões de Kw/h", xlab = "", type = "l")
lines(exp(tendencia.prev), col ="red")
U = exp(prev.electric$pred+prev.electric$se+tendencia.prev[novo.tempo>=1980])
L = exp(prev.electric$pred-prev.electric$se+tendencia.prev[novo.tempo>=1980])
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
grid()