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

z_t = c( 2.1, 2.6, 4.1, 3.9, 6.7, 5.1, 7.8, 9.3, 7.5, 4.1, 2.9, 2.6, 
         2.0, 3.2, 3.7, 4.5, 6.1, 6.5, 8.7, 9.1, 8.1, 4.9, 3.6, 2.0, 
         2.4, 3.3, 3.3, 4.0, 3.6, 6.2, 7.7, 6.8, 5.8, 4.1, 3.0, 1.6, 
         1.9, 3.0, 4.5, 4.2, 4.8, 5.7, 7.1, 4.8, 4.2, 2.3, 2.1, 1.6 )

#postscript("../../WriteUp/Graphics/Chapter9/prob_4_ts_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( z_t, pch=18 )
#dev.off()

#postscript("../../WriteUp/Graphics/Chapter9/prob_4_ts_acf_plot.eps", onefile=FALSE, horizontal=FALSE)
acf( z_t )
#dev.off()


w_t = diff(z_t, lag = 12)

#postscript("../../WriteUp/Graphics/Chapter9/prob_4_wts_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( w_t, pch=18 )
#dev.off()

#postscript("../../WriteUp/Graphics/Chapter9/prob_4_wts_acf_plot.eps", onefile=FALSE, horizontal=FALSE)
acf( w_t )
#dev.off()

# I think the model should be given by:
# 
m1 = arima( w_t, order=c(1,0,0) )
m2 = arima( z_t, order=c(1,0,0), seasonal=list( order=c(0,1,0), period=12 ) )

# The book suggests fitting this model:
#
m3 = arima( z_t, order=c(1,0,0), seasonal=list( order=c(0,1,1), period=12 ) )
m3

# Make the requested forecasts:
preds = predict( m3, n.ahead=24 )

#postscript("../../WriteUp/Graphics/Chapter9/prob_4_zts_predictions.eps", onefile=FALSE, horizontal=FALSE)
N = length(z_t)
plot( 1:N, z_t, type="l", col="black", xlab="discrete time", ylab="z_t", main="time series & predictions" )
t_oos = (N+1):(N+24)
z_oos = preds$pred
lines( c(N,t_oos), c(z_t[N],z_oos), type="l", col="red" )
#dev.off()