Questão 1
phi <- 0.6
theta <- 0.9
set.seed(123)
# ARMA(1,0)
arma10 <- arima.sim(n = 100, list(order = c(1, 0, 0), ar = phi))
# ARMA(0,1)
arma01 <- arima.sim(n = 100, list(order = c(0, 0, 1), ma = theta))
# ARMA(1,1)
arma11 <- arima.sim(n = 100, list(order = c(1, 0, 1), ar = phi, ma = theta))
acf_arma10 <- acf(arma10, plot = FALSE)
acf_arma01 <- acf(arma01, plot = FALSE)
acf_arma11 <- acf(arma11, plot = FALSE)
plot(acf_arma10$acf ~ acf_arma10$lag, type = 'l', col = 'steelblue', lwd = 2, xlab = 'lag', ylab = 'acf')
lines(acf_arma01$acf ~ acf_arma10$lag, type = 'l', col = '#ff5454', lwd = 2)
lines(acf_arma11$acf ~ acf_arma10$lag, type = 'l', col = 'darkgreen', lwd = 2)
legend("topright",
legend = c("ARMA(1,0)", "ARMA(0,1)", "ARMA(1,1)"),
col = c("steelblue", "#ff5454", "darkgreen"),
lwd = 2,
bty = "n")
A ACF do ARMA(1,0) mostra uma queda contínua, a do ARMA(0,1) uma queda mais forte depois do primeiro lag, e a do ARMA(1,1) combina características dos dois modelos.
Questão 2
a)
Podemos ajustar da seguinte maneira:
require(astsa)
require(forecast)
data(cmort)
ajuste <- ar.ols(cmort, order=2,demean = F, intercept = T)
ajuste##
## Call:
## ar.ols(x = cmort, order.max = 2, demean = F, intercept = T)
##
## Coefficients:
## 1 2
## 0.4286 0.4418
##
## Intercept: 11.45 (2.394)
##
## Order selected 2 sigma^2 estimated as 32.32b)
# 4 predições
predicoes <- predict(ajuste,n.ahead = 4)
#Média das predições
pred_med <- predicoes$pred
#EP das predições
pred_ep <- predicoes$se
# IC
IC_inf <- pred_med - 1.96 * pred_ep; IC_sup <- pred_med + 1.96 * pred_ep
# Resultados
knitr::kable(
round(data.frame(Semana = 1:4,
Pred = pred_med,
'IC Inferior' = IC_inf,
'IC Superior' = IC_sup
),2))| Semana | Pred | IC.Inferior | IC.Superior |
|---|---|---|---|
| 1 | 87.60 | 76.46 | 98.74 |
| 2 | 86.76 | 74.64 | 98.89 |
| 3 | 87.34 | 73.35 | 101.32 |
| 4 | 87.21 | 72.33 | 102.10 |
Questão 3
set.seed(123)
# 500 observações em 3 simulações
amostra <- list()
for (i in 1:3) {
amostra[[i]] <- arima.sim(list(ar = 0.9, ma = -0.9), 500)
}
# Plot das séries simuladas
par(mfrow = c(3, 1))
plot(amostra[[1]], main = 'Primeira Simulação')
plot(amostra[[2]], main = 'Segunda Simulação')
plot(amostra[[3]], main = 'Terceira Simulação')
# ACF e PACF das séries
par(mfrow = c(3, 2))
acf(amostra[[1]], main = 'ACF Primeira Amostra'); pacf(amostra[[1]], main = 'PACF Primeira Amostra')
acf(amostra[[2]], main = 'ACF Segunda Amostra'); pacf(amostra[[2]], main = 'PACF Segunda Amostra')
acf(amostra[[3]], main = 'ACF Terceira Amostra'); pacf(amostra[[3]], main = 'PACF Terceira Amostra')
Podemos observar um comportamento muito próximo tanto entre as séries quanto entre as suas respectivas autocorrelações.
# Ajuste ARMA(1,1) aos dados simulados
ajuste1 <- arima(amostra[[1]], order =c(1,0,1))
ajuste1 ##
## Call:
## arima(x = amostra[[1]], order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.6698 -0.7393 0.0319
## s.e. 0.1870 0.1685 0.0341
##
## sigma^2 estimated as 0.9322: log likelihood = -691.94, aic = 1391.88ajuste2 <- arima(amostra[[2]], order =c(1,0,1))
ajuste2 ##
## Call:
## arima(x = amostra[[2]], order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## -0.8053 0.7881 0.0425
## s.e. 0.2691 0.2780 0.0455
##
## sigma^2 estimated as 1.054: log likelihood = -722.66, aic = 1453.32ajuste3 <- arima(amostra[[3]], order =c(1,0,1))
ajuste3 ##
## Call:
## arima(x = amostra[[3]], order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.067 -0.0024 0.0373
## s.e. 0.656 0.6574 0.0470
##
## sigma^2 estimated as 0.9669: log likelihood = -701.04, aic = 1410.09Podemos observar que os coeficientes que o segundo e do primeiro ajuste são os que mais se aproximam dos valores do modelo amostrado.
Questão 4
a)
Podemos começar observando o comportamento da série e vemos que há um padrão de crescimento constante, muito semelhante ao que vemos no Exemplo III.39 do material, mas com alguma sazonalidade.
plot(sales, main = "Vendas (Sales)")
Podemos trabalhar com o trabalhar com o log dos dados e uma diferencaição para observarmos a taxa de crescimento das vendas. Ao fazermos isso vemso que a série apresenta ser um processo bem mais estável.
log_sales <- log(sales)
diff_sales <- diff(log_sales)
plot(diff_sales,xlab="Tempo", ylab=expression(paste(nabla,log,"(",Y[t],")")),
main = 'Taxa de crescimento das vendas (sales)')
Por fim, antes de ajustar o modelo vamos olhar para as funções ACF e PACF amostrais. Ainda não apresenta um padrão satisfatório, onde no PACF podemos observar um lag significativo apenas no primrito valor.
par(mfrow=c(1,2))
acf(diff_sales, main = 'ACF'); pacf(diff_sales, main = 'PACF')
Por fim, como nosso modelo inicial vamos ajustar um modelo \(ARIMA(1,0,1)\).
ajuste1 <- arima(diff_sales, order = c(1,0,1))
ajuste1##
## Call:
## arima(x = diff_sales, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.8400 -0.6149 0.0017
## s.e. 0.0834 0.1182 0.0011
##
## sigma^2 estimated as 3.473e-05: log likelihood = 553.41, aic = -1098.81Analisando os resíduos do modelo podemos observar que i modelo tem um ajuste ben razóavel aos dados.
sarima(diff_sales, 1,0,1)## initial value -5.050724
## iter 2 value -5.080469
## iter 3 value -5.094406
## iter 4 value -5.096106
## iter 5 value -5.097662
## iter 6 value -5.106879
## iter 7 value -5.107101
## iter 8 value -5.111293
## iter 9 value -5.116769
## iter 10 value -5.122205
## iter 11 value -5.127630
## iter 12 value -5.131855
## iter 13 value -5.132116
## iter 14 value -5.132296
## iter 15 value -5.132315
## iter 16 value -5.132333
## iter 17 value -5.132334
## iter 18 value -5.132336
## iter 19 value -5.132336
## iter 19 value -5.132336
## iter 19 value -5.132336
## final value -5.132336
## converged
## initial value -5.132904
## iter 2 value -5.132951
## iter 3 value -5.133070
## iter 4 value -5.133076
## iter 5 value -5.133077
## iter 6 value -5.133078
## iter 7 value -5.133079
## iter 7 value -5.133079
## iter 7 value -5.133079
## final value -5.133079
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 0.8400 0.0834 10.0728 0.0000
## ma1 -0.6149 0.1182 -5.2003 0.0000
## xmean 0.0017 0.0011 1.5210 0.1304
##
## sigma^2 estimated as 3.473379e-05 on 146 degrees of freedom
##
## AIC = -7.374589 AICc = -7.373478 BIC = -7.293946
## 
b)
Podemos exibir o gráfico CCF da seguinte forma:
diff_lead <- diff(log(lead))
ccf(diff_sales, diff_lead, main = 'CCF de diff_sales e diff_lead')
Utilizando lag2.plot() para as séries diferenciadas,
podemos observar que não parece ter uma relação clara entre as duas
séries sem atribuir um lag.
lag2.plot(diff_sales, diff_lead)
c)
diff_lead_lag3 <- lag(diff_lead, -3) # Defasagem de 3
ajuste_arima <- Arima(diff_sales, xreg = diff_lead_lag3, order = c(1, 0, 1))
summary(ajuste_arima)## Series: diff_sales
## Regression with ARIMA(1,0,1) errors
##
## Coefficients:
## ar1 ma1 intercept xreg
## 0.8441 -0.6198 0.0017 -0.0075
## s.e. 0.0815 0.1160 0.0012 0.0161
##
## sigma^2 = 3.564e-05: log likelihood = 553.51
## AIC=-1097.03 AICc=-1096.61 BIC=-1082.01
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4.235582e-05 0.005889205 0.004580906 NaN Inf 0.760782 -0.01055152Questão 5
a)
Podemos observar que o valor médio apresenta um comportamento de queda ao longo dos anos, com um início de estabilização a partir de 1990.
plot(cpg, main = 'Grafico de Ct')
b)
Transfotmando nossa variável em \(log(C_t)\), podemos observar que o modelo linear se ajusta bem aos dados, indicando claramente que a nossa variável original apresenta um comportamento exponencial no período.
log_cpg <- log(cpg)
tempo <- 1980:2008
ajuste <- lm(log_cpg ~ tempo)
plot(tempo, log_cpg, main = "Grafico de Log(Ct) com Regressao Linear",
xlab = "Tempo", ylab = "log(Ct)", pch = 16, col = "blue")
abline(ajuste, col = "red", lwd = 2)
legend("topright", legend = c("Dados log(Ct)", "Regressao Linear"),
col = c("blue", "red"), pch = c(16, NA), lty = c(NA, 1), lwd = c(NA, 2))
c)
Embora no gráfico da alternativa anterior o modelo aparente se ajustar bem aos dados, podemos ver nos gráficos abaixo uma ligeira fuga de normalidade e um padrão claro nos gráficos de resíduos.
par(mfrow = c(2,2))
plot(ajuste)
d)
Podemos ajustar um modelo arima para os dados, utilizando a função
auto.arima().
ajuste_arima <- auto.arima(log_cpg, xreg = tempo)
summary(ajuste_arima)## Series: log_cpg
## Regression with ARIMA(1,0,0) errors
##
## Coefficients:
## ar1 intercept xreg
## 0.8297 1113.0105 -0.5554
## s.e. 0.1190 73.5665 0.0368
##
## sigma^2 = 0.181: log likelihood = -15.37
## AIC=38.73 AICc=40.4 BIC=44.2
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04793293 0.4028656 0.3336757 0.8235877 14.53447 0.5859718
## ACF1
## Training set -0.2191326Também podemos observar o ajuste sobre os dados reais. Dessa vez vemos que o modelo tem um ajuste bem superior ao linear.
previsoes <- fitted(ajuste_arima) # fitted() obtém os valores ajustados
# Plot dos dados reais e do ajuste ARIMA
plot(tempo, log_cpg, main = "Grafico de Log(Ct) com Ajuste ARIMA",
xlab = "Tempo", ylab = "log(Ct)", pch = 16, col = "steelblue")
lines(previsoes, col = 'red4', lwd = 2) # Adiciona o modelo ajustado ao gráfico
legend("topright", legend = c("Dados Reais", "Ajuste ARIMA"),
col = c("steelblue", "red4"), lty = 1, lwd = 2)