Time Series I: Basics

rm(list=ls())
Fact = function(n) if (n == 1) 1 else n * Fact(n - 1)
Fact(5)
data(AirPassengers)
AP = AirPassengers
AP
class(AP)
start(AP)
end(AP)
frequency(AP)
summary(AP)
plot(AP, ylab = "Passengers (1000's)")
layout(1:2)
plot(aggregate(AP))
boxplot(AP ~ cycle(AP))
www = "http://elena.aut.ac.nz/~pcowpert/ts/Maine.dat"
Maine.month = read.table(www, header = TRUE)
attach(Maine.month)
class(Maine.month)
Maine.month.ts = ts(unemploy, start = c(1996, 1), 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)
plot(Maine.Feb)
Maine.Aug = window(Maine.month.ts, start = c(1996,8), freq = TRUE)
plot(Maine.Aug)
Feb.ratio = mean(Maine.Feb) / mean(Maine.month.ts)
Aug.ratio = mean(Maine.Aug) / mean(Maine.month.ts)
Feb.ratio
Aug.ratio
layout(1:1)
www = "http://elena.aut.ac.nz/~pcowpert/ts/USunemp.dat"
US.month = read.table(www, header = T)
attach(US.month)
US.month.ts = ts(USun, start=c(1996,1), end=c(2006,10), freq = 12)
plot(US.month.ts, ylab = "unemployed (%)")
# Multiple time series
www = "http://elena.aut.ac.nz/~pcowpert/ts/cbe.dat"
CBE = read.table(www, header = T)
CBE[1:4, ]
class(CBE)
Elec.ts = ts(CBE[, 3], start = 1958, freq = 12)
Beer.ts = ts(CBE[, 2], start = 1958, freq = 12)
Choc.ts = ts(CBE[, 1], start = 1958, freq = 12)
plot(cbind(Elec.ts, Beer.ts, Choc.ts))
length(AP)
length(Elec.ts)
AP.elec = ts.intersect(AP, Elec.ts)
class(AP.elec)
AP.elec
AP
Elec.ts
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")
layout(1:1)
plot(as.vector(AP), as.vector(Elec),
     xlab = "Air passengers / 1000's",
     ylab = "Electricity production / MWh")
abline(reg = lm(Elec ~ AP))
cor(AP, Elec)
www = "http://elena.aut.ac.nz/~pcowpert/ts/pounds_nz.dat"
Z = read.table(www, header = T)
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)" )

www = "http://elena.aut.ac.nz/~pcowpert/ts/global.dat"
Global = scan(www)
Global.ts = ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.annual = aggregate(Global.ts, FUN = mean)
plot(Global.ts)
plot(Global.annual)
layout(1:1)
New.series = window(Global.ts, start=c(1970, 1), end=c(2005, 12))
New.time = time(New.series)
plot(New.series)
abline(reg=lm(New.series ~ New.time))

# decomposition in R
plot(decompose(Elec.ts))
Elec.decom = decompose(Elec.ts, type = "mult")
plot(Elec.decom)
Trend = Elec.decom$trend
Seasonal = Elec.decom$seasonal
ts.plot(cbind(Trend, Trend * Seasonal), lty = 1:2)

www = "http://elena.aut.ac.nz/~pcowpert/ts/Herald.dat"
Herald.dat = read.table(www, header = T)
attach (Herald.dat)
x = CO; y = Benzoa; n = length(x)
sum((x - mean(x))*(y - mean(y))) / (n - 1)
mean((x - mean(x)) * (y - mean(y)))
cov(x, y)
cov(x,y) / (sd(x)*sd(y))
cor(x,y)
www = "http://elena.aut.ac.nz/~pcowpert/ts/wave.dat"
wave.dat = read.table (www, header=T)
attach(wave.dat)
plot(ts(waveht))
plot(ts(waveht[1:60]))
# acf
acf(waveht)$acf[2]
acf(waveht)$acf[1]
# autocovariance
acf(waveht, type = c("covariance"))$acf[2]
# correlogram
acf(waveht)
acf(AirPassengers)
data(AirPassengers)
AP = AirPassengers
AP.decom = decompose(AP, "multiplicative")
plot(ts(AP.decom$random[7:138]))
acf(AP.decom$random[7:138])
sd(AP[7:138])
sd(AP[7:138] - AP.decom$trend[7:138])
sd(AP.decom$random[7:138])
www = "http://elena.aut.ac.nz/~pcowpert/ts/Fontdsdt.dat"
Fontdsdt.dat = read.table(www, header=T)
attach(Fontdsdt.dat)
plot(ts(adflow), ylab = 'adflow')
acf(adflow, xlab = 'lag (months)', main="")
Build.dat = read.table("http://elena.aut.ac.nz/~pcowpert/ts/ApprovActiv.dat", header=T)
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))

# correlogram and cross-correlogram
acf(ts.union(App.ts, Act.ts))
print(ts.union(App.ts, Act.ts))
print(acf(ts.union(App.ts, Act.ts)))

# Bass model
T79 = 1:10
Tdelt = (1:100) / 10
Sales = c(840,1470,2110,4000, 7590, 10950, 10530, 9470, 7790, 5890)
Cusales = cumsum(Sales)
Bass.nls = nls(Sales ~ M * ( ((P+Q)^2 / P) * exp(-(P+Q) * T79) ) /
                  (1+(Q/P)*exp(-(P+Q)*T79))^2, start = list(M=60630, P=0.03, Q=0.38))
summary(Bass.nls)
Bcoef = coef(Bass.nls)
m = Bcoef[1]
p = Bcoef[2]
q = Bcoef[3]
ngete = exp(-(p+q) * Tdelt)
Bpdf = m * ( (p+q)^2 / p ) * ngete / (1 + (q/p) * ngete)^2
plot(Tdelt, Bpdf, xlab = "Year from 1979",
     ylab = "Sales per year", type='l')
points(T79, Sales)
Bcdf = m * (1 - ngete)/(1 + (q/p)*ngete)
plot(Tdelt, Bcdf, xlab = "Year from 1979",
     ylab = "Cumulative sales", type='l')
points(T79, Cusales)

?HoltWinters
require(graphics)
## Seasonal Holt-Winters
(m <- HoltWinters(co2))
plot(m)
plot(fitted(m))
(m <- HoltWinters(AirPassengers, seasonal = "mult"))
plot(m)
## Non-Seasonal Holt-Winters
x <- uspop + rnorm(uspop, sd = 5)
m <- HoltWinters(x, gamma = FALSE)
plot(m)
## Exponential Smoothing
m2 <- HoltWinters(x, gamma = FALSE, beta = FALSE)
lines(fitted(m2)[,1], col = 3)