# chap11.R

# Exhibit 11.1

win.graph(width=4.875, height=2.5,pointsize=8)
data(airmiles)
plot(log(airmiles),ylab='Log(airmiles)',xlab='Year', )

# Exhibit 11.5
acf(as.vector(diff(diff(window(log(airmiles),end=c(2001,8)),12))),lag.max=48)

# Exhibit 11.6 
air.m1=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1),
period=12),xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)),xreg=data.frame(Dec96=1*(seq(airmiles)==12),
Jan97=1*(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84)),method='ML')
# Additive outliers are incorporated as dummy variables in xreg.
# Transfer function components are incorporated by the xtransf and transfer
# arguments.
# Here, the transfer function consists of two parts omega0*P(t) and 
# omega1/(1-omega2*B)P(t) where the inputs of the two transfer
# functions are identical and equals the dummy variable that is 1 at September
# 2001 (the 69th data point) and zero otherwise.
# xtransf is a matrix whose columns are the input variables.
# transfer is a list consisting of the pair of (MA order, AR order) of each
# transfer function, which in this examples is (0,0) and (1,0).

air.m1

# Exhibit 11.7
plot(log(airmiles),ylab="log(airmiles)")
points(fitted(air.m1))


# Exhibit 11.8
Nine11p=1*(seq(airmiles)==69)
plot(ts(Nine11p*(-0.0949)+
filter(Nine11p,filter=.8139,method='recursive',side=1)*(-0.2715),
frequency=12,start=1996),type='h',ylab='9/11 Effects')
abline(h=0)

# Exhibit 11.9 
set.seed(12345)
y=arima.sim(model=list(ar=.8,ma=.5),n.start=158,n=100)
y[10]
y[10]=10
y=ts(y,freq=1,start=1)
plot(y,type='o')
acf(y)
pacf(y)
eacf(y)
m1=arima(y,order=c(1,0,0))
m1
detectAO(m1)
detectAO(m1, robust=F)
detectIO(m1)
m2=arima(y,order=c(1,0,0),xreg=data.frame(AO=seq(y)==10))
m2
detectAO(m2)
detectIO(m2)
tsdiag(m2)
tsdiag(m1)
m3=arima(y,order=c(1,0,1),xreg=data.frame(AO=seq(y)==10))
detectAO(m3)
detectIO(m3)
tsdiag(m3)
m3

plot(y,type='b') 
arrows(30,6, 11,9.8, length=.1,angle=20)
text(34,6.2, "AO")


# Exhibit 11.10
data(co2)
m1.co2=arima(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
m1.co2
detectAO(m1.co2)
detectIO(m1.co2)
m4.co2=arimax(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),
io=c(57))
m4.co2

# the signs of the MA coefficient estimates appear to be opposite to 
# those in the book, because R uses the plus convention in the MA parameterization.

# Exhibit 11.11
win.graph(width=4.875, height=2.5,pointsize=8)
set.seed(12345)
X=rnorm(105)
Y=zlag(X,2)+.5*rnorm(105)
X=ts(X[-(1:5)],start=1,freq=1)
Y=ts(Y[-(1:5)],start=1,freq=1)
ccf(X,Y,ylab='CCF')


# Exhibit 11.12
phi=seq(0,.95,.15)
rejection=2*(1-pnorm(1.96*sqrt((1-phi^2)/(1+phi^2))))
M=signif(rbind(phi,rejection),2)
rownames(M)=c("phi", "Error Rate")
M


# Exhibit 11.13
set.seed(23457)
correlation.v=NULL
B=1000
n=500
for (i in 1:B) {x=cumsum(arima.sim(model=list(ma=.8),n=n))
y=cumsum(arima.sim(model=list(ma=.8),n=n))
correlation.v=c(correlation.v,ccf(x,y,lag.max=1,plot=F)$acf[2])
}
hist(correlation.v,prob=T,xlab=expression(r[0](X,Y)))

# Exhibit 11.14
data(milk)
data(electricity)
milk.electricity=ts.intersect(milk,log(electricity))
# The ts.intersect function merges serveral time series into a panel of time
# series over the time frame where each series has data.

plot(milk.electricity,yax.flip=T) 
# the option yax.flip=T flips the y-axis for the series alternately so as
# to make the labeling clearer.

# Exhibit 11.15
ccf(as.numeric(milk.electricity[,1]),as.numeric(milk.electricity[,2]),
main='milk & electricity',ylab='CCF')
# The as.numeric function strips the time series attribute from the time series.
# This is done to nullify the default way the ccf function plots the cross-correlations.
# You may want to repeat the command without the as.numeric function to see 
# the default labels of the lags according to the period of the data.
# ccf((milk.electricity[,1]),(milk.electricity[,2]), main='milk & electricity',ylab='CCF')


# Exhibit 11.16
me.dif=ts.intersect(diff(diff(milk,12)),diff(diff(log(electricity),12)))
prewhiten(as.numeric(me.dif[,1]),as.numeric(me.dif[,2]),
,ylab='CCF' )

# Exhibit 11.17
data(bluebird)
plot(bluebird,yax.flip=T)

# Exhibit 11.18
prewhiten(y=diff(bluebird)[,1],x=diff(bluebird)[,2],main="Price & Sales",ylab='CCF')
# As the time series are of unit period, there is no need to apply the as.numeric 
# function.

# Exhibit 11.19
sales=bluebird[,1]
price=bluebird[,2]
chip.m1=lm(sales~price,data=bluebird)
summary(chip.m1)

# Exhibit 11.20
acf(residuals(chip.m1),ci.type='ma')

# Exhibit 11.21
pacf(residuals(chip.m1))

# Exhibit 11.22
eacf(residuals(chip.m1))

# Exhibit 11.23
chip.m2=arima(sales,order=c(1,0,4),xreg=data.frame(price))
chip.m2
# MA(1)& MA(3) estimates are not significant, at 5% level,
# so they are constrained to be zero in the model fit below.
chip.m3=arima(sales,order=c(1,0,4),xreg=data.frame(price),fixed=c(NA,0,NA,0,NA,NA,NA))
# The MA(1) & MA(3) estimates are the second and fourth coefficients listed
# in the model fit chip.m2. They are set to be zero using the fixed option.
# NAs in the fixed option correspond to free parameters. 
chip.m3
# Now, the AR(1) coefficient estimate also seems not significant, so it is
# removed in the next fitted model.
chip.m4=arima(sales,order=c(0,0,4),xreg=data.frame(price),fixed=c(0,NA,0,NA,NA,NA))
chip.m4

# model diagnostic can be done by running the tsdiag command.
tsdiag(chip.m4)


# Exhibit 11.24
data(boardings)
plot(boardings,yax.flip=T)
# The Denver dataset has three time series. Here, we only plot the first 
# two series.

# Exhibit 11.25
m1=arima(boardings[,2],order=c(2,1,0))
prewhiten(x=boardings[,2],y=boardings[,1],x.model=m1)


# Exhibit 11.26
log.boardings=boardings[,1]
log.price=boardings[,2]
boardings.m1=arima(log.boardings,order=c(1,0,0),seasonal=list(order=c(1,0,0),period=12),
xreg=data.frame(log.price))
boardings.m1
detectAO(boardings.m1) 
detectIO(boardings.m1)
# An AO is detected at time point 32, as well as an IO detected at time point 44
# Since the test statistic of the AO has larger magnitude, an AO is added to the
# model fitted below.
boardings.m2=arima(log.boardings,order=c(1,0,3),seasonal=list(order=c(1,0,0),period=12),
xreg=data.frame(log.price,outlier=c(rep(0,31),1,rep(0,36))),
fixed=c(NA,0,0,rep(NA,5)))
boardings.m2
detectAO(boardings.m2)
detectIO(boardings.m2)
# No outliers are detected!
tsdiag(boardings.m2,tol=.15,gof.lag=24)
# Model diagnostics appear to be OK.