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

#
# Seasonal ARIMA Models:
#

# P1 & 2: EPage 278
# 
library("Ecdat")
data(IncomeUK)
consumption = IncomeUK[,2] 

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/consumption_ts_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot(consumption)
grid()
if( save_plots ){ dev.off() }

# Consider consumption on its own:
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/consumption_acfs.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,4))
acf( consumption ) 
acf( diff(consumption,1) )
acf( diff(consumption,4) ) 
acf( diff(diff(consumption,1),4) ) 
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# Consider the logarithm of consumption:
#
lconsumption = log(consumption)
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/log_consumption_acfs.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,4))
acf( lconsumption ) 
acf( diff(lconsumption,1) )
acf( diff(lconsumption,4) ) 
acf( diff(diff(lconsumption,1),4) ) 
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# Fit the model that looks most likely:
# 
consumption_ts_model = arima( consumption, order=c(0,1,0), seasonal=list(order=c(0,1,1), period=4) )
Box.test( residuals(consumption_ts_model), lag=10, type="Ljung" ) 

log_consumption_ts_model = arima( lconsumption, order=c(0,1,0), seasonal=list(order=c(0,1,1), period=4) )
Box.test( residuals(log_consumption_ts_model), lag=10, type="Ljung" ) 

# P3: Epage 278
# 
acf( residuals(log_consumption_ts_model) )

# P4: EPage 278
# 
library(forecast)
auto.arima( lconsumption, ic="bic" )

# P5: EPage 278
# 
consumption_prediction_m_1 = predict( consumption_ts_model, n.ahead=8 ) 
consumption_prediction_m_2 = predict( log_consumption_ts_model, n.ahead=8 )

m_1_pred = consumption_prediction_m_1$pred
m_1_ll = m_1_pred - 2*consumption_prediction_m_1$se
m_1_ul = m_1_pred + 2*consumption_prediction_m_1$se

m_2_pred = exp( consumption_prediction_m_2$pred )
m_2_ll = exp( consumption_prediction_m_2$pred - 2*consumption_prediction_m_2$se )
m_2_ul = exp( consumption_prediction_m_2$pred + 2*consumption_prediction_m_2$se )

plot_y_max = max( c( m_1_ul, m_2_ul ) )
plot_y_min = min( c( m_1_ll, m_2_ll ) )

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/consumption_predictions.eps", onefile=FALSE, horizontal=FALSE) }
plot( m_1_pred, col='red', ylim=c(plot_y_min,plot_y_max) )
points( m_1_ll, type='l', lty=2, col='red' )
points( m_1_ul, type='l', lty=2, col='red' ) 

points( m_2_pred, type='l', col='green' )
points( m_2_ll, type='l', lty=2, col='green' )
points( m_2_ul, type='l', lty=2, col='green' ) 

grid()
if( save_plots ){ dev.off() }


#
# VAR Models:
#

data(Tbrate,package="Ecdat")
# r = the 91-day Treasury bill rate
# y = the log of real GDP
# pi = the inflation rate
#
del_dat = diff(Tbrate)
var1 = ar(del_dat,order.max=4,aic=T)
var1

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/var_residual_acf.eps", onefile=FALSE, horizontal=FALSE) }
acf(var1$resid[-1,])
if( save_plots ){ dev.off() }

# Predict the next time points 
r = c(-1.41,-0.48,0.66)
y = c(-0.019420,0.015147,0.003303)
pi = c(2.31,-1.01,0.31)

df = data.frame( r, y, pi )
newdata = ts( df, start=c(1997, 1), end=c(1997,3), frequency=4 ) 
predict(var1, newdata=newdata )


#
# Long-Memory Processes
#

# P8: EPage 279
# 
data(Mishkin,package="Ecdat")
cpi = as.vector(Mishkin[,5])
DiffSqrtCpi = diff(sqrt(cpi))

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/diff_sqrt_cpi_ts_N_acf.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
plot( DiffSqrtCpi, type='l' )
acf( DiffSqrtCpi )
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

library("fracdiff")
fit.frac = fracdiff(DiffSqrtCpi,nar=0,nma=0) # estimate the amount of fractional differencing to apply
fit.frac$d
fdiff = diffseries(DiffSqrtCpi,fit.frac$d)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter10/frac_differenced_diff_sqrt_cpi_acf.eps", onefile=FALSE, horizontal=FALSE) }
acf(fdiff)
if( save_plots ){ dev.off() }

# This gives a ARIMA(0,1,1) model: 
auto.arima(fdiff,ic="aic")

# This gives a ARIMA(0,0,0) model: 
auto.arima(fdiff,d=0,ic="aic")


#
# Model-Based Bootstrapping of an ARIMA Process
#

library(AER)
library(forecast)
data("FrozenJuice")
price = FrozenJuice[,1]
plot(price)
auto.arima(price,ic="bic")

# Generate bootstrap samples 
n = length(price)
sink("priceBootstrap.txt")
set.seed(1998852)
for( iter in 1:10 ){
  eps = rnorm(n+20)
  y = rep(0,n+20)
  for( t in 3:(n+20) ){
    y[t] = 0.2825 * y[t-1] + 0.057 * y[t-2] + eps[t]
  }
  y = y[101:n+20]
  y = cumsum(y)
  y = ts(y,frequency=12)
  fit = auto.arima(y,d=1,D=0,ic="bic")
  print(fit)
}
sink()


set.seed(1998852)
niter = 250
estimates = matrix(0,nrow=niter,ncol=2)
for( iter in 1:niter){
  eps = rnorm(n+20)
  y = rep(0,n+20)
  for( t in 3:(n+20) ){
    y[t] = 0.2825 * y[t-1] + 0.057 * y[t-2] + eps[t]
  }
  y = y[101:n+20]
  y = cumsum(y)
  t = ts(y,frequency=12)
  fit = arima(y,order=c(2,1,0))
  estimates[iter,] = fit$coef
}
# Compute the bias, standard deviations, and MSE  of the given estimate of the two AR(2) parameters: 
bias = c( 0.2825, 0.057 ) - apply( estimates, 2, mean )
std_devs = apply( estimates, 2, sd )
mses = std_devs / sqrt(dim(estimates)[1])
t(data.frame( bias, std_devs, mses ))