去年从系里借了一本名为 Learning Introductory Time Series with R 的书，觉得很有用，开始很认真地做笔记，后来越来越难，学得步履维艰，最终还是半途而废了。开头几节我亲手敲下的代码笔记，本来想整理一下分开发到博客上，拖来拖去，到现在也没有实现。今天整理电脑又看见了它，干脆全部放在这里。有些事，越拖越没有意义。
后来一查，书里的代码早就被作者全部放到网上了……
http://tur-www1.massey.ac.nz/~pscowper/ts/#Contents
[sourcecode language=”r”]
Notes for Learning Introductory Time Series with R
Peng Zhao
Example 1.1: AirPassengers. Page 4.
data(AirPassengers).
AP <- AirPassengers
class(AP)
start(AP);end(AP);frequency(AP)
plot(AP,ylab=“Passenger (1000’s)“)
layout(1:2)
plot(aggregate(AP)) # Peng: aggregate() returns an aggregated annual series.
boxplot(AP ~ cycle(AP)) # Peng: cycle() returns an seasonal series.
Example 1.2: Maine. Page 7.
www <- “http://www.massey.ac.nz/~pscowper/ts/Maine.dat"
Maine.month <- read.table(www,header = TRUE)
class(Maine.month)
attach(Maine.month) # Peng: attach() and detach() are smart ideas to call Maine.month$unemploy as unemploy for short.
Maine.month.ts <- ts(unemploy, start = c(1996,1), end = c(2001,2), freq=12)
Maine.annual.ts <- aggregate(Maine.month.ts)/12
layout(1:2)
plot(Maine.month.ts, ylab = “unemployed (%)”)
plot(Maine.annual.ts, ylab = “unemployed (%)”)
Maine.Feb <- window(Maine.month.ts, start = c(1996,2), freq = TRUE) # Peng: window() can extract a subset of a time series.
Maine.Aug <- window(Maine.month.ts, start = c(1996,8), freq = TRUE)
Feb.ratio <- mean(Maine.Feb) / mean(Maine.month.ts)
Aug.ratio <- mean(Maine.Aug) / mean(Maine.month.ts)
Example 1.3: Multiple. Page 10.
www <- “http://www.massey.ac.nz/~pscowper/ts/cbe.dat"
CBE <- read.table(www, header = TRUE)
class(CBE)
Elec.ts <- ts(CBE[,3], start = 1958, freq = 12) # Peng: ts() transfer a data frame into a time series.
Beer.ts <- ts(CBE[,2], start = 1958, freq = 12)
Choc.ts <- ts(CBE[,1], start = 1958, freq = 12)
Mult.ts <- cbind(Elec.ts, Beer.ts, Choc.ts) # Peng: cbind() means column bind.
plot(Mult.ts)
AP.elec <- ts.intersect(AP, Elec.ts) # Peng: ts.intersect() produces the intersection(交集) from different time series.
str(AP.elec)
summary(AP.elec)
start(AP.elec);end(AP.elec)
AP <- AP.elec$AP; Elec <- AP.elec$Elec.ts # Peng: error. why? the following is correct.
AP <- AP.elec[,1]; Elec <- AP.elec[,2]
layout(1:2)
plot(AP, main = “”, ylab = “Air passengers / 1000’s”)
plot(Elec, main = “”, ylab = “Electricity production / MkWh”)
plot(as.vector(AP), as.vector(Elec), xlab = “Air passengers / 1000’s”, ylab = “Electricity producton / MWh”) # Peng: why vector?
abline(reg = lm(Elec ~ AP)) # Peng: linear regression.
Example 1.4: Quarterly. Page 14.
www <- “http://www.massey.ac.nz/~pscowper/ts/pounds_nz.dat"
Z <- read.table(www, header = TRUE)
Z.ts <- ts(Z, st = 1991, fr = 4)
plot(Z.ts, xlab = “time / years”, ylab = “Quarterly exchange rate in $NZ / pound”)
Z.92.96 <- window(Z.ts, start = c(1992,1), end = c(1996,1))
Z.96.98 <- window(Z.ts, start = c(1996,1), end = c(1998,1))
layout(1:2)
plot(Z.92.96, ylab = “Exchange rate in $NZ/pound”, xlab = “Time (years)”)
plot(Z.96.98, ylab = “Exchange rate in $NZ/pound”, xlab = “Time (years)”)
Example 1.5: Temperature. Page 16.
www <- “http://www.massey.ac.nz/~pscowper/ts/global.dat"
Global <- scan(www) # Peng: why not read.table() but scan()?
Global.ts <- ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.annual <- aggregate(Global.ts, FUN = mean) # Peng: aggregate() can also produce an annually mean series.
plot(Global.ts)
plot(Global.annual)
abline(reg=lm(Global.ts ~ time(Global.ts)))
New.series <- window(Global.ts, start = c(1970,1),end=c(2005,12))
New.time <- time(New.series) # Peng: time() extract the time stamps.
plot(New.series)
abline(reg=lm(New.series ~ New.time))
Example 1.6: Decomposition. Page 22.
plot(Elec.ts)
Elec.decompose <- decompose(Elec.ts) # Peng: decompose() returns a trend series, a seasonal series and residuals.
plot(Elec.decompose)
Elec.decompose <- decompose(Elec.ts, type = “mult”) # Peng: different decomposing method.
plot(Elec.decompose)
Example 2.1: Covariance. Page 29.
www <- “http://www.massey.ac.nz/~pscowper/ts/Herald.dat"
Herald.dat <- read.table(www, header = TRUE)
attach(Herald.dat)
x <- CO; y <- Benzoa; n <- length(x)
sum((x – mean(x)) * (y – mean(y)))/(n – 1) # Peng: calculation of covariance.
mean((x – mean(x)) * (y – mean(y)))
cov(x,y) # Peng: covariance(协方差). variance(方差) when x and y are the same variable.
sd(x) == (sum((x – mean(x)) ^ 2) / (n -1)) ^ 0.5 # Peng: standard deviation is the square root of variance.
cov(x,y) / (sd(x) * sd(y)) # Peng: calculation of correlation.
cor(x,y) # Peng: correlation(相关系数).
Example 2.2: Autocorrelation. Page 34.
www <- “http://www.massey.ac.nz/~pscowper/ts/wave.dat"
layout(1:2)
plot(ts(waveht)); plot(ts(waveht[1:60]))
waveht.acf <- acf(waveht) # Peng: autocorrelation(自相关) function.
waveht.acf.ck <- acf(waveht)$acf
acf(waveht)$acf[2]# Peng: when k = 1.
waveht.acf
length(waveht.acf)
acf(waveht, type = c(“covariance”))$acf[2]
plot(waveht[1:396], waveht[2:397])
Example 2.3: Correlogram. Page 38.
data(AirPassengers)
AP <- AirPassengers
AP.decom <- decompose(AP, “multiplicative”)
plot(ts(AP.decom$random[7:138]))
acf(AP.decom$random[7:138])
Example 2.4: Font Reservoir. Page 40.
www <- “http://www.massey.ac.nz/~pscowper/ts/Fontdsdt.dat"
Fontdsdt.dat <- read.table(www, header = TRUE)
attach(Fontdsdt.dat)
layout(1:2)
plot(ts(adflow), ylab = ‘adflow’)
acf(adflow, xlab = “lag (months)”, main = “”)
Example 3.1: Approvals. Page 46.
www <- “http://www.massey.ac.nz/~pscowper/ts/ApprovActiv.dat"
Build.dat <- read.table(www, header = TRUE)
attach(Build.dat)
App.ts <- ts(Approvals,start = c(1996,1),freq=4)
Act.ts <- ts(Activity, start = c(1996,1),freq=4)
ts.plot(App.ts, Act.ts, lty = c(1,3)) # Peng: ts.plot() plots two time series in one figure.
Build.union <- ts.union(App.ts, Act.ts) # Peng: ts.union() binds time series with a common frequency, padding with “NA”
to the union of their time coverages.
acf(Build.union) # Peng: correlograms and cross-correlograms(互相关)
app.ran <- decompose(App.ts)$random
app.ran.ts <- window(app.ran, start = c(1996,3))
act.ran <- decompose(Act.ts)$random
act.ran.ts <- window(act.ran, start = c(1996,3))
acf(ts.union(app.ran.ts,act.ran.ts))
ccf(app.ran.ts, act.ran.ts)
print(acf(ts.union(app.ran.ts,act.ran.ts)))
[/sourcecode]