Problem 1

Uniform white noise

n <- 1000
W <- runif(n, -0.5, 0.5)
plot(W, type = "l", main = "Uniform white noise")

Autocovariance and autocorrelation for h = 0, 1, …, 20:

acf(W, lag.max = 20, type = "covariance", main = "gamma_W, h <= 20")

acf(W, lag.max = 20, main = "rho_W, h <= 20")

Same for h = 0, 1, …, 100:

acf(W, lag.max = 100, type = "covariance", main = "gamma_W, h <= 100")

acf(W, lag.max = 100, main = "rho_W, h <= 100")

The true values are gamma_W(0) = 1/12 and gamma_W(h) = 0 for h not equal to 0. The estimates fluctuate around zero with no visible pattern.

A few sample autocorrelations cross the confidence bands. This is what should happen rather than a sign of a problem: the bands are 95 percent bands, so about 1 lag in 20 crosses by chance, and over 100 lags around five crossings are expected under white noise.

Box-Pierce and Ljung-Box on the white noise

The rule of thumb gives k no larger than n/5 = 200.

for (k in c(5, 10, 20, 50)) {
  print(Box.test(W, lag = k, type = "Box-Pierce"))
  print(Box.test(W, lag = k, type = "Ljung-Box"))
}
## 
##  Box-Pierce test
## 
## data:  W
## X-squared = 1.2446, df = 5, p-value = 0.9405
## 
## 
##  Box-Ljung test
## 
## data:  W
## X-squared = 1.2486, df = 5, p-value = 0.9401
## 
## 
##  Box-Pierce test
## 
## data:  W
## X-squared = 9.258, df = 10, p-value = 0.5078
## 
## 
##  Box-Ljung test
## 
## data:  W
## X-squared = 9.345, df = 10, p-value = 0.4997
## 
## 
##  Box-Pierce test
## 
## data:  W
## X-squared = 22.921, df = 20, p-value = 0.2927
## 
## 
##  Box-Ljung test
## 
## data:  W
## X-squared = 23.25, df = 20, p-value = 0.2767
## 
## 
##  Box-Pierce test
## 
## data:  W
## X-squared = 49.837, df = 50, p-value = 0.4799
## 
## 
##  Box-Ljung test
## 
## data:  W
## X-squared = 51.22, df = 50, p-value = 0.4256

All p-values lie between 0.28 and 0.94, so there is no evidence against white noise at any of these lags. The two tests give almost identical statistics, which is expected: they differ only by the factor (n + 2) / (n - h), and with n = 1000 that factor is close to 1.

MA(2) process

X_t = W_t + W_(t-1) + W_(t-2)

X <- W[3:n] + W[2:(n - 1)] + W[1:(n - 2)]
plot(X, type = "l", main = "MA(2) process")

acf(X, lag.max = 20, type = "covariance", main = "gamma_X, h <= 20")

acf(X, lag.max = 20, main = "rho_X, h <= 20")

acf(X, lag.max = 100, main = "rho_X, h <= 100")

Comparison with the true functions

With sigma^2 = 1/12 the true autocovariances are gamma_X(0) = 3 sigma^2, gamma_X(1) = 2 sigma^2, gamma_X(2) = sigma^2 and 0 for h > 2, so the true autocorrelations are 1, 2/3, 1/3, 0, 0, …

s2 <- 1 / 12
h <- 0:5
gamma_true <- c(3, 2, 1, 0, 0, 0) * s2
gamma_hat <- drop(acf(X, lag.max = 5, type = "covariance", plot = FALSE)$acf)
rho_hat <- drop(acf(X, lag.max = 5, plot = FALSE)$acf)
data.frame(h, gamma_hat, gamma_true, rho_hat, rho_true = gamma_true / gamma_true[1])
##   h     gamma_hat gamma_true      rho_hat  rho_true
## 1 0  0.2606872730 0.25000000  1.000000000 1.0000000
## 2 1  0.1779762418 0.16666667  0.682719335 0.6666667
## 3 2  0.0908648368 0.08333333  0.348558776 0.3333333
## 4 3  0.0061118300 0.00000000  0.023445065 0.0000000
## 5 4  0.0005744443 0.00000000  0.002203576 0.0000000
## 6 5 -0.0036813992 0.00000000 -0.014121898 0.0000000

The estimates are close to the true values. The sample autocorrelations are 0.683 at h = 1 against a true 2/3, and 0.349 at h = 2 against a true 1/3, and for h of 3 or more they are within about 0.02 of zero. The estimated variance gamma_hat(0) = 0.261 sits slightly above the true 3 sigma^2 = 0.25, which is ordinary sampling error at this sample size.

Box-Pierce and Ljung-Box on the MA(2)

for (k in c(5, 10, 20, 50)) {
  print(Box.test(X, lag = k, type = "Box-Pierce"))
  print(Box.test(X, lag = k, type = "Ljung-Box"))
}
## 
##  Box-Pierce test
## 
## data:  X
## X-squared = 587.18, df = 5, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  X
## X-squared = 589.07, df = 5, p-value < 2.2e-16
## 
## 
##  Box-Pierce test
## 
## data:  X
## X-squared = 595.18, df = 10, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  X
## X-squared = 597.14, df = 10, p-value < 2.2e-16
## 
## 
##  Box-Pierce test
## 
## data:  X
## X-squared = 601.56, df = 20, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  X
## X-squared = 603.65, df = 20, p-value < 2.2e-16
## 
## 
##  Box-Pierce test
## 
## data:  X
## X-squared = 639.62, df = 50, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  X
## X-squared = 643.06, df = 50, p-value < 2.2e-16

Every p-value is below 2.2e-16, so there is overwhelming evidence against white noise.

The statistic barely increases with k, from 587 at k = 5 to 640 at k = 50, because almost all of the dependence sits in the first two lags. This is the cut-off property of an MA(2): adding lags beyond h = 2 adds degrees of freedom without adding signal.

Problem 2

Import the data

p_long  <- scan("sp500_31Dec59-30Dec16.txt")
p_short <- scan("sp500_2009_07_21-2017_07_18.txt")

c(length(p_long), length(p_short))
## [1] 14349  2013

Daily returns

r_long  <- diff(log(p_long))
r_short <- diff(log(p_short))

par(mfrow = c(2, 1))
plot(r_long, type = "l", main = "Returns 1959-2016")
plot(r_short, type = "l", main = "Returns 2009-2017")

par(mfrow = c(1, 1))

If the log prices are a random walk with drift, the returns are white noise plus a constant, so the returns should show no autocorrelation.

acf(r_long, lag.max = 50, main = "rho of returns 1959-2016")

acf(r_short, lag.max = 50, main = "rho of returns 2009-2017")

Box.test(r_long, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  r_long
## X-squared = 87.742, df = 20, p-value = 1.835e-10
Box.test(r_short, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  r_short
## X-squared = 44.191, df = 20, p-value = 0.001419

Both series reject. For 1959 to 2016 the statistic is 87.74 on 20 degrees of freedom with a p-value of 1.8e-10, and for 2009 to 2017 it is 44.19 with a p-value of 0.0014. The 5 percent critical value of a chi-squared distribution on 20 degrees of freedom is 31.4, so both statistics are comfortably beyond it. The returns are therefore not white noise, and the log prices are not a pure random walk in either period. This is the conclusion Lo and MacKinlay reached.

Two qualifications are worth attaching to that.

The individual autocorrelations are very small, well under 0.05 in absolute value. With more than 14000 observations in the long series, even negligible dependence attains significance, so this is evidence against the strict random walk rather than evidence of predictability worth trading on.

The Ljung-Box test also assumes constant variance under the null, and the return plots show obvious volatility clustering, with calm stretches and violent ones. Under that kind of heteroskedasticity the test over-rejects. So part of the rejection may reflect the fact that returns are dependent through their variance rather than correlated in their mean.

Drift estimator and drift line

The drift estimator is the mean of the log returns, which equals (y_n - y_1) / (n - 1) for the log price series y. The drift line is y_1 + mu t.

drift_plot <- function(p, main) {
  y <- log(p)
  n <- length(y)
  t <- 0:(n - 1)
  mu <- mean(diff(y))
  line <- y[1] + mu * t
  plot(t, y, type = "l", main = main, ylab = "log price")
  lines(t, line, col = "red", lwd = 2)
  list(mu = mu, line = line)
}

d_long  <- drift_plot(p_long,  "log SP500 1959-2016 with drift line")

d_short <- drift_plot(p_short, "log SP500 2009-2017 with drift line")

c(mu_long = d_long$mu, mu_short = d_short$mu)
##      mu_long     mu_short 
## 0.0002523836 0.0004706228

The estimated drift is about 0.00025 per day for the long series and 0.00048 for the short one, so roughly 6 percent and 12 percent per year.

Exponential of the drift line

plot(p_long, type = "l", main = "SP500 1959-2016 with exponential trend")
lines(exp(d_long$line), col = "red", lwd = 2)

plot(p_short, type = "l", main = "SP500 2009-2017 with exponential trend")
lines(exp(d_short$line), col = "red", lwd = 2)

Does this give a better fit? Yes, in the sense that matters. Constant drift in logs means proportional growth in levels, so the exponential curve reproduces the shape of the price series in a way no straight line on the untransformed data could.

It is not a close fit in either sample, and for the same reason in both. The drift estimator depends only on the first and last observations, so the curve is pinned at the two end points and is unconstrained in between. For 1959 to 2016 the prices sit below the curve through most of the middle of the sample. For 2009 to 2017 they sit above it from roughly the 200th to the 1500th trading day, by as much as 150 index points.

Problem 3

Subtract the exponential trend from the data

plot(p_long - exp(d_long$line), type = "l",
     main = "SP500 1959-2016 minus exponential trend")
abline(h = 0, col = "grey")

plot(p_short - exp(d_short$line), type = "l",
     main = "SP500 2009-2017 minus exponential trend")
abline(h = 0, col = "grey")

Neither series is stationary, though for somewhat different reasons.

For 1959 to 2016 the deviations are almost imperceptible before about 1980 and very large after 2000, so the variability grows with the level. For 2009 to 2017 the deviations run from roughly -250 to +200 and move in long slow swings rather than fluctuating about a fixed level, so the mean is not constant either.

Subtract the linear trend from the log transformed data

plot(log(p_long) - d_long$line, type = "l",
     main = "log SP500 1959-2016 minus linear trend")
abline(h = 0, col = "grey")

plot(log(p_short) - d_short$line, type = "l",
     main = "log SP500 2009-2017 minus linear trend")
abline(h = 0, col = "grey")

This series is much more regular than anything preceding it, and the gain is in the variability. The fluctuations now have roughly the same amplitude at the start of the sample as at the end, whereas in the previous plot they grew by orders of magnitude. That is the effect of the log transformation, which turns proportional variation into additive variation.

The mean is a separate question. The 1959 to 2016 series wanders between about -0.8 and +0.6 in slow excursions lasting decades, and the 2009 to 2017 series does the same on its own scale. This is what should be expected, since subtracting a deterministic line from a random walk leaves a random walk. The series is close to stationary in the sense of having stabilised variance, which is what makes it look regular, but the level still drifts.