I'd like to perform vector autoregression on a two variable system. I know that the signals $x$ and $y$ have a time lag of > 100 time points, and thus any fit with that many time lag parameters is likely to be bad. While I could take the data at a lower frequency, the data is produced at one given frequency. Striding the data (taking every $s$-th observation) is possible, but it seems like it would be throwing out quite a bit of potentially useful data, and would make the analysis far more susceptible to the noise - in my experience, if I stride and then vary the time point I start from just a little, I get very different results.

I believe the right thing to do is obviously to fit this equation:

$$x_t = x_{t-s} + x_{t-2s} + x_{t-3s} + ... + x_{t-ps} + y_{t-s} + y_{t-2s} + y_{t-3s} + ... + y_{t-ps}$$

However, I can't seem to find existing code in R to do that. VAR and ar don't seem to have options to do this, as far as I can find. Does anyone have any suggestions, or do I need to code this up from scratch? Is there a name for this type of model that I'm unaware of?

I've thought about applying smoothing with the window size of $s$ and then striding, and when I do that, I get a significant result irregardless of if I vary the time point I begin at. To me, this indicates that there is certainly a significant time lag contribution (using the Wald test), but I think using the above equation is more rigorous.

EDIT: I wrote R code that allows you do to OLS fitting of vector autoregression when the lag is t - i*s. It functionally gives you everying ar() or VAR() does, but has the added flexibility of the s parameter. I'm working it into my own workflow, so it's a bit optimized, but if anyone needs it, comment and I'll send it over in a more workable form.

  • $\begingroup$ I've made some formatting changes and edits to clarify your mathematics. Please check it says what you intended. $\endgroup$ – Glen_b Oct 7 '14 at 23:04
  • $\begingroup$ Hi @Mike, I'm looking at some data and have lag terms with some variability. I'd love to hear more details about your approach or take a look at your code if you've still got it. $\endgroup$ – 7yl4r May 4 '16 at 15:47
  • $\begingroup$ I'd try regularizing the solution using the Lasso. Just estimate each equation separately using, for example, glmnet, throwing all the available lags. $\endgroup$ – mbiron May 5 '16 at 21:13

I just had time to go over this more carefully. I know this is 1 year and a half old, but it's still an interesting (and unanswered) question.

Here is a toy example of my approach to solving this using a regularized regression estimator like Lasso. The only important parameters in this approach are:

  1. min_lag ($L_{min} = 100$): the assumed minimum separation between both series.
  2. nb_lags ($\Delta L = 20$): how many lags to scan (to include in the data matrix for regression).

There is also sigma_e ($\sigma_{\varepsilon} = 0.1$), which is the amount of noise in the series we will simulate.

Using these parameters, the data matrix will include:

  1. Recent lags for x and y, from 1 to nb_lags.
  2. Older lags for x and y, from min_lag to min_lag + nb_lags - 1.

This amounts to 80 regressors for each equation.

On the other hand, the process that generates the data consists of the equations (four coefficients):

$x_t = 0.9x_{t-1} + 0.1y_{t-110} + \varepsilon_t^x$

$y_t = 0.9y_{t-1} + 0.1x_{t-110} + \varepsilon_t^y$

where $\{\varepsilon_t^x, \varepsilon_t^y\}_{t\in\{1...N\}}\overset{iid}{\sim} N(0, \sigma_{\varepsilon}^2)$, which we try to discover by fitting the Lasso to

$x_t = \beta_0^x + \sum_{l=1}^{\Delta L}{\beta_l^x x_{t-l}} + \sum_{l=L_{min}}^{L_{min}+\Delta L-1}{\beta_l^x x_{t-l}} + \sum_{l=L_{min}}^{L_{min}+\Delta L-1}{\gamma_l^x y_{t-l}} + u_t^x$

$y_t = \beta_0^y + \sum_{l=1}^{\Delta L}{\beta_l^y y_{t-l}} + \sum_{l=L_{min}}^{L_{min}+\Delta L-1}{\beta_l^y y_{t-l}} + \sum_{l=L_{min}}^{L_{min}+\Delta L-1}{\gamma_l^y x_{t-l}} + u_t^y$

Note that I've left out the recent lags of the opposed variable in each equation because we know beforehand that they are not relevant. I still leave in the old lags of the same variable, though.

Using the code shown below, the estimated coefficients are non-zero only for the correct lags (1 and 110 in both regressions) and with magnitudes

  • $\hat{\beta_1^x} = 0.86$
  • $\hat{\gamma_{110}^x} = 0.09$
  • $\hat{\beta_1^y} = 0.86$
  • $\hat{\gamma_{110}^y} = 0.05$

which is pretty close IMHO to the real coefficients (you would expect them to be lower than the real value because of the regularization imposed). I imagine including more data could yield a more accurate answer, but for me what's amazing is the fact that the Lasso can pick up the right lags in both series.

Hope this helps.

# seed

# create data
N = 1000
x = numeric(N)
y = numeric(N)

xy_lag = 110 # position of the lag between series
sigma_e = 0.1

# fill first data randomly
x[1:xy_lag] = rnorm(xy_lag, 0, sigma_e)
y[1:xy_lag] = rnorm(xy_lag, 0, sigma_e)

# construct the series
for (t in (xy_lag+1):N) {
  x[t] = 0.9 * x[t - 1] + 0.1 * y[t - xy_lag] + rnorm(1, 0, sigma_e)
  y[t] = 0.9 * y[t - 1] + 0.1 * x[t - xy_lag] + rnorm(1, 0, sigma_e)

# drop warm-up data
x = x[(xy_lag+1):N]
y = y[(xy_lag+1):N]


# create data matrix for regression
min_lag = 100 # assumed minimum separation
nb_lags = 20 # number of lags to scan

# util
multi_lag = function(x, k1, k2){
  return(sapply((k1):(k2), function(n){dplyr::lag(x, n)}))

data_mat = cbind(multi_lag(x, 1, nb_lags), # recent lags of x
                 multi_lag(x, min_lag, min_lag + nb_lags - 1), # old lags of x
                 multi_lag(y, 1, nb_lags), # recent lags of y
                 multi_lag(y, min_lag, min_lag + nb_lags - 1)) # old lags of y
colnames(data_mat) = c(paste0("x_l", 1:nb_lags),
                       paste0("x_l", min_lag:(min_lag + nb_lags - 1)),
                       paste0("y_l", 1:nb_lags),
                       paste0("y_l", min_lag:(min_lag + nb_lags - 1)))

# remove first rows (NA's)
data_mat = data_mat[-(1:(min_lag + nb_lags - 1)),]
x = x[-(1:(min_lag + nb_lags - 1))]
y = y[-(1:(min_lag + nb_lags - 1))]

# fit lasso

fit_x_cv = cv.glmnet(x = data_mat[, -((2*nb_lags+1):(3*nb_lags))], # don't need recent lags of y
                     y = x)
coef(fit_x_cv) # correct lags, coefficients are close!

fit_y_cv = cv.glmnet(x = data_mat[, -(1:nb_lags)], # don't need recent lags of x
                     y = y)
coef(fit_y_cv) # correct lags, coefficients are close!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.