Modelling relation between two persistent AR(1) processes? Is my approach reasonable and how to adjust standard errors?

Question

I want to run a simulation in which I want to find out whether there is a relation between the independent variable $x_t$ and the dependent variable $y_t$. I.e., in the following regression I want to find out if $\beta$ is signficantly different from zero:

$$ y_t = \beta x_t + \epsilon_t. $$

In empirical data, $y_t$ is stationary, but close to a random walk.I want to model it as an AR(1) process. Of course, I could just model $x_t$ as an AR(1) process and if $\beta$ is different from zero, $y_t$ would end up as such a process as well. However, I also want to consider the case in which $\beta=0$ and I don't want to confuse the reader with changing descriptions for different scenarios (i.e., "if $\beta=1$ I simulate $y_t$ from the above equation, if $\beta=0$, I simulate $y_t$ as an AR(1) process.")

I came up with the following solution to this dilemma. Why not model both $y_t$ and $x_t$ as AR(1) processes with potentially correlated errors. That is:

$$ y_{t+1} = \tau_y y_t + u_{t+1}, \quad 0 < \tau_y < 1, \quad u_{t+1} \sim N(0, \sigma_u^2) $$ and $$ x_{t+1} = \tau_x x_t + v_{t+1} , \quad 0 < \tau_x < 1, \quad v_{t+1} \sim N(0, \sigma_{v}^2), $$

where

$$ cov\left(\begin{bmatrix} v_{t+1} \\ u_{t+1} \end{bmatrix}, \begin{bmatrix} v_{t+1} & u_{t+1} \end{bmatrix}\right) = \begin{bmatrix} \sigma_v^2 & \sigma_{v} \sigma_u \rho_{u,v} \\ \sigma_{v} \sigma_u \rho_{u,v} & \sigma_u^2 \end{bmatrix} $$

Since $\beta = \frac{Cov(y, x)}{\sigma_x^2}$ and

$$ \begin{aligned} Cov(y, x) &= E[y_{t+1}x_{t+1}] \\ &= E[(\tau_y y_t + u_{t+1})(\tau_x x_t + v_{t+1})] \\ &= \tau_y \tau_x E[y_t x_t] + E[u_{t+1} v_{t+1}] \\ &= \tau_y \tau_x Cov(y, x) + \sigma_{u} \sigma_{v} \rho_{v, u} \\ \end{aligned} $$

we get by rearranging

$$ \beta = \frac{\sigma_u}{\sigma_{v}}\rho_{v, u} \frac{1 - \tau_x^2}{1-\tau_x \tau_y}. $$

Fair enough, this is a little bit complicated, but it allows me to control $\beta$ without changing the structure of either $y_t$ or $x_t$ and both those processes are AR(1) in this setup.

I basically have two questions:

1. Is this a reasonable approach given the requirement that preferably both $x_t$ and $y_t$ should be AR(1) processes. The interpretation of the relation should really be just like in a simple OLS setup, i.e. has $x_t$ an impact on $y_t$ (I'm not concerned about causality here, just relation). Would you consider my setup still as a reasonable way of modelling it? (I just don't make a direct link, but use the error structure for that. I don't see a problem with this approach. If $x_t$ and $y_t$ are related via correlated shocks, so be it.)

2. I already simulated this and it works fine. However, I also want to obtain correct confidence intervals for each run. That is, I want to run the regression in the first equation, get $\widehat{\beta}$ and also standard errors that should be reasonable. Now I noticed that the errors are highly autocorrelated and this autocorrelation seems to be identical to $\tau_y$. However, I could not formalize it. So it would be great to know how I would have to adjust the standard errors. (It would be extra awesome if I could get a hint how to implement that in R.)

   set.seed(123)
   library(MASS)
   ### Set start values
   nrT <- 1e5
   burnin <- 1e3
   sd_y_shock  <- 1
   sd_x_shock  <- 1
   corr_x_y    <- 0.5
   tau_y  <- 0.8
   tau_x  <- 0.1
   ### Simulate the correlated shocks
   shocks <- mvrnorm(nrT + burnin, 
                     mu = c(0,0), 
                     Sigma = matrix(c(sd_y_shock^2, corr_x_y * sd_y_shock * sd_x_shock, 
                                      corr_x_y * sd_y_shock * sd_x_shock, sd_x_shock^2), 
                                    nrow=2))
   
   vec_y <- arima.sim(list(order = c(1, 0, 0), ar = tau_y),
                       n = nrT + burnin, 
                       innov = shocks[, 1]) 
   
   vec_x <- arima.sim(list(order = c(1, 0, 0), ar = tau_x),
                      n = nrT + burnin, 
                      innov = shocks[,2])
   ### Check that formula derived above is correct; the two should be similar
   sd_y_shock/sd_x_shock * corr_x_y * (1 - tau_x^2)/(1 - tau_x * tau_y) 
   coef(lm(vec_y ~ vec_x))[2]
   ### Plot ACF
   # Note that, independent from corr_x_y and tau_x, the autocorrelation structure seems
   # always to be tau_y
   acf(lm(vec_y ~ vec_x)$resid)
   

Answer 1

Given the information in the question I do not see any problems with such simulation set-up. You are modelling VAR(1) process and then explore the contemporaneous relationship between its components. Perfectly reasonable, although I would double check that this conforms to your empirical model.

As for the second question, the best way would be to figure out the exact covariance structure of the stationary process $\varepsilon_t=y_t-\beta x_t$ and then use GLS instead of OLS with the covariance matrix exactly specified. This way you would get the efficient estimates of $\beta$ and subsequently the smallest possible standard errors. Another way is to use autocorrelation-robust standard errors. In R you can calculate them using package sandwich:

library(lmtest)
library(sandwich)
mod <- lm(vec_y ~ vec_x)
vv <- vcovHAC(mod)
coeftest(mod,vcov=vv)


  t test of coefficients:

              Estimate Std. Error  t value Pr(>|t|)    
(Intercept) -0.0093830  0.0111428  -0.8421   0.3998    
vec_x        0.5333422  0.0052502 101.5857   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Here I used function coeftest from the package lmtest to get the usual summary output. Note that for large sample sizes the calculation of robust covariance matrix can take up some time.