#============================================================
#ARMA Model Fitting and ARMA Forecasting
#============================================================
#Examples from Z&W.
#US and Canadian Interest Rate Differentials
#====================================================

#Getting the Data and Putting in Labels.

uscn.id = 100*(lexrates.dat[,"USCNF"]-lexrates.dat[,"USCNS"])
colIds(uscn.id) = "USCNID"
uscn.id@title = "US/CA 30 day interest rate differential"

#Looking at the Picture

par(mfrow=c(2,1))
plot(uscn.id,reference.grid=F)
abline(h=0)
tmp = acf(uscn.id)


#Finding Out about arima.mle

args(arima.mle)

#Examples of ways to specfy the model and "starting values"

#   mod.list = list(order=c(1,0,1))
#   mod.list = list(order=c(1,0,1),ar=0.75,ma=0)
#   mod.list = list(ar=c(0.75,-0.25),ma=c(0,0))

#Caveat: arima.mle expects a "demeaned" series.
#Thus, either first demean, or (as we see later) set xreg=1

uscn.id.dm = uscn.id - mean(uscn.id)

#Our model selection

arma11.mod = list(ar=0.75,ma=0)

#Get the fit

arma11.fit = arima.mle(uscn.id.dm,model=arma11.mod)

#Look at what we have created.

names(arma11.fit)


#Or Look all at onece
arma11.fit

# Get standard errors for the estimators Standard
# by taking the square roots of the diagonal elements of 
# the variance-covariance matrix

std.errs = sqrt(diag(arma11.fit$var.coef))
names(std.errs) = colIds(arma11.fit$var.coef)
std.errs

#Note that is is pretty cool.

#Now, a second view of the same problem, but we use xreg
#When we take xreg=1 in essense we tell S-Plus
#Our series may not be "demeaned", so you take a look an tell us.
#Call any mean you find "the intercept"

arma11.fit2 = arima.mle(uscn.id,model=arma11.mod,xreg=1)


arma11.fit2


#Now lets take a look at our fits.
#It turns out that the plot function 
#applied to arma11.fit gives the basic plot
#and 3 diagnostic plots

plot(arma11.fit)

#Finally, let's produce monthly forecasts for
#the demeaned interest rate differential from July 1996 through June 1997

fcst.dates = timeSeq("7/1/1996", "6/1/1997", by="months", format="%b %Y")
uscn.id.dm.fcst = arima.forecast(uscn.id.dm, n=12, model=arma11.fit$model, future.positions=fcst.dates)


#Let's take a look at the stucture of our answer

names(uscn.id.dm.fcst)


#Forecasts and 95% forecast confidence intervals are produced by

smpl = positions(uscn.id.dm) >= timeDate("6/1/1995")

plot(
uscn.id.dm[smpl,],
uscn.id.dm.fcst$mean,
uscn.id.dm.fcst$mean+2*uscn.id.dm.fcst$std.err,
uscn.id.dm.fcst$mean-2*uscn.id.dm.fcst$std.err,
plot.args=list(lty=c(1,2,3,3))
)