Time Series Analysis with Applications in R, Second Edition Jonathan D. Cryer

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

123

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

关键词：Applications Time Series Application Analysis Jonathan Series Robert

Time Series Analysis with Applications in R, Second Edition Jonathan D. Cryer an.rar (6.82 MB, 需要: 1 个论坛币)

2014-1-6 06:14:43 上传

Time Series Analysis with Applications in R, Second Edition Jonathan D. Cryer and Kung-Sik Chan
需要: 1 个论坛币  [购买]

http://homepage.stat.uiowa.edu/~kchan/TSA.htm
http://cran.r-project.org/web/packages/TSA/TSA.pdf

# Exhibit 2.1
win.graph(width=4.875, height=2.5,pointsize=8)
# rwalk contains a simulated random walk
data(rwalk)
plot(rwalk,type='o',ylab='Random Walk')

# R code for simulating a random walk with, say 60, iid standard normal errors
n=60
set.seed(12345) # intialize the random number so that the simulation can be
# reproducible.
sim.random.walk=ts(cumsum(rnorm(n)),freq=1,start=1)
plot(sim.random.walk,type='o',ylab='Another Random Walk')

# Exhibit 3.1
# time(rwalk) yields a time series of the time epoches when the random walk was sampled.
data(rwalk)
model1=lm(rwalk~time(rwalk))
summary(model1)

# Exhibit 3.2
win.graph(width=4.875, height=2.5,pointsize=8)
# rwalk contains a simulated random walk
plot(rwalk,type='o',ylab='y')
abline(model1) # add the fitted least squares line

# Exhibit 3.3
# season(tempdub) creates a vector of the month index of the data as a factor
data(tempdub)
month.=season(tempdub) # the period sign is included to make the printout from
# the commands two line below clearer; ditto below.
model2=lm(tempdub~month.-1) # -1 removes the intercept term
summary(model2)

# Exhibit 3.4
model3=lm(tempdub~month.) # intercept is automatically included so one month (Jan) is dropped
summary(model3)

# Exhibit 3.5
# first creates the first pair of harmonic functions and then fit the model
har.=harmonic(tempdub,1)
model4=lm(tempdub~har.)
summary(model4)

# Exhibit 3.6
win.graph(width=4.875, height=2.5,pointsize=8)
plot(ts(fitted(model4),freq=12,start=c(1964,1)),ylab='Temperature',type='l',
ylim=range(c(fitted(model4),tempdub))) # the ylim option ensures that the
# y axis has a range that fits the raw data and the fitted values
points(tempdub)

# Exhibit 3.7
data(rwalk)
model1=lm(rwalk~time(rwalk))
summary(model1)

# Exhibit 3.8
plot(y=rstudent(model3),x=as.vector(time(tempdub)),xlab='Time',
ylab='Standardized Residuals',type='o')

# Exhibit 3.9
plot(y=rstudent(model3),x=as.vector(time(tempdub)),xlab='Time',
ylab='Standardized Residuals',type='l')
points(y=rstudent(model3),x=as.vector(time(tempdub)),
pch=as.vector(season(tempdub)))

# Exhibit 3.10
plot(y=rstudent(model3),x=as.vector(fitted(model3)),xlab='Fitted Trend Values',
ylab='Standardized Residuals',type="n")
points(y=rstudent(model3),x=as.vector(fitted(model3)),
pch=as.vector(season(tempdub)))

# Exhibit 3.11
hist(rstudent(model3),xlab='Standardized Residuals',main='')

# Exhibit 3.12
win.graph(width=3, height=3,pointsize=8)
qqnorm(rstudent(model3),main='')

# Exhibit 3.13
win.graph(width=4.875, height=3,pointsize=8)
acf(rstudent(model3),main='')

# Exhibit 3.14
plot(y=rstudent(model1),x=as.vector(time(rwalk)),ylab='Standardized Residuals',
xlab='Time',type='o')

# Exhibit 3.15
win.graph(width=4.875, height=3,pointsize=8)
plot(y=rstudent(model1),x=fitted(model1),ylab='Standardized Residuals',
xlab='Fitted Trend Values',type='p')

# Exhibit 3.16
acf(rstudent(model1),main='')

# Exhibit 4.2
win.graph(width=4.875, height=3,pointsize=8)
data(ma1.2.s)
plot(ma1.2.s,ylab=expression(Y[t]),type='o')

# An MA(1) series with MA coefficient equal to -0.9 and of length n=100 can be simulated by the following command
set.seed(12345) # initialize the seed of the random number generator so that
# the simulations can be reproduced.
y=arima.sim(model=list(ma=-c(-0.9)),n=100)
# Note that R uses the plus convention in the model formula so the additional minus sign.  

# Exhibit 4.3
win.graph(width=3, height=3,pointsize=8)
plot(y=ma1.2.s,x=zlag(ma1.2.s),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p')

# Exhibit 4.4
plot(y=ma1.2.s,x=zlag(ma1.2.s,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p')


# Exhibit 4.5
win.graph(width=4.875, height=3,pointsize=8)
data(ma1.1.s)
plot(ma1.1.s,ylab=expression(Y[t]),type='o')

# An MA(1) series with ma coefficient equal to 0.9 and of length n=100 can be simulated by the following command
y=arima.sim(model=list(MA=-c(0.9)),n=100)
# Note that R uses the plus convention in the MA model formula so the additional minus sign.  


# Exhibit 4.6
win.graph(width=3, height=3,pointsize=8)
plot(y=ma1.1.s,x=zlag(ma1.1.s),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p')

# Exhibit 4.7
plot(y=ma1.1.s,x=zlag(ma1.1.s,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p')

# Exhibit 4.8
win.graph(width=4.875, height=3,pointsize=8)
data(ma2.s)
plot(ma2.s,ylab=expression(Y[t]),type='o')

# An MA(2) series with MA coefficients equal to 1 and -0.6 and of length n=100 can be simulated by the following command
y=arima.sim(model=list(ma=-c(1, -0.6)),n=100)
# Note that R uses the plus convention in the MA model formula so the additional minus sign.  

# Exhibit 4.9
win.graph(width=3, height=3,pointsize=8)
plot(y=ma2.s,x=zlag(ma2.s),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p')

# Exhibit 4.10
plot(y=ma2.s,x=zlag(ma2.s,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p')

# Exhibit 4.11
plot(y=ma2.s,x=zlag(ma2.s,3),ylab=expression(Y[t]),xlab=expression(Y[t-3]),type='p')

# Exhibit 4.13
win.graph(width=4.875, height=3,pointsize=8)
data(ar1.s)
plot(ar1.s,ylab=expression(Y[t]),type='o')

# An AR(1) series with AR coefficient equal to 0.9 and of length n=100 can be simulated by the following command
y=arima.sim(model=list(ar=c(0.9)),n=100)
# Note that the R convention for the AR model formula is same as the book, so  NO additional minus sign.  

# Exhibit 4.14
win.graph(width=3, height=3,pointsize=8)
plot(y=ar1.s,x=zlag(ar1.s),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p')

# Exhibit 4.15
plot(y=ar1.s,x=zlag(ar1.s,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p')

# Exhibit 4.16
plot(y=ar1.s,x=zlag(ar1.s,3),ylab=expression(Y[t]),xlab=expression(Y[t-3]),type='p')

# Exhibit 4.19
win.graph(width=4.875, height=3,pointsize=8)
data(ar2.s)
plot(ar2.s,ylab=expression(Y[t]),type='o')

# An AR(2) series with AR coefficients equal to 1.5 and -0.75 and of length n=100 can be simulated by the following command
y=arima.sim(model=list(ar=c(1.5,-0.75)),n=100)
# Note that the R convention for the AR model formula is same as the book, so  NO additional minus sign.

# Chap5.R

# Exhibit 5.1
win.graph(width=4.875, height=3,pointsize=8)
data(oil.price)
plot(oil.price, ylab='Price per Barrel',type='l')

# Exhibit 5.3
data(explode.s)
plot(explode.s,ylab=expression(y[t]),type='o')

# Exhibit 5.4
plot(diff(log(oil.price)),ylab='Change in Log(Price)',type='l')

# Exhibit 5.5
data(ima22.s)
plot(ima22.s,ylab="IMA(2,2) Simulation",type='o')

# Exhibit 5.6
plot(diff(ima22.s),ylab='First Difference',type='o')

# Exhibit 5.7
plot(diff(ima22.s,difference=2),ylab='Differenced Twice',type='o')

# Note that plot(diff(ima22.s,2),ylab='Differenced Twice',type='o') will draw a wrong figure because the second argument is the lag not the times of differencing. That is, diff(ima22.s,2) is the series of ima22.s(t)-ima22.s(t-2).


# Exhibit 5.8
data(electricity)
plot(electricity)

# Exhibit 5.9
plot(log(electricity),ylab='Log(electricity)') # without specifying the y-label the default y-label is "electricity" rather than "log(electricity)"!

# Exhibit 5.10
plot(diff(log(electricity)),ylab='Difference of Log(electricity)')

# Exhibit 5.11
win.graph(width=3, height=3,pointsize=8)
BoxCox.ar(electricity)
# In case the function fails, check wheher all data are positive. If not, shift all data by a fixed constant. If the function still fails,
# try setting method='ols', as an alternative to the default method of 'mle'.

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