Vector autoregression with interval lag terms in R? 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.
 A: 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:


*

*min_lag ($L_{min} = 100$): the assumed minimum separation between both series.

*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:


*

*Recent lags for x and y, from 1 to nb_lags.

*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
set.seed(1313)

# 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]

plot(x)
plot(y)

# 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
library(glmnet)

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!

# END

