#Scripts are given below for all Figures and Tables in McLeod and Zhang (2008b).
#

#Figure 1. Plot of lynx time series using plot.ts
plot(lynx)

#Figure 2. Plot of lynx series using TimeSeriesPlot
TimeSeriesPlot(lynx, type="o", pch=16, ylab="# pelts", main="Lynx Trappings")

#Figure 3. Trellis plot for Ninemile series
graphics.off() #clear previous graphics
data(Ninemile)
print(TimeSeriesPlot(Ninemile, SubLength=200))

#Figure 4. Partial autocorrelation plot of lynx series 
graphics.off() #clear previous graphics
PacfPlot(log(lynx))

## Not run: #takes some time for all these examples
#Figure 5. Using SelectModel to select the best subset ARz or ARp and
#          comparing BIC and UBIC subset selection.
#
graphics.off() #clear previous graphics
layout(matrix(1:4,ncol=2),respect=TRUE)
ansBICp<-SelectModel(log(lynx),lag.max=15,Criterion="BIC", ARModel="ARp", Best=3)
ansUBICp<-SelectModel(log(lynx),lag.max=15, ARModel="ARp", Best=3)
ansBICz<-SelectModel(log(lynx),lag.max=15,Criterion="BIC", ARModel="ARz", Best=3)
ansUBICz<-SelectModel(log(lynx),lag.max=15, ARModel="ARz", Best=3)
par(mfg=c(1,1))
plot(ansBICp)
par(mfg=c(2,1))
plot(ansUBICp)
par(mfg=c(1,2))
plot(ansBICz)
par(mfg=c(2,2))
plot(ansUBICz)

#Figure 6. Logged spectral density function fitted to square-root of monthly
#          sunspot series using the non-subset AR and subset ARz.
#          AIC and BIC are used for the AR while BIC and UBIC are used
#          for the ARz. Takes about 115 seconds on 3.6 GHz Pentium PC.
graphics.off() #clear previous graphics
layout(matrix(1:4,ncol=2),respect=TRUE)
z<-sqrt(sunspots)
P<-200
pAIC<-SelectModel(z, lag.max=P, ARModel="AR", Best=1, Criterion="AIC")
ARAIC<-FitAR(z, pAIC)
par(mfg=c(1,1))
sdfplot(ARAIC)
title(main="AIC Order Selection")
pBIC<-SelectModel(z, lag.max=P, ARModel="AR", Best=1, Criterion="BIC")
ARBIC<-FitAR(z, pBIC)
par(mfg=c(1,2))
sdfplot(ARBIC)
title(main="BIC Order Selection")
SunspotMonthARzBIC<-SelectModel(z,lag.max=P, ARModel="ARz", Best=1, Criterion="BIC")
ARzBIC<-FitAR(z, SunspotMonthARzBIC)
par(mfg=c(2,1))
sdfplot(ARzBIC)
title(main="BIC Subset Selection")
SunspotMonthARzUBIC<-SelectModel(z,lag.max=P, ARModel="ARz", Best=1)
ARzUBIC<-FitAR(z, SunspotMonthARzUBIC)
par(mfg=c(2,2))
sdfplot(ARzUBIC)
title(main="UBIC Subset Selection")

#Table 3.
#First part of table: AR(1) and AR(2).
#Only timings for GetFitAR and FitAR since the R function ar produces too many
#  warnings and an error message as noted in McLeod and Zhang (2008b, p.12). 
#The ar function with mle option is not recommended.
start.time<-proc.time()
set.seed(661177723)
NREP<-100 #takes about 156 sec
NREP<-10 #takes about 16 sec
ns<-c(50,100,200,500,1000)
ps<-c(1,2) #AR(p), p=1,2
tmsA<-matrix(numeric(4*length(ns)*length(ps)),ncol=4)
ICOUNT<-0
for (IP in 1:length(ps)){
p<-ps[IP]
for (ISIM in 1:length(ns)){
    ICOUNT<-ICOUNT+1
    n<-ns[ISIM]
    ptm <- proc.time()
    for (i in 1:NREP){
        phi<-PacfToAR(runif(p, min=-1, max =1))
        z<-SimulateGaussianAR(phi,n)
        phiHat<-try(GetFitAR(z,p,MeanValue=mean(z))$phiHat)
        }
    t1<-(proc.time() - ptm)[1]
#
    ptm <- proc.time()
    for (i in 1:NREP){
        phi<-PacfToAR(runif(p, min=-1, max =1))
        z<-SimulateGaussianAR(phi,n)
        phiHat<-try(FitAR(z,p,MeanMLEQ=TRUE)$phiHat)
        }
    t2<-(proc.time() - ptm)[1]
#
    ptm <- proc.time()
    for (i in 1:NREP){
        phi<-PacfToAR(runif(p, min=-1, max =1))
        z<-SimulateGaussianAR(phi,n)
        #uncomment this line and next two lines for ar timings -- expect lots of
        #   warnings and an error message!!
        #phiHat<-try(ar(z,aic=FALSE,order.max=p,method="mle")$ar)
        #delete this line and the next one
        phiHat<-NA
    }
    #uncomment this line for ar timings
         #t3<-(proc.time() - ptm)[1]
         t3<-NA #delete this line for ar timings

        tmsA[ICOUNT,]<-c(n,t1,t2,t3)
  } 
} 
rnames<-c(rep("AR(1)", length(ns)),rep("AR(2)", length(ns)) )
cnames<-c("n", "GetFitAR", "FitAR", "ar")  
dimnames(tmsA)<-list(rnames,cnames)
tmsA[,-1]<-round(tmsA[,-1]/NREP,2)
end.time<-proc.time()
total.time<-(end.time-start.time)[1]

#Second part of table: AR(20) and AR(40).
#NOTE: ar is not recommended with method="mle" produces numerous warnings
#      and also takes a long time!
        start.time<-proc.time()
        set.seed(661177723)
        NREP<-100 #takes 7.5 hours
        NREP<-10 #takes 45 minutes
        ns<-c(1000,2000,5000)
        ps<-c(20,40)
        tmsB<-matrix(numeric(4*length(ns)*length(ps)),ncol=4)
        ICOUNT<-0
        for (IP in 1:length(ps)){
        p<-ps[IP]
        phi<-PacfToAR(0.8/(1:p))
        for (ISIM in 1:length(ns)){
            ICOUNT<-ICOUNT+1
            n<-ns[ISIM]
            ptm <- proc.time()
            for (i in 1:NREP){
                z<-SimulateGaussianAR(phi,n)
                phiHat<-try(GetFitAR(z,p,MeanValue=mean(z))$phiHat)
                }
            t1<-(proc.time() - ptm)[1]
            ptm <- proc.time()
            for (i in 1:NREP){
                z<-SimulateGaussianAR(phi,n)
                phiHat<-try(FitAR(z,p,MeanMLEQ=TRUE)$phiHat)
                }
            t2<-(proc.time() - ptm)[1]
            ptm <- proc.time()
            for (i in 1:NREP){
                z<-SimulateGaussianAR(phi,n)
                phiHat<-try(ar(z,aic=FALSE,order.max=p,method="mle")$ar)
                }
            t3<-(proc.time() - ptm)[1]
            tmsB[ICOUNT,]<-c(n,t1,t2,t3)
        } 
        }   
        rnames<-c( rep("AR(20)", length(ns)), rep("AR(40)", length(ns)) )
        cnames<-c("n", "GetFitAR", "FitAR", "ar")  
        dimnames(tmsB)<-list(rnames,cnames)
        tmsB[,-1] <- round(tmsB[,-1]/NREP,2)
        end.time<-proc.time()
        total.time<-(end.time-start.time)[1]

#Figure 7. Comparing Box-Cox analyses using FitAR and MASS
library(MASS)
graphics.off() #clear previous graphics
layout(matrix(c(1,2,1,2),ncol=2))
pvec<-c(1,2,4,10,11)
out<-FitAR(lynx, ARModel="ARp", pvec)
BoxCox(out)
PMAX<-max(pvec)
Xy <- embed(lynx, PMAX + 1)
y <- Xy[, 1]
X <- (Xy[, -1])[, pvec] #pvec != 1
outlm<-lm(y~X)
boxcox(outlm,lambda=seq(0.0,0.6,0.05))

#Figure 8
graphics.off() #clear previous graphics
BoxCox(AirPassengers) #takes about 30 sec

#Figure 9
graphics.off() #clear previous graphics
data(rivers)
BoxCox(rivers)
title(sub="Length of 141 North American Rivers")

#Figure 10
graphics.off() #clear previous graphics
data(USTobacco)
TimeSeriesPlot(USTobacco, aspect=1)

#Figure 11
graphics.off() #clear previous graphics
data(USTobacco)
outUST<-arima(USTobacco, c(0,1,1))
BoxCox(outUST)

#Figure 12. Basic diagnostic plots for ARp fitted to the log lynx series
graphics.off() #clear previous graphics
out<-FitAR(log(lynx), ARModel="ARp", c(1,2,4,10,11))
plot(out, terse=TRUE)

#Figure 13. RSF plot for ARp fitted to log lynx series
graphics.off() #clear previous graphics
out<-FitAR(log(lynx), ARModel="ARp", c(1,2,4,10,11))
rfs(out)

#Table 6. Comparison of bootstrap and large-sample sd
#Use bootstrap to compute standard errors of parameters
#takes about 34 seconds on a 3.6 GHz PC
ptm <- proc.time() #user time
set.seed(2491781) #for reproducibility
R<-100  #number of bootstrap iterations
p<-c(1,2,4,7,10,11)
ans<-FitAR(log(lynx),p)
out<-Boot(ans, R)
fn<-function(z) FitAR(z,p)$zetaHat
sdBoot<-sqrt(diag(var(t(apply(out,fn,MARGIN=2)))))
sdLargeSample<-coef(ans)[,2][1:6]
sd<-matrix(c(sdBoot,sdLargeSample),ncol=2)
dimnames(sd)<-list(names(sdLargeSample),c("Bootstrap","LargeSample"))
ptm<-(proc.time()-ptm)[1]
sd


## End(Not run)