STAT 520 R Code: Chapter 6

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

# Author: Joshua M. Tebbs #

# Date: 20 Dec 2009 #

# Update: 25 Jul 2011 #

# Purpose: STAT 520 R code #

# CHAPTER 6 #

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

# Example 6.2.

# Monte Carlo simulation

# MA(1) with theta = -0.7

# Simulate sampling distribution of r_1, r_2, r_5, and r_10

# n = 200 (sample size)

# M = 2000 (number of MC simulations)

# Page 141

M<-2000

r.1<-rep(0,M)

r.2<-rep(0,M)

r.5<-rep(0,M)

r.10<-rep(0,M)

for(i in 1:M){ts <- arima.sim(list(order = c(0,0,1), ma = 0.7), n =

200)

temp<-as.numeric(acf(ts,plot=FALSE)[[1]])

r.1[i]<-temp[2]

r.2[i]<-temp[3]

r.5[i]<-temp[6]

r.10[i]<-temp[11]}

par(mfrow=c(2,2))

hist(r.1,xlab=expression(r[1]),main="")

hist(r.2,xlab=expression(r[2]),main="")

hist(r.5,xlab=expression(r[5]),main="")

hist(r.10,xlab=expression(r[10]),main="")

# Example 6.3.

# MA(1) and MA(2) simulations and ACFs

# MA margin of error bounds are used

# Page 142

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

ma2.sim <- arima.sim(list(order = c(0,0,2), ma = c(-0.5,0.5)), n = 100)

par(mfrow=c(2,2))

plot(ma1.sim,ylab=expression(Y[t]),main="MA(1) process",xlab="Time",type="o")

acf(ma1.sim,ci.type='ma',main="Sample ACF")

plot(ma2.sim,ylab=expression(Y[t]),main="MA(2) process",xlab="Time",type="o")

acf(ma2.sim,ci.type='ma',main="Sample ACF")

# Example 6.4.

# AR(1) and AR(2) population ACF/PACF

# Uses ARMAacf function

# ARMAacf function includes the k=0 lag for ACF

# Use y = y[2:21] to remove k=0 lag from ARMAacf output; only for ACF

# Not needed for PACF

# Page 149

par(mfrow=c(2,2))

y = ARMAacf(ar = c(0.9), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(0.9), lag.max = 20, pacf=TRUE)

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Partial autocorrelation", main = "Population PACF")

abline(h = 0)

y = ARMAacf(ar = c(-0.5,0.25), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ar = c(-0.5,0.25), lag.max = 20, pacf=TRUE)

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Partial autocorrelation", main = "Population PACF")

abline(h = 0)

# Example 6.4.

# AR(1) and AR(2) simulated ACFs and PACFs

# Page 150

ar1.sim <- arima.sim(list(order = c(1,0,0), ar = c(0.9)), n = 150)

ar2.sim <- arima.sim(list(order = c(2,0,0), ar = c(-0.5,0.25)), n = 150)

par(mfrow=c(3,2))

plot(ar1.sim,ylab=expression(Y[t]),main="AR(1) process",xlab="Time",type="o")

plot(ar2.sim,ylab=expression(Y[t]),main="AR(2) process",xlab="Time",type="o")

acf(ar1.sim,main="AR(1) sample ACF")

acf(ar2.sim,main="AR(2) sample ACF")

pacf(ar1.sim,main="AR(1) sample PACF")

pacf(ar2.sim,main="AR(2) sample PACF")

# Example 6.5.

# MA(1) and MA(2) population ACF/PACF

# Uses ARMAacf function

# ARMAacf function includes the k=0 lag for ACF

# Use y = y[2:21] to remove k=0 lag from ARMAacf output; only for ACF

# Not needed for PACF

# Page 151

par(mfrow=c(2,2))

y = ARMAacf(ma = c(-0.9), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ma = c(-0.9), lag.max = 20, pacf=TRUE)

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Partial autocorrelation", main = "Population PACF")

abline(h = 0)

y = ARMAacf(ma = c(0.5,-0.25), lag.max = 20)

y = y[2:21]

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Autocorrelation", main = "Population ACF")

abline(h = 0)

y = ARMAacf(ma = c(0.5,-0.25), lag.max = 20, pacf=TRUE)

plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",

ylab = "Partial autocorrelation", main = "Population PACF")

abline(h = 0)

# Example 6.5.

# MA(1) and MA(2) simulated ACFs and PACFs

# Page 152

ma1.sim <- arima.sim(list(order = c(0,0,1), ma = c(-0.9)), n = 150)

ma2.sim <- arima.sim(list(order = c(0,0,2), ma = c(0.5,-0.25)), n = 150)

par(mfrow=c(3,2))

plot(ma1.sim,ylab=expression(Y[t]),main="MA(1) process",xlab="Time",type="o")

plot(ma2.sim,ylab=expression(Y[t]),main="MA(2) process",xlab="Time",type="o")

acf(ma1.sim,main="MA(1) sample ACF")

acf(ma2.sim,main="MA(2) sample ACF")

pacf(ma1.sim,main="MA(1) sample PACF")

pacf(ma2.sim,main="MA(2) sample PACF")

# Example 6.6.

# ARMAacf function to compute theoretical ACF and PACF

# AR(2) and MA(2) calculations

# Page 154

round(ARMAacf(ar = c(0.6,-0.4), lag.max = 10),digits=3)

round(ARMAacf(ar = c(0.6,-0.4), lag.max = 10, pacf=TRUE),digits=3)

round(ARMAacf(ma = c(0.6,-0.4), lag.max = 10),digits=3)

round(ARMAacf(ma = c(0.6,-0.4), lag.max = 10, pacf=TRUE),digits=3)

# Example 6.7.

# Sample EACF calculation

# Page 162-163

#ARMA(1,1)

arma11.sim <- arima.sim(list(order = c(1,0,1), ar = c(0.6), ma = c(0.8)), n = 200)

eacf(arma11.sim)

#ARMA(2,2)

arma22.sim <- arima.sim(list(order = c(2,0,2), ar = c(0.5,-0.5), ma = c(0.8,-0.2)), n = 200)

eacf(arma22.sim)

#ARMA(3,3)

arma33.sim <- arima.sim(list(order = c(3,0,3), ar = c(0.8,0.8,-0.9), ma = c(-0.9,0.8,-0.2)), n = 200)

eacf(arma33.sim)

# Example 6.8.

# Global temperature and LA rainfall data

# Time series plots

# Plots were constructed separately

# Page 168

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

plot(globaltemps,ylab="Global temperature deviations",xlab="Year",type="o")

data(larain)

plot(larain,ylab="LA rainfall amounts",xlab="Year",type="o")

# Example 6.8.

# Global temperature and LA rainfall data

# Augmented Dickey-Fuller unit root test output

# The user must install uroot package first!!

# Page 168-169

# Global temperature data

library(uroot)

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

# Run this command to get best value of k (here, k=3)

ar(diff(globaltemps))

ADF.test(globaltemps,selectlags=list(mode=c(1,2,3),Pmax=3),itsd=c(1,0,0))

# LA rainfall data

library(uroot)

data(larain)

# Run this command to get best value of k (here, k=4)

ar(diff(larain))

ADF.test(larain,selectlags=list(mode=c(1,2,3,4),Pmax=4),itsd=c(1,0,0))

# Example 6.9.

# Ventilation data

# Plots were created separately

# Page 171

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

plot(ventilation,ylab="Ventilation (L/min)",xlab="Observation time",type="o")

plot(diff(ventilation),ylab="Ventilation differences",xlab="Observation time",type="o")

# Perform Dickey-Fuller unit root test

library(uroot)

ar(diff(ventilation))

# Last command suggests k=4

ADF.test(ventilation,selectlags=list(mode=c(1,2,3,4),Pmax=4),itsd=c(1,0,0))

# Last command gives p-value > 0.10 (original series is not stationary)

# Therefore, use BIC on data differences

# ARMA best subsets

# Arranged according to BIC

# Using first differences

# Page 172

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

plot(res)