Ajustar um modelo \(ARIMA(p,d,q)\) para os dados da temperatura global em globtemp realizando todos os diagnósticos necessários. Depois de decidir sobre um modelo apropriado, faça uma previsão com limites de confiança para os próximos 10 anos. Comente. 

library(astsa)
par(mfrow=c(1,1), mar=c(3,2,3,1)+.5, mgp=c(1.6,.6,0))
data(globtemp)
plot(globtemp, type="o", xlab="Tempo", ylab="Desvios Globais de Temperatura", pch=19, , ylim = c(-0.45,0.85), 
     main="Desvios de temperatura global média anual entre 1880 e 2015\n em graus centígrados")
grid()

Observamos estabilidade nos desvios globais de temperatura desde 1880 até 1940, ou seja, desde começos da Primeira Revolução Industrial até o começo da Segunda Guerra Mundial. Durante a guerra houve um pico de alta na temperatura e depois estabilzou-se ao nível de 1940 até 1980, quando aconteceu a Segunda Revolução Industrial, com o começo da industrialização no leste da Asia, Oceania e América Latina. Desde então apresenta crescimento contínuo.

Pode-se então sugerir uma regressão segmentada ou uma deferenciação dos dados: 

library(segmented)
tempo = time(globtemp)
Segunda.guerra = rep(0,length(tempo))
Segunda.guerra[60] = Segunda.guerra[61] = Segunda.guerra[62] = Segunda.guerra[63] = Segunda.guerra[64] = Segunda.guerra[65] = Segunda.guerra[66] = 1 
ajuste1a = lm(globtemp ~ tempo+Segunda.guerra, na.action=NULL)
summary(ajuste1a)
## 
## Call:
## lm(formula = globtemp ~ tempo + Segunda.guerra, na.action = NULL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32676 -0.11175 -0.02114  0.11620  0.38571 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.364e+01  6.642e-01  -20.54   <2e-16 ***
## tempo           7.011e-03  3.409e-04   20.56   <2e-16 ***
## Segunda.guerra  1.418e-01  6.057e-02    2.34   0.0208 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.156 on 133 degrees of freedom
## Multiple R-squared:  0.7619, Adjusted R-squared:  0.7584 
## F-statistic: 212.8 on 2 and 133 DF,  p-value: < 2.2e-16
ajuste2a = segmented(ajuste1a, seg.Z = ~ tempo, psi = c(1940,1980))
summary(ajuste2a)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = ajuste1a, seg.Z = ~tempo, psi = c(1940, 1980))
## 
## Estimated Break-Point(s):
##                 Est. St.Err
## psi1.tempo 1909.000  3.911
## psi2.tempo 1974.389  3.053
## 
## Meaningful coefficients of the linear terms:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.349651   3.633313   2.298   0.0232 *  
## tempo          -0.004541   0.001918  -2.368   0.0194 *  
## Segunda.guerra  0.246528   0.036379   6.777 3.95e-10 ***
## U1.tempo        0.010180   0.002010   5.065       NA    
## U2.tempo        0.011967   0.001342   8.916       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09092 on 129 degrees of freedom
## Multiple R-Squared: 0.9216,  Adjusted R-squared: 0.9179 
## 
## Boot restarting based on 7 samples. Last fit:
## Convergence attained in 2 iterations (rel. change -7.955e-08)
slope(ajuste2a)
## $tempo
##              Est.    St.Err. t value  CI(95%).l   CI(95%).u
## slope1 -0.0045406 0.00191780 -2.3676 -0.0083350 -0.00074617
## slope2  0.0056391 0.00060107  9.3818  0.0044499  0.00682840
## slope3  0.0176060 0.00120000 14.6710  0.0152320  0.01998100
intercept(ajuste2a)
## $tempo
##                Est.
## intercept1   8.3497
## intercept2 -11.0830
## intercept3 -34.7110
par(mfrow=c(1,1), mar=c(3,2,1,1)+.5, mgp=c(1.6,.6,0))
plot.segmented(ajuste2a, lwd=2, col = "red", ylim = c(-0.45,0.85))
points(tempo, globtemp, pch=19, col="black", type="b")
grid()

Percebemos que nossas observações não foram corretas nas datas de mudanças de comportamento nem sequer na observação de estabilidade na temperatura nos primeiros anos de registro. O ajuste da regressão segmentada indica diminuição dos desvios de temperatura até 1908, um primeiro incremento do desvio de temperatura até 1975 e depois um incremento mais acentuado dos desvios de temperatura

Série sem tendência

Resposta = globtemp - fitted(ajuste2a)
par(mfrow=c(1,1), mar=c(3,2,3,1)+.5, mgp=c(1.6,.6,0))
plot(Resposta, lwd=2, col = "black", type = "b", pch = 19)
abline(h=0, col="red", lwd=3, lty=2)
grid()

Observemos as funções de autocorrelação e de autocorrelação parcial. 

acf2(Resposta)

##      [,1]  [,2]  [,3] [,4]  [,5]  [,6]  [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.23 -0.10 -0.10 0.02 -0.12 -0.03  0.02 0.04 -0.11  0.01 -0.06 -0.11 -0.21
## PACF 0.23 -0.16 -0.04 0.05 -0.17  0.05 -0.01 0.00 -0.11  0.06 -0.12 -0.09 -0.17
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## ACF  -0.04  0.01 -0.01  0.05  0.10  0.01  0.06  0.03  0.03
## PACF -0.02 -0.04 -0.07  0.07 -0.01  0.00  0.07 -0.02  0.01

Estas funções sugerem fortemente um modelo autoregressivo de ordem 1. Por esse motivo procedemos a ajustarmos três modelos autoregressivos de até ordem 3, para depois utilizarmos algum critério de escolha de modelos. 

modelo1 = sarima(Resposta, 1,0,0, details = FALSE); modelo1a = arima(Resposta, order=c(1,0,0))
modelo2 = sarima(Resposta, 2,0,0, details = FALSE)
modelo3 = sarima(Resposta, 3,0,0, details = FALSE)

Procedemos à escolha de modelos. 

modelo1$AIC; modelo1$AICc; modelo1$BIC
## [1] -2.0218
## [1] -2.021136
## [1] -1.95755
modelo2$AIC; modelo2$AICc; modelo2$BIC
## [1] -2.034557
## [1] -2.03322
## [1] -1.948891
modelo3$AIC; modelo3$AICc; modelo3$BIC
## [1] -2.021487
## [1] -2.019242
## [1] -1.914404

Selecionamos pelo BIC o modelo1, autoregressivo de ordem 1. 

modelo1$ttable
##       Estimate     SE t.value p.value
## ar1     0.2324 0.0835  2.7848  0.0061
## xmean   0.0002 0.0096  0.0210  0.9833

Resíduos

sarima(Resposta, 1,0,0)
## initial  value -2.420603 
## iter   2 value -2.448494
## iter   3 value -2.448499
## iter   4 value -2.448500
## iter   5 value -2.448500
## iter   5 value -2.448500
## iter   5 value -2.448500
## final  value -2.448500 
## converged
## initial  value -2.451895 
## iter   2 value -2.451897
## iter   3 value -2.451897
## iter   4 value -2.451897
## iter   4 value -2.451897
## iter   4 value -2.451897
## final  value -2.451897 
## converged

## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xmean, include.mean = FALSE, transform.pars = trans, 
##     fixed = fixed, optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1   xmean
##       0.2324  0.0002
## s.e.  0.0835  0.0096
## 
## sigma^2 estimated as 0.007415:  log likelihood = 140.48,  aic = -274.96
## 
## $degrees_of_freedom
## [1] 134
## 
## $ttable
##       Estimate     SE t.value p.value
## ar1     0.2324 0.0835  2.7848  0.0061
## xmean   0.0002 0.0096  0.0210  0.9833
## 
## $AIC
## [1] -2.0218
## 
## $AICc
## [1] -2.021136
## 
## $BIC
## [1] -1.95755

PREVISÃO

Utilizamos como referência o Exemplo III.25 

fore = predict(modelo1a, n.ahead = 10) # utilizamos o resultado do ajuste com o comando arima
par(mfrow = c(1,1),mar=c(4,3,1,1),mgp=c(1.6,.6,0), pch=19)
ts.plot(Resposta, fore$pred, col=1:2, xlim=c(1975,2025), lwd=2, ylab="Desvios de temperatura sem tendência", xlab="Tempo")
U = fore$pred+fore$se; L = fore$pred-fore$se
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
lines(fore$pred, type="p", col=2)
grid()

Para transformar os dados à escala original devemos prever a tendência, obtida com a regressão segmentada, para os próximos 10 anos.

Estimativa da reta de regressão segmentada. 

reta.pred = intercept(ajuste2a)$tempo[3]+slope(ajuste2a)$tempo[3]*time(fore$pred)

Obtemos agora os valores preditos na escala original dos dados. 

par(mfrow = c(1,1),mar=c(4,3,1,1),mgp=c(1.6,.6,0), pch=19)
ts.plot(globtemp,fore$pred+reta.pred, col=1:2, xlim=c(1880,2025), lwd=2, ylab="Desvios de temperatura sem tendência", xlab="Tempo")
U = fore$pred+reta.pred+fore$se; L = fore$pred+reta.pred-fore$se
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
lines(fore$pred+reta.pred, type="p", col=2)
grid()

Poderíamos utilizar primeira diferença nesta situação. 

par(mfrow=c(1,1), mar=c(3,2,3,1)+.5, mgp=c(1.6,.6,0))
plot(diff(globtemp), lwd=2, col = "black", type="b", pch=19)
abline(h=0, col="red", lwd=3, lty=2)
grid()

e procedermos propondo e ajustando um novo modelo adequado: ARIMA(1,1,0), modelo autoregressivo na primeira diferença. 

modelo4 = sarima(globtemp, 1,1,0, details = FALSE); modelo4a = arima(globtemp, order=c(1,1,0))

ou ARIMA(1,0,0), modelo autoregressivo aplicado à primeira diferença. 

modelo5 = sarima(diff(globtemp), 1,0,0, details = FALSE); modelo5a = arima(diff(globtemp), order=c(1,0,0))

é o mesmo! Acreditamos que estimar a tendência e depois ajustar o modelo fornece mais informação.