R exercises from Shumway and Stoffer (2011 - 4th edition)

# The data we are going to use are in the astsa package, which was previously installed
library(astsa)

Exercise 1.1: Compare earthquakes and explosions

ts.plot(EQ5, EXP6, ylim = c(-.5,.5), main="Earthquakes vs. Explosions", 
        gpars = list(col = c("darkgreen", "brown1")))
legend("bottomleft", legend=c("Earthquakes", "Explosions"),
       col=c("darkgreen", "brown1"), lty=1:2, cex=0.5)

As in the example in the book, we denote the two phases or arrivals in the graph by P (t = 1, …, 1024) and S (t= 1025, …, 2048). The distinguishing features of the time series are the rough ampitude ratios of the first phase P and the second phase S. We can see clearly in the graph above that in both phases the amplitude ratios are larger for the earthquake series than for the explosion series, but remarkably at the phase S.

Exercise 1.2 - Signal-plus-noise model

(a) Consider \(x_t = s_t + w_t\), for t = 1,…,200, where

\[ s_{t} = \begin{cases} 0, & \mbox{if t = 1,...,100} \\ 10* \exp{\frac{-(t-100)}{20}} \cos(\frac{2* \pi *t}{4}), & \mbox{if t = 101,...,200} \end{cases} \]

#Generate the data
set.seed(1)
s = c(rep(0,100), 10*exp(-(1:100)/20)*cos(2*pi*101:200/4)) 
x = s + rnorm(200, 0 ,1) 
plot.ts(x)

(b) Consider \(x_t = s_t + w_t\), for t = 1,…,200, where

\[ s_{t} = \begin{cases} 0, & \mbox{if t = 1,...,100} \\ 10* \exp{\frac{-(t-100)}{200}} \cos(\frac{2* \pi *t}{4}), & \mbox{if t = 101,...,200} \end{cases} \]

#Generate the data
set.seed(2)
s_2 = c(rep(0,100), 10*exp(-(1:100)/200)*cos(2*pi*101:200/4)) #generate 100 first obs with value 0 and 100 other observations applying the formula given.
x_2 = s_2 + rnorm(200, 0 ,1) #generate the series summing s and the Gaussian white noise
plot.ts(x_2)

(c) Compare the two series.

We observe that the generated series and the ones from exercise 1.1 are somewhat similiar, specially because there exist very clear phases in both cases. Hence, it seems that both the earthquake and the explosion series data could have been derived from a generating process of the sort of the current exercise. Moreover, the explosion series is more similar to the series in item (a), where the signal modulator is larger (since the denominator is larger), while the earthquake series is similar to that in item b.

# Comparioson for the signal modulators
signal_a = c(exp(-(1:100)/20))
signal_b = c(exp(-(1:100)/200))

ts_a <- ts(signal_a, start=1, end=100, frequency =  1)
ts_b <- ts(signal_b, start=1, end=100, frequency =  1)

ts.plot(ts_a, ts_b, main="Comparison of signal modulators", 
        gpars = list(col = c("cornflowerblue", "darkblue")))
legend("bottomleft", legend=c("Signal a", "Signal b"),
       col=c("cornflowerblue", "darkblue"), lty=1:2, cex=0.7)

We can see that the signal modulator for the series in (a) dissipates faster (approaches zero faster) than the other series. Since this term is multiplicative, this nature of modulators leads the series in (a) to have a smaller amplitude in the second phase, while the amplitude in (b) is more persistent at higher levels because of the slowly decaying modulator.

1.3 Different series with 100 obs

(a) Generate n = 100 obs from the autoregressive model \(x_t = -.9x_{t-2} + w_t\). Apply the moving average filter. Comment on the behavior of xt and how applying the moving average filter changes that behavior.

# Autorregression
set.seed(3)
w = rnorm(150,0,1)
x_a = filter(w, filter=c(0,-.9), method="recursive")[-(1:50)] # as in ex. 1.10

# Moving average filter
v_a = filter(x_a, rep(1/4, 4), sides = 1)

# Plot
par(mfrow=c(1,1))
plot(x_a, xlab = "Time", ylab = "Series", lty = 1, type="l")
lines(v_a, pch = 18, col = "blue", type = "l", 
      lty = 2, lwd = 1)
legend("bottomleft", legend = c("Autoregression", "Moving average"),
       col = c("black", "blue"), lty = 1:2, cex = 0.7)

We can see that the moving average filter smooths out a good portion of the oscillations.

(b) Repeat (a) with \(x_t=\cos(\frac{\pi*t}{4})\)

# Generating the series
set.seed(4)
x_b <- c(cos(2*pi*1:100/4))

# Moving average filter
v_b = filter(x_b, rep(1/4, 4), sides = 1)

# Plot
par(mfrow=c(1,1))
plot(x_b, xlab = "Time", ylab = "Series", lty = 1, type="l")
lines(v_b, pch = 18, col = "blue", type = "l", 
      lty = 2, lwd = 1)
legend("bottomleft", legend = c("Deterministic", "Moving average"),
       col = c("black", "blue"), lty = 1:2, cex = 0.7)

Now, we have a completely deterministic sequence with constant oscillatorry pattern. The moving average filter, then, smooths completely the series, generating a horizontal line in its mean value (0).

(c) Repeat (b), but add the white noise

set.seed(5)
w_1 = w[-(1:50)]
x_c = x_b + w_1
# moving average filter
v_c = filter(x_c, rep(1/4, 4), sides = 1)

# plot
plot(x_c, xlab = "Time", ylab = "Series", lty = 1, type="l")
lines(v_c, pch = 18, col = "blue", type = "l", 
      lty = 2, lwd = 1)
legend("bottomleft", legend = c("Deterministic + WN", "Moving average"),
       col = c("black", "blue"), lty = 1:2, cex = 0.7)

(d) Compare and contrast (a)-(c); i.e., how does the moving average change each series.

The moving average filter is a smoothing operator. In the case of a stationary autorregressive process, the moving average filter approximates the graph to the horizontal line close to the mean. In the case of a stationary completely deterministic process, the moving average filter generates a completely flat line which coincides with the series mean. When a stochastic term is added, the moving average series recalls the graph of a random walk process (not stationary).

1.8 Consider the random walk with drift model

\[ x_t = \delta + x_{t-1} + w_t \] for t = 1,2,…, with \(x_0 = 0\), where \(w_t\) is white noise with with variance \(\sigma_w^2\).

(a) Show that the model can be written as \(x_t = \delta t + \sum_{k=1}^{t} w_k\)

By the formula given, we have \[ (t=1): x_1 = \delta + x_0 + w_1 \] \[ (t=2): x_2 = \delta + x_1 + w_2 \] \[ (t=3): x_3 = \delta + x_2 + w_3 \] \[ . \] \[ . \] \[ . \] \[ (t=T): x_T = \delta + x_{T-1} + w_T \]

Substituting iteralively every \(x_s\) in each line \(\forall s\), we have \[ (t=1): x_1 = \delta + x_0 + w_1 \] \[ (t=2): x_2 = 2* \delta + x_0 + w_1 + w_2 \] \[ (t=3): x_3 = 3* \delta + x_0 + w_1 + w_2 + w_3 \] \[ . \] \[ . \] \[ . \] \[ (t=T): x_T = T* \delta + x_0 + w_1 + w_2 + w_3 +...+ w_T \] Now, considering \(x_0=0\), we can generalize

\[ x_t = t* \delta +\sum_{k=1}^{t} w_k \]

(b) Find the mean function and the autocovariance function of \(x_t\).

The mean function can be expressed as \[ \mathbb{E}[x_t]= \mathbb{E}[ \delta t + \sum_{k=1}^{t} w_k] = \delta t + \mathbb{E}[\sum_{k=1}^{t} w_k] = \delta t + \sum_{k=1}^{t} \mathbb{E}[w_k] = \delta t \] Assuming that \(w_t\) is white noise with mean zero.

The autocovariance function is \[ Cov(x_t, x_{t-l}) = \mathbb{E}[(x_t - \mu_t)(x_{t-l}-\mu_{t-l})] = \mathbb{E}[\sum_{k=1}^{t} w_k \sum_{k=1}^{t-l} w_k] = \min \{t,t-l\} \sigma_w^2 = (t-l) \sigma_w^2 \] Since a white noise process is not serially correlated, i.e., \(Cov(w_t, w_{t-l}) = 0\) for \(l \neq 0\).

(c) Argue that \(x_t\) is not stationary.

By definition, a sequence \(\{r_t\}\) is weakly stationary if its mean function is constant and if the autocovariance function is a function of the time lapse between observations. Here, the mean and the autocovariance functions are time-dependent, which implies that the random walk with drift model cannot be weakly stationary (nor strictly stationary, as consequence).

Indeed, for a simulated random walk with drift model, we can see clearly the stochastic trend, while the estimated ACF shows a very slowly decaying pattern.

set.seed(1234)
n     = 10000
w     = rnorm(n,0,0.0637)
x     = rep(0,n)
x[1]  = 0
for (t in 2:n){
  x[t]  = 0.0103 + x[t-1] + w[t]
}
par(mfrow=c(1,2))
ts.plot(x[1:1000],ylab="",ylim=range(x[1:1000]))
acf(x,main="",lag=100)

(d) Show \(\rho_x(t-1,t) = \sqrt{\frac{t-1}{t}} \rightarrow \infty\) as \(t \rightarrow \infty\). What is the implication of this result?

Take the coefficient of autocorrelation \[ \rho_x(t-1,t) = \frac{Cov(x_{t-1},x_t)}{\sqrt{Var(x_{t-1}) Var(x_t)}} = \frac{(t-1) \sigma_w^2}{\sqrt{(t-1) \sigma_w^2 t \sigma_w^2}} = \frac{\sigma_w^2 t}{\sigma_w^2 \sqrt{(t-1)t}} = \sqrt{\frac{(t-1)^2}{t(t-1)}}= \sqrt{\frac{t-1}{t}} \] Taking the limit \[ \lim_{t\to\infty} \sqrt{\frac{t-1}{t}} = \lim_{t\to\infty} \sqrt{ 1 -\frac{1}{t}} = \sqrt{1} =1 \] This means that two consecutive realizations of a random walk process as t approaches infinity (or t large enough) are perfectly correlated, i.e., \(x_t\) and \(x_{t-1}\) are almost the same thing.

(e) Suggest a transformation to make the series stationary, and prove that the transformed series is stationary.

We can make the series stationary by taking the first difference: \[ \Delta x_t = x_t - x_{t-1} = \delta + w_t \] Which is a function of a constant and a stationary term.

To see that the transformed series is also stationary, we check the mean and the autocovariance functions:

  1. Mean function \[ \mathbb{E}[\Delta x_t] = \mathbb{E}[\delta + w_t] = \delta \]

  2. Autocovariance function \[ Cov(\Delta x_t, \Delta x_{t-l}) = \mathbb{E}[(\Delta x_t - \delta)(\Delta x_{t-l} - \delta)] = \mathbb{E}[w_t w_{t-l}] = 0, \forall l \neq 0 \] \[ Var(\Delta x_t) =\mathbb{E}[w_t^2] = \sigma_w^2 \] Thus, the mean and the autocovariance are constants, which implies that the series taken in the first difference is stationary.

In the simulated series, we have

par(mfrow=c(1,2))
ts.plot(diff(x[1:1000]),ylab="",ylim=range(x[1:1000]))
acf(diff(x),main="",lag=100)

1.19 Suppose \(x_t = \mu + w_t + \theta w_{t-1}\), where \(w_t ~ wn(0,\sigma_w^2)\)

(a) Show that the mean function is

\[ \mathbb{E}(x_t) = \mu \]

We have, for \(w_t\) white noise with mean zero, \[ \mathbb{E}(x_t) = \mathbb{E}(\mu + w_t + \theta w_{t-1}) = \mu + \mathbb{E}(w_t) + \theta \mathbb{E}(w_{t-1}) = \mu \]

(b) Show that the autocovariance function of \(x_t\) is given by \(\gamma_x(0) = \sigma_w^2 ( 1+\theta^2)\), \(\gamma_x(\pm 1) = \sigma_w^2 \theta\), and \(\gamma_x(h) = 0\) otherwise.

The general autocovariance function is \[ Cov(x_t,x_{t-l}) = \mathbb{E}[(w_t + \theta w_{t-1})(w_{t-l}+ \theta w_{t-l-1})] = \mathbb{E}[w_t w_{t-l}] + \theta \mathbb{E}[w_{t-1} w_{t-l}] + \theta \mathbb{E}[w_t w_{t-l-1}] + \theta^2 \mathbb{E}[w_{t-1} w_{t-l-1}] = \gamma_x(l) \] If \(l = 0\), then \[ \gamma_x(0) = \mathbb{E}[w_t^2] + \theta \mathbb{E}[w_{t-1} w_t] + \theta \mathbb{E}[ w_t w_{t-1}] +\theta^2 \mathbb{E}[w_{t-1} w_{t-1}] = \sigma_w^2 + \theta * 0 + \theta * 0 + \theta^2 \sigma_w^2 = \sigma_w^2 (1+ \theta^2) \] If \(l = \pm 1\), then \[ \gamma_x(1) = \mathbb{E}[w_t w_{t-1}] + \theta \mathbb{E}[w_{t-1}^2] + \theta \mathbb{E}[ w_t w_{t-2}] +\theta^2 \mathbb{E}[w_{t-1} w_{t-2}] = 0 + \theta * \sigma_w^2 + \theta * 0 + \theta^2 *0 = \sigma_w^2 \theta \] \[ \gamma_x(-1) = \mathbb{E}[w_t w_{t+1}] + \theta \mathbb{E}[w_{t-1} w_{t+1}] + \theta \mathbb{E}[ w_t^2] +\theta^2 \mathbb{E}[w_{t-1} w_{t}] = 0 + \theta * 0 + \theta * \sigma_w^2 + \theta^2 *0 = \sigma_w^2 \theta \] If \(l \neq \pm 1, 0\), then \[ \gamma_x(l) = \mathbb{E}[w_t w_{t-l}] + \theta \mathbb{E}[w_{t-1} w_{t-l}] + \theta \mathbb{E}[w_t w_{t-l-1}] + \theta^2 \mathbb{E}[w_{t-1} w_{t-l-1}] =0 + \theta * 0 + \theta * 0 + \theta^2 *0 = 0 \]

(c) Show that \(x_t\) is stationary for all values of \(\theta \in \mathbb{R}\).

We note that the first two moments are time invariant (\(\mathbb{E}[x_t] = \mu\) and \(Var(x_t) = \gamma_x(0) = \sigma_w^2 (1+\theta)\)). Hence, for any value of parameter \(\theta\), the series is weakly stationary.

(d) Calculate \(Var(\bar{x})\) for estimating \(\mu\) when (i) \(\theta =1\), (ii) \(\theta = 0\), and (iii) \(\theta = -1\).

We have \[ Var(\bar{x}) = Var (\frac{1}{n} \sum_{t=1}^{n} x_t) = \] \[ \frac{1}{n^2} Cov(\sum_{t=1}^{n}x_t, \sum_{s=1}^{n}x_s) = \frac{1}{n^2} (n \gamma_x(0) + (n-1) \gamma_x(1) + (n-2) \gamma_x(2) +...+ \gamma_x(n-1) + (n-1) \gamma_x(-1)+ (n-2) \gamma_x (-2) + ...+ \gamma_x(1-n)) = \] \[ \frac{1}{n^2}(n \sigma_w^2 (1+ \theta^2) + (n-1) \sigma_w^2 \theta + (n-1) \sigma_w^2 \theta) = \] \[ \frac{1}{n^2}(n \sigma_w^2 (1+ \theta^2) + 2(n-1) \sigma_w^2 \theta) \] (i) If \(\theta = 1\), then \[ Var(\bar{x}) = \frac{2 \sigma_w^2}{n} ( 1 + \frac{n-1}{n}) \] (ii) If \(\theta = 0\), then \[ Var(\bar{x}) = \frac{ \sigma_w^2}{n} \] (iii) If \(\theta = -1\), then \[ Var(\bar{x}) = \frac{2 \sigma_w^2}{n} ( 1 - \frac{n-1}{n}) \]

(e) In time series, the sample size n is typically large, so that \(\frac{n-1}{n} \approx 1\). With this as a consideration, comment on the results of part (d); in particular, how does the accuracy in the estimate of the mean \(\mu\) change for the three different cases?

From basic statistics, we know that \(\bar{x}\) is a consistent estimator of \(\mu\) for a series \(\{x_t\}\), but its accuracy depends on \(var(\bar{x})\). Considering that \(\frac{n-1}{n} \approx 1\), we have

  1. For \(\theta = 1\), \(Var(\bar{x}) = \frac{4 \sigma_w^2}{n} \rightarrow 0\) as \(n \rightarrow \infty\).
  2. For \(\theta = 0\), \(Var(\bar{x}) = \frac{\sigma_w^2}{n} \rightarrow 0\) as \(n \rightarrow \infty\).
  3. For \(\theta = -1\), \(Var(\bar{x}) = 0, \forall n\).

Hence, in the third case, which corresponds to the case of an MA(1) process, the estimate of the mean is always accurate, whilst in the others, the estimate is accurate only asymptotically (if n is very large), with the variance of the second case being smaller than that of the first for each n.

1.20 ACF for Gaussian White Noise series

(a) Simulate a series of n = 500 Gaussian white noise observations and compute the sample ACF, \(\hat{\rho}(h)\), to lag 20. Compare the sample ACF you obtain to the actual ACF, \(\rho(h)\).

set.seed(6)
# Generating 500 obs from Gaussian WN
u_a = rnorm(500,0,1)

# Plotting the sample autocorrelation function
par(mfrow=c(1,2))
ts.plot(u_a,ylab="")
acf(u_a, 20, ylab = 'ACF(u_a)')

We can see that the estimated ACF is statistically 0 in any lag other than the first, implying that the sample is not autocorrelated, and that the first plotted observation is the equivalent to the variance of the Gaussian white noise, which is 1. This is consistent with the theoretical ACF, which builds up on the autocovariance function \[ \gamma_u(h) = cov (u_{t+h}, u_t) = \begin{cases} \sigma_u^2, & \mbox{h=0} \\ 0, & \mbox{h} \neq 0 \end{cases} \] where \(\sigma_u^2 = 1\), in this case.

(b) Repeat part (a) using only n = 50.

set.seed(7)
# Generating 50 obs from Gaussian WN
u_b = rnorm(50,0,1)

# Plotting the sample autocorrelation function
par(mfrow=c(1,2))
ts.plot(u_b,ylab="")
acf(u_b, 20, ylab = 'ACF(u_b)')

The estimated ACF now shows an oscilatory pattern, but the estimated coefficients are still statistically insignificant, although they are closer to the confidence interval limits as compared to the series in item (a). We can infer that the larger the sample, the smaller the confidence intervals for the estimated ACF, and the more it is consistent with the theoretical ACF.

1.21 ACF for Moving Average series

(a) Simulate a series of n = 500 moving average observations as in Example 1.9 and compute the sample ACF, \(\hat{\rho}(h)\), to lag 20. Compare the sample ACF you obtain to the actual ACF, \(\rho (h)\).

# Generating a moving average series as in example 1.9
w_2a = rnorm(500,0,1) 
v_2a = filter(w_2a, sides=2, filter=rep(1/3,3)) # moving average
#Dropping the first and the last observations one because they are NA 
v_2a <- v_2a[-(1:1)]
v_2a <- v_2a[-(499:499)]

# Plotting the estimated ACF
par(mfrow=c(1,2))
ts.plot(v_2a,ylab="")
acf(v_2a, 20, ylab = 'ACF(Moving Average of WN)')

We constructed a three-point moving average process upon a Gaussian white noise series. The autocorrelation graph shows that the series is truncated in lag 2, but other lags (e.g. 3) are also statistically different than zero due to sampling. The theoretical ACF for a moving average process suggests that the graph be truncated at lag 2 as well, while the others lags should indicate 0 autocorrelation.

(b) Repeat part (a) using only n = 50. How does n affect the results?

# simulating the series again for n = 50
w_2b = rnorm(50,0,1) 
v_2b = filter(w_2b, sides=2, filter=rep(1/3,3)) # moving average
#Dropping the first obs and the last one because they are NA 
v_2b <- v_2b[-(1:1)]
v_2b <- v_2b[-(49:49)]

# Plotting the estimated ACF and saving
par(mfrow=c(1,2))
ts.plot(v_2b,ylab="")
acf.2b = acf(v_2b, 20, ylab = 'ACF(Moving Average of WN)')

Now, considering a smaller sample, we can see that the confidence intervals for the estimated ACF are larger than before. As a consequence, the series is truncated at lag 1, and some other lags are statistically different than zero in an oscillatory pattern.