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

set.seed(0)

library(forecast)

# Exercise 5
# 
library(Ecdat)
data(IncomeUK)
income = IncomeUK[,1]

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/ex_5_income_ts_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot(income)
if( save_plots ){ dev.off() }

# Plot the ACF of some differences: 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/ex_5_income_acf_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,4))
acf(income)
acf(diff(income))
acf(diff(income,4))
acf(diff(diff(income,4)))
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }
#

# What do we get if we look for the best fitting model using auto.arima:
auto.arima( income, ic="aic" )

auto.arima( income, ic="bic" )

arima( income, order=c(0,1,0), seasonal=list( order=c(0,1,0), period=4 ) )

# Exercise 6
#

# Part (a): 
library(AER)
library(forecast)
data(USMacroG)
unemp = USMacroG[,9]

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/ex_6_unemp_ts_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot(unemp)
grid()
if( save_plots ){ dev.off() }

# Plot the ACF of some differences: 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/ex_6_unemp_acf_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,4))
acf(unemp)
acf(diff(unemp))
acf(diff(unemp,4))
acf(diff(diff(unemp,4)))
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }
#

fit = auto.arima(unemp,ic="bic")

# Verify that our residuals have no autocorrelation: 
acf(fit$residuals)

# Part (b):
# 
# Perform bootstrap on this ARIMA(1,1,0)(0,0,2)[4] model (the model found above):
#
# Generate bootstrap samples 
n = length(unemp)
sink("unempBootstrap.txt")
set.seed(1998852)
for( iter in 1:100 ){
  y = simulate.Arima( fit ) # use the "simulate.Arima" command in the forecast library/package
  fit = auto.arima(y,ic="bic")
  print(fit)
}
sink()


# Exercise 7
#
# In Ecdat we have three columns: 
# 
# r = 91-day treasury bill rate
# y = the log or real GDP
# pi = the inflation rate
#
library(forecast)
data(Tbrate,package="Ecdat")
pi = Tbrate[,3]
auto.arima( pi, ic="aic" ) 
auto.arima( pi, ic="bic" )