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

set.seed(66)

require(graphics)

# lets generate an example AR(3) time series (with some arbitary coefficients):
nsamples = 100
x = arima.sim( n=nsamples, list( ar=c(0.9,-0.5,+0.2) ), sd = sqrt(2.0) )
postscript(file="../../WriteUp/Graphics/Chapter2/fig_2_5_a.eps", onefile=FALSE, horizontal=FALSE)
plot(1:nsamples,x,'o',main="AR(3) Series",xlab="index",ylab="ts value")
dev.off()

# construct various ARIMA(p,0,0) model on this data for different values of p
# Note: p=0 is the estimate of the mean only and is a 1 parameter model
# Q: why is not a two parameter model since we have to estimate the variance sigma^2 also? 
#
numOfParams = c()
allAIC = c()
for( p in 0:10 ){
  fit = arima( x, order=c(p,0,0) )
  AIC = nsamples * log( sum( fit$residuals^2 )/nsamples ) + 2 * (p+1)
  numOfParams = c(numOfParams,p+1)
  allAIC = cbind(allAIC,AIC)
}

postscript(file="../../WriteUp/Graphics/Chapter2/fig_2_5_b.eps", onefile=FALSE, horizontal=FALSE)
plot( numOfParams, allAIC, 'o', main="AIC Plot", xlab="Number of Parameters", ylab="AIC Value" )
dev.off()