Descrição

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

Estudo descritivo 

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.


Logaritmizando a séries original 

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).


Estimação paramétrica

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)

Modelo 

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()


Previsão nos dados originais

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()