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

set.seed(0)

# Exercise 1
# 
library(Ecdat)
data(CRSPday)
crsp = CRSPday[,7]

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/ex_1_acf_crsp.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
acf(crsp)
acf(as.numeric(crsp))
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# Exercise 2
#
arima(crsp,order=c(1,0,0))
arima(crsp,order=c(2,0,0))

arima(crsp,order=c(0,0,1)) # just to see how this model stacks up
#

# Exercise 14
#
library(Ecdat)
data(Mishkin)
tb1 = log(Mishkin[,3])

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

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/ex_14_tb1_acf_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
acf(tb1)
acf(diff(tb1))
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# b
#
# Note: My results are obtained using the forecast library version 5.9:
#
# You can determine your version of the forecast library by calling sessionInfo()
# 
library(forecast)
auto.arima( tb1, max.P=0, max.Q=0, ic="aic" ) 

# Fit the model we found:
#
aic_bm_fit = arima( tb1, order=c(5,1,3) )
bic_bm_fit = arima( tb1, order=c(0,1,1) )

# Plot the acf of the residuals of this model:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/ex_14_tb1_residual_acf_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
acf( residuals( aic_bm_fit ) )
acf( residuals( bic_bm_fit ) )
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }