# Chap6.R

# Exhibit 6.5
data(ma1.1.s)
win.graph(width=4.875, height=3,pointsize=8)
acf(ma1.1.s,xaxp=c(0,20,10)) # if the arument xaxp is omitted, the
# tick marks will be generated according to the default convention. Run the command ?par to learn more the xaxp argment.  
acf(ma1.1.s) # see the new tickmarks

# So how was xaxp determined? First run the command without xaxp to find out how many lags are on the x-axis. Suppose there are 20 lags. Then we specify xaxp=c(0,20,10), i.e. the two extreme tickmarks are 0 and 20, and we want 10 tickmarks in between. How to set the xaxp argument if we want 1 tickmark for each lag?

# Exhibit 6.6
acf(ma1.1.s,ci.type='ma',xaxp=c(0,20,10))

# Exhibit 6.7
data(ma1.2.s)
acf(ma1.2.s,xaxp=c(0,20,10))

# Exhibit 6.8
data(ma2.s)
acf(ma2.s,xaxp=c(0,20,10))

# Exhibit 6.9
acf(ma2.s,ci.type='ma',xaxp=c(0,20,10))

# Exhibit 6.10
data(ar1.s)
acf(ar1.s,xaxp=c(0,20,10))


# Exhibit 6.11
pacf(ar1.s,xaxp=c(0,20,10))

# Exhibit 6.12
data(ar2.s)
acf(ar2.s,xaxp=c(0,20,10))

# Exhibit 6.13
pacf(ar2.s,xaxp=c(0,20,10))

# Exhibit 6.14
data(arma11.s)
plot(arma11.s, type='b',ylab=expression(y[t]))

# Exhibit 6.15
acf(arma11.s,xaxp=c(0,20,10))

# Exhibit 6.16
pacf(arma11.s,xaxp=c(0,20,10))

# Exhibit 6.17
eacf(arma11.s)

# If desirable,  a title can be added to the EACF table by the following command.

cat('Exhibit 6.17\n');eacf(arma11.s)

# Exhibit 6.18
data(oil.price)
acf(as.vector(oil.price),main='Sample ACF of the Oil Price Time Series',xaxp=c(0,24,12))
# The as.vector function strips the ts attaribute of oil.price in order to prevent the plotting convention for seasonal time series. The main argument supplies the main heading of the figure.Try the following command to appreciate the effects of applying the as.vector function.

acf(oil.price,main='Sample ACF of the Oil Price Time Series')
# The tickmark 1 on the x-axis now refers to a lag of 1 period with each period containing 12 months, so it is lag 12!

# Exhibit 6.19
acf(diff(as.vector(log(oil.price))),main='Sample ACF of the Difference of the Oil Price Time Series',xaxp=c(0,24,12))

# Exhibit 6.20
data(rwalk)
acf(diff(rwalk,difference=2),ci.type='ma',xaxp=c(0,18,9))

# Exhibit 6.21
acf(diff(rwalk),ci.type='ma',xaxp=c(0,18,9))

library(uroot)
# Carry out the Dickey-Fuller unit root tests, Find out the AR order for the differenced series
ar(diff(rwalk)) # order 8 is indicated by AIC

library(uroot)
ADF.test(rwalk,selectlags=list(mode=c(1,2,3,4,5,6,7,8),Pmax=8),itsd=c(1,0,0))
# which is equivalent to ADF.test(rwalk,selectlags=list(mode=c(1,2,3,4,5,6,7,8)),itsd=c(1,0,0)). The argument mode can be either a numerical vector, or "aic","bic" and "signf". If it is a numerical vector, the test assumes the vector corresponds to the lags of the response that must be included in the test, and the Pmax option is ignored.However, if mode is specified as "aic", the function uses the AIC to select lags up to and including Pmax to be included in the test. Specifying mode to "bic" or "signf" changes the selection criterion to BIC or significance criterion. Note "signf" not "signif" is the valid option.

# In comparison, setting the true order to be zero for the first difference ADF.test(rwalk,selectlags=list(Pmax=0),itsd=c(1,0,0)). Repeat the test with the alternative have an intercept term and a linear trend Run ?uroot to learn more about the selectlags option and the itsd option.

ADF.test(rwalk,selectlags=list(mode=c(1,2,3,4,5,6,7,8),Pmax=8),itsd=c(1,1,0))
# Similarly, a simplified command is ADF.test(rwalk,selectlags=list(mode=c(1,2,3,4,5,6,7,8)),itsd=c(1,1,0)). Repeat the preceding test but now setting the true order of the first difference as zero.
ADF.test(rwalk,selectlags=list(Pmax=0),itsd=c(1,1,0))


# Exhibit 6.22
set.seed(92397)
test=arima.sim(model=list(ar=c(rep(0,11),.8),ma=c(rep(0,11),0.7)),n=120)
res=armasubsets(y=test,nar=14,nma=14,y.name='test',ar.method='ols')
plot(res)
# A title may be added to the plot, by adding the main="..." option in the above plot function. However, the title may crash with other labels. A better approach is to add a title by the following command.
title(main="Exhibit 6.22", line=6)
# The option line=6 means put the title 6 lines from the edge of the plot region. You can experiment with this option to fine tune the placement of the label.

# Exhibit 6.23
data(larain)
win.graph(width=3, height=3,pointsize=8)
qqnorm(log(larain))
qqline(log(larain))

# Exhibit 6.24
win.graph(width=4.875, height=3,pointsize=8)
acf(log(larain),xaxp=c(0,20,10)) # note the main heading now includes the data name.

# Exhibit 6.25
data(color)
acf(color,ci.type='ma')

# Exhibit 6.26
pacf(color)

# Exhibit 6.27  
win.graph(width=4, height=4,pointsize=8)
data(hare)
bxh=BoxCox.ar(hare)
bxh$mle # the mle of the power parameter
bxh$ci # corresponding 95% C.I.

# Exhibit 6.28
acf(hare^.5)

# Exhibit 6.29
pacf(hare^.5)

# Exhibit 6.30
eacf(diff(log(oil.price)))

# Exhibit 6.31
res=armasubsets(y=diff(log(oil.price)),nar=7,nma=7,
y.name='test', ar.method='ols')
plot(res)

# Exhibit 6.32
acf(as.vector(diff(log(oil.price))),xaxp=c(0,22,11))

# Exhibit 6.33
pacf(as.vector(diff(log(oil.price))),xaxp=c(0,22,11))

# Chap7.R

# Below is a function that computes the method of moments estimator of
# the MA(1) coefficient of an MA(1) model.
estimate.ma1.mom=function(x){r=acf(x,plot=F)$acf[1]; if (abs(r)<0.5)
return((-1+sqrt(1-4*r^2))/(2*r)) else return(NA)}

# Exhibit 7.1
data(ma1.2.s)
estimate.ma1.mom(ma1.2.s)

data(ma1.1.s)
estimate.ma1.mom(ma1.1.s)

set.seed(1234)
ma1.3.s=arima.sim(list(ma=c(.9)),n=60)
estimate.ma1.mom(ma1.3.s)

ma1.4.s=arima.sim(list(ma=c(-0.5)),n=60)
estimate.ma1.mom(ma1.4.s)

arima(ma1.4.s,order=c(0,0,1),method='CSS',include.mean=F)
data(ar1.s)
ar(ar1.s,order.max=1,AIC=F,method='yw')
data(ar1.2.s)
ar(ar1.2.s,order.max=1,AIC=F,method='yw')
data(ar2.s)
ar(ar2.s,order.max=2,AIC=F,method='yw')


# Exhibit 7.4
data(ar1.s)
ar(ar1.s,order.max=1,AIC=F,method='yw') # method of moments
ar(ar1.s,order.max=1,AIC=F,method='ols') # conditional sum of squares
ar(ar1.s,order.max=1,AIC=F,method='mle') # maximum likelihood
# The AIC option is set to be False otherwise the function will choose
# the AR order by minimizing AIC, so that zero order might be chosen.

data(ar1.2.s)
ar(ar1.2.s,order.max=1,AIC=F,method='yw') # method of moments
ar(ar1.2.s,order.max=1,AIC=F,method='ols') # conditional sum of squares
ar(ar1.2.s,order.max=1,AIC=F,method='mle') # maximum likelihood

# Exhibit 7.5
data(ar2.s)
ar(ar2.s,order.max=2,AIC=F,method='yw') # method of moments
ar(ar2.s,order.max=2,AIC=F,method='ols') # conditional sum of squares
ar(ar2.s,order.max=2,AIC=F,method='mle') # maximum likelihood

# Exhibit 7.6
data(arma11.s)
arima(arma11.s, order=c(1,0,1),method='CSS') # conditional sum of squares
arima(arma11.s, order=c(1,0,1),method='ML') # maximum likelihood
#
# Recall that R uses the plus convention whereas our book uses the minus
# convention in the specification of the MA part, i.e. R specifies an
# ARMA(1,1) model as z_t=theta_0+phi*z_{t-1}+e_t+theta_1*e_{t-1}
# versus our convention
# z_t=theta_0+phi*z_{t-1}+e_t-theta_1*e_{t-1}

# Exhibit 7.7
data(color)
ar(color,order.max=1,AIC=F,method='yw') # method of moments
ar(color,order.max=1,AIC=F,method='ols') # conditional sum of squares
ar(color,order.max=1,AIC=F,method='mle') # maximum likelihood


# Exhibit 7.8
data(hare)
arima(sqrt(hare),order=c(3,0,0))

# Exhibit 7.9
data(oil.price)
arima(log(oil.price),order=c(0,1,1),method='CSS') # conditional sum of squares
arima(log(oil.price),order=c(0,1,1),method='ML') # maximum likelihood

# Exhibit 7.10

res=arima(sqrt(hare),order=c(3,0,0),include.mean=T)
set.seed(12345)
coefm.cond.norm=arima.boot(res,cond.boot=T,is.normal=T,B=1000,init=sqrt(hare))
signif(apply(coefm.cond.norm,2,function(x){quantile(x,c(.025,.975),na.rm=T)}),3)
# Method I

coefm.cond.replace=arima.boot(res,cond.boot=T,is.normal=F,B=1000,init=sqrt(hare))
signif(apply(coefm.cond.replace,2,function(x){quantile(x,c(.025,.975),na.rm=T)}),3)
# Method II

coefm.norm=arima.boot(res,cond.boot=F,is.normal=T,ntrans=100,B=1000,init=sqrt(hare))
signif(apply(coefm.norm,2,function(x){quantile(x,c(.025,.975),na.rm=T)}),3)
# Method III

coefm.replace=arima.boot(res,cond.boot=F,is.normal=F,ntrans=100,B=1000,init=sqrt(hare))
signif(apply(coefm.replace,2,function(x){quantile(x,c(.025,.975),na.rm=T)}),3)
# Method IV

# Some bootstrap series may be explosive which will be discarded. To see
# the number of usable bootstrap series, run the command

dim(coefm.replace)
# the output should be
# [1] 952   5
# i.e. we have only 952 usable (i.e. finite) bootstrap time series even though
# we simulate 1000 series.

# The theoretical confidence intervals were computed by the output in Exhibit 7.8.


# Compute the quasi-period of the bootstrap series based on the method of
# stationary bootstrap with the errors drawn from the residuals with replacement.

period.replace=apply(coefm.replace,1,function(x){
roots=polyroot(c(1,-x[1:3]))
# find the complex root with smalles magnitude
min1=1.e+9
rootc=NA
for (root in roots) {
if( abs(Im(root))<1e-10) next
if (Mod(root)< min1) {min1=Mod(root); rootc=root}
}
if(is.na(rootc)) period=NA else period=2*pi/abs(Arg(rootc))
period
})

sum(is.na(period.replace)) # number of bootstap series that do not admit a  well-defined quasi-period.

quantile(period.replace, c(.025,.975),na.rm=T)

# Exhibit 7.11
win.graph(width=3.9,height=3.8,pointsize=8)
hist(period.replace,prob=T,main="",xlab="quasi-period",axes=F,xlim=c(5,16))
axis(2)
axis(1,c(4,6,8,10,12,14,16),c(4,6,8,10,12,14,NA))


# Exhibit 7.12

win.graph(width=3,height=3,pointsize=8)
qqnorm(period.replace,main="") #Normal Q-Q Plot for the Bootstrap Quasi-period Estimates")
qqline(period.replace)

# chap8.R

# Exhibit 8.1
win.graph(width=4.875, height=3,pointsize=8)
data(color)
m1.color=arima(color,order=c(1,0,0))
m1.color
plot(rstandard(m1.color),ylab='Standardized residuals',type='b')
abline(h=0)

# Exhibit 8.2
data(hare)
m1.hare=arima(sqrt(hare),order=c(3,0,0))
m1.hare # the AR(2) coefficient is not significant; it is second in the list of coefficients.
m2.hare=arima(sqrt(hare),order=c(3,0,0),fixed=c(NA,0,NA,NA)) # fixed the AR(2) coefficient to be 0 via the fixed argument.
m2.hare
# Note that the intercept term is actually the mean in the centered form of the ARMA model, i.e. if y(t)=sqrt(hare)-intercept, then the model is y(t)=0.919*y(t-1)-0.5313*y(t-3)+e(t), So the "ture" intercept equals 5.6889*(1-0.919+0.5313)=3.483, as stated in the book!
plot(rstandard(m2.hare),ylab='Standardized residuals',type='b')
abline(h=0)

# Exhibit 8.3
data(oil.price)
m1.oil=arima(log(oil.price),order=c(0,1,1))
plot(rstandard(m1.oil),ylab='Standardized residuals',type='l')
abline(h=0)

# Exhibit 8.4
win.graph(width=3, height=3,pointsize=8)
qqnorm(residuals(m1.color))
qqline(residuals(m1.color))

# Exhibit 8.5
qqnorm(residuals(m1.hare))
qqline(residuals(m1.hare))

# Exhibit 8.6
qqnorm(residuals(m1.oil))
qqline(residuals(m1.oil))

# Exhibit 8.9
win.graph(width=4.875, height=3,pointsize=8)
acf(residuals(m1.color),main='Sample ACF of Residuals from AR(1) Model for Color')

# Exhibit 8.10
acf(residuals(arima(sqrt(hare),order=c(2,0,0))),main='Sample ACF of Residuals from AR(2) Model for Hare')

# Exhibit 8.11
acf(residuals(m1.color),plot=F)$acf
signif(acf(residuals(m1.color),plot=F)$acf[1:6],2)# to display the first 6 acf to 2 significant digits.

# Exhibit 8.12
win.graph(width=4.875, height=4.5)
tsdiag(m1.color,gof=15,omit.initial=F) # the tsdiag function is modified from that in the stats package of R.

# Exhibit 8.13
m1.color

# Exhibit 8.14
m2.color=arima(color,order=c(2,0,0))
m2.color

# Exhibit 8.15
m3.color=arima(color,order=c(1,0,1))
m3.color

# Exhibit 8.16
m4.color=arima(color,order=c(2,0,1))
m4.color

# chap9.R

library(TSA)

# Exhibit 9.1
data(color)
m1.color=arima(color,order=c(1,0,0))
m1.color

# Exhibit 9.2
# append 2 years of missing values to the tempdub data as we want to forecast the temperature for two years.
data(tempdub)

tempdub1=ts(c(tempdub,rep(NA,24)),start=start(tempdub),freq=frequency(tempdub))

# creates the first pair of harmonic functions and then fit the model
har.=harmonic(tempdub,1)
m5.tempdub=arima(tempdub,order=c(0,0,0),xreg=har.)
m5.tempdub
# The result is same as that from the fit using lm function.
har.=harmonic(tempdub,1)
model4=lm(tempdub~har.)
summary(model4)


# create the harmonic functions over the period of forecast.
newhar.=harmonic(ts(rep(1,24), start=c(1976,1),freq=12),1)
# Compute and plot the forecasts.
win.graph(width=4.875, height=3,pointsize=8)
plot(m5.tempdub,n.ahead=24,n1=c(1972,1),newxreg=newhar.,type='b',ylab='Temperature',xlab='Year')

# Exhibit 9.3
data(color)
m1.color=arima(color,order=c(1,0,0))
plot(m1.color,n.ahead=12,type='b', xlab='Time', ylab='Color Property')
# add the horizontal line at the estimated mean ("intercept")
abline(h=coef(m1.color)[names(coef(m1.color))=='intercept'])


# Exhibit 9.4
data(hare)
m1.hare=arima(sqrt(hare),order=c(3,0,0))
plot(m1.hare, n.ahead=25,type='b',xlab='Year',ylab='Sqrt(hare)')
abline(h=coef(m1.hare)[names(coef(m1.hare))=='intercept'])

# chap10.R

# Exhibit 10.1
library(TSA)
data(co2)
# The stats pacakage has another dataset of the same name 'co2', so
# make sure that you have loaded the TSA package before loading the co2 data.
win.graph(width=4.875, height=3,pointsize=8)
plot(co2,ylab='CO2',main='Monthly Carbon Dioxide Levels at Alert, NWT, Canada')


# Exhibit 10.2
plot(window(co2,start=c(2000,1)),main='Carbon Dioxide Levels with Monthly Symbols', ylab='CO2')
# specify the y-label as "CO2", otherwise it will use "window(co2,start=c(2000,1))"
Month=c("J","F","M","A","M","J","J","A","S","O","N","D")
points(window(co2,start=c(2000,1)),pch=Month)

# Exhibit 10.3

# First divide the plotting area into a 1 by 2 matrix of plots, i.e. one row of two figures.
win.graph(width=4.875, height=3,pointsize=8)
par(mfrow=c(1,2))
# The ARMAacf function computes the theoretical acf of an ARMA model. Note that the seasonal MA part (1+0.5B)(1+0.8B**12)=(1+0.5B+0.8B**12+0.4B**13).
plot(y=ARMAacf(ma=c(0.5,rep(0,10),0.8,0.4),lag.max=13)[-1],x=1:13,type='h',
xlab='Lag k',ylab=expression(rho[k]),axes=F,ylim=c(0,0.6))
points(y=ARMAacf(ma=c(0.5,rep(0,10),0.8,0.4),lag.max=13)[-1],x=1:13,pch=20)
abline(h=0)
axis(1,at=1:13, labels=c(1,NA,3,NA,5,NA,7,NA,9,NA,11,NA,13))
axis(2)
text(x=7,y=.5,labels=expression(list(theta==-0.5,Theta==-0.8)))

plot(y=ARMAacf(ma=c(-0.5,rep(0,10),0.8,-0.4),lag.max=13)[-1],x=1:13,type='h',
xlab='Lag k',ylab=expression(rho[k]),axes=F)
points(y=ARMAacf(ma=c(-0.5,rep(0,10),0.8,-0.4),lag.max=13)[-1],x=1:13,pch=20)
abline(h=0)
axis(1,at=1:13, labels=c(1,NA,3,NA,5,NA,7,NA,9,NA,11,NA,13))
axis(2)
text(x=7,y=.35,labels=expression(list(theta==0.5,Theta==-0.8)))


# Exhibit 10.4
plot(y=ARMAacf(ar=c(rep(0,11), 0.75), ma=0.4,lag.max=61)[-1],x=1:61,type='h',
xlab='Lag k',ylab=expression(rho[k]),axes=F,ylim=c(-0.4,.8),xlim=c(0,61))
points(y=ARMAacf(ar=c(rep(0,11), 0.75), ma=0.4,lag.max=61)[c(1,11,12,13,23,24,25,35,36,37,47,48,49,59,60,61)+1],
x=c(1,11,12,13,23,24,25,35,36,37,47,48,49,59,60,61),pch=20)
abline(h=0)
axis(1,at=c(0,1,12,24,36,48,60,61),labels=c(NA,1,12,24,36,48,60,NA))
axis(2, at=c(-0.4,0.0,0.4,0.8))
text(x=40,y=.8,labels=expression(list(Phi==0.75,theta==-0.4)))

plot(y=ARMAacf(ar=c(rep(0,11), 0.75), ma=-0.4,lag.max=61)[-1],x=1:61,type='h',
xlab='Lag k',ylab=expression(rho[k]),axes=F,ylim=c(-0.4,.8))
points(y=ARMAacf(ar=c(rep(0,11), 0.75), ma=-0.4,lag.max=61)[c(1,11,12,13,23,24,25,35,36,37,47,48,49,59,60,61)+1],
x=c(1,11,12,13,23,24,25,35,36,37,47,48,49,59,60,61),pch=20)
abline(h=0)
axis(1,at=c(0,1,12,24,36,48,60,61),labels=c(NA,1,12,24,36,48,60,NA))
axis(2, at=c(-0.4,0.0,0.4,0.8))
text(x=40,y=.8,labels=expression(list(Phi==0.75,theta==0.4)))

# Exhibit 10.5
par(mfrow=c(1,1))
acf(as.vector(co2),lag.max=36,
main=expression(Sample~~ACF~~of~~CO[2]~~Levels))

# Exhibit 10.6
plot(diff(co2),main=expression(Time~~Series~~Plot~~of~~the~~First~~Differences~~of~~
CO[2]~~Levels), ylab=expression(First~~Difference~~of~~CO[2]))

# Exhibit 10.7
acf(as.vector(diff(co2)),lag.max=36,main=expression(Sample~~ACF~~of~~the~~First~~Differences~~of~~CO[2]~~Levels))

# Exhibit 10.8
plot(diff(diff(co2),lag=12),main=expression(Time~~Series~~Plot~~of~~the~~First~~and~~Seasonal~~Differences~~of~~CO[2]~~Levels),
ylab=expression(First~~and~~Seasonal~~Difference~~of~~C~O[2]))

# Exhibit 10.9
acf(as.vector(diff(diff(co2),lag=12)),lag.max=36,ci.type='ma',
main=expression(Sample~~ACF~~of~~the~~First~~and~~Seasonal~~Differences~~of~~CO[2]~~Levels))

# Exhibit 10.10
# Do remember that in the book the MA parameterization uses the minus convention but R uses# the positive convention, lest you find the estimates from R to be different from the reported values # in the book!
m1.co2=arima(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
m1.co2

# Exhibit 10.11
# first thirteen residuals are omitted from the plot.
plot(window(rstandard(m1.co2),start=c(1995,2)),ylab='Standardized Residuals', type='o',
main=expression(Residuals~~from~~the~~ARIMA(list(0,1,1))%*%(list(0,1,1))[12]~~Model))
abline(h=0)

# Exhibit 10.12
acf(as.vector(window(rstandard(m1.co2),start=c(1995,2))),lag.max=36,
main=expression(ACF~~of~~Residuals~~from~~the~~ARIMA(list(0,1,1))%*%(list(0,1,1))[12]~~Model))

# The figures in the above two exhibits can also be obtained as the first two sub-plots from the # following command. The plotting convention is slightly different for the first sub-plot.
win.graph(width=4.875, height=5,pointsize=8)
# by default, the first 13 residuals are ommited; please use ?tsdiag to learn more about the function.
tsdiag(m1.co2, gof.lag=36)

# Exhibit 10.13
hist(window(rstandard(m1.co2),start=c(1995,2)),xlab='Standardized Residuals',
main=expression(Residuals~~from~~the~~ARIMA(list(0,1,1))%*%(list(0,1,1))[12]~~Model))

# Exhibit 10.14
win.graph(width=4, height=4,pointsize=8)
qqnorm(window(rstandard(m1.co2),start=c(1995,2)),main=expression(Normal~~Q-Q~~Plot))
title(main=expression(Residuals~~from~~the~~ARIMA(list(0,1,1))%*%(list(0,1,1))[12]~~Model),
line=3)
qqline(window(rstandard(m1.co2),start=c(1995,2)))

# Exhibit 10.15
m2.co2=arima(co2,order=c(0,1,2),seasonal=list(order=c(0,1,1),period=12))
m2.co2

# Exhibit 10.16
win.graph(width=4.875, height=3,pointsize=8)
plot(m1.co2,n1=c(2003,1),n.ahead=24,col='red',xlab='Year',type='o',
ylab=expression(CO[2]~~Levels),
main=expression(Forecasts~~and~~Forecast~~Limits~~'for'~~the~~CO[2]~~Model))
# Note that for is a reserved word in R so it has to be enclosed in quotation marks.

# Exhibit 10.17
plot(m1.co2,n1=c(2004,1),n.ahead=48,col='red',xlab='Year',type='o',
ylab=expression(CO[2]~~Levels),
main=expression(Long~~Term~~Forecasts~~'for'~~the~~CO[2]~~Model))

Thank you.

Thank you!