jj = ts(scan("/mydata/jj.dat"), start=1960, frequency=4)
plot(jj, ylab="Quarterly Earnings per Share")
Example 1.2
globtemp = ts(scan("/mydata/globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900) # use data from 1990
plot(gtemp, type="o", ylab="Global Temperature Deviations")
Example 1.3
plot(speech <- ts(scan("/mydata/speech.dat")), ylab="speech")
Example 1.4
plot(nyse <- ts(scan("/mydata/nyse.dat")), ylab="NYSE Returns")
Example 1.5
soi = ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(2,1), mar=c(3, 4, 3, 2)) # sets up the graphics - ?par for info
plot(soi, ylab="", xlab="", main="Southern Oscillation Index")
plot(rec, ylab="", xlab="", main="Recruitment")
Example 1.6
fmri = read.table("/mydata/fmri.dat")
par(mfrow=c(2,1), mar=c(5, 4, 3, 2))
ts.plot(fmri[,2:5],lty=1:4,ylab="BOLD",xlab="Time (1 pt = 2 sec)",main="Cortex",col=1:4)
ts.plot(fmri[,6:9],lty=1:4,ylab="BOLD",xlab="",main="Thalamus & Cerebellum",col=1:4)
# col=1:4 adds color but if you print monochrome you may want to remove it
Example 1.7
x = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
par(mfrow=c(2,1))
plot.ts(x[,1], main="Earthquake", ylab="EQ5")
plot.ts(x[,2], main="Explosion", ylab="EXP6")
Example 1.9
w = rnorm(500,0,1) # 500 N(0,1) variates
v = filter(w, sides=2, rep(1,3)/3) # moving average
par(mfrow=c(2,1))
plot.ts(w, main="white noise")
plot.ts(v, ylim=c(-3,3), main="moving average")
# now try this:
ts.plot(w, v, lty=2:1, col=1:2, lwd=1:2)
Example 1.10
w = rnorm(550,0,1) # 50 extra to avoid startup problems - the [-(1:50)] on the next line...
x = filter(w, filter=c(1,-.9), method="recursive")[-(1:50)] #... removes the first 50 values
plot.ts(x)
Example 1.11
set.seed(154)
w = rnorm(200,0,1)
x = cumsum(w)
wd = w +.2
xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55))
lines(x)
lines(.2*(1:200), lty="dashed")
Example 1.12
c = 2*cos(2*pi*(1:500)/50 + .6*pi)
w = rnorm(500,0,1)
par(mfrow=c(3,1), mar=c(3,2,2,1), cex.main=1.5)
plot.ts(c, main=expression(x[t]==2*cos(2*pi*t/50+.6*pi)))
plot.ts(c+w, main=expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,1)))
plot.ts(c+5*w, main=expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,25)))
Example 1.24
speech = scan("/mydata/speech.dat")
acf(speech,250)
Example 1.25*
soi = ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(3,1))
acf(soi, 50)
acf(rec, 50)
ccf(soi, rec, 50)
* You don't have to read the data in again if you already did it for Example 1.5 and you've saved your workspace. This goes for the rest of the examples. Type objects() to see a list of your objects.
top
CHAPTER 2
Example 2.1
globtemp = ts(scan("/mydata/globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900) # use years 1990 to 1997
fit = lm(gtemp~time(gtemp), na.action=NULL) # regress gtemp on time
#-- note: na.action=NULL is used to retain time series attributes - ?lm for info
summary(fit) # view the results
plot(gtemp, type="o", ylab="Global Temperature Deviation") # plot the series
abline(fit) # add regression line to the plot
Example 2.2
mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
temp = ts(scan("/mydata/temp.dat"), start=1970, frequency=52)
part = ts(scan("/mydata/part.dat"), start=1970, frequency=52)
par(mfrow=c(3,1), mar=c(3,4,3,2))
plot(mort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(temp, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
x11() # open another graphic device
pairs(cbind(mort, temp, part)) # scatterplot matrix
temp = temp-mean(temp)
temp2 = temp^2
trend = time(mort) # time
fit = lm(mort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit) # regression results
anova(lm(mort~cbind(trend, temp, temp2, part))) # pretty ANOVA table
num = length(mort) # sample size
AIC(fit)/num - log(2*pi) # AIC (corresponds to def 2.1)
AIC(fit, k=log(num))/num - log(2*pi) # BIC (corresponds to def 2.3)
(AICc = log(sum(resid(fit)^2)/num)+(num+5)/(num-5-2)) # AICc (as in def 2.2)
Examples 2.3 and 2.4
globtemp = ts(scan("/mydata/globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900) # use years 1990 to 1997
fit = lm(gtemp~time(gtemp), na.action=NULL) # regress gtemp on time
par(mfrow=c(2,1))
plot(resid(fit), type="o", main="detrended")
plot(diff(gtemp), type="o", main="first difference") # this and below will do example 2.4
x11() # open a new graphic device
par(mfrow=c(3,1)) # plot ACFs
acf(gtemp, 48)
acf(resid(fit), 48)
acf(diff(gtemp), 48)
Example 2.6
You can download a function called lag.plot2 that will do the second part of this example over here.
soi = ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
lag.plot(soi, lags=12, layout=c(3,4), diag=FALSE)
par(mfrow=c(3,3), mar=c(2.5, 4, 4, 1)) # set up plot area
for(h in 0:8){ # loop through lags 0-8
plot(lag(soi,-h),rec, main=paste("soi(t-",h,")",sep=""), ylab="rec(t)",xlab="")
}
Example 2.8
x = 2*cos(2*pi*(1:500)/50 + .6*pi) + rnorm(500,0,5)
I = abs(fft(x)/sqrt(500))^2 # the periodogram
P = (4/500)*I # the scaled periodogram
f = 0:250/500
plot(f, P[1:251], type="l", xlab="frequency", ylab=" ")
abline(v=seq(0,.5,.02), lty="dotted")
# Note: the dominant period is 50 points or 1/50 = .02 cycles per point
Example 2.10
mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
ma5 = filter(mort, sides=2, rep(1,5)/5)
ma53 = filter(mort, sides=2, rep(1,53)/53)
plot(mort, type="p", xlab="week", ylab="mortality")
lines(ma5)
lines(ma53)
Example 2.11
wk = time(mort) - mean(time(mort)) # week (centered) - wk is basically t/52 ...
wk2 = wk^2 # ... because frequency=52 for mort
wk3 = wk^3
c = cos(2*pi*wk)
s = sin(2*pi*wk)
reg1 = lm(mort~wk + wk2 + wk3, na.action=NULL)
reg2 = lm(mort~wk + wk2 + wk3 + c + s, na.action=NULL)
plot(mort, type="p", ylab="mortality")
lines(fitted(reg1))
lines(fitted(reg2))
# here's the original code where mort isn't a ts object
mort = scan("/mydata/cmort.dat")
t = 1:length(mort)
t2 = t^2
t3 = t^3
c = cos(2*pi*t/52)
s = sin(2*pi*t/52)
fit1 = lm(mort~t + t2 + t3)
fit2 = lm(mort~t + t2 + t3 + c + s)
plot(t, mort)
lines(fit1$fit)
lines(fit2$fit)
Example 2.12
plot(mort, type="p", ylab="mortality")
lines(ksmooth(time(mort), mort, "normal", bandwidth=10/52)) # or try bandwidth=5/52
lines(ksmooth(time(mort), mort, "normal", bandwidth=2))
# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(ksmooth(t, mort, "normal", bandwidth=10))
lines(ksmooth(t, mort, "normal", bandwidth=104))
# Red indicates a change; Fig. 2.15 shows a bandwidth of 10.
Example 2.13
par(mfrow=c(2,1))
plot(mort,type="p", ylab="mortality", main="nearest neighbor")
lines(supsmu(time(mort), mort, span=.5))
lines(supsmu(time(mort), mort, span=.01))
plot(mort, type="p", ylab="mortality", main="lowess")
lines(lowess(mort, f=.02))
lines(lowess(mort, f=2/3))
Example 2.14
plot(mort, type="p", ylab="mortality")
lines(smooth.spline(time(mort), mort))
lines(smooth.spline(time(mort), mort, spar=1))
# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(smooth.spline(t, mort, spar=.0000001))
lines(smooth.spline(t, mort, spar=.1))
Example 2.15
par(mfrow=c(2,1))
plot(temp, mort, main="lowess")
lines(lowess(temp,mort))
plo
评论0