Stat 6555 Homework 4
Enrique Rodriguez
May 24, 2016
TODO: 3.36, 3.37
31
In Example 3.39, we presented the diagnostics for the MA(2) fit to the GNP growth rate series. Using that example as a guide, complete the diagnostics for the AR(1) fit.
Solution
First we can visualize the time series data,
library(astsa)
#plot the time series data
plot(gnp)
#plot the acf of the data
acf(gnp, lag.max = 55)
The trend is hiding too much of the underlying signal, so we will take a log and difference.
gnp_diff <- diff(log(gnp), 1)
plot(gnp_diff)
acf2(gnp_diff)
## ACF PACF
## [1,] 0.35 0.35
## [2,] 0.19 0.08
## [3,] -0.01 -0.11
## [4,] -0.12 -0.12
## [5,] -0.17 -0.09
## [6,] -0.11 0.01
## [7,] -0.09 -0.03
## [8,] -0.04 -0.02
## [9,] 0.04 0.05
## [10,] 0.05 0.01
## [11,] 0.03 -0.03
## [12,] -0.12 -0.17
## [13,] -0.13 -0.06
## [14,] -0.10 0.02
## [15,] -0.11 -0.06
## [16,] 0.05 0.10
## [17,] 0.07 0.00
## [18,] 0.10 0.02
## [19,] 0.06 -0.04
## [20,] 0.07 0.01
## [21,] -0.09 -0.11
## [22,] -0.05 0.03
## [23,] -0.10 -0.03
## [24,] -0.05 0.00
## [25,] 0.00 0.01Model Fit and Diagnoses
sarima(gnp_diff, 1, 0, 0) # AR(1)## initial value -4.589567
## iter 2 value -4.654150
## iter 3 value -4.654150
## iter 4 value -4.654151
## iter 4 value -4.654151
## iter 4 value -4.654151
## final value -4.654151
## converged
## initial value -4.655919
## iter 2 value -4.655921
## iter 3 value -4.655922
## iter 4 value -4.655922
## iter 5 value -4.655922
## iter 5 value -4.655922
## iter 5 value -4.655922
## final value -4.655922
## 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, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 xmean
## 0.3467 0.0083
## s.e. 0.0627 0.0010
##
## sigma^2 estimated as 9.03e-05: log likelihood = 718.61, aic = -1431.22
##
## $AIC
## [1] -8.294403
##
## $AICc
## [1] -8.284898
##
## $BIC
## [1] -9.263748Standardized Resioduals: From this plot we can see that there is no pattern to make us believe this is something other than white noise. It should be noted that towards the mid 80’s to 2000’s it appears that the variance is smaller than before.
ACF of Residuals: The ACF of the residuals appear like what we would hope for, white noise. Very much like the standardized residuals, there appears to be nothing that would make us believe this is something other than white noise.
Normal QQ Plot: The qqplot does show some departure from normality on the tails of the residuals. As was noted by the author this is mainly due to outliers that occured in the 50’s and the 80
32
Crude oil prices in dollars per barrel are in oil; see Appendix R for more details. Fit an ARIMA(p, d, q) model to the growth rate performing all necessary diagnostics. Comment.
Solution
First lets visualize the time series.
plot.ts(oil)
We now look at the ACF and PACF to see if we can infer any neccesary transformations or differencing
acf(oil)
pacf(oil)
We can see that the ACF decays away very slowly, this is and the time series plot are indicative that we should take a log and difference.
oil_d <- diff(oil)
plot(oil_d)
We can see that the trend has been eliminated but it appears that the variance is non constant, since the tail end of the series appears to change. We can alliviate this problem by taking the log.
oil_dl <- diff(log(oil), 1)
plot(oil_dl)
At this point the series looks much better but it appears that we may have some outliers, namely and around 2009?
The next step is to establish some initially orders for the ARIMA.
par(mfrow = c(2,1))
acf(oil_dl, 300)
pacf(oil_dl, 300)
par(mfrow = c(1,1))From the above we can see that both the ACF and PACF seem to cutoff after about lag 1. This suggest, we should fit the model ARIMA(1,1,1).
# tried a buinch around the same cutoff point, this one is good
fit.1 <- sarima(oil_dl, 1.8,0,1, details = FALSE)
Diagnoses
Standardized Residuals: There appears to be some outliers in the data, but overall the residuals tend to look like white noise.
ACF of Residuals: Our first intuition about the residuals above is confirmed when we inspect the ACF of these. This plot stringly suggest that we have white noise left over, which is good.
Normal QQ Plot: The resdiuals tend to diverge from the theoretical quantiles at the tails, but overall a strong fit can be seen.
33
Fit an ARIMA(p, d, q) model to the global temperature data gtemp performing all of the necessary diagnostics. After deciding on an appropriate model, forecast (with limits) the next 10 years. Comment.
Solution
First lets see the time series
plot.ts(gtemp)
We can see a clear trend, let’s look at the ACF to confirm this,
acf(gtemp)
The slow decay of the ACF suggest us to take a difference, we do that now.
gtemp_d <- diff(gtemp)
plot.ts(gtemp_d)
Looks much better.
The next step is to plot the ACF and PACF of this differenced gtemp, and establish orders for the ARIMA model.
par(mfrow = c(2,1))
acf(gtemp_d)
pacf(gtemp_d)
par(mfrow = c(1,1))From the above we can see that the ACF appears to cut off at lag 2 and the PACF tails off. This is suggesting the model MA(2) for the differenced temp, or ARIMA(0, 1,2) for the gtemp. Let use this and fit the model.
temp.fit <- sarima(gtemp, 0, 1, 2, details = FALSE)
Diagnoses
Standardized Residuals: They have to apparent pattern, and so we can conclude tha they resemble white noise.
ACF of Residuals: Much like the Residuals Vs Time plot the ACF of these appears to resemble the acf of white noise and so we are happy.
QQ Plot: The qqplot is a great fit, even towards the tails, the divergence from the theoretical quantiles is minimal.
Forecast
sarima.for(gtemp, 10, 0, 1, 2)
## $pred
## Time Series:
## Start = 2010
## End = 2019
## Frequency = 1
## [1] 0.5623494 0.5543659 0.5608301 0.5672943 0.5737586 0.5802228 0.5866870
## [8] 0.5931513 0.5996155 0.6060797
##
## $se
## Time Series:
## Start = 2010
## End = 2019
## Frequency = 1
## [1] 0.09544714 0.10614790 0.10982195 0.11337700 0.11682392 0.12017201
## [7] 0.12342932 0.12660285 0.12969874 0.1327224536.
Plot (or sketch) the ACF of the seasonal \(ARIMA(0, 1) \times (1, 0)_{12}\) model with \(\phi = .8\) and \(\theta = .5\).
Solution
From results in the book we have that the ACF for this mixed model is given by,
\[\rho(12h + 1) = \frac{\theta}{1 + \theta^2}\phi^h \;\; h = 0, 1, 2,...\]
h <- 1:100 # lag vals
acf_mm <- function(p, t, h) {
((t)/(1 + t^2))*p^h
}
plot(h, acf_mm(.8, .5, h), xlab ="Lags", ylab = "ACF", type = "l")
37
Fit a seasonal ARIMA model of your choice to the unemployment data (unemp) displayed in Figure 3.21. Use the estimated model to forecast the next 12 months.