R: FORECASTING AND TIME SERIES # 1

# Example 1.1.

# Page 2

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

plot(globaltemps,ylab="Global temperature deviations",xlab="Year",type="o")

# Example 1.2.

# Page 3

data(milk)

plot(milk,ylab="Amount of milk produced",xlab="Year",type="o")

# Example 1.3.

# Page 4

data(CREF)

plot(CREF,ylab="CREF stock values",type="o")

# Example 1.4.

# Page 5

homeruns <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\homeruns.txt"),start=1909)

plot(homeruns,ylab="Number of homeruns",xlab="Year",type="o")

# Example 1.5.

# Page 6

earthquake <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\earthquake.txt"),start=1900)

plot(earthquake,ylab="Number of earthquakes (7.0 or greater)",xlab="Year",type="o")

# Example 1.6.

# Page 7

enrollment <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\enrollment.txt"),start=1954)

plot(enrollment,ylab="USC Columbia fall enrollment",xlab="Year",type="o")

# Example 1.7.

# Page 8

data(star)

plot(star,ylab="Star brightness",type="o")

# Example 1.8.

# Page 9

data(airmiles)

plot(airmiles,ylab="Airline miles",xlab="Year",type="o")

# Example 1.9.

# Page 10

sp500 <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\sp500.txt"))

plot(sp500,ylab="SP500 Index",xlab="Time",type="o")

# Example 1.10.

# Page 11

ventilation <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\ventilation.txt"))

plot(ventilation,ylab="Ventilation (L/min)",xlab="Observation time",type="o")

# Example 1.11.

# Page 12

exchangerate <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\exchangerate.txt"),start=1980,frequency=52)

plot(exchangerate,ylab="British pounds",xlab="Year",type="o")

# Example 1.12.

# Page 13

data(oil.price)

plot(oil.price,ylab="Oil prices",xlab="Year",type="o")

# Example 1.13.

# Page 14

data(larain)

plot(larain,ylab="LA rainfall amounts",xlab="Year",type="o")

# Example 1.14.

# Page 15

brick <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\brick.txt"),start=1956,freq=4)

plot(brick,ylab="Australian clay brick production (in millions)",xlab="Time",type="o")

# Example 1.15.

# Page 16

supremecourt <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\supremecourt.txt"),start=1926)

plot(supremecourt,ylab="Percentage granted review",xlab="Time",type="o")

# Example 1.8.

# Add special plotting symbols for months

# Page 19

data(airmiles)

plot(airmiles,ylab="Airline miles",xlab="Year",type='l')

points(y=airmiles,x=time(airmiles),pch=as.vector(season(airmiles)),cex=1)

R: FORECASTING AND TIME SERIES # 2

# Example 2.1.

# Page 30

white.noise <- rnorm(150,0,1)

plot(white.noise,ylab="Simulated white noise process",xlab="Time",type="o")

# Example 2.2.

# Uses white noise process from Example 2.1.

# Page 33

random.walk <- white.noise*0

for(i in 1:length(white.noise)){

      random.walk[i]<-sum(white.noise[1:i])

      }

plot(random.walk,ylab="Simulated random walk process",xlab="Time",type="o")

# Example 2.3.

# Uses white noise process from Example 2.1.

# Page 36

moving.average <- filter(x = white.noise, filter = rep(x = 1/3, times = 3), method = "convolution", sides = 1)

plot(moving.average,ylab="Simulated moving average process",xlab="Time",type="o")

# Example 2.4.

# Autoregressive model simulation

# Page 37

autoregressive <- arima.sim(model = list(ar = c(0.75)), n = 150, rand.gen = rnorm, sd = 1)

plot(autoregressive,ylab="Simulated autoregressive process",xlab="Time",type="o")

# Example 2.5.

# Sinusoidal process

# Page 38

mean.function <- 2*sin(2*pi*1:156/52+0.6*pi)

w <- rnorm(156,0,1)

par(mfrow=c(2,2))

plot(mean.function,ylab="",xlab="Time",type="l")

plot(mean.function+w,ylab="",xlab="Time",type="o")

plot(mean.function+2*w,ylab="",xlab="Time",type="o")

plot(mean.function+4*w,ylab="",xlab="Time",type="o")

R: FORECASTING AND TIME SERIES # 3

# Example 3.4.

# Global temperature data from 1900

# Page 53

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

plot(globaltemps.1900,ylab="Global temperature deviations (since 1900)",xlab="Year",type="o")

# Example 3.4.

# Global temperature data from 1900

# Straight line model fit with output

# Page 54

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

summary(fit)

plot(globaltemps.1900,ylab="Global temperature deviations (since 1900)",xlab="Year",type="o")

abline(fit)

# Example 3.4.

# Global temperature data from 1900

# Residuals from straight line model fit

# Page 55

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

plot(resid(fit),ylab="Residuals",xlab="Year",type="o")

# Example 3.4.

# Global temperature data from 1900

# Plot of first data differences

# Page 57

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

plot(diff(globaltemps.1900),ylab="Global temperature deviation differences",xlab="Year",type="o")

# Example 3.5.

# Gold price data

# Page 59

data(gold)

plot(gold,ylab="Price",xlab="Time",type="o")

# Example 3.5.

# Gold price data

# Quadratic model fit with output

# Page 59-60

data(gold)

t <- time(gold)

t.sq <- t^2

fit <- lm(gold ~ t + t.sq)

summary(fit)

plot(gold,ylab="Price",xlab="Time",type="o")

curve(expr = fit$coefficients[1] +

             fit$coefficients[2]*x +

             fit$coefficients[3]*x^2, col = "black", 

             lty = "solid", lwd = 1, add = TRUE)

# Example 3.5.

# Gold price data

# Residuals from quadratic model fit

# Page 61

data(gold)

t <- time(gold)

t.sq <- t^2

fit <- lm(gold ~ t + t.sq)

plot(resid(fit),ylab="Residuals",xlab="Time",type="o")

# Example 3.6.

# Beer sales data

# Monthly plotting symbols added

# Regression output included

# Page 62-63

data(beersales)

b.sales<-window(beersales,start=1980)

plot(b.sales,ylab="Sales",xlab="Year",type='l')

points(y=b.sales,x=time(b.sales),pch=as.vector(season(b.sales)))

month.<-season(b.sales)

fit<-lm(b.sales ~ month.-1)

summary(fit)

# Example 3.6.

# Beer sales data

# Residuals from seasonal means model fit

# Page 64

data(beersales)

b.sales<-window(beersales,start=1980)

fit<-lm(b.sales ~ month.-1)

plot(resid(fit),ylab="Residuals",xlab="Time",type="o")

# Example 3.6.

# Beer sales data

# Cosine trend model fit and output

# Page 66-67

data(beersales)

b.sales<-window(beersales,start=1980)

har. <- harmonic(b.sales,1)

fit <- lm(b.sales~har.)

summary(fit)

plot(ts(fitted(fit),freq=12,start=c(1980,1)),ylab="Sales",type='l',ylim=range(c(fitted(fit),b.sales)))

points(b.sales)

# Example 3.6.

# Beer sales data

# Residuals from cosine trend model fit

data(beersales)

# Page 68

b.sales<-window(beersales,start=1980)

har. <- harmonic(b.sales,1)

fit <- lm(b.sales~har.)

plot(resid(fit),ylab="Residuals",xlab="Time",type="o")

# Example 3.4.

# Global temperature data from 1900

# Standardised residuals from straight line model fit

# Histogram and qq plot

# Figures were constructed separately

# Page 72

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

hist(rstudent(fit),main="Histogram of standardized residuals",xlab="Standardized residuals")

qqnorm(rstudent(fit),main="QQ plot of standardized residuals")

# Example 3.4.

# Global temperature data from 1900

# Shapiro-Wilk test for normality

# Page 73

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

shapiro.test(rstudent(fit))

# Example 3.4.

# Global temperature data from 1900

# Standardised residuals from straight line model fit

# Horizontal line added at 0

# Page 74

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

plot(rstudent(fit),ylab="Residuals",xlab="Year",type="o")

abline(h=0)

# Example 3.4.

# Global temperature data from 1900

# Runs test for independence on standardised residuals

# Page 75

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

runs(rstudent(fit))

# Example 3.4.

# Global temperature data from 1900

# Calculate sample ACF for standardised residuals

# Page 78

globaltemps <- ts(read.table(file = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\globaltemps.txt"),start=1856)

globaltemps.1900 <- window(globaltemps,start=1900)

fit <- lm(globaltemps.1900~time(globaltemps.1900))

acf(rstudent(fit),main="Sample ACF for standardized residuals")

# Simulation exercise

# Simulate white noise processes and plot ACFs

# Page 79

w.n.1 <- rnorm(100,0,1)

w.n.2 <- rnorm(100,0,1)

par(mfrow=c(2,2))

plot(w.n.1,ylab="White noise process.1",xlab="Time",type="o")

acf(w.n.1,main="Sample ACF")

plot(w.n.2,ylab="White noise process.2",xlab="Time",type="o")

acf(w.n.2,main="Sample ACF")

R: FORECASTING AND TIME SERIES # 4

# Example 4.1.

# MA(1) simulation

# Important! R's convention is to use positive thetas for MA models (so we have to negate)

# E.g., ma = 0.9 means theta = -0.9.

# Page 83

ma.sim <- arima.sim(list(order = c(0,0,1), ma = 0.9), n = 100)

par(mfrow=c(2,2))

plot(ma.sim,ylab=expression(Y[t]),xlab="Time",type="o",main="MA(1) simulation")

acf(ma.sim,main="Sample ACF")

plot(y=ma.sim,x=zlag(ma.sim,1),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p',main="Lag 1 scatterplot")

plot(y=ma.sim,x=zlag(ma.sim,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p',main="Lag 2 scatterplot")

# Example 4.1.

# MA(1) simulation

# Important! R's convention is to use positive thetas for MA models (so we have to negate)

# E.g., ma = -0.9 means theta = 0.9.

# Page 84

ma.sim <- arima.sim(list(order = c(0,0,1), ma = -0.9), n = 100)

par(mfrow=c(2,2))

plot(ma.sim,ylab=expression(Y[t]),xlab="Time",type="o",main="MA(1) simulation")

acf(ma.sim,main="Sample ACF")

plot(y=ma.sim,x=zlag(ma.sim,1),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p',main="Lag 1 scatterplot")

plot(y=ma.sim,x=zlag(ma.sim,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p',main="Lag 2 scatterplot")

# Example 4.2.

# MA(2) simulation

# Important! R's convention is to use positive thetas for MA models (so we have to negate)

# Page 87

ma.sim <- arima.sim(list(order = c(0,0,2), ma = c(-0.9,0.7)), n = 100)

par(mfrow=c(2,2))

plot(ma.sim,ylab=expression(Y[t]),xlab="Time",type="o",main="MA(2) simulation")

acf(ma.sim,main="Sample ACF")

plot(y=ma.sim,x=zlag(ma.sim,1),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p',main="Lag 1 scatterplot")

plot(y=ma.sim,x=zlag(ma.sim,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p',main="Lag 2 scatterplot")

# Example 4.3.

# True AR(1) autocorrelation functions (out to k=20 lags)

# Uses ARMAacf function

# ARMAacf function includes the k=0 lag

# Use y = y[2:21] to remove k=0 lag from ARMAacf output

# Page 91

par(mfrow=c(2,2))

y = ARMAacf(ar = c(0.9), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(-0.9), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(0.5), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(-0.5), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

# Example 4.3.

# AR(1) simulations

# Page 92

par(mfrow=c(2,2))

ar.sim.1 <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = 100)

ar.sim.2 <- arima.sim(list(order = c(1,0,0), ar = -0.9), n = 100)

ar.sim.3 <- arima.sim(list(order = c(1,0,0), ar = 0.5), n = 100)

ar.sim.4 <- arima.sim(list(order = c(1,0,0), ar = -0.5), n = 100)

plot(ar.sim.1,ylab=expression(Y[t]),xlab="Time",type="o")

plot(ar.sim.2,ylab=expression(Y[t]),xlab="Time",type="o")

plot(ar.sim.3,ylab=expression(Y[t]),xlab="Time",type="o")

plot(ar.sim.4,ylab=expression(Y[t]),xlab="Time",type="o")

# AR(1) sample ACFs

# Page 93

par(mfrow=c(2,2))

acf(ar.sim.1,main="Sample ACF")

acf(ar.sim.2,main="Sample ACF")

acf(ar.sim.3,main="Sample ACF")

acf(ar.sim.4,main="Sample ACF")

# Figure 4.7.

# AR(2) stationarity region

# Page 97

par(xaxs = "i", yaxs = "i")

plot(x = -2, y = -1, xlim = c(-2,2), ylim = c(-1,1),

    type = "n", frame.plot = FALSE, xlab = expression(phi[1]), ylab = expression(phi[2]))

#abline() draws y = mx + b      

abline(a = 1, b = -1) #Draw line of phi2 = 1 – phi1

abline(a = 1, b =  1) #Draw line of phi2 = 1 + phi1

#Plot the phi1 and phi2 values

points(x = 0, y = 0.5, pch = 1, col = "red")

points(x = 0, y = -0.5, pch = 2, col = "darkgreen")

points(x = -0.2, y = -0.5, pch = 2, col = "darkgreen")

points(x = -1, y = -0.5, pch = 2, col = "darkgreen")

points(x = -1.8, y = -0.9, pch = 2, col = "darkgreen")

points(x = 0.5, y = 0.25, pch = 1, col = "red")

points(x = 1.8, y = 0.9, pch = 3, col = "blue")

points(x = -1.2, y = 0.8, pch = 3, col = "blue")

legend(locator(1), legend = c("Real roots",

      "Complex roots", "Outside stationarity region"), pch =

      c(1,2,3), col = c("red", "darkgreen", "blue"), cex =

      0.75, bty = "n")

# Example 4.4.

# True AR(2) autocorrelation functions (out to k=20 lags)

# Uses ARMAacf function

# ARMAacf function includes the k=0 lag

# Use y = y[2:21] to remove k=0 lag from ARMAacf output

# Page 100

par(mfrow=c(2,2))

y = ARMAacf(ar = c(0.5,-0.5), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(1.1,-0.3), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(-0.5,0.25), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(1,-0.5), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

# Example 4.4.

# AR(2) simulations

# Page 101

par(mfrow=c(2,2))

ar2.sim.1 <- arima.sim(list(order = c(2,0,0), ar = c(0.5,-0.5)), n = 100)

ar2.sim.2 <- arima.sim(list(order = c(2,0,0), ar = c(1.1,-0.3)), n = 100)

ar2.sim.3 <- arima.sim(list(order = c(2,0,0), ar = c(-0.5,-0.25)), n = 100)

ar2.sim.4 <- arima.sim(list(order = c(2,0,0), ar = c(1,-0.5)), n = 100)

plot(ar2.sim.1,ylab=expression(Y[t]),xlab="Time",type="o")

plot(ar2.sim.2,ylab=expression(Y[t]),xlab="Time",type="o")

plot(ar2.sim.3,ylab=expression(Y[t]),xlab="Time",type="o")

plot(ar2.sim.4,ylab=expression(Y[t]),xlab="Time",type="o")

# AR(2) sample ACFs

# Page 102

par(mfrow=c(2,2))

acf(ar2.sim.1,main="Sample ACF")

acf(ar2.sim.2,main="Sample ACF")

acf(ar2.sim.3,main="Sample ACF")

acf(ar2.sim.4,main="Sample ACF")

# Figure 4.11.

# True ARMA(1,1) autocorrelation functions (out to k=20 lags)

# Uses ARMAacf function

# ARMAacf function includes the k=0 lag

# Use y = y[2:21] to remove k=0 lag from ARMAacf output

# Page 112

par(mfrow=c(2,2))

y = ARMAacf(ar = 0.9, ma = 0.25, lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = -0.9, ma = 0.25, lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = 0.5, ma = 0.25, lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = -0.5, ma = 0.25, lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

      ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)