##4.3.

Consider the following two models: \[(i) \;\;\; x_t = .80x_{t-1} - .15x_{t-2} + w_t - .30w_{t-1}\] \[(ii) \;\;\; x_t = x_{t-1} - .50x_{t-2} + w_t - w_{t-1}\]

  1. Using Example 4.10 as a guide, check the models for parameter redundancy. If a model has redundancy, find the reduced form of the model.

We rewrite the model in the following form:

\[ i)\;\;\; x_t(1-.8B+.15B^2) = w_t(1-.30B) \] \[ ii)\;\; \; x_t(1-B-.50B^2) = w_t(1-B) \]

and use polyroot to find the roots of both sides.

#for model i
ARi = c(1, -.8, .15) # original AR coefs on the left
polyroot(ARi)
## [1] 2.000000-0i 3.333333+0i
MAi = c(1, -.3) # original MA coefs on the right
polyroot(MAi)
## [1] 3.333333+0i
#for model ii
ARii = c(1, -1, .5) # original AR coefs on the left
polyroot(ARii)
## [1] 1+1i 1-1i
MAii = c(1, -1) # original MA coefs on the right
polyroot(MAii)
## [1] 1+0i

Model i is overparametrized, while model ii is already in its simplest form.

The work for the simplified model of i is the following:

for AR: \[x_t(1−0.80B+0.15B^2) = x_t(B−2)(B−3.333333)\]

for MA:

\[w_t(1−0.30B)=(B−3.333333)w_t\]

returning both AR and MA terms equal to each other:

\[x_t(B−2)(B−3.333333) = (B−3.333333)w_t\]

we then divide by the common term, \((B-3.333333)\) and we have:

\[x_t(B−2) = w_t \]

Expand and rearrange:

\[ x_tB - 2x_t = w_t\] \[ x_{t-1} - 2x_t = w_t\] solve for \(x_t\) :

\[x_t = .5x_{t-1} +.5w_t\]

  1. A way to tell if an ARMA model is causal is to examine the roots of AR term \(\phi(B)\) to see if there are no roots less than or equal to one in magnitude. Likewise, to determine invertibility of a model, the roots of the MA term \(\theta(B)\) must not be less than or equal to one in magnitude. Use Example 4.11 as a guide to determine if the reduced (if appropriate) models (i) and (ii), are causal and/or invertible.

model i, reduced model AR(1):

\(|\phi|\) = .5

\(|\theta|\) = .5

The model is causal because \(|\phi| < 1\)

The model is invertible because $theta < 1 $

model ii, given:

\[(ii) \;\;\; x_t = x_{t-1} - .50x_{t-2} + w_t - w_{t-1}\] rewrite using backshift operator:

AR part:

\[ 1-B+.5B^2 = (B- (1+i))(B-(1-i)) \] \[|a+bi| = \sqrt{a^2+b^2}\\ |1+i| = \sqrt{1^2+1^2} = \sqrt{2} \\ |1-i| = \sqrt{1^2+1^2} = \sqrt{2} \\ \sqrt{2} > 1 \] From the AR part the model is causal because the roots are > 1

MA part:

\[ 1 - B = 0 \] Where the root for the MA process is exactly one and is therefore not invertible.

  1. In Example 4.3 and Example 4.12, we used ARMAtoMA and ARMAtoAR to exhibit some of the coefficients of the causal [\(MA(\infty)\)] and invertible [\(AR(\infty)\)] representations of a model. If the model is in fact causal or invertible, the coefficients must converge to zero fast. For each of the reduced (if appropriate) models (i) and (ii), find the first 50 coefficients and comment.

reduced model i

library(astsa)
## Warning: package 'astsa' was built under R version 4.1.3
round(ARMAtoMA(ar = .5, ma = .5, lag.max = 50),2)
##  [1] 1.00 0.50 0.25 0.12 0.06 0.03 0.02 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [16] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [31] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [46] 0.00 0.00 0.00 0.00 0.00
round(ARMAtoAR(ar = .5, ma = .5, lag.max = 50),2)
##  [1] -1.00  0.50 -0.25  0.12 -0.06  0.03 -0.02  0.01  0.00  0.00  0.00  0.00
## [13]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [25]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [37]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [49]  0.00  0.00

for model i, we can say that the earlier work was correct, the reduced model is causal and invertible, as both coefficients quickly converge to 0

model ii

round( ARMAtoMA(ar=c(1,-.5), ma=c(-1), 50), 2)
##  [1]  0.00 -0.50 -0.50 -0.25  0.00  0.12  0.12  0.06  0.00 -0.03 -0.03 -0.02
## [13]  0.00  0.01  0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [25]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [37]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [49]  0.00  0.00
round( ARMAtoAR(ar=c(1,-.5), ma=c(-1), 50), 2)
##  [1] 0.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
## [20] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
## [39] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

For model ii, we can see from the AR representation that it is not invertible because it does not converge to 1, but it is causal because the MA representation does converge.

#4.4. (a) Compare the theoretical ACF and PACF of an ARMA(1, 1), an ARMA(1, 0), and an ARMA(0, 1) series by plotting the ACFs and PACFs of the three series for \[\phi = .6 \] and \[ \theta = .9 \]. Comment on the capability of the ACF and PACF to determine the order of the models. Hint: See the code for Example 4.18.

ARMA(1,1)

phi <- 0.6
theta <- 0.9
acf11 = ARMAacf(ar=phi, ma=theta, 24)[-1]
pacf11 = ARMAacf(ar=phi, ma=theta, 24, pacf=TRUE)
par(mfrow=1:2)
tsplot(acf11, type="h", xlab="lag", ylim=c(-.8,1))
abline(h=0)
tsplot(pacf11, type="h", xlab="lag", ylim=c(-.8,1))
abline(h=0)

ARMA(1,0)

acf10 = ARMAacf(ar=phi, ma=0, 24)[-1]
pacf10 = ARMAacf(ar=phi, ma=0, 24, pacf=TRUE)
par(mfrow=1:2)
tsplot(acf10, type="h", xlab="lag", ylim=c(-1,1))
abline(h=0)
tsplot(pacf10, type="h", xlab="lag", ylim=c(-1,1))
abline(h=0)

ARMA(0,1)

acf01 = ARMAacf(ar=0, ma=theta, 24)[-1]
pacf01 = ARMAacf(ar=0, ma=theta, 24, pacf=TRUE)
par(mfrow=1:2)
tsplot(acf01, type="h", xlab="lag", ylim=c(-1,1))
abline(h=0)
tsplot(pacf01, type="h", xlab="lag", ylim=c(-1,1))
abline(h=0)

ARMA(1,1) tails off in both ACF and PACF ARMA(1,0) ACF tails off and PACF cuts off past lag 1 ARMA(0,1) ACF cuts off and PACF tails off

And table 4.1 says that they should behave that way. ARMA(1,0) is really an AR(1), whose ACF tails off and PACF cuts off and ARMA(0,1) is really an MA(1) whose ACF and PACF behavior cuts off and tails off respectively. Given this, all models work pretty well, which leads me to say that perhaps another measure is needed to better select the model, such as AIC, BIC.

  1. Use arima.sim to generate n = 100 observations from each of the three models discussed in (a). Compute the sample ACFs and PACFs for each model and compare it to the theoretical values. How do the results compare with the general results given in Table 4.1?
n <- 100
arma11_100 <- arima.sim(model = list(ar = phi, ma = theta), n = n)
arma10_100 <- arima.sim(model = list(ar = phi, ma = 0), n = n)    
arma01_100 <- arima.sim(model = list(ar = 0, ma = theta), n= n)
par(mfrow=1:2)
acf(arma11_100)
pacf(arma11_100)

par(mfrow=1:2)
acf(arma10_100)
pacf(arma10_100)

par(mfrow=1:2)
acf(arma01_100)
pacf(arma01_100)

in comparison to the theoretical evaluations in a), The behavior of the acf and pacf of the simulated models is completely different for the AR and MA models, while ARMA remained constant. Solely from these observations, it looks like ARMA is the better model for this problem.

  1. Repeat (b) but with n = 500. Comment.
n <- 500
arma11_500 <- arima.sim(model = list(ar = phi, ma = theta), n = n)
arma10_500 <- arima.sim(model = list(ar = phi, ma = 0), n = n)    
arma01_500 <- arima.sim(model = list(ar = 0, ma = theta), n= n)
par(mfrow = c(1,2))
acf(arma11_500)
pacf(arma11_500)

par(mfrow = c(1,2))
acf(arma10_500)
pacf(arma10_500)

par(mfrow = c(1,2))
acf(arma01_500)
pacf(arma01_500)

My answer in this question doesn’t change from the answer for b. The behaviors of MA and AR for b) and c) remained consistent, while ARMA was consistent across all trials. I would still say that ARMA(1,1) is the model to choose.

#4.5.

Let \(c_t\) be the cardiovascular mortality series (cmort) discussed in Example 3.5 and let \(x_t = \nabla c_t\) be the differenced data.

  1. Plot \(x_t\) and compare it to the actual data plotted in Figure 3.2. Why does differencing seem reasonable in this case?
data(cmort)
c_t <- cmort
x_t <- diff(cmort, lag = 1)
tsplot(c_t, col = "blue")

tsplot(x_t, col = "red")

The differentiation makes makes \(x_t\) look more stationary and it helps remove the cyclic effects imprinted in \(c_t\)

  1. Calculate and plot the sample ACF and PACF of \(x_t\) and using Table 4.1, argue that an AR(1) is appropriate for \(x_t\)
par(mfrow = c(1,2))
acf(x_t, lag.max = 5)
pacf(x_t, lag.max = 5)

AR(1) is an appropiate model as the ACF tails off rather quickly and the PACF looks like it cuts off past lag 1, or in this case, lag p.

  1. Fit an AR(1) to xt using maximum likelihood (basically unconditional least squares) as in Section 4.3. The easiest way to do this is to use sarima from astsa. Comment on the significance of the regression parameter estimates of the model. Wha◘?
sarima(xdata = x_t, p = 1, d = 0, q = 0)
## initial  value 1.908706 
## iter   2 value 1.760342
## iter   3 value 1.760341
## iter   4 value 1.760341
## iter   4 value 1.760341
## iter   4 value 1.760341
## final  value 1.760341 
## converged
## initial  value 1.760656 
## iter   2 value 1.760656
## iter   3 value 1.760656
## iter   3 value 1.760656
## iter   3 value 1.760656
## final  value 1.760656 
## converged

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = xmean, include.mean = FALSE, transform.pars = trans, fixed = fixed, 
##     optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1    xmean
##       -0.5064  -0.0263
## s.e.   0.0383   0.1715
## 
## sigma^2 estimated as 33.81:  log likelihood = -1612.05,  aic = 3230.11
## 
## $degrees_of_freedom
## [1] 505
## 
## $ttable
##       Estimate     SE  t.value p.value
## ar1    -0.5064 0.0383 -13.2233  0.0000
## xmean  -0.0263 0.1715  -0.1533  0.8782
## 
## $AIC
## [1] 6.371023
## 
## $AICc
## [1] 6.37107
## 
## $BIC
## [1] 6.396044

\(\sigma^2_w = 33.81\)

  1. Examine the residuals and comment on whether or not you think the residuals are white.

The residuals look like white noise. We can see this in the standardized residual and QQ plot. In the Standardized residual, it looks like it has a mean of 0, and the QQ plot tells us they’re normal, so the residuals are white noise.