############################
# Author: Joshua M. Tebbs  #
# Date: 20 Dec 2009        #
# Update: 25 Jul 2011      #
# Purpose: STAT 520 R code #
# CHAPTER 8                #
############################

# Example 8.1
# Gota River discharge data
# Residual analysis
# Page 212-213
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
par(mfrow=c(2,2))
plot(gota,ylab="Water discharge rate",xlab="Year",type="o")
plot(rstandard(gota.ma1.fit),xlab="Time",ylab="Standardised residuals",type='o')
abline(h=0)
hist(rstandard(gota.ma1.fit),xlab="Standardised residuals",main="")
qqnorm(rstandard(gota.ma1.fit),main="")
qqline(rstandard(gota.ma1.fit))
# Shapiro-Wilk and runs tests
shapiro.test(rstandard(gota.ma1.fit))
runs(rstandard(gota.ma1.fit))

# Example 8.2
# Lake Huron elevation data
# Residual analysis
# Page 213-215
huron = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/huron.txt"),start=1880)
huron.ar1.fit = arima(huron,order=c(1,0,0),method='ML')
huron.ar1.fit
par(mfrow=c(2,2))
plot(huron,ylab="Elevation level (in feet)",xlab="Year",type="o")
plot(rstandard(huron.ar1.fit),xlab="Time",ylab="Standardised residuals",type='o')
abline(h=0)
hist(rstandard(huron.ar1.fit),xlab="Standardised residuals",main="")
qqnorm(rstandard(huron.ar1.fit),main="")
qqline(rstandard(huron.ar1.fit))
# Shapiro-Wilk and runs tests
shapiro.test(rstandard(huron.ar1.fit))
runs(rstandard(huron.ar1.fit))

# Example 8.3
# Gota River discharge data
# Sample ACF for residuals
# Page 218-219
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')
acf(residuals(gota.ma1.fit),main="Sample ACF for MA(1) residuals")
acf(residuals(gota.ma1.fit),plot=F,lag.max=10)
acf(residuals(arima(gota,order=c(0,0,1))),plot=F,lag.max=10)

# Example 8.4
# Gota River discharge data
# Ljung-Box test for model adequacy
# tsdiag output
# Page 221-222
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')
Box.test(residuals(gota.ma1.fit),lag = 10,type="Ljung-Box",fitdf=1)
tsdiag(gota.ma1.fit,gof=20,omit.initial=F)

# Example 8.5.
# Earthquake data
# tsdiag output
# Page 223
earthquake = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/earthquake.txt"),start=1900)
# Fit ARMA(1,1) model to square-root transformed data
arima(sqrt(earthquake),order=c(1,0,1),method='ML') # maximum likelihood
sqrt.earthquake.arma11.fit = arima(sqrt(earthquake),order=c(1,0,1),method='ML') # maximum likelihood
tsdiag(sqrt.earthquake.arma11.fit,gof=20,omit.initial=F)
# Shapiro-Wilk and runs tests
shapiro.test(rstandard(sqrt.earthquake.arma11.fit))
runs(rstandard(sqrt.earthquake.arma11.fit))

# Example 8.6.
# Supreme Court data
# tsdiag output
# Page 224
supremecourt = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/supremecourt.txt"),start=1926)
# Fit ARIMA(0,1,1) model to the log-transformed data
log.supremecourt.ima11.fit = arima(log(supremecourt),order=c(0,1,1),method='ML') # ML
tsdiag(log.supremecourt.ima11.fit,gof=20,omit.initial=F)
# Shapiro-Wilk and runs tests
shapiro.test(rstandard(log.supremecourt.ima11.fit))
runs(rstandard(log.supremecourt.ima11.fit))

# Example 8.7
# Oil price data
# Residual analysis for log-transformed IMA(1,1) model
# Page 225
data(oil.price)
ima11.log.oil.fit = arima(log(oil.price),order=c(0,1,1),method='ML')
# Page 226
par(mfrow=c(2,2))
plot(diff(log(oil.price)),ylab="1st differences of logarithms",xlab="Year",type="o")
plot(rstandard(ima11.log.oil.fit),xlab="Time",ylab="Standardised residuals",type='o')
abline(h=0)
hist(rstandard(ima11.log.oil.fit),xlab="Standardised residuals",main="")
qqnorm(rstandard(ima11.log.oil.fit),main="")
qqline(rstandard(ima11.log.oil.fit))
## Shapiro-Wilk and runs tests
shapiro.test(rstandard(ima11.log.oil.fit))
runs(rstandard(ima11.log.oil.fit))
# Page 227
tsdiag(ima11.log.oil.fit,gof=20,omit.initial=F)

# Example 8.8
# Gota River discharge data
# Overfitting
# Page 229
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')
# Two overfit models
gota.ma2.overfit = arima(gota,order=c(0,0,2),method='ML')
gota.arma11.overfit = arima(gota,order=c(1,0,1),method='ML')
gota.ma1.fit
gota.ma2.overfit
gota.arma11.overfit