######################################################
# Stat 485
#
# Sarah Rathwell - 301084687
#
# Objective - Investigate ARMA models
######################################################
rm(list=ls())
options(width=200)
library(astsa)
######################################################
cat(" Investigate ARMA models\n\n",
"Last modification:",date(),'\n')
## Investigate ARMA models
##
## Last modification: Mon Nov 24 10:10:21 2014
######################################################
#
# 3.9
# simulate 100 obs for three models: ARMA(1,1), ARMA(0,1),
# ARMA(1,0); theta = .9, phi = .6
#
ARMA11 <- arima.sim(model=list(ar=.6,ma=.9),n=100)
ARMA01 <- arima.sim(model=list(ma=.9),n=100)
ARMA10 <- arima.sim(model=list(ar=.6),n=100)
#
# compute ACF for all three sims
#
ARMA11.acf <- acf(ARMA11,plot=F)
ARMA01.acf <- acf(ARMA01,plot=F)
ARMA10.acf <- acf(ARMA10,plot=F)
#
# compute theoretical ACFs for all models
#
acf.11 <- ARMAacf(ar=.6,ma=.9,lag.max=20)
acf.01 <- ARMAacf(ma=.9,lag.max=20)
acf.10 <- ARMAacf(ar=.6,lag.max=20)
#
# plot against each other
#
plot(ARMA11.acf)
lines(ARMA11.acf$lag,acf.11)
plot(ARMA01.acf)
lines(ARMA01.acf$lag,acf.01)
plot(ARMA10.acf)
lines(ARMA10.acf$lag,acf.10)
#
# compute sample PACF and plot against sample ACF
#
ARMA11.pacf <- pacf(ARMA11,plot=F)
ARMA01.pacf <- pacf(ARMA01,plot=F)
ARMA10.pacf <- pacf(ARMA10,plot=F)
plot(ARMA11.acf$acf[-1],type='l',lty=2,ylim=c(-.8,1),ylab='ACF_PACF',
xlab='lag',main="ACF vs PACF for Series ARMA11")
abline(h=0)
abline(h=c(.2,-.2),lty=3)
lines(ARMA11.acf$lag[-1],ARMA11.pacf$acf)
legend(15,-.25,c("ACF","PACF"),lty=c(2,1))
plot(ARMA01.acf$acf[-1],type='l',lty=2,ylim=c(-.8,1),ylab='ACF_PACF',
xlab='lag',main="ACF vs PACF for Series ARMA01")
abline(h=0)
abline(h=c(.2,-.2),lty=3)
lines(ARMA01.acf$lag[-1],ARMA01.pacf$acf)
legend(15,-.25,c("ACF","PACF"),lty=c(2,1))
plot(ARMA10.acf$acf[-1],type='l',lty=2,ylim=c(-.8,1),ylab='ACF_PACF',
xlab='lag',main="ACF vs PACF for Series ARMA10")
abline(h=0)
abline(h=c(.2,-.2),lty=3)
lines(ARMA10.acf$lag[-1],ARMA10.pacf$acf)
legend(15,-.25,c("ACF","PACF"),lty=c(2,1))
#
# We can see that the sample PACF for the ARMA11 and MA1
# 'tail off' around lag 2, where the sample PACF for AR1 cuts off
# after lag p=1. Similarly, the sample ACF for the ARMA11 and
# appear to tail off slightly while the MA1 cutts off after lag
# 1. Note however, for AR1/MA1, the tailing off happens early (lag
# 1-2), which makes sense since we had set p=1; q=1.
#
######################################################
# 3.10
# (a)
# Use cmort data set to fit AR2 regression
data(cmort)
str(cmort)
## Time-Series [1:508] from 1970 to 1980: 97.8 104.6 94.4 98 95.8 ...
acf2(cmort,20)
## ACF PACF
## [1,] 0.77 0.77
## [2,] 0.77 0.44
## [3,] 0.68 0.03
## [4,] 0.65 0.03
## [5,] 0.58 -0.01
## [6,] 0.53 -0.05
## [7,] 0.48 -0.02
## [8,] 0.41 -0.05
## [9,] 0.39 0.05
## [10,] 0.32 -0.08
## [11,] 0.28 -0.03
## [12,] 0.23 0.00
## [13,] 0.18 -0.06
## [14,] 0.13 -0.06
## [15,] 0.11 0.05
## [16,] 0.04 -0.10
## [17,] 0.01 -0.04
## [18,] -0.03 -0.01
## [19,] -0.07 -0.03
## [20,] -0.08 0.03
cmort.fit <- ar.ols(cmort,order=2,demean=F,intercept=T)
cmort.fit
##
## Call:
## ar.ols(x = cmort, order.max = 2, demean = F, intercept = T)
##
## Coefficients:
## 1 2
## 0.4286 0.4418
##
## Intercept: 11.45 (2.394)
##
## Order selected 2 sigma^2 estimated as 32.32
cmort.fit$asy.se.coef # SEs
## $x.mean
## [1] 2.393673
##
## $ar
## [1] 0.03979433 0.03976163
#
# Our estimates are Phi0=11.45 se(2.394), Phi1=0.43 (se=.04),
# Phi2=0.44 (se=.04), and sigma^2=32.32.
#
# Est. model --> Xt = 11.45 + 0.43*X(t-1) + 0.44*X(t-2)
#
# (b)
# Find 4 week forcasts and 95% PI's for model (a)
#
cmort.pred <-predict(cmort.fit,n.ahead=4)
#
# list the 95% PI for each increase
#
upper <-c(cmort.pred$pred+2*cmort.pred$se)
lower <-c(cmort.pred$pred-2*cmort.pred$se)
list(m1=c(lower[1],upper[1]),m2=c(lower[2],upper[2]),
m3=c(lower[3],upper[3]),m4=c(lower[4],upper[4]))
## $m1
## [1] 76.23017 98.96956
##
## $m2
## [1] 74.39354 99.13344
##
## $m3
## [1] 73.06868 101.60559
##
## $m4
## [1] 72.02679 102.40021
######################################################
# 3.18
# (a)
# compare LM and YW coef estimates
cmort.yw <-ar.yw(cmort,order=2)
list(LMests=c(cmort.fit$ar),YWests=c(cmort.yw$ar))
## $LMests
## [1] 0.4285906 0.4417874
##
## $YWests
## [1] 0.4339481 0.4375768
#
# Our Yule-Walker coeffs estimates are Phi1=0.43, Phi2=0.44,
# and sigma^2=32.8
# The estimates are nearly identical for both lm and yw.
#
# (b)
# compare var of estimators
list(LMse=cmort.fit$asy.se.coef$ar,YWse=sqrt(diag(cmort.yw$asy.var.coef)))
## $LMse
## [1] 0.03979433 0.03976163
##
## $YWse
## [1] 0.04001303 0.04001303
#
# Again, the SE's of our estimated coeff's are almost identical
######################################################