Time Series III
#SOCY7113arima.r
library(TSA)
library(forecast)
library(expsmooth)
# Let's look at two time series. In this data frame,
# Y is the monthly unemployment rate for women, 20+ yrs,
# and X is the same rate for men, 20+.
USunemp<-read.csv("http://www.courseserve.info/files/us-unempl-gender.csv")
plot(USunemp$Month,USunemp$Y)
plot(USunemp$Month,USunemp$X)
# We can ask whether or not the unemployment rate for men predicts
# the unemployment rate for women. We might use regular linear
# regression:
summary(lm(USunemp$Y~USunemp$X))
# But because we have time series, the estimated effect of X on Y
# is incorrect because we have violated the assumption that the
# observations are independent. This is the problem facing a
# researcher testing causal models of time series data.
# We can use some of the techniques we've looked at recently to
# make adjustments to the series in order to get a more accurate
# test of the relationship. (This is the Box & Jenkins model.)
# We can use moving average and autoregressive functions to find
# the best fit, and then use the residuals in our causal model.
# The residuals represent the series after the non-stationarity
# and autocorrelative components are removed.
acf(USunemp$Y)
acf(USunemp$X)
source("http://www.courseserve.info/files/SOCY7113ma3.r")
usy.ma3<-ma3(USunemp$Y)
usy.ma3
usx.ma3<-ma3(USunemp$X)
usx.ma3
# We can use a new function arima() to fit a model with
# moving average and autoregressive components.
# The function takes the form: arima(p,d,q) where
# p is the autoregressive period, q is the moving average
# period, and d is the differencing lag.
# We fit a model to the independent variable, X, and use
# it on Y as well. First, we'll vary p:
usx.100<-arima(USunemp$X, order=c(1,0,0))
usx.200<-arima(USunemp$X, order=c(2,0,0))
usx.300<-arima(USunemp$X, order=c(3,0,0))
usx.400<-arima(USunemp$X, order=c(4,0,0))
usx.500<-arima(USunemp$X, order=c(5,0,0))
# We'll use the goodness-of-fit measure, AIC, to test which
# model fits best:
AIC(usx.100,usx.200,usx.300,usx.400,usx.500)
# The model with p=4 fits best (smallest AIC).
# Now we'll vary q:
usx.401<-arima(USunemp$X, order=c(4,0,1))
usx.402<-arima(USunemp$X, order=c(4,0,2))
usx.403<-arima(USunemp$X, order=c(4,0,3))
usx.404<-arima(USunemp$X, order=c(4,0,4))
usx.405<-arima(USunemp$X, order=c(4,0,5))
AIC(usx.401,usx.402,usx.403,usx.404,usx.405)
# You could continue testing moving averages to find the
# minimum. It appears that q=3 fits best.
# Finally, we'll fit d:
usx.413<-arima(USunemp$X, order=c(4,1,3))
usx.423<-arima(USunemp$X, order=c(4,2,3))
usx.433<-arima(USunemp$X, order=c(4,3,3))
usx.443<-arima(USunemp$X, order=c(4,4,3))
usx.453<-arima(USunemp$X, order=c(4,5,3))
AIC(usx.413,usx.423,usx.433,usx.443,usx.453)
# The best fit is with d=2. So, our model for X
# is p=4, d=2, q=3 -- arima(4,2,3) -- which
# we named usx.423
# We take the residuals:
usx.r<-residuals(usx.423)
# Now we fit the same model, arima(4,2,3), to our
# dependent variable, Y:
usy.423<-arima(USunemp$Y, order=c(4,2,3))
# We take the residuals:
usy.r<-residuals(usy.423)
# Now, at last, we can use linear least-squares regression
# to test the causal relationship:
summary(lm(usy.r~usx.r))
# Compare this to the uncorrected model:
summary(lm(USunemp$Y~USunemp$X))
# NOTE: I've discovered that the forecast library has not
# been updated for R 3.0.1, which is the version of the base
# that we've got installed in the lab. So we won't be able to
# use the auto.arima() function.
# In fact, there is an easier way to fit the model, using the
# auto.arima() function in the forecast package. This function
# has the advantage of also iteratively testing seasonality.
usx.auto<-auto.arima(USunemp$X)
# Once you've autofit the model for X, use the arima() function
# with the same values of p, d, and q, on Y. Do not autofit Y!
# You must use the same model fit for both X and Y. If your autofit
# model for X was (3,0,0), your command for Y would be:
# usy.auto<-arima(USunemp$Y, order=c(3,0,0))
usx.autor<-residuals(usx.auto)
usy.autor<-residuals(usy.auto)
summary(lm(usy.autor~usx.autor))