library(astsa)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(pander)
Part A
rm(list = ls())
# Loading the data
data(cmort, package="astsa")
# Autoregressive Model (AR(2)) Fitting
max.lag <- 52 # Define maximum lag for ACF/PACF
invisible(acf2(cmort, max.lag = max.lag)) # Display ACF and PACF plots
cmort.ar2 <- ar.ols(cmort, order=2, demean=FALSE, intercept=TRUE) # Fit AR(2) model
pander(data.frame(mean=cmort.ar2$asy.se.coef$x.mean), caption="Mean of Estimates") # Displaying the mean of estimates
pander(data.frame(se.estimates=cmort.ar2$asy.se.coef$ar), caption="Standard Errors of the Estimates") # Displaying the standard errors
Standard Errors of the Estimates
| 0.03979 |
| 0.03976 |
# ARIMA(2,0,0) Model Fitting and Diagnostics
cmort.arima.fit <- arima(cmort, order=c(2,0,0)) # Fit ARIMA(2,0,0) model
print(cmort.arima.fit) # Printing ARIMA model summary
##
## Call:
## arima(x = cmort, order = c(2, 0, 0))
##
## Coefficients:
## ar1 ar2 intercept
## 0.4301 0.4424 88.8538
## s.e. 0.0397 0.0398 1.9407
##
## sigma^2 estimated as 32.37: log likelihood = -1604.71, aic = 3217.43
# Simulating an ARIMA process based on the fitted AR(2) model's coefficients
cmort.ar2.sim <- arima.sim(list(order=c(2,0,0), ar=c(0.4280, 0.4412)), sd = sqrt(32.32), n.start=52*5, n=500) + mean(cmort)
# Plot the original cmort data and the simulated ARIMA process for comparison
plot.zoo(cbind(ts(cmort), ts(cmort.ar2.sim)), plot.type = "single", col = c("red", "blue"))
legend("topleft", c("cmort", "AR(2) Simulated"), fill = c("red", "blue"))
# Perforiming ACF analysis on the simulated series to check autocorrelation
acf(cmort.ar2.sim, lag.max = max.lag, main="ACF of Simulated Series")
cmort.fit <- sarima(cmort, 2, 0, 0)
## initial value 2.300097
## iter 2 value 2.148325
## iter 3 value 1.756362
## iter 4 value 1.751728
## iter 5 value 1.737837
## iter 6 value 1.737835
## iter 7 value 1.737834
## iter 8 value 1.737832
## iter 9 value 1.737831
## iter 10 value 1.737830
## iter 11 value 1.737829
## iter 12 value 1.737824
## iter 13 value 1.737818
## iter 14 value 1.737810
## iter 15 value 1.737805
## iter 16 value 1.737804
## iter 17 value 1.737804
## iter 18 value 1.737804
## iter 18 value 1.737804
## iter 18 value 1.737804
## final value 1.737804
## converged
## initial value 1.740035
## iter 2 value 1.740030
## iter 3 value 1.740016
## iter 4 value 1.740012
## iter 5 value 1.740006
## iter 6 value 1.740000
## iter 7 value 1.739986
## iter 8 value 1.739968
## iter 9 value 1.739959
## iter 10 value 1.739956
## iter 11 value 1.739954
## iter 12 value 1.739953
## iter 13 value 1.739952
## iter 14 value 1.739952
## iter 15 value 1.739952
## iter 16 value 1.739951
## iter 17 value 1.739951
## iter 18 value 1.739951
## iter 19 value 1.739950
## iter 20 value 1.739950
## iter 21 value 1.739949
## iter 22 value 1.739949
## iter 23 value 1.739949
## iter 24 value 1.739948
## iter 24 value 1.739948
## iter 24 value 1.739948
## final value 1.739948
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 0.4301 0.0397 10.8404 0
## ar2 0.4424 0.0398 11.1245 0
## xmean 88.8538 1.9407 45.7839 0
##
## sigma^2 estimated as 32.3709 on 505 degrees of freedom
##
## AIC = 6.333521 AICc = 6.333615 BIC = 6.366832
##
print(cmort.fit) # Display the SARIMA model fit summary
## $fit
##
## Call:
## 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 ar2 xmean
## 0.4301 0.4424 88.8538
## s.e. 0.0397 0.0398 1.9407
##
## sigma^2 estimated as 32.37: log likelihood = -1604.71, aic = 3217.43
##
## $degrees_of_freedom
## [1] 505
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.4301 0.0397 10.8404 0
## ar2 0.4424 0.0398 11.1245 0
## xmean 88.8538 1.9407 45.7839 0
##
## $ICs
## AIC AICc BIC
## 6.333521 6.333615 6.366832
Part B
ols.fit = ar.ols(cmort, order=2, demean=FALSE, intercept=TRUE)
foreca = predict(ols.fit, n.ahead=24,interval="prediction")
ts.plot(cmort, foreca$pred, col=1:2, ylab="Mortality")
lines(foreca$pred, type="p", col=2)
lines(foreca$pred+1.96*foreca$se, lty="dashed", col=4)
lines(foreca$pred-1.96*foreca$se, lty="dashed", col=4)
foreca = predict(ols.fit, n.ahead=4,interval="prediction")
predicted <- foreca$pred
right.pi <- foreca$pred+1.96*foreca$se
left.pi <- foreca$pred-1.96*foreca$se
pi.width <- 1.96*foreca$se
pander(data.frame(left.pi=left.pi, predicted=predicted,right.pi=right.pi, pi.width=pi.width),caption="Predicted with the limits of prediction interval.")
Predicted with the limits of prediction interval.
| 76.46 |
87.6 |
98.74 |
11.14 |
| 74.64 |
86.76 |
98.89 |
12.12 |
| 73.35 |
87.34 |
101.3 |
13.98 |
| 72.33 |
87.21 |
102.1 |
14.88 |