Everything you see in a box below is R code. You can copy-and-paste any line (or multiple lines) into R. Make sure the data files are in the mydata directory (or change the code accordingly). R code for Chapter 6 has its own page and we also have some code for Chapter 7. We have improved, corrected or streamlined some of the code below, so there are a few cases where the code here differs slightly from what is in the text.
Disclaimer: Sometimes I use t for time, e.g., t = 1984:2001. This is considered bad R practice because t is used for matrix transpose, i.e., t(A) is the transpose of matrix A. If you don't want to be bad, you can replace our t in the code below with something else, BUT NOT time because that's reserved, too. Also, when I can, I prefer to use the simple assignment character (=) to the more complex assignment characters (<- or ->); this is also considered bad R practice. I like to be bad sometimes, but you've been warned ... so don't blame me if you turn to a life of crime.
jj = ts(scan("/mydata/jj.dat"), start=1960, frequency=4)
plot(jj, ylab="Quarterly Earnings per Share")
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")
plot(speech <- ts(scan("/mydata/speech.dat")), ylab="speech")
plot(nyse <- ts(scan("/mydata/nyse.dat")), ylab="NYSE Returns")
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")
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
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")
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)
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)
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")
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)))
speech = scan("/mydata/speech.dat")
acf(speech,250)
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.
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
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)
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)
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=F)
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="")
}
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
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)
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)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.
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))
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))
par(mfrow=c(2,1)) plot(temp, mort, main="lowess") lines(lowess(temp,mort)) plot(temp, mort, main="smoothing splines") lines(smooth.spline(temp,mort))
par(mfrow=c(2,1)) plot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x",main=(expression(AR(1)~~~phi==+.9))) # ~ is a space and == is equal plot(arima.sim(list(order=c(1,0,0), ar=-.9), n=100), ylab="x",main=(expression(AR(1)~~~phi==-.9)))
par(mfrow=c(2,1)) plot(arima.sim(list(order=c(0,0,1), ma=.5), n=100), ylab="x",main=(expression(MA(1)~~~theta==+.5))) plot(arima.sim(list(order=c(0,0,1), ma=-.5), n=100), ylab="x",main=(expression(MA(1)~~~theta==-.5)))
set.seed(5) ar2 = arima.sim(list(order=c(2,0,0), ar=c(1.5,-.75)), n = 144) t = (1:144)/12 plot(t, ar2, type="l", xlab="Time (one unit = 12 points)") abline(v=0:12, lty="dotted")
ARMAtoMA(ar=.9, ma=.5, 50) # for a list plot(ARMAtoMA(ar=.9, ma=.5, 50)) # for a graph
ar2.acf = ARMAacf(ar=c(1.5,-.75), ma=0, 24) ar2.pacf = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=T) par(mfrow=c(1,2)) plot(ar2.acf, type="h", xlab="lag") abline(h=0) plot(ar2.pacf, type="h", xlab="lag") abline(h=0)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(2,1)); acf(rec, 48); pacf(rec, 48)
(fit = ar.ols(rec,order=2,demean=F,intercept=T)) # estimates
fit$asy.se # standard errors
rec.yw = ar.yw(rec, order=2) rec.yw$x.mean rec.yw$ar sqrt(diag(rec.yw$asy.var.coef)) rec.yw$var.pred rec.pr = predict(rec.yw, n.ahead=24) U = rec.pr$pred + rec.pr$se L = rec.pr$pred - rec.pr$se minx = min(rec,L) maxx = max(rec,U) ts.plot(rec, fore$pred, xlim=c(1980,1990), ylim=c(minx,maxx)) lines(rec.pr$pred, col="red", type="o") lines(U, col="blue", lty="dashed") lines(L, col="blue", lty="dashed")
rec.mle = ar.mle(rec, order=2) rec.mle$x.mean rec.mle$ar sqrt(diag(rec.mle$asy.var.coef)) rec.mle$var.pred
x = scan("/mydata/ar1boot.dat") # see note below
m = mean(x) # estimate of mu
fit = ar.yw(x, order=1)
phi = fit$ar # estimate of phi
nboot = 200 # number of bootstrap replicates
resids = fit$resid[-1] # the first resid is NA (so remove it)
x.star = x # initialize x.star
phi.star = matrix(0, nboot, 1)
for (i in 1:nboot) {
resid.star = sample(resids, replace=TRUE)
for (t in 1:99){
x.star[t+1] = m + phi*(x.star[t]-m) + resid.star[t]
}
phi.star[i] = ar.yw(x.star, order=1)$ar
}
hist(phi.star, prob=T) # histogram of bootstrapped φs
# If you want to simulate your own data, use the following commands:
e = rexp(150, rate = .5)
u = runif(150,-1,1)
de = e*sign(u)
x = 50 + arima.sim(n = 100, list(ar =.95), innov = de, n.start = 50)
gnp96 = read.table("/mydata/gnp96.dat")
gnp = ts(gnp96[,2], start=1947, frequency=4)
plot(gnp)
acf(gnp, 50)
gnpgr = diff(log(gnp)) # growth rate
plot.ts(gnpgr)
par(mfrow=c(2,1))
acf(gnpgr, 24)
pacf(gnpgr, 24)
# Fit ARIMA and view results:
(gnpgr.ar = arima(gnpgr, order = c(1, 0, 0))) # potential problem here (see R Issue 1)
(gnpgr.ma = arima(gnpgr, order = c(0, 0, 2)))
ARMAtoMA(ar=.35, ma=0, 10) # prints psi-weights
tsdiag(gnpgr.ma, gof.lag=20)
varve = scan("/mydata/varve.dat")
(varve.ma = arima(log(varve), order = c(0, 1, 1)))
tsdiag(varve.ma)
(varve.arma = arima(log(varve), order = c(1, 1, 1)))
tsdiag(varve.arma, gof.lag=20)
n = length(gnpgr) # sample size
kma = length(gnpgr.ma$coef) # number of parameters in ma model
sma=gnpgr.ma$sigma2 # mle of sigma^2
kar = length(gnpgr.ar$coef) # number of parameters in ar model
sar=gnpgr.ar$sigma2 # mle of sigma^2
# AIC Returned Value
log(sma) + (n+2*kma)/n # MA2 -8.298
log(sar) + (n+2*kar)/n # AR1 -8.294
# AICc
log(sma) + (n+kma)/(n-kma-2) # MA2 -8.288
log(sar) + (n+kar)/(n-kar-2) # AR1 -8.285
# BIC
log(sma) + kma*log(n)/n # MA2 -9.252
log(sar) + kar*log(n)/n # AR1 -9.264
## Note: You can use sarima.R to get the above
## (remember to source the code first):
sarima(gnpgr,0,0,2) # for the MA(2)
sarima(gnpgr,1,0,0) # for the AR(1)
phi = c(rep(0,11),.8) acf = ARMAacf(ar=phi, ma=-.5, 50) pacf = ARMAacf(ar=phi, ma=-.5, 50, pacf=T) par(mfrow=c(1,2)) plot(acf, type="h", xlab="lag") abline(h=0) plot(pacf, type="h", xlab="lag") abline(h=0)
prod=ts(scan("/mydata/prod.dat"), start=1948, frequency=12)
par(mfrow=c(2,1)) # (P)ACF of data
acf(prod, 48)
pacf(prod, 48)
par(mfrow=c(2,1)) # (P)ACF of d1 data
acf(diff(prod), 48)
pacf(diff(prod), 48)
par(mfrow=c(2,1)) # (P)ACF of d1-d12 data
acf(diff(diff(prod),12), 48)
pacf(diff(diff(prod),12), 48)
### fit model (iii)
prod.fit3 = arima(prod, order=c(1,1,1), seasonal=list(order=c(2,1,1), period=12))
prod.fit3 # to view the results
tsdiag(prod.fit3, gof.lag=48) # diagnostics
### forecasts for the final model
prod.pr = predict(prod.fit3, n.ahead=12)
U = prod.pr$pred + 2*prod.pr$se
L = prod.pr$pred - 2*prod.pr$se
ts.plot(prod,prodpr$pred, col=1:2, type="o", ylim=c(105,175), xlim=c(1975,1980))
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
## Note: This can be done using acf2.R, sarima.R
## and sarima.for.R (over here) as follows
## (remember to source the code first):
prod=scan("/mydata/prod.dat")
acf2(prod,48)
acf2(diff(prod), 48)
acf2(diff(diff(prod),12), 48)
### fit model (iii)
sarima(prod,1,1,1,2,1,1,12)
### forecasts for the final model
sarima.for(prod,12,1,1,1,2,1,1,12)
##
## -- See the note for this example in the errata.
t = 1:100 x1 = 2*cos(2*pi*t*6/100) + 3*sin(2*pi*t*6/100) x2 = 4*cos(2*pi*t*10/100) + 5*sin(2*pi*t*10/100) x3 = 6*cos(2*pi*t*40/100) + 7*sin(2*pi*t*40/100) x = x1 + x2 + x3 par(mfrow=c(2,2)) plot.ts(x1, ylim=c(-10,10), main = expression(omega==6/100~~~A^2==13)) plot.ts(x2, ylim=c(-10,10), main = expression(omega==10/100~~~A^2==41)) plot.ts(x3, ylim=c(-10,10), main = expression(omega==40/100~~~A^2==85)) plot.ts(x, ylim=c(-16,16),main="sum")
P = abs(2*fft(x)/100)^2 f = 0:50/100 plot(f, P[1:51], type="o", xlab="frequency", ylab="periodogram")
x = c(1,2,3,2,1) c1 = cos(2*pi*1:5*1/5) s1 = sin(2*pi*1:5*1/5) c2 = cos(2*pi*1:5*2/5) s2 = sin(2*pi*1:5*2/5) omega1 = cbind(c1, s1) omega2 = cbind(c2, s2) anova(lm(x~omega1+omega2)) # ANOVA Table abs(fft(x))^2/5 # the periodogram (as a comparison)
soi = scan("/mydata/soi.dat")
rec = scan("/mydata/recruit.dat")
par(mfrow=c(2,1))
soi.per = spec.pgram(soi, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=1/48, lty="dotted")
rec.per = spec.pgram(rec, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=1/48, lty="dotted")
soi.per$spec[40] # soi pgram at freq 1/12 = 40/480
soi.per$spec[10] # soi pgram at freq 1/48 = 10/480
# -- conf intervals -- # returned value:
U = qchisq(.025,2) # 0.05063562
L = qchisq(.975,2) # 7.377759
2*soi.per$spec[10]/L # 0.1747835
2*soi.per$spec[10]/U # 25.46648
2*soi.per$spec[40]/L # 3.162688
2*soi.per$spec[40]/U # 460.813
# -- replace soi with rec above to get recruit values
#######################################################################################
###- Repeat example using "frequency=12" for monthly data and note the differences -###
###- basically, "everything" is in multiples of Δ = 1/12 -###
#######################################################################################
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))
soi.per = spec.pgram(soi, taper=0, log="no")
abline(v=1, lty="dotted")
abline(v=1/4, lty="dotted")
rec.per = spec.pgram(rec, taper=0, log="no")
abline(v=1, lty="dotted")
abline(v=1/4, lty="dotted")
soi.per$spec[40] # soi pgram at freq 1/12 = 40/480
soi.per$spec[10] # soi pgram at freq 1/48 = 10/480
# -- conf intervals -- # returned value:
U = qchisq(.025,2) # 0.05063562
L = qchisq(.975,2) # 7.377759
2*soi.per$spec[10]/L # 0.01456530
2*soi.per$spec[10]/U # 2.122207
2*soi.per$spec[40]/L # 0.2635573
2*soi.per$spec[40]/U # 38.40108
par(mfrow=c(2,1))
k = kernel("daniell",4)
soi.ave = spec.pgram(soi, k, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=2/12, lty="dotted")
abline(v=3/12, lty="dotted")
abline(v=1/48, lty="dotted")
# -- Repeat 5 lines above using rec in place of soi
soi.ave$bandwidth # reported bandwidth
soi.ave$bandwidth*sqrt(12) # Bw
df = soi.ave$df # df = 16.9875 (returned values)
U = qchisq(.025,df) # U = 7.555916
L = qchisq(.975,df) # L = 30.17425
soi.ave$spec[10] # 0.5942431
soi.ave$spec[40] # 1.428959
# -- intervals --
df*soi.ave$spec[10]/L # 0.334547
df*soi.ave$spec[10]/U # 1.336000
df*soi.ave$spec[40]/L # 0.8044755
df*soi.ave$spec[40]/U # 3.212641
# -- repeat above commands with soi replaced by rec
kernel("modified.daniell", c(3,3)) # for a list
plot(kernel("modified.daniell", c(3,3))) # for a graph
par(mfrow=c(2,1))
k = kernel("modified.daniell", c(3,3))
soi.smo = spec.pgram(soi, k, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=1/48, lty="dotted")
# -- Repeat above 3 lines with rec replacing soi
df = soi.smo$df # df=17.42618
Lh = 1/sum(k[-k$m:k$m]^2) # Lh=9.232413
Bw = Lh/480 # Bw=0.01923419
# An easier way to obtain soi.smo (you can add log="no"):
soi.smo = spectrum(soi, spans=c(7,7), taper=0)
par(mfrow=c(3,1)) spectrum(soi, spans=c(7,7), taper=0, main="No Taper") abline(v=1/12,lty="dashed") abline(v=1/48,lty="dashed") spectrum(soi, spans=c(7,7), main="10% Taper") abline(v=1/12,lty="dashed") abline(v=1/48,lty="dashed") spectrum(soi, spans=c(7,7), taper=.5, main="50% Taper") abline(v=1/12,lty="dashed") abline(v=1/48,lty="dashed")
x = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
eqP = x[1:1024, 1]; eqS = x[1025:2048, 1]
exP = x[1:1024, 2]; exS = x[1025:2048, 2]
par(mfrow=c(2,2))
eqPs=spectrum(eqP, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.04))
eqSs=spectrum(eqS, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.4))
exPs=spectrum(exP, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.04))
exSs=spectrum(exS, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.4))
exSs$df
x = ts(cbind(soi,rec))
s = spec.pgram(x, kernel("daniell",9), taper=0)
s$df # df = 35.8625
f = qf(.999, 2, s$df-2) # f = 8.529792
c = f/(18+f) # c = 0.3188779
plot(s, plot.type = "coh", ci.lty = 2)
abline(h = c)
par(mfrow=c(3,1))
plot.ts(soi) # the data
plot.ts(diff(soi)) # first difference
k = kernel("modified.daniell", 6) #-- 12 month filter
soif = kernapply(soi,k)
plot.ts(soif)
x11() # open new graphics device
spectrum(soif, spans=9, log="no") #-- spectral analysis
abline(v=1/52, lty="dotted")
windows()
w = seq(0,.5, length=1000) #-- frequency response
FR = abs(1-exp(2i*pi*w))^2
plot(w, FR, type="l")
spec.ar(soi, log="no") # plot min AIC spectrum
abline(v=1/52, lty="dotted") # locate El Niño period
abline(v=1/12, lty="dotted") # locate yearly period
soi.ar = ar(soi, order.max=30) # obtain AICs **(see below)
plot(0:30, soi.ar$aic, type="l") # plot AICs
soi.ar # view results
#** Here's a nicer display with the three information criteria:
n = length(soi)
AIC = matrix(0,30,1) -> AICc -> BIC -> sig2
for (k in 1:30){
sigma2=ar(soi, order=k, aic=F, method="burg", var.method=2)$var.pred
#-- ?ar details suggest the above is the mle of σ2
BIC[k]=log(sigma2)+(k*log(n)/n)
AICc[k]=log(sigma2)+((n+k)/(n-k-2))
AIC[k]=log(sigma2)+((n+2*k)/n)
sig2[k]=sigma2
}
IC=cbind(AIC, AICc, BIC+1) # BIC+1 reduces white space in the plot
ts.plot(IC,type="o",xlab="order",ylab="AIC / BIC",main="Information Criteria",col=c(1,2,4))
text(12.5,-1.55,"AIC", col=1)
text(14.7,-1.48,"AICc", col=2)
text(15.5,-1.39,"BIC", col=4)
abline(v=1:30, lty=3)
eqexp = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
ex = eqexp[,2] # the explosion series
## -- dynamic spectral analysis -- ##
nobs = length(ex) # number of observations
wsize = 256 # window size
overlap = 128 # overlap
ovr = wsize-overlap
nseg = floor(nobs/ovr)-1; # number of segments
krnl = kernel("daniell", c(1,1)) # kernel
ex.spec = matrix(0,wsize/2,nseg)
for (k in 1:nseg) {
a = ovr*(k-1)+1
b = wsize+ovr*(k-1)
ex.spec[,k]=spectrum(ex[a:b],krnl,taper=.5,plot=F)$spec
}
## --plot results -- ##
x = seq(0, .5, len = nrow(ex.spec))
y = seq(0, ovr*nseg, len = ncol(ex.spec))
persp(x, y, ex.spec, zlab="Power", xlab="frequency",
ylab="time", ticktype = "detailed", theta=25, d=2)
# -- uses the waveslim package --#
library(waveslim)
eqexp = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
eq = eqexp[,1] # the earthquake series
eq.dwt = dwt(eq, n.levels=6)
# -- plot the dwt and calculate TP -- #
TP = matrix(0,7,1)
par(mfcol=c(7,1), pty="m", mar=c(3,4,2,2))
for(i in 1:6){
plot.ts(up.sample(eq.dwt[[i]], 2^i), type="h", axes=F,
ylab=names(eq.dwt)[i])
abline(h=0)
axis(side=2)
TP[i]=sum(eq.dwt[[i]]^2)
}
plot.ts(up.sample(eq.dwt[[7]], 2^6), type="h", axes=F,
ylab=names(eq.dwt)[7])
abline(h=0)
axis(side=2)
axis(side=1)
TP[7]=sum(eq.dwt[[7]]^2)
TP/sum(eq^2) # the energy distribution
# To plot the ordinates of the wavelet transforms on the
# same scale, include a command like ylim=c(-1.5,1.5) in
# each plot.ts() command.
# -- uses the waveslim package --# # -- assumes eq is retained from previous example --# library(waveslim) eq.dwt = dwt(eq, n.levels=6) eq.trsh = universal.thresh(eq.dwt, hard=F) eq.smo = idwt(eq.trsh) par(mfrow=c(3,1)) plot.ts(eq, ylab="Earthquake", ylim=c(-.5,.5)) plot.ts(eq.smo,ylab="Smoothed Earthquake ",ylim=c(-.5,.5)) plot.ts(eq-eq.smo, ylab="Noise", ylim=c(-.5,.5))
Example 4.23
# -- uses the fracdiff package --# library(fracdiff) lvarve = log(varve)-mean(log(varve)) varve.fd = fracdiff(lvarve, nar=0, nma=0, M=30) varve.fd$d varve.fd$stderror.dpq
# -- uses the tseries package --#
library(tseries)
gnp96 = read.table("/mydata/gnp96.dat")
gnpr = diff(log(gnp96[,2]))
gnpr.ar = ar.mle(gnpr, order=1) # recall phi1 =.347
y = gnpr.ar$resid[-1] # first resid is NA
arch.y = garch(y, order=c(0,1))
summary(arch.y)
NB: The code for Examples 5.3 and 5.4 are slightly
different than what is in text because the
tseries package changed. Also, it may be possible
to fit simultaneous [ARMA mean and GARCH variance] models using garchFit()
from the package fSeries. Unfortunately, I couldn't get it to work with this example (not that
I tried very hard).
# -- uses the tseries package --#
library(tseries)
nyse = scan("/mydata/nyse.dat")
nyse.g = garch(nyse, order=c(1,1))
summary(nyse.g)
u = predict(nyse.g)
plot(800:1000, nyse[800:1000], type="l", xlab="Time", ylab="NYSE Returns")
lines(u[,1], col="blue", lty="dashed")
lines(u[,2], col="blue", lty="dashed")
library(nlme)
mort = scan("/mydata/cmort.dat")
temp = scan("/mydata/temp.dat")
part = scan("/mydata/part.dat")
temp = temp - mean(temp)
temp2 = temp^2
trend = 1:length(mort)
# -- lm() fit in Example 2.2 indicates AR(2) for the residuals
fit.gls = gls(mort~trend + temp + temp2 + part, correlation=corARMA(p=2), method="ML")
# -- the results (abbreviated)
summary(fit.gls)
Correlation Structure: ARMA(2,0)
Parameter estimate(s):
Phi1 Phi2
0.3848530 0.4326282
Coefficients: Value Std.Error t-value p-value
(Intercept) 87.63747 2.7327956 32.06880 0.0000
trend -0.02915 0.0081705 -3.56781 0.0004
temp -0.01881 0.0430264 -0.43715 0.6622
temp2 0.01542 0.0020289 7.60244 0.0000
part 0.15437 0.0248071 6.22297 0.0000
Residual standard error: 7.699336
Degrees of freedom: 508 total; 503 residual
# -- check noise for whiteness (* see note below)
w = filter(residuals(fit.gls), filter=c(1,-.3848530, -.4326282), sides=1)
w = w[3:508] # first two are NA
Box.test(w, 20, type="Ljung") # Ljung-Box-Pierce Statistic
X-squared = 26.793, df = 20, p-value = 0.1412
pchisq(26.793, 18, lower=F) # the p-value (they are resids from an ar2 fit)
0.08296072
* The errors from the gls fit, say e(t), are estimated to be of the form
e(t) = .3848530 e(t-1) + .4326282 e(t-2) + w(t). The filter statement gets at the w(t), which should be
white.
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)
dmort = resid(lm(mort~time(mort))) # detrend series first
dtemp = resid(lm(temp~time(temp)))
dpart = resid(lm(part~time(part)))
lag.max = 10
BIC = matrix(NA,lag.max,1)
fit = vector("list",lag.max)
n = length(mort)
k = 3 # number of series
for (p in 1:lag.max){
fit[[p]] = ar(cbind(dmort,dtemp,dpart), order=p)
BIC[p] = log(det(fit[[p]]$var.pred)) + (k^2*p*log(n)/n)
}
BIC # Example 5:10 - view the BICs, note p=2 is min, and
fit[[2]] # view the parameter estimates for the min BIC model- VAR(2)
#-- for Example 5.9, look at fit[[1]]
#-- if you look at
fit[[lag.max]]$aic
#-- you will see all the AICs and you'll notice the min is at p=8