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

# Author: Joshua M. Tebbs  #

# Date: 20 Dec 2009        #

# Update: 25 Jul 2011      #

# Purpose: STAT 520 R code #

# CHAPTER 7                #

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

# Example 7.1.

# Page 180-181

# Function that computes MOM estimate in MA(1) model

estimate.ma1.mom=function(x){

      r=acf(x,plot=F)$acf[1]; if (abs(r)<0.5) return((-1+sqrt(1-4*r^2))/(2*r))

      else return(NA)

      }

# Number of MC simulations

M=2000

ma.mom = rep(0,M)

wn.var.mom = rep(0,M)

for (i in 1:M){

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

      ma.mom[i] <- estimate.ma1.mom(ma.sim)

wn.var.mom[i] <- var(ma.sim)/(1+(estimate.ma1.mom(ma.sim))^2)

      }

# Create historgrams (separately)

hist(ma.mom,xlab=expression(theta),main="")

hist(wn.var.mom,xlab=expression(sigma[e]^2),main="")

# Count how many times (out of M=2000) MOM estimate does not exist

M-length(ma.mom[ma.mom>1])

# Example 7.2

# Lake Huron elevation data

# Page 182

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

plot(huron,ylab="Elevation level (in feet)",xlab="Year",type="o")

# Sample ACF and PACF

# Plots are constructed separately

# Page 183

acf(huron,main="Sample ACF")

pacf(huron,main="Sample PACF")

# Compute sample statistics

# First two sample autocorrelations

acf(huron,lag.max=2,plot=FALSE)

# Sample mean and variance

mean(huron)

var(huron)

# Fit AR(1) and AR(2) models automatically to Huron data

ar(huron,order.max=1,AIC=F,method='yw') # method of moments

ar(huron,order.max=2,AIC=F,method='yw') # method of moments

# Example 7.3

# Gota river discharge data

# Page 190

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

plot(gota,ylab="Water discharge rate",xlab="Year",type="o")

# Sample ACF and PACF

# Plots are constructed separately

# Page 191

acf(gota,main="Sample ACF")

pacf(gota,main="Sample PACF")

# First sample autocorrelation

acf(gota,lag.max=1,plot=FALSE)

# Sample mean and variance

mean(gota)

var(gota)

# Fit MA(1) model using CLS

arima(gota,order=c(0,0,1),method='CSS') # conditional least squares

# Example 7.4

# Lake Huron elevation data

# CLS estimation

# Page 192-193

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

# Fit AR(1) and AR(2) models using CLS

arima(huron,order=c(1,0,0),method='CSS') # conditional least squares

arima(huron,order=c(2,0,0),method='CSS') # conditional least squares

# ARMA best subsets

# Arranged according to BIC

# This figure is not shown in the notes

res=armasubsets(y=huron,nar=6,nma=6,y.name='huron',ar.method='ols')

plot(res)

# Example 7.5

# Bovine blood sugar data

# Page 195

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

plot(cows,ylab="Blood sugar level (mg/100ml blood)",xlab="Days",type="o")

# Sample ACF and PACF

# Plots are constructed separately

# Page 196

acf(cows,main="Sample ACF")

pacf(cows,main="Sample PACF")

# Fit ARMA(1,1) model using CLS

arima(cows,order=c(1,0,1),method='CSS') # conditional least squares

# This figure is not shown in the notes

res=armasubsets(y=cows,nar=6,nma=6,y.name='cows',ar.method='ols')

plot(res)

# Example 7.6

# Gota river discharge data

# ML estimation

# Page 202

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

arima(gota,order=c(0,0,1),method='ML') # maximum likelihood

# Example 7.7.

# Earthquake data

# Page 203-204

# Plots were constructed separately

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

BoxCox.ar(earthquake)

# Square root transformed series

# Sample ACF/PACF and armasubsets output

# Page 205

par(mfrow=c(2,2))

plot(sqrt(earthquake),ylab="No. of earthquakes (Square-root scale)",xlab="Year",type="o")

res=armasubsets(y=sqrt(earthquake),nar=6,nma=6,y.name='sqrt(e.qu.)',ar.method='ols')

plot(res)

acf(sqrt(earthquake),main="Sample ACF")

pacf(sqrt(earthquake),main="Sample PACF")

# Fit ARMA(1,1) model to square-root transformed data

arima(sqrt(earthquake),order=c(1,0,1),method='ML') # maximum likelihood

# Example 7.8.

# Supreme Court data

# Page 206-207

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

par(mfrow=c(2,2))

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

BoxCox.ar(supremecourt)

plot(log(supremecourt),ylab="Percentage granted review (Log scale)",xlab="Year",type="o")

plot(diff(log(supremecourt)),ylab="Difference of logarithms",xlab="Year",type="o")

# These results are not shown in the notes

acf(diff(log(supremecourt)))

pacf(diff(log(supremecourt)))

eacf(diff(log(supremecourt)))

res=armasubsets(y=diff(log(supremecourt)),nar=6,nma=6,y.name='LD(SC)',ar.method='ols')

plot(res)

# Fit ARIMA(0,1,1) model to the log-transformed data

arima(log(supremecourt),order=c(0,1,1),method='ML') # ML