Lista de Exercícios - Séries Temporais
Victor Viana
27/01/2021
Neste documento estará contido todas as respostas solicitadas pelo prof. Adriano Lorena na 1ª parte da disciplina de Séries Temporais da Pós graduação do CIn - UFPE.
Capítulo 3 - Exercício 7
Letra a)
Na letra a, foi solicitado que plotemos a série original, a série média anual e o boxplot que resuma os valores observados para cada temporada.
www <- "/cloud/project/global.dat"
Global <- scan(www)
Global.ts <- ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.annual <- aggregate(Global.ts, FUN = mean)
par(mfrow=c(1,3))
plot(Global.ts)
plot(Global.annual)
boxplot(Global.ts ~ cycle(Global.ts))
Podemos observar que a série “Global” possui uma tendência, como pode ser visto na brusca subida do meio pro fim da série. Quando transformamos a série na média anual ficou mais claro de visualizar a tendência. Já quanto aos boxplots para os 12 meses, podemos observar que os meses de Junho, Julho, Agosto e Setembro possuem uma menor variabailidade dos dados perante os demais meses. Vale ressaltar que em alguns meses, são encontrados outliers, ou seja, pontos discrepantes perante os demais contidos naquele mês.
Letra b)
Na letra b, foi pedido para decompor a série nas componentes de tendência, efeito sazonal e resíduos e plotar a série decomposta. Além disso, produzir um gráfico da tendência com um efeito sazonal sobreposto.
Global.month.ts <- ts(Global, st = c(1856, 1), end = c(2005, 12),
frequency = 12)
Global.decom <- decompose(Global.month.ts, type = "mult")
trend <- Global.decom$trend
seasonal <- Global.decom$seasonal
plot(cbind(trend,trend*seasonal),lty=1:2)
Letra c)
Aqui, foi pedido o correlograma dos resíduos da questão 7b e explicá-lo.
Global.decom <- decompose(Global.ts, type = "mult")
acf(Global.decom$random[7:138])
pacf(Global.decom$random[7:138])
Em que poucos lags romperam a barreira, mas o lag 12 se tornou significativamente diferente de zero.
Letra d)
Foi pedido para aplicar um modelo Holt-Winters.
Global.hw <- HoltWinters(Global.month.ts, seasonal = 'multiplicative')
Global.hw## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Global.month.ts, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.7616377
## beta : 0.5693481
## gamma: 0.01225672
##
## Coefficients:
## [,1]
## a 0.1878636
## b -0.1255702
## s1 1.9764446
## s2 0.4762221
## s3 0.6218828
## s4 0.9229188
## s5 0.8111828
## s6 0.4177652
## s7 0.5730167
## s8 0.6151843
## s9 0.8981872
## s10 1.2871224
## s11 1.6120281
## s12 1.4695998plot(Global.hw)
O método Holt-Winters foi escolhido porque generaliza a equação de suavização exponencial simples e leva em consideração os efeitos sazonais e tendências. Ajustando os dados ao modelo de Holt-Winters, obtivemos estimativas de parâmetros conforme acima. Uma vez que os dados exibem alguma tendência sistemática e sazonal, é apropriado ajustar os dados usando o modelo Holt-Winters.
Letra e)
Nesta letra foi solicitado que projetemos os valores para os anos 2005–2010 com o modelo proposto acima. E quando é válido utilizar essa abordagem.
#library(forecast)
Global.month.ts <- ts(Global, st = c(1856, 1), end = c(2004, 12),
frequency = 12)
Global.hw <- HoltWinters(Global.month.ts, seasonal = 'multiplicative')
Global.hw## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Global.month.ts, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.6555837
## beta : 0.3804906
## gamma: 0.02359869
##
## Coefficients:
## [,1]
## a 0.26457179
## b -0.06271146
## s1 1.03227109
## s2 0.95958733
## s3 0.79937517
## s4 0.70475737
## s5 0.81129585
## s6 0.87817614
## s7 0.79620693
## s8 1.00713448
## s9 1.04654626
## s10 1.08560872
## s11 1.48118138
## s12 1.96875606global.predict <- predict(Global.hw,n.head=5*12,prediction.interval = TRUE)
global.predict## fit upr lwr
## Jan 2005 0.2083746 0.5018777 -0.08512853plot(fitted(Global.hw))
ts.plot(Global.month.ts,global.predict,lty=1:2)
Vale lembrar que a previsão é baseada no efeito de tendência durante o qual o modelo foi ajustado. É válido se a tendência e o efeito sazonal de 2006 a 2010 forem os mesmos que em anos anteriores.
Capítulo 4 - Exercício 2
Simule séries temporais de comprimento 100 de um modelo AR (1) com \(\alpha\) igual a −0,9, −0,5, 0,5 e 0,9.
set.seed(1)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.9 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
for (t in 2:100) x[t] <- -0.9 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
for (t in 2:100) x[t] <- 0.5 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
for (t in 2:100) x[t] <- -0.5 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
Capítulo 5 - Exercício 1
Letra a
Time = 1:100
w = rnorm(100, sd=25)
z = w
for (i in 2:100) z[i] = 0.5*z[i-1] + w[i]
x = z
x = 70 + 2*Time - 3*Time^2 + z
plot(x, type='l')
Letra b
Ajuste uma tendência quadrática para a série {xt}. Dê os coeficientes do modelo ajustado.
x.lm = lm(x ~ Time + I(Time^2))
coef(x.lm)## (Intercept) Time I(Time^2)
## 75.797847 1.466587 -2.994285Letra c
Encontre um intervalo de confiança de 95% para os parâmetros do modelo quadrático e comente
confint(x.lm)## 2.5 % 97.5 %
## (Intercept) 59.3115698 92.284124
## Time 0.7131194 2.220055
## I(Time^2) -3.0015131 -2.987058Ou seja, com 95% de confiança, as previsões feitas para o tempo estará entre 0.71 e 2.22.
Letra d
Trace o correlograma dos resíduos e comente.
acf(resid(x.lm))
O correlograma dos resíduos (obtido a partir do código acima) apresenta um valor signi cativo no lag 1, devido à autocorrelação introduzida usando um processo AR simulado de erros residuais.
Letra e
Forneça os erros padrão das estimativas dos parâmetros e comente.
acf(resid(x.lm))
Não estou conseguindo pôr o código dessa letra no R Markdown, mas irei colocar a análise baseada no código do console R que fiz. O GLS tem erros padrão mais amplos para as estimativas dos parâmetros. Isso é mais aparente em um intervalo de confiança.
Capítulo 5 - Exercício 3
Essa questão é baseada na série de produção elétrica (1958-190).
Letra a
Dê duas razões pelas quais uma transformação logarítmica pode ser apropriada para a série de eletricidade.
Os dados têm uma variação crescente. Se o desvio padrão for aproximadamente proporcional à média (coeficiente de variação constante), os logaritmos dos dados terão um desvio padrão constante (e, portanto, variância constante). Isso torna a modelagem mais fácil. Além disso, simular o logaritmo de uma variável garante que todos os valores simulados da própria variável são positivos.
Letra b
Ajustar um modelo de indicador sazonal com uma tendência quadrática para o logaritmo (natural) da série. Use a regressão stepwise para selecionar o melhor modelo com base no AIC.
cbe <- read.delim("/cloud/project/cbe.dat")
cbe[1:2,]## choc beer elec
## 1 1451 96.3 1497
## 2 2037 84.4 1463attach(cbe)
length(elec)## [1] 396length(elec)/12## [1] 33imth = rep(1:12,33)
T1 = 1:396
T2 = T1^2
elec.lm = lm(log(elec)~T1+T2+factor(imth))
step(elec.lm)## Start: AIC=-2717.56
## log(elec) ~ T1 + T2 + factor(imth)
##
## Df Sum of Sq RSS AIC
## <none> 0.3860 -2717.6
## - factor(imth) 11 2.4917 2.8777 -1944.1
## - T2 1 2.5629 2.9489 -1914.4
## - T1 1 20.3972 20.7833 -1141.1##
## Call:
## lm(formula = log(elec) ~ T1 + T2 + factor(imth))
##
## Coefficients:
## (Intercept) T1 T2 factor(imth)2 factor(imth)3
## 7.271e+00 7.960e-03 -6.883e-06 -1.991e-02 6.598e-02
## factor(imth)4 factor(imth)5 factor(imth)6 factor(imth)7 factor(imth)8
## 3.288e-02 1.462e-01 1.777e-01 2.375e-01 1.994e-01
## factor(imth)9 factor(imth)10 factor(imth)11 factor(imth)12
## 1.074e-01 9.044e-02 4.278e-02 2.350e-02Letra c
Ajuste um modelo harmônico com uma tendência quadrática para o logaritmo da série. Use a regressão stepwise para selecionar o melhor modelo com base no AIC.
c1=cos(2*pi*T1/12); s1=sin(2*pi*T1/12)
c2=cos(2*2*pi*T1/12); s2=sin(2*2*pi*T1/12)
c3=cos(3*2*pi*T1/12); s3=sin(3*2*pi*T1/12)
c4=cos(4*2*pi*T1/12); s4=sin(4*2*pi*T1/12)
c5=cos(5*2*pi*T1/12); s5=sin(5*2*pi*T1/12)
c6=cos(6*2*pi*T1/12)
elec.harm = lm(log(elec)~T1+T2+c1+s1+c2+s2+c3+s3+c4+s4+c5+s5+c6)
step(elec.harm)## Start: AIC=-2717.56
## log(elec) ~ T1 + T2 + c1 + s1 + c2 + s2 + c3 + s3 + c4 + s4 +
## c5 + s5 + c6
##
## Df Sum of Sq RSS AIC
## - s4 1 0.0002 0.3863 -2719.3
## - c3 1 0.0003 0.3863 -2719.2
## - c4 1 0.0005 0.3866 -2719.0
## <none> 0.3860 -2717.6
## - c5 1 0.0199 0.4059 -2699.7
## - c6 1 0.0252 0.4113 -2694.5
## - c2 1 0.0443 0.4303 -2676.6
## - s2 1 0.0451 0.4312 -2675.8
## - s3 1 0.0472 0.4332 -2673.9
## - s5 1 0.0949 0.4810 -2632.5
## - s1 1 0.6664 1.0525 -2322.4
## - c1 1 1.5472 1.9332 -2081.6
## - T2 1 2.5629 2.9489 -1914.4
## - T1 1 20.3972 20.7833 -1141.1
##
## Step: AIC=-2719.31
## log(elec) ~ T1 + T2 + c1 + s1 + c2 + s2 + c3 + s3 + c4 + c5 +
## s5 + c6
##
## Df Sum of Sq RSS AIC
## - c3 1 0.0003 0.3866 -2721.0
## - c4 1 0.0005 0.3868 -2720.8
## <none> 0.3863 -2719.3
## - c5 1 0.0199 0.4062 -2701.4
## - c6 1 0.0252 0.4115 -2696.2
## - c2 1 0.0443 0.4306 -2678.3
## - s2 1 0.0451 0.4314 -2677.5
## - s3 1 0.0472 0.4335 -2675.7
## - s5 1 0.0949 0.4812 -2634.3
## - s1 1 0.6664 1.0527 -2324.3
## - c1 1 1.5472 1.9335 -2083.6
## - T2 1 2.5629 2.9492 -1916.3
## - T1 1 20.3974 20.7837 -1143.1
##
## Step: AIC=-2720.99
## log(elec) ~ T1 + T2 + c1 + s1 + c2 + s2 + s3 + c4 + c5 + s5 +
## c6
##
## Df Sum of Sq RSS AIC
## - c4 1 0.0005 0.3871 -2722.4
## <none> 0.3866 -2721.0
## - c5 1 0.0199 0.4065 -2703.1
## - c6 1 0.0252 0.4118 -2697.9
## - c2 1 0.0443 0.4309 -2680.0
## - s2 1 0.0452 0.4317 -2679.2
## - s3 1 0.0472 0.4338 -2677.4
## - s5 1 0.0949 0.4815 -2636.0
## - s1 1 0.6664 1.0530 -2326.2
## - c1 1 1.5472 1.9338 -2085.5
## - T2 1 2.5629 2.9495 -1918.3
## - T1 1 20.3977 20.7842 -1145.1
##
## Step: AIC=-2722.43
## log(elec) ~ T1 + T2 + c1 + s1 + c2 + s2 + s3 + c5 + s5 + c6
##
## Df Sum of Sq RSS AIC
## <none> 0.3871 -2722.4
## - c5 1 0.0199 0.4070 -2704.6
## - c6 1 0.0252 0.4124 -2699.4
## - c2 1 0.0443 0.4314 -2681.5
## - s2 1 0.0452 0.4323 -2680.8
## - s3 1 0.0472 0.4343 -2678.9
## - s5 1 0.0949 0.4821 -2637.6
## - s1 1 0.6664 1.0536 -2328.0
## - c1 1 1.5472 1.9343 -2087.4
## - T2 1 2.5629 2.9500 -1920.2
## - T1 1 20.3980 20.7852 -1147.1##
## Call:
## lm(formula = log(elec) ~ T1 + T2 + c1 + s1 + c2 + s2 + s3 + c5 +
## s5 + c6)
##
## Coefficients:
## (Intercept) T1 T2 c1 s1 c2
## 7.363e+00 7.961e-03 -6.883e-06 -8.840e-02 -5.803e-02 1.496e-02
## s2 s3 c5 s5 c6
## 1.510e-02 -1.544e-02 1.003e-02 2.190e-02 -7.985e-03Letra d
Plote o correlograma e o correlograma parcial dos resíduos do modelo geral mais adequado e comente sobre os gráficos.
elec.best = lm(log(elec)~T1+T2+c1+s1+c2+s2+s3+c5+s5+c6)O melhor modelo adquirido possuiu um AIC de -2722. O correlograma indica que ainda há alguns efeitos sazonais presentes. O correlograma parcial é difícil de interpretar na presença de efeitos sazonais existentes. No entanto, o termo de defasagem 1 alto, seguido por valores baixos sugere que um processo AR (1) com termos sazonais para a série de erros poderia ser uma boa tentativa.
Letra e
Ajuste um modelo AR aos resíduos do modelo de melhor ajuste. Forneça a ordem do modelo AR de melhor ajuste e os parâmetros estimados do modelo.
ar(resid(elec.best))##
## Call:
## ar(x = resid(elec.best))
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.3812 0.2187 0.0853 0.0274 -0.0004 0.0350 -0.0241 0.0071
## 9 10 11 12 13 14 15 16
## 0.1094 0.0317 0.1208 0.1715 -0.0889 -0.0664 -0.0458 -0.0896
## 17 18 19 20 21 22 23
## 0.0211 0.0038 -0.1251 0.0162 -0.0539 -0.0161 0.1609
##
## Order selected 23 sigma^2 estimated as 0.0003933Letra f
Trace o correlograma dos resíduos do modelo AR e comente.
fit.ar <- ar(resid(elec.best))
acf(fit.ar$res[-(1:23)])
O correlograma indica que os resíduos do modelo de AR criado são aproximadamente RUÍDO BRANCO. Observe que os primeiros 23 valores tiveram que ser removidos, uma vez que um modelo AR (23) foi ajustado, então os primeiros 23 resíduos foram NA.
Letra h
Use o modelo de melhor ajuste para prever a produção de eletricidade para os anos de 1991 a 2000.
new.t = 397:(397+119)
new.c1 = cos(2*pi*new.t/12); new.s1 = sin(2*pi*new.t/12)
new.c2 = cos(2*2*pi*new.t/12); new.s2 = sin(2*2*pi*new.t/12)
new.c3 = cos(3*2*pi*new.t/12); new.s3 = sin(3*2*pi*new.t/12)
new.c5 = cos(5*2*pi*new.t/12); new.s5 = sin(5*2*pi*new.t/12)
new.c6 = cos(6*2*pi*new.t/12);
new.t2 = new.t^2
new.dat = data.frame(T1=new.t, T2=new.t2, c1=new.c1, s1=new.s1, c2=new.c2, s2=new.s2, s3=new.s3, c5=new.c5, s5=new.s5, c6=new.c6)
ar.pred = predict(ar(resid(elec.best)), n.ahead=120)
log.pred = predict(elec.best, new.dat)
elec.pred = exp(log.pred + ar.pred$pred + 0.5*0.0003933)
elec.ts = ts(elec, st=1958, fr=12)
elec.pred.ts = ts(elec.pred, st=1991, fr=12)
ts.plot(elec.ts, elec.pred.ts, lty=1:2)
Os valores previstos estão em elec.pred e foram ajustados usando um fator de correção de 1/2\(\sigma\)^2 para contabilizar o viés devido à obtenção de registros. Finalmente, os valores previstos são adicionados ao gráfico de tempo da série original. Ao visualizar este gráfico, é evidente que as previsões não são particularmente boas, uma vez que não seguem as tendências, e que melhores previsões provavelmente seriam obtidas usando um procedimento de Holt-Winters multiplicativo.
Capítulo 6 - Exercício 4
Letra a
Nesse exercício, foi pedido para ajustar um modelo de regressão adequado para a série de passageiros aéreos. Comente o correlograma dos resíduos do modelo de regressão ajustado.
Ajustamos um modelo OLS com tendência quadrática para a transformação Log do conjunto de dados AP. Como pode ser visto abaixo.
data(AirPassengers)
AP <- AirPassengers
AP## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432Time <- 1:length(AP)
Imth <- cycle(AP)
AP.lm <- lm(log(AP) ~ Time + I(Time^2) + factor(Imth))
coef(AP.lm)## (Intercept) Time I(Time^2) factor(Imth)2 factor(Imth)3
## 4.651379e+00 1.318368e-02 -2.148187e-05 -2.226964e-02 1.077856e-01
## factor(Imth)4 factor(Imth)5 factor(Imth)6 factor(Imth)7 factor(Imth)8
## 7.638788e-02 7.392931e-02 1.960325e-01 2.999749e-01 2.907230e-01
## factor(Imth)9 factor(Imth)10 factor(Imth)11 factor(Imth)12
## 1.461743e-01 8.144975e-03 -1.354009e-01 -2.132107e-02acf(resid(AP.lm))
O gráfico ACF dos resíduos do modelo OLS ajustado. O gráfico mostra uma correlação decrescente, mas com muitas defasagens significativas.
Letra b
Ajuste um modelo ARMA (p, q) para valores de p e q não superiores a 2 para a série residual do modelo de regressão ajustado. Escolha o modelo de melhor ajuste com base no AIC e comente seu correlograma.
Com base no AIC, o modelo ARMA de melhor ajuste é o ARMA (1, 0, 2), com valor AIC de -560,39.
Correlograma
A maioria das defasagens estão abaixo do nível significativo. Isso indica que o modelo é potencialmente um bom ajuste. As defasagens significativas podem ser ocasionadas devido a algum efeito sazonal que o modelo ARMA não consegue mapear.
Letra c
Aqui foi pedido uma previsão do número de passageiros viajando na companhia aérea em 1961.
Previsão