#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = F

# In Ecdat we have three columns: 
# 
# r = 91-day treasury bill rate
# y = the log or real GDP
# pi = the inflation rate
#
data(Tbrate,package="Ecdat")

library(tseries)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/rlab_tbrate_ts_plots.eps", onefile=FALSE, horizontal=FALSE) }
plot(Tbrate)
if( save_plots ){ dev.off() }

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/rlab_tbrate_acf.eps", onefile=FALSE, horizontal=FALSE) }
acf(Tbrate)
if( save_plots ){ dev.off() }

adf.test(Tbrate[,1])
adf.test(Tbrate[,2])
adf.test(Tbrate[,3])

# Lets take the first difference and see if that gives us stationary timeseries: 
diff_rate = diff(Tbrate)

pairs(diff_rate) # scatterplot matrix

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/rlab_tbrate_diff_ts_plots.eps", onefile=FALSE, horizontal=FALSE) }
plot(diff_rate)
if( save_plots ){ dev.off() }

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/rlab_tbrate_diff_acf_plots.eps", onefile=FALSE, horizontal=FALSE) }
acf(diff_rate)
if( save_plots ){ dev.off() }

adf.test(diff_rate[,1])
adf.test(diff_rate[,2])
adf.test(diff_rate[,3])

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/rlab_tbrate_mean_diff.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,1))
boxplot(diff_rate[,1] ~ cycle(diff_rate))
if( save_plots ){ dev.off() }

library(forecast)
auto.arima( Tbrate[,1], max.P=0, max.Q=0, ic="aic" )
#auto.arima( Tbrate[,1], max.P=0, max.Q=0, ic="bic" )

# Refit our model:
#
fit1 = arima( Tbrate[,1], order=c(0,1,1) )
acf( residuals( fit1 ) )
Box.test( residuals(fit1), lag=10, type="Ljung" )

resid2 = residuals(fit1)^2
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/rlab_tbrate_resid2_acf_plot.eps", onefile=FALSE, horizontal=FALSE) }
acf(resid2)
if( save_plots ){ dev.off() }

Box.test(resid2, lag=10, type="Ljung")

# Next lets forecast some ARIMA models using R:
#
data(Tbrate,package="Ecdat")
# r = the 91-day Treasury bill rate
# y = log of real GDP
# pi = the inflation rate
# fit the nonseasonal ARIMA model found by auto.arima
library(forecast)
pi = Tbrate[,3]
auto.arima( pi, max.P=0, max.Q=0, ic="bic" )
fit = arima( pi, order=c(1,1,1) )
forecasts = predict(fit,36)

if( save_plots ){ postscript(file="../../WriteUp/Graphics/Chapter9/rlab_forecast_plots.eps", onefile=FALSE, horizontal=FALSE) }
plot(pi,xlim=c(1980,2006),ylim=c(-7,12))
lines(seq(from=1997,by=0.25,length=36),forecasts$pred,col="red")
lines(seq(from=1997,by=0.25,length=36),forecasts$pred + 1.96*forecasts$se,col="blue")
lines(seq(from=1997,by=0.25,length=36),forecasts$pred - 1.96*forecasts$se,col="blue")
grid()
if( save_plots ){ dev.off() }