bpfaff / vars

Multivariate Time Series Models: VAR, SVAR and SVEC
40 stars 23 forks source link

SVAR irf error with Amat called "A" + some feature requests #11

Open SebKrantz opened 3 years ago

SebKrantz commented 3 years ago

There is an error in SVAR estimation when matrices are called "A" or "B", which is presumably due to the depreciated A() and B() functions:

library(vars)
data(Canada)
var.2c <- VAR(Canada, p = 2, type = "const")
A <- diag(4)
diag(A) <- NA
A[2, 1] <- NA
A[4, 1] <- NA
## Estimation method scoring
mod <- SVAR(x = var.2c, estmethod = "scoring", Amat = A, Bmat = NULL,
            max.iter = 100, maxls = 1000, conv.crit = 1.0e-8) 
irf(mod)
#> Error in SVAR(x = varboot, estmethod = "scoring", Amat = A, Bmat = NULL, : 
#> No parameters provided for optimisation, i.e.
#> neither 'Amat' nor 'Bmat' do contain NA elements.
## Estimation method direct
mod2 <- SVAR(x = var.2c, estmethod = "direct", Amat = A, Bmat = NULL,
             hessian = TRUE, method="BFGS") 
irf(mod2)
#> Error in SVAR(x = varboot, estmethod = "direct", Amat = A, Bmat = NULL, : 
#> No parameters provided for optimisation, i.e.
#> neither 'Amat' nor 'Bmat' do contain NA elements.

Created on 2021-10-02 by the reprex package (v0.3.0)

Then, I would like to request asymptotic computation of IRF's to be added as a feature, since bootstrapping can take a while with some models and is also succeptible to stoachastic singularity etc. in small samples. I also find myself writing my own IRF plotting functions a lot of times, the main reason being that comparing IRF's on the same y-axis scale just doesn't make sense for some models, and I'm more interested in the shape than seeing the order of magnitude of the response. So the ability to put ylim = NULL or ylim = NA in plot.irf and get individual scales is very important for these plot functions to be useful to me. (See here for a earlier request of this feature by @MatthieuStigler ). A nother request I would have regarding IRFs is an option to scale them so that the response of variable to its own shock in period 0 is 1, which facilitates inpterpretation of IRFs in a log-level VAR as elasticities. Eviews has an add-in called sirf to do this. So far I have used the following function:

unitize_irf <- function(x) {
  nam <- x$impulse
  if(length(x$Lower)) {
    x[1:3] <- collapse::t_list(setNames(lapply(nam, function(i) {
      irfi <- x$irf[[i]]
      m <- irfi[1, i]
      list(irf = irfi / m, 
           Lower = x$Lower[[i]] / m, 
           Upper = x$Upper[[i]] / m)
    }), nam))
  } else {
    x$irf <- setNames(lapply(nam, function(i) {
      irfi <- x$irf[[i]]
      irfi / irfi[1, i]
    }), nam)
  }
  return(x)
}

Finally, it would be cool if the optimization function could be passed as an argument to SVAR etc. allowing the user to pass alternative (and more powerful) optimization functions (like optimx etc.). Thanks! (I'm happy to implement and send PRs for some of these requests, but would first like to hear your thoughts).