## (I) difference-stationary noise
X <- read.csv("airlines.csv",header=FALSE)
X <- as.numeric(X[,1])
plot(X)

plot(log(X),type="l")
plot(Z <- diff(log(X),12),type="l")

## subtract mean
mu <- mean(Z)
Z <- Z-mu

fitted <- ar.yw( Z,order.max=5)
## addn new observations at the end 
addn <- 30
newZ <- filter( rnorm(addn)*sqrt(fitted$var.pred), filter= fitted$ar, method="recursive", init=Z[length(Z)-((length(fitted$ar)-1):0)]) 

## make variance and mean equal to observations
## combina observations and simulation
allZ <- c(Z,newZ) + mu 
## plot new time series
plot(allZ,type="l")
abline(v=length(Z))

Extended <- diffinv( allZ, 12, xi=log(X)[1:12])
## plot simulated and observed time-series in original space
plot( Extended)
abline(v=length(Z))

## repeat a 100 times and plot all simulations
for (sim in 1:20){

    fitted <- ar.yw( Z,order.max=5)
    addn <- 50
    newZ <- filter( rnorm(addn)*sqrt(fitted$var.pred), filter= fitted$ar, method="recursive", init=Z[length(Z)-((length(fitted$ar)-1):0)])
    allZ <- c(Z,newZ) +mu 

    Ext <- diffinv( allZ, 12, xi=log(X)[1:12])
    if(sim==1) plot( exp(Ext),type="l", col=rgb(0.2,0.2,0.2,0.2),lwd=1,log="y") else lines( exp(Ext), col=rgb(0.2,0.2,0.2,0.2),lwd=1)
}



## (II) seasonal / stationary noise
X <- read.csv("airlines.csv",header=FALSE)
X <- as.numeric(X[,1])
plot(X)


plot(log(X),type="l")
time <- 1:length(X)
months <- rep(1:12, ceiling(length(time)/12))[1:length(time)]
## first remove linear time trend, then monthly effect (better: remove jointly)
Z <- residuals( ll <- lm(log(X) ~ time))
monthsmean <- aggregate( Z, by=list(mo=months),FUN="mean")
Z <- Z - monthsmean[ match(months,monthsmean[,"mo"]) ,2]
plot(Z ,type="l")

fitted <- ar.yw( Z,order.max=5)
addn <- 30
newZ <- filter( rnorm(addn)*sqrt(fitted$var.pred), filter= fitted$ar, method="recursive", init=Z[length(Z)-((length(fitted$ar)-1):0)])
allZ <- c(Z ,newZ)
plot(allZ,type="l")
abline(v=length(Z))
allZ <- allZ + monthsmean[match(rep(1:12,ceiling(length(allZ)/12))[1:length(allZ)],monthsmean[,"mo"]),2] + coef(ll)[1] + coef(ll)[2]* as.numeric(1:length(allZ))
plot(allZ)

for (sim in 1:20){

    fitted <- ar.yw( Z,order.max=5)
    addn <- 50
    newZ <- filter( rnorm(addn)*sqrt(fitted$var.pred), filter= fitted$ar, method="recursive", init=Z[length(Z)-((length(fitted$ar)-1):0)])
    allZ <- c(Z,newZ)
    newmonths <- rep(1:12, ceiling( (length(allZ))/12))[1:length(allZ)]
    allZ <- allZ + monthsmean[match(newmonths,monthsmean[,"mo"]),2]

    allZ <- allZ + coef(ll)[1] + coef(ll)[2]* as.numeric(1:length(allZ))
    Ext <- (allZ)
    if(sim==1) plot( exp(Ext),type="l", col=rgb(0.2,0.2,0.2,0.2),lwd=1,log="y") else lines( exp(Ext), col=rgb(0.2,0.2,0.2,0.2),lwd=1)
    
}