R: Constructing VAR lag structure

Question

I am searching for a function or code, that can create the lag structure of a VAR(p) model, since I need it to run some other functions on it. I've written a function that creates the lag structure matrix of one time series in the univariate case, but I've no idea to write a function for the multivariate case. Or are there any implementations of it in R ?

A short reproducible example would be 3 dimensional VAR of:

Y <- matrix(1:30, ncol=3)

and assume I want construct lag matrix of an 3 dimensional VAR(2) model

EDIT

The final goal is to estimate an 5-dimnensional VAR(20) lasso Regression with the package glmnet

Answer 1

I am searching for a function or code, that can create the lag structure of a VAR(p) model, ...

Here is one way

n <- 1000
p <- 3
set.seed(1)
Y <- matrix(rnorm(n * p), ncol = p)
X <- cbind(lag_1 = rbind(NA, head(Y, -1)), lag_2 = rbind(NA, NA, head(Y, -2)))
head(X)
#R         [,1]    [,2]     [,3]    [,4]    [,5]    [,6]
#R [1,]      NA      NA       NA      NA      NA      NA
#R [2,] -0.6265  1.1350 -0.88615      NA      NA      NA
#R [3,]  0.1836  1.1119 -1.92225 -0.6265  1.1350 -0.8861
#R [4,] -0.8356 -0.8708  1.61970  0.1836  1.1119 -1.9223
#R [5,]  1.5953  0.2107  0.51927 -0.8356 -0.8708  1.6197
#R [6,]  0.3295  0.0694 -0.05585  1.5953  0.2107  0.5193
head(Y)
#R         [,1]    [,2]     [,3]
#R [1,] -0.6265  1.1350 -0.88615
#R [2,]  0.1836  1.1119 -1.92225
#R [3,] -0.8356 -0.8708  1.61970
#R [4,]  1.5953  0.2107  0.51927
#R [5,]  0.3295  0.0694 -0.05585
#R [6,] -0.8205 -1.6626  0.69642

You can use it to e.g., perform conditional log-likelihood estimation as follows

colnames(X) <- c(paste0("X_lag_1", 1:3), paste0("X_lag_2", 1:3))
lm.fit(x = X[-(1:2), ], y = Y[-(1:2), ])$coefficients
#R               [,1]       [,2]     [,3]
#R X_lag_11 -0.041859  0.0048129 -0.02624
#R X_lag_12 -0.013648  0.0005279 -0.02306
#R X_lag_13  0.037641 -0.0087508  0.05377
#R X_lag_21 -0.033324  0.0637967 -0.05455
#R X_lag_22 -0.007617  0.0384764 -0.08435
#R X_lag_23  0.006812 -0.0420907 -0.02983

---

Update

and you can make a function that does this like

n <- 1000
p <- 2
set.seed(1)
Y <- matrix(rnorm(n * p), ncol = p)

lag_series <- function(Y, max_lag)
  do.call(cbind, lapply(1:max_lag, function(i)
    do.call(rbind, c(as.list(rep(NA, i)), list(head(Y, -i))))))

head(lag_series(Y, 1))
#R            [,1]        [,2]
#R [1,]         NA          NA
#R [2,] -0.6264538  1.13496509
#R [3,]  0.1836433  1.11193185
#R [4,] -0.8356286 -0.87077763
#R [5,]  1.5952808  0.21073159
#R [6,]  0.3295078  0.06939565
head(lag_series(Y, 2))
#R          [,1]        [,2]       [,3]       [,4]
#R [1,]         NA          NA         NA         NA
#R [2,] -0.6264538  1.13496509         NA         NA
#R [3,]  0.1836433  1.11193185 -0.6264538  1.1349651
#R [4,] -0.8356286 -0.87077763  0.1836433  1.1119318
#R [5,]  1.5952808  0.21073159 -0.8356286 -0.8707776
#R [6,]  0.3295078  0.06939565  1.5952808  0.2107316
head(lag_series(Y, 3))
#R             [,1]        [,2]       [,3]       [,4]       [,5]       [,6]
#R [1,]         NA          NA         NA         NA         NA         NA
#R [2,] -0.6264538  1.13496509         NA         NA         NA         NA
#R [3,]  0.1836433  1.11193185 -0.6264538  1.1349651         NA         NA
#R [4,] -0.8356286 -0.87077763  0.1836433  1.1119318 -0.6264538  1.1349651
#R [5,]  1.5952808  0.21073159 -0.8356286 -0.8707776  0.1836433  1.1119318
#R [6,]  0.3295078  0.06939565  1.5952808  0.2107316 -0.8356286 -0.8707776

#...