#APTS Statistical Modelling 2012 #Practical 1 - time series and model selection #function from the worksheet 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") } } #function to plot *IC plot.ic <- function(fit, new=T, sd=0.1, pen="AIC", pl=T,ptitle="") { # code to plot AIC against order of AR model fitted p<-(0:fit$order.max) n<-fit$n.used tl<- -(fit$aic-2-2*p) ic<-switch(pen, BIC = p*log(n)-tl, AICc = n*(n+p)/(n-p-2)-tl, 2*p+2-tl ) ic<-ic-min(ic) if(pl) { if (new) { plot((1:length(fit$aic))-1,ic,type="l",xlab="Order",ylab=pen,ylim=c(0,50),main=ptitle) } else { lines((1:length(fit$aic))-1,ic,type="l") points(rnorm(1,as.double(which.min(ic))-1,sd=sd),rnorm(1,sd=sd),pch=16,col="red") } } return(ic=ic) } # 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) sim.plot<-function(n,R,pen="AIC",slist=list(ar=c(0.5,0.1)),ptitle="") { # first dataset to get things started y <- sim.y(n,slist) fit <- ar(y) #split.screen(c(2,2)) #screen(1) plot.ic(fit,pen=pen,ptitle=ptitle) # we will store the orders chosen using AIC, BIC, and AICC IC.order <- NULL # Now make R replicates, plot the corresponding *IC curves for (i in 1:R ) { fit <- ar( sim.y(n, model=list(ar=c(0.5,0.1))),max.order=19 ) plot.ic(fit, new=F,pen=pen) IC.order <- c(IC.order, as.double(which.min(plot.ic(fit, pen=pen, pl=F)))-1) } # tabulate the order of the chosen model table(IC.order) } #(a) n <-c(25,50,100,500) # length of time series R <- 100 # number of replicates par(mfrow=c(2,2)) sim.plot(n[1],R,ptitle="n=25") sim.plot(n[2],R,ptitle="n=50") sim.plot(n[3],R,ptitle="n=100") sim.plot(n[4],R,ptitle="n=500") #(b) n <-c(25,50,100,500) # length of time series R <- 100 # number of replicates par(mfrow=c(2,2)) sim.plot(n[1],R,ptitle="n=25",slist=list(ma=0.9)) sim.plot(n[2],R,ptitle="n=50",slist=list(ma=0.9)) sim.plot(n[3],R,ptitle="n=100",slist=list(ma=0.9)) sim.plot(n[4],R,ptitle="n=500",slist=list(ma=0.9)) #(c) n <-c(25,50,100,500) # length of time series R <- 100 # number of replicates par(mfrow=c(2,2)) sim.plot(n[1],R,ptitle="n=25",pen="BIC") sim.plot(n[2],R,ptitle="n=50",pen="BIC") sim.plot(n[3],R,ptitle="n=100",pen="BIC") sim.plot(n[4],R,ptitle="n=500",pen="BIC") par(mfrow=c(2,2)) sim.plot(n[1],R,ptitle="n=25",pen="AICc") sim.plot(n[2],R,ptitle="n=50",pen="AICc") sim.plot(n[3],R,ptitle="n=100",pen="AICc") sim.plot(n[4],R,ptitle="n=500",pen="AICc") par(mfrow=c(2,2)) sim.plot(n[1],R,ptitle="n=25",pen="BIC",slist=list(ma=0.9)) sim.plot(n[2],R,ptitle="n=50",pen="BIC",slist=list(ma=0.9)) sim.plot(n[3],R,ptitle="n=100",pen="BIC",slist=list(ma=0.9)) sim.plot(n[4],R,ptitle="n=500",pen="BIC",slist=list(ma=0.9)) par(mfrow=c(2,2)) sim.plot(n[1],R,ptitle="n=25",pen="AICc",slist=list(ma=0.9)) sim.plot(n[2],R,ptitle="n=50",pen="AICc",slist=list(ma=0.9)) sim.plot(n[3],R,ptitle="n=100",pen="AICc",slist=list(ma=0.9)) sim.plot(n[4],R,ptitle="n=500",pen="AICc",slist=list(ma=0.9))