#install.packages("astsa") require(astsa) #help(package = astsa) # Original series plot(prodn, main = "frb production series") win.graph() resid1 = ts(lm(prodn~time(prodn))$resid) acf2(resid1, 48) # Note that with 'acf2' there is no "lag 0" acf # First difference win.graph() Dprod = diff(prodn) # No need now to detrend - differencing will also do this. plot(Dprod, main = "First difference of frb production series") win.graph() acf2(Dprod, 48) # Seasonal difference as well win.graph() D12D1prod = diff(Dprod, 12) plot(D12D1prod, main = "D12D1 series") win.graph() acf2(D12D1prod, 48) # fit.R = arima(prodn, order = c(1,1,1), seasonal = list(order = c(1,1,0), period = 12)) # R function fit.SS= sarima(prodn, p=1,d=1,q=1, P=1,D=1,Q=0, S=12) fit.SS names(fit.SS) # Put D=0 to see the difference between the fits: print(arima(prodn, order = c(1,1,1), seasonal = list(order = c(1,0,0), period = 12))) print(sarima(prodn, p=1,d=1,q=1, P=1,D=0,Q=0, S=12, details = F)$fit) # Try several models: out = matrix(ncol = 10) d = 1 D = 1 S = 12 Q = 1 for (p in 1:2) { for (q in 0:2) { for (P in 1:2) { fits = sarima(prodn, p, d, q, P, D, Q, S) out = rbind(out, c(p, d, q, P, D, Q, S, fits$AIC, fits$AICc, fits$BIC)) }}} out = out[-1,] colnames(out) = c("p", "d", "q", "P", "D", "Q", "S", "AIC", "AICc","BIC") out out = round(out,3) out[out[,8] == min(out[,8])] # min AIC out[out[,9] == min(out[,9])] # min AICc out[out[,10] == min(out[,10])] # min BIC # Examine the candidates # (i) 1,1,0, 2,1,1, 12 Chosen by BIC # (ii) 2,1,0, 2,1,1, 12 Chosen by AIC and AICc # (iii) 1,1,1, 2,1,1, 12 Chosen by S&S # more closely: win.graph() p=1;q=0;P=2;Q=1 # Chosen by BIC fit1 = sarima(prodn, p,1,q, P,1,Q, S=12, details = F) title(sub = paste("p,q,P,Q =", p,q,P,Q), outer = T, line = -1) win.graph() p=2;q=0;P=2;Q=1 # Chosen by AIC fit2 = sarima(prodn, p,1,q, P,1,Q, S=12, details = F) title(sub = paste("p,q,P,Q =", p,q,P,Q), outer = T, line = -1) win.graph() p=1;q=1;P=2;Q=1 # Chosen by S&S fit3 = sarima(prodn, p,1,q, P,1,Q, S=12, details = F) title(sub = paste("p,q,P,Q =", p,q,P,Q), outer = T, line = -1) stdres1 = fit1$fit$residuals/sqrt(fit1$fit$sigma2) # Standardized residuals stdres2 = fit2$fit$residuals/sqrt(fit2$fit$sigma2) stdres3 = fit3$fit$residuals/sqrt(fit3$fit$sigma2) win.graph() par(mfrow=c(3,1)) plot(stdres1, lag(D12D1prod,-1), ylab = "res/std.dev(resid)") plot(stdres2, lag(D12D1prod,-1), ylab = "res/std.dev(resid)") plot(stdres3, lag(D12D1prod,-1), ylab = "res/std.dev(resid)") # Model 2 - arima(2,1,0, 2,1,1, 12) - has now been chosen (tentatively): win.graph() fit = fit2 print(fit) acf2(fit2$fit$residuals) stdres = stdres2 # Test normality par(mfrow=c(2,1)) hist(stdres) qqnorm(stdres) shapiro.test(stdres) # Look at some obviously poor fits: win.graph() poorfit1 = sarima(prodn, 0,1,0, 2,1,1, S = 12) win.graph() acf2(poorfit1$fit$residuals) win.graph() poorfit2 = sarima(prodn, 2,1,0, 0,1,0, S=12) win.graph() acf2(poorfit2$fit$residuals) # Forecasts for the final model win.graph() prod.pr = sarima.for(prodn, n.ahead = 12, 2,1,0, 2,1,1, 12) abline(v=1979-.05)