library(ts)

load("serie9.RData")

##1. SARIMA-Modelle

##number of airline passeneges
##monthly number of departing airline passengers in the U.S. for the years 1949-1960.  
 
airl <- ts(scan("http://stat.ethz.ch/Teaching/Datasets/airline.dat"))

plot(airl)

#differencing
plot(diff.ts(airl,lag=1,differences=1))

#--> not successful!

#with seasonal differences (seasonal period 12)
plot(diff.ts(diff.ts(airl,lag=12,differences=2),lag=1,differences=1))

xdiff <- diff.ts(diff.ts(airl,lag=12,differences=2),lag=1,differences=1)

par(mfrow=c(1,2))
acf(xdiff)
acf(xdiff,type="partial")

aic.res <- matrix(0,ncol=3,nrow=5)
for(i in 1:5)
  {
    for(j in 1:2)
      {
        aic.res[i,j] <- arima(airl,order=c(j,1,0),seasonal=list(order=c(i,2,0),period=12))$aic
        print(c(i,j,aic.res[i,j]))
        
      }
  }

aic.res

#best value for SARIMA
fit <- arima(airl,order=c(1,1,0),seasonal=list(order=c(5,2,0),period=12))
#Fit fuer SARIMA((1,1,0),(5,2,0)_12)

par(mfrow=c(1,2))
acf(fit$resid)
acf(fit$resid,type="partial")

par(mfrow=c(2,2))
acf(fit$resid,ylim=c(-0.5,1))
acf(fit$resid,type="partial",ylim=c(-0.5,0.3))
acf(xdiff,ylim=c(-0.5,1))
acf(xdiff,type="partial",ylim=c(-0.5,0.3))

par(mfrow=c(1,1))
qqnorm(fit$resid)

#improvement
fit2 <- arima(airl,order=c(1,1,1),seasonal=list(order=c(1,2,1),period=12))
#fit for SARIMA((1,1,1),(1,2,1)_12)

fit2$aic
[1] 948.3984

fit$aic
[1] 951.7593

par(mfrow=c(1,2))
acf(fit2$resid)
acf(fit2$resid,type="partial")

par(mfrow=c(1,1))
qqnorm(fit2$resid)

lag.plot(fit2$resid,set.lags=1:4)

#prediction
pred <- predict(fit2,n.ahead=12,se=T)

#data, prediction, prediction interval
plot(1:156,1:156,type="n",xlab="months",ylab="number of passengers",ylim=c(100,690))
lines(1:144,airl)
lines(145:156,pred$pred,col=2)
lines(145:156,pred$pred + 1.96 * pred$se,col=3)
lines(145:156,pred$pred - 1.96 * pred$se,col=3)


##2. with log-Transformation
plot(log(airl))

#differencing
plot(diff.ts(log(airl),lag=1,differences=1))

plot(diff.ts(diff.ts(log(airl),lag=12,differences=1),lag=1,differences=1))

xdiff <- diff.ts(diff.ts(log(airl),lag=12,differences=1),lag=1,differences=1)

par(mfrow=c(1,2))
acf(xdiff)
acf(xdiff,type="partial")

aic.res2 <- matrix(0,ncol=4,nrow=4)
for(i in 1:4)
  {
    for(j in 1:4)
      {
        aic.res2[i,j] <- arima(log(airl),order=c(j,1,0),seasonal=list(order=c(i,1,0),period=12))$aic
        print(c(i,j,aic.res2[i,j]))
        
      }
  }

aic.res2

#improvement
fit3 <- arima(log(airl),order=c(2,1,1),seasonal=list(order=c(1,1,1),period=12))
#Fit for SARIMA((2,1,1),(1,1,1)_12)

fit3$aic
[1] -480.4212

par(mfrow=c(1,2))
acf(fit3$resid)
acf(fit3$resid,type="partial")

par(mfrow=c(2,2))
acf(fit3$resid,ylim=c(-0.5,1))
acf(fit3$resid,type="partial",ylim=c(-0.5,0.3))
acf(xdiff,ylim=c(-0.5,1))
acf(xdiff,type="partial",ylim=c(-0.5,0.3))


par(mfrow=c(1,1))
qqnorm(fit3$resid)

lag.plot(fit3$resid,set.lags=1:4)


#prediction
pred3 <- predict(fit3,n.ahead=12,se=T)

#data, prediction, prediction interval
par(mfrow=c(1,1))
plot(1:156,1:156,type="n",xlab="months",ylab="number of passengers",ylim=c(100,690))
lines(1:144,airl)
lines(145:156,exp(pred3$pred),col=2)
lines(145:156,exp(pred3$pred + 1.96 * pred3$se),col=3)
lines(145:156,exp(pred3$pred - 1.96 * pred3$se),col=3)

#comparision with and without log-transformation
par(mfrow=c(2,1))
plot(1:156,1:156,type="n",xlab="months",ylab="number of passengers",ylim=c(100,690))
lines(1:144,airl)
lines(145:156,pred$pred,col=2)
lines(145:156,pred$pred + 1.96 * pred$se,col=3)
lines(145:156,pred$pred - 1.96 * pred$se,col=3)

plot(1:156,1:156,type="n",xlab="months",ylab="number of passengers",ylim=c(100,690))
lines(1:144,airl)
lines(145:156,exp(pred3$pred),col=2)
lines(145:156,exp(pred3$pred + 1.96 * pred3$se),col=3)
lines(145:156,exp(pred3$pred - 1.96 * pred3$se),col=3)

par(mfrow=c(1,1))

plot(pred$pred,ylim=c(350,750))
lines(pred$pred + 1.96 * pred$se,col=3)
lines(pred$pred - 1.96 * pred$se,col=3)

lines(exp(pred3$pred),col=2)
lines(exp(pred3$pred + 1.96 * pred3$se),col=4)
lines(exp(pred3$pred - 1.96 * pred3$se),col=4)

legend(145,730,c("without log","without log","log","log"),lty=rep(1,4),col=c(1,3,2,4))

save(airl,file="serie9.RData")



library(ts)

##2. spectral density of ARMA processes

spec.arma <- function(ar = NULL, ma = NULL, sigma = 1.0)
  {
      ## Purpose: plot of spectral density (on log-scle) for ARMA processes
      ## -------------------------------------------------------------------------
      ## Arguments: AR- and MA-coefficients, sigma
      ## -------------------------------------------------------------------------
      ## Author: Peter Buehlmann, Date: 18 Dec 2000, 11:56
    lambda <- seq(0,pi,0.01)
    Phi <- rep(1,length(lambda))
    Theta <- rep(1,length(lambda))
    if(length(ar)) {
        for(j in 1:length(ar)){
              Phi <- Phi - ar[j]*exp(-(1i*lambda))^j
            }
      }
    if(length(ma)) {
        for(j in 1:length(ma)){
              Theta <- Theta + ma[j]*exp(-(1i*lambda))^j
            }
      }
    plot(lambda,log((sigma*abs(Theta)/abs(Phi))^2/(2*pi)),type="l",main="ARMA spectrum",xlab="frequency",ylab="log(spectral density)")
  }

##Spektraldichte von ARMA Prozess: variiere die Parameter
spec.arma(ar=c(0.99),ma=seq(0.5,0.3))
spec.arma(ar=c(0.8,-0.5))

##sunspot data
data(sunspots)

par(mfrow=c(1,1))
ts.plot(sunspots)

fit <- ar.yw(ts(sunspots))
fit

#AR-Spektrum fuer sunspot Daten
spec.ar(fit)

#roots for sunspot data
#Nullstellen von AR- und MA-Polynom

ns.arma <- function(ar = NULL,ma = NULL)
  {
      ## Purpose: roots of AR- and MA-polynomial
      ## -------------------------------------------------------------------------
      ## Arguments: AR- and MA-coefficients
      ## -------------------------------------------------------------------------
      ## Author: Peter Buehlmann, Date: 18 Dec 2000, 11:05
    list(ar=abs(polyroot(c(1,- ar))), ma=abs(polyroot(c(1,ma))))
  }

ns.arma(ar=fit$ar)

#Resultat:
> $ar
 [1] 1.089403 1.099563 1.081002 1.015015 1.091106 1.081002 1.080131 1.015015
 [9] 1.116835 1.080267 1.095476 1.159490 1.080131 1.080267 1.091106 1.089745
[17] 1.095476 1.104369 1.116835 1.115337 1.104369 1.099563 1.115337 1.089745
[25] 1.089403 1.159490 1.854561 1.854561

$ma
numeric(0)

#plot of the roots
plot(polyroot(c(1,-fit$ar)),ylim=c(-1.4,1.4))
points(exp(-(1i)*seq(0,2*pi,0.01)),type="l")

#frequencies where the first two roots are approximately equal to exp(i \lambda)
#with frequencies at index [4],[8]

freq <- Arg(polyroot(c(1,-fit$ar)))[c(4,8)]

spec.ar(fit)
lines(rep(freq[2]/(2*pi),49901),seq(100,50000,1),col=2)


#quasi-period:
2*pi/freq[2]
[1] 128.4377
#approx. 10.7 years