# Pacotes
require(astsa)
require(Rssa)
require(ggplot2)
Questão 1
a)
data("USAccDeaths")
mod <- ssa(USAccDeaths, L=24)
plot(log(mod$sigma)) # Log dos 24 valores singulares
lines(log(mod$sigma))
# Outra forma:
# log_mod <- ssa(log(USAccDeaths), L=24)
# plot(log_mod, type = "values")
b) e c)
rec <- reconstruct(mod, groups = list(Trend = 1, `2_12` = 2:12, `2_5` = 2:5, `7_12` = 7:12))
plot(rec, add.residuals = F, add.original = TRUE, plot.method = "xyplot",
superpose = TRUE, auto.key = list(columns = 2))
d)
for3 <- forecast(mod, groups = list(1:6),
method = "recurrent", interval = "confidence",
only.intervals = FALSE,
len = 6, R = 100, level = 0.95)
plot(for3, include = 24, shadecols = "green", type = "l",
main = "Confidence intervals")
e)
mae <- mean(abs(USAccDeaths - (rec$Trend + rec$`2_12`))) ; mae
## [1] 95.81257
mape <- mean(abs((USAccDeaths - (rec$Trend + rec$`2_12`)) / USAccDeaths)) * 100 ; mape
## [1] 1.106093
Questão 2 (Exercício 9)
data("sunspotz")
# Sem suavização
# sunspotz.per <- mvspec(sunspotz, log="no")
# Com suavização - Kernel de Daniell
sunspotz.ave <- mvspec(sunspotz, kernel('daniell', 4), log='no')
## Bandwidth: 0.038
## Degrees of Freedom: 17.21
abline(v=c(.09,.19), lty=2) # Interessante notar que o período de 11 anos (1/11 = 0.09) é um ciclo conhecido na atividade solar
# Encontrando os índices dos picos
target <- .09
index <- which.min(abs(sunspotz.ave$freq - target)) ; index # Índice 22
## [1] 22
target <- .19
index <- which.min(abs(sunspotz.ave$freq - target)) ; index # Índice 46
## [1] 46
# Intervalo de confiança
df <- sunspotz.ave$df; df
## [1] 17.2125
L <- qchisq(.975, df); L
## [1] 30.4756
U <- qchisq(.025, df); U
## [1] 7.70502
sunspotz.ave$spec[22]
## [1] 11879.22
df*sunspotz.ave$spec[22]/L
## [1] 6709.336
df*sunspotz.ave$spec[22]/U
## [1] 26537.38
# Contribuição significativa para a variação observada.
# Intervalo amplo, indicando alguma variabilidade na estimativa.
sunspotz.ave$spec[46]
## [1] 681.6344
df*sunspotz.ave$spec[46]/L
## [1] 384.9844
df*sunspotz.ave$spec[46]/U
## [1] 1522.726
# Potência bem menor.
# Intervalo mais estreito, indicando uma estimativa mais precisa.
# Escala log(10)
sunspotz.ave <- mvspec(sunspotz, kernel('daniell',4), log = "y")
## Bandwidth: 0.038
## Degrees of Freedom: 17.21
abline(v=c(.09,.19), lty=2)
Questão 3 (Exercício 36)
a)
data("climhyd")
climhyd$Inflow <- log10(climhyd$Inflow)
climhyd$Precip <- sqrt(climhyd$Precip)
# Não consegui plotar as coerências quadradas sem os intervalos de confiança,
# então fiz um gráfico para cada variável
sr_Precip <- mvspec(cbind(climhyd$Inflow, climhyd$Precip), kernel("daniell",9), plot=FALSE)
sr_Temp <- mvspec(cbind(climhyd$Inflow, climhyd$Temp) , kernel("daniell",9), plot=FALSE)
sr_DewPt <- mvspec(cbind(climhyd$Inflow, climhyd$DewPt) , kernel("daniell",9), plot=FALSE)
sr_CldCvr <- mvspec(cbind(climhyd$Inflow, climhyd$CldCvr), kernel("daniell",9), plot=FALSE)
sr_WndSpd <- mvspec(cbind(climhyd$Inflow, climhyd$WndSpd), kernel("daniell",9), plot=FALSE)
par(mfrow=c(2,3))
plot(sr_Precip, plot.type = "coh", ci.lty = 2, ylim=c(0,1)) ; grid()
plot(sr_Temp, plot.type = "coh", ci.lty = 2, ylim=c(0,1)) ; grid()
plot(sr_DewPt, plot.type = "coh", ci.lty = 2, ylim=c(0,1)) ; grid()
plot(sr_CldCvr, plot.type = "coh", ci.lty = 2, ylim=c(0,1)) ; grid()
plot(sr_WndSpd, plot.type = "coh", ci.lty = 2, ylim=c(0,1)) ; grid()
b)
ajuste <- LagReg(climhyd$Precip, climhyd$Inflow, L=15, M=32)
## INPUT: climhyd$Precip OUTPUT: climhyd$Inflow L = 15 M = 32
##
## The coefficients beta(0), beta(1), beta(2) ... beta(M/2-1) are
##
## 0.02049094 0.009540118 0.007590805 0.003181126 0.002572341 0.00236956
## 0.0003186685 -0.002572417 -0.0007618132 -0.0007784056 -0.0008056257
## 0.0002694471 -0.0007644752 -0.003599572 -0.001470321 -0.001502055
##
##
## The coefficients beta(0), beta(-1), beta(-2) ... beta(-M/2+1) are
##
## 0.02049094 -0.0009190138 -0.001661455 -0.00278144 -0.002333444 -0.003982277
## -0.002226545 -0.001995739 -0.001160625 -0.002320983 -0.002110916 -0.002258649
## -0.001951727 -0.003696525 -0.0006236561 -0.002546621
## The positive lags, at which the coefficients are large
## in absolute value, and the coefficients themselves, are:
## lag s beta(s)
## [1,] 0 0.0204909428
## [2,] 1 0.0095401180
## [3,] 2 0.0075908052
## [4,] 3 0.0031811258
## [5,] 4 0.0025723411
## [6,] 5 0.0023695598
## [7,] 6 0.0003186685
## [8,] 7 -0.0025724173
## [9,] 8 -0.0007618132
## [10,] 9 -0.0007784056
## [11,] 10 -0.0008056257
## [12,] 11 0.0002694471
## [13,] 12 -0.0007644752
## [14,] 13 -0.0035995723
## [15,] 14 -0.0014703215
## [16,] 15 -0.0015020548
##
## The prediction equation is
## climhyd$Inflow(t) = alpha + sum_s[ beta(s)*climhyd$Precip(t-s) ], where alpha = 1.962931
## MSE = 0.01665116
