############################ # 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 = "http://people.stat.sc.edu/hansont/stat520/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") # Tim added: also try "default" prediction plots, e.g. plot(huron.ar1.fit,n.ahead=20,xlab="year",ylab="feet") plot(huron.ar1.fit,n.ahead=20,n1=1950,xlab="year",ylab="feet") # Example 9.4 # Gota River discharge data (1807-1956) # MA(1) forecasting # Page 246 gota = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/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") # Tim added: the "default" prediction plot is fine... plot(gota.ma1.fit,n.ahead=10,col='red',type='b',pch=16,n1=1890,ylab="Water discharge rate",xlab="Year") # Example 9.5 # Bovine blood sugar data (176 observations) # ARMA(1,1) forecasting # Page 251 cows= ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/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") # Tim added: Naive approach is "quick and dirty" plot(ima11.log.oil.fit,n.ahead=12,transform=exp,pch=16,n1=c(1998,1),ylab="Oil prices",xlab="Year") # Example 9.7 # USC Fall enrollment data (1954-2010) # ARI(1,1) forecasting # Page 256 enrollment = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/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))