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 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érie 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,1970,1975))
summary(seg.electric)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = modelo1, seg.Z = ~tempo, psi = c(1963, 1970, 
##     1975))
## 
## Estimated Break-Point(s):
##                 Est. St.Err
## psi1.tempo 1965.029  0.785
## psi2.tempo 1966.513  0.762
## psi3.tempo 1971.979  0.646
## 
## Meaningful coefficients of the linear terms:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -152.58351    5.76318 -26.476   <2e-16 ***
## tempo          0.08250    0.00294  28.058   <2e-16 ***
## U1.tempo      -0.06132    0.05141  -1.193       NA    
## U2.tempo       0.06453    0.05187   1.244       NA    
## U3.tempo      -0.04691    0.00855  -5.486       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09414 on 292 degrees of freedom
## Multiple R-Squared: 0.9655,  Adjusted R-squared: 0.9647 
## 
## Boot restarting based on 10 samples. Last fit:
## Convergence attained in 3 iterations (rel. change 4.9113e-06)
slope(seg.electric)
## $tempo
##            Est.   St.Err.  t value CI(95%).l CI(95%).u
## slope1 0.082502 0.0029404 28.05800  0.076715  0.088289
## slope2 0.021187 0.0513250  0.41279 -0.079828  0.122200
## slope3 0.085719 0.0074688 11.47700  0.071019  0.100420
## slope4 0.038812 0.0041609  9.32790  0.030623  0.047002
intercept(seg.electric)
## $tempo
##                Est.
## intercept1 -152.580
## intercept2  -32.097
## intercept3 -159.000
## intercept4  -66.502
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.38 -0.66 -0.77 -0.68 -0.42 -0.02  0.38  0.63  0.83
## [13]  0.63  0.38 -0.01 -0.41 -0.66 -0.77 -0.67 -0.42 -0.04  0.34  0.59  0.79
## [25]  0.61  0.39  0.00 -0.37 -0.62 -0.72 -0.63 -0.40 -0.05  0.33  0.56  0.76
## [37]  0.59  0.39  0.02 -0.33 -0.57 -0.67 -0.60 -0.38 -0.05  0.31  0.54  0.74
lag1.plot(log.electric3, 9)