plot.aic <- function(fit, new=T, sd=0.1)
{  # code to plot AIC against order of AR model fitted
	if (new) plot((1:length(fit$aic))-1,fit$aic,type="l",xlab="Order",ylab="AIC") else 
	         lines((1:length(fit$aic))-1,fit$aic,type="l")
    points(rnorm(1,fit$order,sd=sd),rnorm(1,sd=sd),pch=16,col="red") 
}

# generates data from an autoregressive process, by default of order 1
sim.y <- function(n, model=list(ar=c(0.9)))  arima.sim(model=model, n)

n <- 25  # length of time series
R <- 1000  # number of replicates

# first dataset to get things started
y <- sim.y(n,list(ar=c(0.5,0.1)))
fit <- ar(y,order.max=19)
plot.aic(fit)

# we will store the orders chosen using AIC, BIC, and AICC
AIC.order <- NULL
BIC.order <- NULL
AICC.order <- NULL

# Now make R replicates, plot the corresponding AIC curves
for (i in 1:R ) 
{
	fit <- ar( sim.y(n, model=list(ar=c(0.5,0.1))), order.max=19 )
	plot.aic(fit, new=F)
	AIC.order <- c(AIC.order, fit$order)
# The next two lines should be uncommented and modified to give the
#  optimal orders when BIC and AICC are used for order selection
#	BIC.order <- c(BIC.order, NA)    
#	AICC.order <- c(AICC.order, NA)
}

# tabulate the order of the chosen model
table(AIC.order)