############################
# 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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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 = "http://people.stat.sc.edu/hansont/stat520/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') # ML

# Example 7.8.
# Supreme Court data
# Page 206-207
supremecourt = ts(read.table(file = "http://people.stat.sc.edu/hansont/stat520/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