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