Section 5.6 Homework for Time Series Analysis

Section 5.6 homework for Time Series Analysis

Below are partial answers to the homework problems.

5.11 Using the set up given for the problem in Shumway and Stoffer, complete the following:

a)Examine plots of St and Lt over time.

> load(file = "C:\\chris\\unl\\Dropbox\\NEW\\STAT_time_series\\tsa3.rda")

> S<-sales

> L<-lead

> head(S)

[1] 200.1 199.5 199.4 198.9 199.0 200.2

> head(L)

[1] 10.01 10.07 10.32 9.75 10.33 10.13

> ######################################################################

> # Initial plots

> win.graph(width = 8, height = 6, pointsize = 10)

> par(mfrow = c(2,1))

> plot(x = S, ylab = expression(S[t]), xlab = "t", type = "l", col = "red", main

= "Plot of sales.dat")

> points(x = S, pch = 20, col = "blue")

> abline(v = c(0, 50, 100, 150, 200), lty = "dotted", col = "gray")

> abline(h = c(200, 220, 240, 260), lty = "dotted", col = "gray")

> plot(x = L, ylab = expression(L[t]), xlab = "t", type = "l", col = "red", main

= "Plot of lead.dat")

> points(x = L, pch = 20, col = "blue")

> abline(v = c(0, 50, 100, 150, 200), lty = "dotted", col = "gray")

> abline(h = 10:14, lty = "dotted", col = "gray")

> par(mfrow = c(1,1))

Both series are non-stationary in the mean.

b)Examine plots of (1-B)St and (1-B)Lt over time.

> par(mfrow = c(2,1))

> plot(x = diff(x = S, lag = 1, differences = 1), ylab = expression(S[t] –

S[t-1]), xlab = "t", type = "l", col = "red", main = "Plot of first

differenced sales.dat")

> points(x = diff(x = S, lag = 1, differences = 1), pch = 20, col = "blue")

> abline(v = c(0, 50, 100, 150, 200), lty = "dotted", col = "gray")

> abline(h = c(-2, 0, 2, 4), lty = "dotted", col = "gray")

> plot(x = diff(x = L, lag = 1, differences = 1), ylab = expression(L[t] –

L[t-1]), xlab = "t", type = "l", col = "red",

main = "Plot of first differenced lead.dat")

> points(x = diff(x = L, lag = 1, differences = 1), pch = 20, col = "blue")

> abline(v = c(0, 50, 100, 150, 200), lty = "dotted", col = "gray")

> abline(h = c(-0.5, 0, 0.5), lty = "dotted", col = "gray")

> par(mfrow = c(1,1))

Both series now appear to be stationary in the mean.

c)Fit a simple linear regression model assuming independent error terms as shown in #5.11. State the model. Be careful with obtaining the “lead” part of the model. Below is how I got the data set into the correct form before using lm() to fit the model.

> set1<-ts.intersect(delta.S = diff(x = S, lag = 1, differences = 1),

delta.B3.L = lag(diff(x = L, lag = 1, differences = 1), k = -3))

> set1<-ts.intersect(delta.S = diff(x = S, lag = 1, differences = 1),

delta.B3.L = lag(diff(x = L, lag = 1, differences = 1), k = -3))

> mod.fit<-lm(formula = delta.S ~ delta.B3.L, data = set1)

> summary(mod.fit)

Call:

lm(formula = delta.S ~ delta.B3.L, data = set1)

Residuals:

Min 1Q Median 3Q Max

-2.18876 -0.65502 -0.07291 0.60347 2.84546

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 0.35538 0.08301 4.281 3.38e-05 ***

delta.B3.L 3.33733 0.26190 12.743 < 2e-16 ***

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1 on 144 degrees of freedom

Multiple R-Squared: 0.53, Adjusted R-squared: 0.5267

F-statistic: 162.4 on 1 and 144 DF, p-value: < 2.2e-16

The estimated model is (1 – B) = 0.3554 + 3.3373(1-B)Lt-3 where the error term is assumed to be independent N(0,2) with 2 estimated to be 1.

d)Examine the residuals for part c’s model fit. What ARMA model would be appropriate for the residuals? Explain.

> #Examine residuals

> par(mfrow = c(1,2))

> acf(x = mod.fit$residuals, type = "correlation", lag.max = 20, ylim = c(-1,1),

main = "Estimated ACF plot for residuals")

> pacf(x = mod.fit$residuals, lag.max = 20, ylim = c(-1,1), xlab = "h", main =

"Estimated ACF plot for residuals")

> par(mfrow = c(1,1))

An AR(1) model would be appropriate for the residuals because the estimated PACF cuts off to 0 after lag 1 and the estimated ACF tails off to 0.

e)Fit the same regression as in part c, but using the ARMA model from part d for the residuals. Show how to use both the gls() and arima() functions. State the model.

> #Need time indicator

> set2<-data.frame(set1, time.ind = 1:nrow(set1))

> mod.fit<-gls(model = delta.S ~ delta.B3.L, data = set2, correlation =

corARMA(form = ~ time.ind, p = 1))

> summary(mod.fit)

Generalized least squares fit by REML

Model: delta.S ~ delta.B3.L

Data: set2

AIC BIC logLik

349.0892 360.9685 -170.5446

Correlation Structure: AR(1)

Formula: ~time.ind

Parameter estimate(s):

Phi

0.6541633

Coefficients:

Value Std.Error t-value p-value

(Intercept) 0.3622145 0.1825202 1.984517 0.0491

delta.B3.L 2.7853560 0.1426505 19.525731 0.0000

Correlation:

(Intr)

delta.B3.L -0.02

Standardized residuals:

Min Q1 Med Q3 Max

-2.1446477 -0.6881815 -0.1018380 0.6098130 2.6662434

Residual standard error: 1.021178

Degrees of freedom: 146 total; 144 residual

The estimated model is (1 – B) = 0.3622 + 2.7854(1-B)Lt-3 where the error term has an AR(1) model of (1 – 0.6542B)xt = wt where wt ~ independent N(0,2) with 2 estimated to be 1.02.

> mod.fit.SS<-arima(x = set2$delta.S, order = c(1, 0, 0), xreg = set2$delta.B3.L)

> mod.fit.SS

Call:

arima(x = set2$delta.S, order = c(1, 0, 0), xreg = set2$delta.B3.L)

Coefficients:

ar1 intercept set2$delta.B3.L

0.6451 0.3624 2.7876

s.e. 0.0628 0.1767 0.1432

sigma^2 estimated as 0.5884: log likelihood = -168.72, aic = 345.43

The estimated model is (1 – B) = 0.3624 + 2.7876(1-B)Lt-3 where the error term has an AR(1) model of (1 – 0.6451B)xt = wt where wt ~ independent N(0,2) with 2 estimated to be 1.02.

f)Examine the appropriate residuals from the fit of the model in part e that account for the ARMA error term. Does there need to be any adjustment to the ARMA model chosen for the residuals? Explain.

> w.resid2<-residuals(object = mod.fit, type = "normalized")

> head(w.resid2)

[1] -0.4204317 0.5467604 -0.6047486 -0.4307466 1.1300007 -0.5522774

> par(mfrow = c(1,2))

> acf(x = w.resid2, lag.max = 20, type = "correlation", main = "Estimated ACF for

normalized residuals", xlim = c(1,20), ylim = c(-1,1))

> pacf(x = w.resid2, lag.max = 20, main = "Estimated PACF for normalized

residuals", xlim = c(1,20), ylim = c(-1,1))

> par(mfrow = c(1,1))

The estimated ACF and PACF plot show almost no significant values. The PACF has one marginally significant value at lag 9. Given these results, the AR(2) model for the error term has accounted for well the dependence in the error structure.

From using arima(), we obtain a similar conclusion:

> examine.mod(mod.fit.obj = mod.fit.SS, mod.name = "Reg. model with AR(1) error",

max.lag = 20)

$z

ar1 intercept set2$delta.B3.L

10.266335 2.050868 19.461046

$p.value

ar1 intercept set2$delta.B3.L

0.00000000 0.04027983 0.00000000

1