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))