# 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)