Simulating Long-Horizon Return Regressions with an AR(1) Predictor

# Beta estimates of Regressions on AR(1) Process


# Beta estimates of Regressions on AR(1) Process












I am currently working through the paper The Myth of Long-Horizon Predictability [1] and I got stuck in reproducing the empirical results in Section 1.4.

It is my understanding that time series of length 75 years (with each year corresponding to one sample) are simulated under the assumption of no predictability and with an AR(1) process on the predictor $X_t$, with \begin{align} R_{t,t+1} &= \epsilon_{t,t+1} \\ X_{t+1} &= \rho X_t + u_{t,t+1} \end{align} where $R_{t,t+1}$ is the response variable, consisting of annual returns from year $t$ to $t+1$, $X_t$ is the predictor variable at time $t$ and $\epsilon_{t,t+1} \sim N(0,\sigma_{\epsilon}^2)$ as well as $u_{t,t+1} \sim N(0,\sigma_u^2)$.

The text is explicit in using the labels $\sigma_{u}$ and $\sigma_{\epsilon}$ for the standard deviations of $u_{t,t+1}$ and $\epsilon_{t,t+1}$. Additionally, the text is using the label $\sigma_{\epsilon u}$ for the correlation between $u_{t,t+1}$ and $\epsilon_{t,t+1}$.

The simulations involve 100,000 replications each.

In footnote 8, the values used in the simulations are listed: $\rho = 0.953, \sigma_{\epsilon}=0.202, \sigma_{u}=0.154$ and $\sigma_{\epsilon u} = -0.712$.

Five different regressions are run:

\begin{align} R_{t,t+1} &= \alpha_1 + \beta_1 X_t + \epsilon_{t,t+1} \\ \vdots \\ R_{t,t+j} &= \alpha_j + \beta_j X_t + \epsilon_{t,t+j} \\ \vdots \\ R_{t,t+J} &= \alpha_J + \beta_J X_t + \epsilon_{t,t+J} \end{align}

where $j\in\{1,\dots,J\}$ with $J=5$. The return $R_{t,t+j}$ is the return from $t$ to $t+j$, the returns are thus overlapping for $j>2$.

What I am stuck with is the following: My mean for $\beta_1$ across the simulations is $-0.1610587$, whereas in the text it seems to be $0.055$ (Table 1, Horizon 1, Column Mean). A negative value for $\beta_1$ makes sense to me, as $u$ and $\epsilon$ are negatively correlated by construction. Am I missing something about the setup of the simulation?

EDIT: Due to the kind advice of a friend, I have modified the generateSeries function below in the following way: Instead of taking the first 75 samples, 1000+75 samples of $X_t$ are generated and the last 75 samples are taken to ensure that the process has converged. This indeed solves some issues, but the problem of the question posed here still remains.

Below is my R code for the simulations and the regressions.

```
library(MASS)
library(Matrix)
library(data.table)

library(doParallel)
registerDoParallel(cores = 10)

generateSeries <- function(samples = 75, sdX = 0.154, sdY = 0.202, corrXY = -0.712, rho = 0.953) {
  buffer <- 1000

  varX = sdX^2
  varY = sdY^2
  covXY = corrXY * sdX * sdY

  Sigma = matrix(c(varX, covXY, covXY, varY), nrow=2)

  data = as.data.frame(mvrnorm(n=samples+buffer, mu=c(0, 0), Sigma=Sigma, empirical=TRUE))
  setDT(data)

  u <- data$V1
  epsilon <- data$V2

  X0 = 0
  X <- vector()

  X[[1]] <- rho * X0 + u[[1]]

  for (i in 2:(samples+buffer)) {
    X[[i]] <- rho * X[[i-1]] + u[[i]]
  }

  timeSeries <- data.table(R = tail(epsilon, samples), X = tail(X, samples))

  timeSeries
}

M <- 100000

beta.estimators <- foreach(i = 1:M, .combine = 'rbind') %dopar% {
  beta.estimators.instance <- vector()

  single.instance <- generateSeries()
  for (j in 1:5) {
    single.instance[, R.j := frollsum(R, n = j, align = "left")]

    regression.model <- lm(R.j ~ X, data = na.omit(single.instance))

    beta.estimators.instance[j] <- regression.model$coefficients[2]
  }

  beta.estimators.instance
}

beta.estimators <- as.data.frame(beta.estimators)
setDT(beta.estimators)
lapply(beta.estimators, mean)
lapply(beta.estimators, function(x) cor(beta.estimators$V1, x))
```

[1] The Myth of Long-Horizon Predictability, Jacob Boudoukh, Matthew Richardson and Robert F. Whitelaw, The Review of Financial Studies, Vol. 21, No. 4 (Jul., 2008), pp. 1577-1605