############################

# Author: Joshua M. Tebbs  #

# Date: 20 Dec 2009        #

# Update: 25 Jul 2011      #

# Purpose: STAT 520 R code #

# CHAPTER 9                #

############################

# Examples 9.1 and 9.2

# Global temperature and beer sales data

# Fitted trend models added

# See Chapter 3 notes for R code

# Page 234-235

# Example 9.3

# Lake Huron elevation data (1880-2006)

# AR(1) forecasting

# Page 241

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

# Fit AR(1) model using ML

huron.ar1.fit <- arima(huron,order=c(1,0,0),method='ML')

huron.ar1.fit

# Obtain MMSE forecasts

huron.ar1.predict <- predict(huron.ar1.fit,n.ahead=20)

round(huron.ar1.predict$pred,3)

round(huron.ar1.predict$se,3)

# Create lower and upper prediction interval bounds

lower.pi<-huron.ar1.predict$pred-qnorm(0.975,0,1)*huron.ar1.predict$se

upper.pi<-huron.ar1.predict$pred+qnorm(0.975,0,1)*huron.ar1.predict$se

# Display prediction intervals (2007-2026)

data.frame(Year=c(2007:2026),lower.pi,upper.pi)

# Original series starts at Year = 1880

# Note: Argument n1=1940 starts plot at 1940

# Note: Argument pch=16 produces a small black circle (for MMSE forecasts)

plot(huron.ar1.fit,n.ahead=20,col='red',type='b',pch=16,n1=1940,ylab="Elevation level (in feet)",xlab="Year")

# Put horizontal line at ML estimate of overall mean

abline(h=coef(huron.ar1.fit)[names(coef(huron.ar1.fit))=='intercept'])

# Put prediction interval lines on plot (darker than default)

lines(y=lower.pi,x=2007:2026,lwd=2,col="red",lty="dashed")

lines(y=upper.pi,x=2007:2026,lwd=2,col="red",lty="dashed")

# Example 9.4

# Gota River discharge data (1807-1956)

# MA(1) forecasting

# Page 246

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

gota.ma1.fit = arima(gota,order=c(0,0,1),method='ML')

gota.ma1.fit

# Obtain MMSE forecasts

gota.ma1.predict <- predict(gota.ma1.fit,n.ahead=10)

round(gota.ma1.predict$pred,3)

round(gota.ma1.predict$se,3)

# Create lower and upper prediction interval bounds

lower.pi<-gota.ma1.predict$pred-qnorm(0.975,0,1)*gota.ma1.predict$se

upper.pi<-gota.ma1.predict$pred+qnorm(0.975,0,1)*gota.ma1.predict$se

# Display prediction intervals (1957-1966)

data.frame(Year=c(1957:1966),lower.pi,upper.pi)

# Original series starts at Year = 1807

# Note: Argument n1=1890 starts plot at 1890

# Note: Argument pch=16 produces a small black circle (for MMSE forecasts)

plot(gota.ma1.fit,n.ahead=10,col='red',type='b',pch=16,n1=1890,ylab="Water discharge rate",xlab="Year")

# Put horizontal line at ML estimate of overall mean

abline(h=coef(gota.ma1.fit)[names(coef(gota.ma1.fit))=='intercept'])

# Put prediction interval lines on plot (darker than default)

lines(y=lower.pi,x=1957:1966,lwd=2,col="red",lty="dashed")

lines(y=upper.pi,x=1957:1966,lwd=2,col="red",lty="dashed")

# Example 9.5

# Bovine blood sugar data (176 observations)

# ARMA(1,1) forecasting

# Page 251

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

cows.arma11.fit = arima(cows,order=c(1,0,1),method='ML')

cows.arma11.fit

# Obtain MMSE forecasts

cows.arma11.predict <- predict(cows.arma11.fit,n.ahead=10)

round(cows.arma11.predict$pred,3)

round(cows.arma11.predict$se,3)

# Create lower and upper prediction interval bounds

lower.pi<-cows.arma11.predict$pred-qnorm(0.975,0,1)*cows.arma11.predict$se

upper.pi<-cows.arma11.predict$pred+qnorm(0.975,0,1)*cows.arma11.predict$se

# Display prediction intervals (Days 177-186)

data.frame(Day=c(177:186),lower.pi,upper.pi)

# Original series starts at Day = 1

# Note: Argument n1=80 starts plot at Day = 81

# Note: Argument pch=16 produces a small black circle (for MMSE forecasts)

plot(cows.arma11.fit,n.ahead=10,col='red',type='b',pch=16,n1=81,ylab="Blood sugar level (mg/100ml blood)",xlab="Day")

# Put horizontal line at ML estimate of overall mean

abline(h=coef(cows.arma11.fit)[names(coef(cows.arma11.fit))=='intercept'])

# Put prediction interval lines on plot (darker than default)

lines(y=lower.pi,x=177:186,lwd=2,col="red",lty="dashed")

lines(y=upper.pi,x=177:186,lwd=2,col="red",lty="dashed")

# Example 9.6

# Oil price data (1/86-1/06)

# IMA(1,1) forecasting

# Log transformed

# Page 255

data(oil.price)

ima11.log.oil.fit = arima(log(oil.price),order=c(0,1,1),method='ML')

ima11.log.oil.fit

# Obtain MMSE forecasts (on log scale)

ima11.log.oil.predict <- predict(ima11.log.oil.fit,n.ahead=12)

round(ima11.log.oil.predict$pred,3)

round(ima11.log.oil.predict$se,3)

# Create lower and upper prediction interval bounds

lower.pi<-ima11.log.oil.predict$pred-qnorm(0.975,0,1)*ima11.log.oil.predict$se

upper.pi<-ima11.log.oil.predict$pred+qnorm(0.975,0,1)*ima11.log.oil.predict$se

# Display prediction intervals (12 months ahead)

year.temp = c(2006.083,2006.166,2006.250,2006.333,2006.416,2006.500,2006.583,2006.666,2006.750,2006.833,2006.916,2007)

data.frame(Month=year.temp,lower.pi,upper.pi)

# Original series starts at Month = Jan 1986

# Note: Argument n1=c(1998,1) starts plot at Time = Jan 1998

# Note: Argument pch=16 produces a small black circle (MMSE prediction)

plot(ima11.log.oil.fit,n.ahead=12,col='red',type='b',pch=16,n1=c(1998,1),ylab="Oil prices (on log scale)",xlab="Year")

# Put prediction interval lines on plot (darker than default)

lines(y=lower.pi,x=year.temp,lwd=2,col="red",lty="dashed")

lines(y=upper.pi,x=year.temp,lwd=2,col="red",lty="dashed")

# Example 9.7

# USC Fall enrollment data (1954-2010)

# ARI(1,1) forecasting

# Page 256

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

enrollment.ari11.fit = arima(enrollment,order=c(1,1,0),method='ML')

enrollment.ari11.fit

# Obtain predictions

enrollment.ari11.predict <- predict(enrollment.ari11.fit,n.ahead=10)

round(enrollment.ari11.predict$pred,3)

round(enrollment.ari11.predict$se,3)

# Create lower and upper prediction interval bounds

lower.pi<-enrollment.ari11.predict$pred-qnorm(0.975,0,1)*enrollment.ari11.predict$se

upper.pi<-enrollment.ari11.predict$pred+qnorm(0.975,0,1)*enrollment.ari11.predict$se

# Display prediction intervals (10 years ahead)

data.frame(Year=c(2011:2020),lower.pi,upper.pi)

# Original series starts at 1954

# Note: Argument n1=1974 starts plot at 1974

# Note: Argument pch=16 produces a small black circle (MMSE prediction)

plot(enrollment.ari11.fit,n.ahead=10,col='red',type='b',pch=16,n1=1974,ylab="USC Columbia fall enrollment",xlab="Year")

# Put prediction interval lines on plot (darker than default)

lines(y=lower.pi,x=c(2011:2020),lwd=2,col="red",lty="dashed")

lines(y=upper.pi,x=c(2011:2020),lwd=2,col="red",lty="dashed")

# Example 9.8

# Global temperature data

# Linear trend model fit

# See Chapter 3 notes for R code

# Page 260

# Example 9.9

# Lake Huron elevation data (1880-2006)

# Prediction limits

# See prediction limit code in Example 9.3

# Page 262

# Example 9.10

# Oil price data (1/86-1/06)

# IMA(1,1) forecasting

# Log transformed

# Prediction limits

# Page 265

data(oil.price)

ima11.log.oil.fit = arima(log(oil.price),order=c(0,1,1),method='ML')

# MMSE forecasts on log scale

ima11.log.oil.predict <- predict(ima11.log.oil.fit,n.ahead=12)

round(ima11.log.oil.predict$pred,3)

round(ima11.log.oil.predict$se,3)

# MMSE forecasts back-transformed (to original scale)

oil.price.predict <- round(exp(ima11.log.oil.predict$pred + (1/2)*(ima11.log.oil.predict$se)^2),3)

oil.price.predict

# Compute prediction intervals (on original scale)

lower.pi<-ima11.log.oil.predict$pred-qnorm(0.975,0,1)*ima11.log.oil.predict$se

upper.pi<-ima11.log.oil.predict$pred+qnorm(0.975,0,1)*ima11.log.oil.predict$se

year.temp = c(2006.083,2006.166,2006.250,2006.333,2006.416,2006.500,2006.583,2006.666,2006.750,2006.833,2006.916,2007)

data.frame(Month=year.temp,lower.pi=exp(lower.pi),upper.pi=exp(upper.pi))