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

# 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 = "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

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 = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\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 = "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')

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 = "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')

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 = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\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 = "C:\\Users\\Tebbs\\My Documents\\texfiles\\Classes\\USC\\stat520\\f11\\data\\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 = "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')

# 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