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.
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=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="") }
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=TRUE) 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=FALSE,intercept=TRUE)) # 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=TRUE) # 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=TRUE) 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)
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=FALSE, 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=FALSE)$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=FALSE, 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=FALSE, 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=FALSE) 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=FALSE) # 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