# the following commands simulate different length series
# from an AR(4) process with coefs ar4.real

ar4.real <- c(2.7607, -3.8106,  2.6535, -0.9258)
ar4.64 <- arima.sim(64,model=list(ar=ar4.real))
ar4.256 <-  arima.sim(256,model=list(ar=ar4.real))
ar4.1024 <- arima.sim(1024,model=list(ar=ar4.real))

# get estimates of coefs using different methods:

ar4.64.YW <- lev.dur(ar4.64,4)
ar4.64.YWtap <- lev.dur.tap(ar4.64,4)
ar4.64.Fls <- Fls(ar4.64,4)
ar4.64.Bls <- Bls(ar4.64,4)
ar4.64.FBls <- FBls(ar4.64,4)

# best way to compare is to look at theoretical spectra:

spec.64.YW <- spec.ar(ar4.64.YW)
spec.64.YWtap <- spec.ar(ar4.64.YWtap)
spec.64.Fls  <- spec.ar(ar4.64.Fls)
spec.64.Bls  <- spec.ar(ar4.64.Bls)
spec.64.FBls  <- spec.ar(ar4.64.FBls)



spec.true_spec.ar(list(ar=ar4.real,var.pred=1,order=4))
spec.range <- range(spec.64.YW$spec,spec.64.Fls$spec)
spec.range <- range(spec.range,spec.64.Bls$spec,spec.64.FBls$spec)
spec.range <- range(spec.range,spec.64.YWtap$spec,spec.true$spec)

plot(spec.true$freq,spec.true$spec,ylim=range(spec.range),type="l")
lines(spec.64.YW$freq,spec.64.YW$spec,lty=2)
title("comparison of true and YW spectra")

plot(spec.true$freq,spec.true$spec,ylim=range(spec.range),type="l")
lines(spec.64.YWtap$freq,spec.64.YWtap$spec,lty=2)
title("comparison of true and tapered YW spectra")

plot(spec.true$freq,spec.true$spec,ylim=range(spec.range),type="l")
lines(spec.64.Fls$freq,spec.64.Fls$spec,lty=2)
title("comparison of true and Forward LS spectra")

plot(spec.true$freq,spec.true$spec,ylim=range(spec.range),type="l")
lines(spec.64.Bls$freq,spec.64.Bls$spec,lty=2)
title("comparison of true and Backward LS spectra")

plot(spec.true$freq,spec.true$spec,ylim=range(spec.range),type="l")
lines(spec.64.FBls$freq,spec.64.FBls$spec,lty=2)
title("comparison of true and Forward/Backward LS spectra")


# fit some YW AR(8) models to the data sets and plot
# note we are over parameterising as the true model is AR(4)
# produce similar plots for other methods
YW8.plot(ar4.64, ar4.256, ar4.1024, ar4.real)


# Do some predictions
# First fit the models excluding npred points 
# (where npred is the number of points you which to predict)

ar.data <- ar4.64
npred <- 10
n <- 64

ar.fit<- lev.dur(ar.data[1:(n-npred)],4)
ar.pred <- arima.forecast(ar.data[1:(n-npred)], n=npred,model=list(order=c(ar.fit$order,0,0),ar=ar.fit$ar,ndiff=0))

# plot the predictions

plot((n-npred+1):n,ar.data[(n-npred+1):n], ylim=range(ar.data,ar.pred$mean), type="l", xlab="Time", ylab="y")
lines((n-npred+1):n,ar.pred$mean,lty=2)
legend(n-npred+1,max(ar.data), legend=c("true","forecast"), lty=1:2, bty="n")


# If you want to plot 10 points (for example) before the forecast horizon:

plot((n-npred-10+1):n,ar.data[(n-npred-10+1):n], ylim=range(ar.data,ar.pred$mean), type="l", xlab="Time", ylab="y")
lines((n-npred+1):n,ar.pred$mean,lty=2)
legend(n-npred-10+1,max(ar.data), legend=c("true","forecast"), lty=1:2, bty="n")

#  To compare different models (fit model orders 1 to 20 and keep a track of 
# the residual sum of squares)

res.sqq <- vector("numeric")
for(p in 1:20){
 ar.fit<- lev.dur(ar.data[1:(n-npred)],p)
 ar.pred <- arima.forecast(ar.data[1:(n-npred)], n=npred,model=list(order=c(ar.fit$order,0,0),ar=ar.fit$ar,ndiff=0))
 res.sqq[p] <- sum((ar.data[(n-npred+1):n] - ar.pred$mean)^2)
}

# find the value of p which minimises the residual sum of squares

final.p <- (1:20)[res.sqq == min(res.sqq)]

plot(1:20,res.sqq,type="l", ylas="RSS",xlab="p")
abline(v=final.p)

# fit and plot the "best" p

ar.fit<- lev.dur(ar.data[1:(n-npred)],final.p)
ar.pred <- arima.forecast(ar.data[1:(n-npred)], n=npred,model=list(order=c(ar.fit$order,0,0),ar=ar.fit$ar,ndiff=0))

plot((n-npred+1):n,ar.data[(n-npred+1):n], ylim=range(ar.data,ar.pred$mean), type="l", xlab="Time", ylab="y")
lines((n-npred+1):n,ar.pred$mean,lty=2)
legend(n-npred+1,max(ar.data), legend=c("true","forecast"), lty=1:2, bty="n")
title(paste("Fitted model is AR(",format(final.p),") \n Chosen to minimise RSS", sep=""))