Enhancement / feature request - identification strategy

[msh855]

Hi, Thank you for this library, it is great so far.

However, I think this library has way more potentials if there is some guidance on how to implement standard identification strategies.

I am more familiar with ECB's BEAR tool box, and I think the bsvars library can become way more user friendly if the identification functions get more simplified. For example, I am struggling to understand how I can implement some standard identification strategies like:

- Cholesky decomposition
- Sign Restrictions
- Zero restrictions

Also, can the library handle anticipated vs unanticipated shocks as in Structural scenario analysis with SVARs ?

[donotdespair]

Hey Moustafa @msh855!

Thank you so much for your comment and encouragement! I really appreciate your support and feedback. You included a big ask, but many of these points are planned for implementation. So, here we go:

1. The next points to develop is to make the package user-friendly with summary and plot functions. More on that here
2. The plan also includes more resources such as instructive examples in the website, a vignette, and so on! It's all coming up!
3. I do not use the Cholesky decomposition as the package estimates the structural model explicitly. But the lower-triangular SVAR can be easily set up following all the examples provided in the package, such as this one:

   
   library(bsvars)

upload data

data(us_fiscal_lsuw) # upload dependent variables data(us_fiscal_ex) # upload exogenous variables

specify the model and set seed

specification = specify_bsvar$new(us_fiscal_lsuw, p = 4, exogenous = us_fiscal_ex) set.seed(123)

run the burn-in

burn_in = estimate(specification, 10)

estimate the model

posterior = estimate(burn_in, 10, thin = 2)

compute impulse responses 2 years ahead

irf = compute_impulse_responses(posterior, horizon = 8)

The code will estimate a lower-triangular SVAR with exogenous variables. Executing line:
```{r}
specification  = specify_bsvar$new(us_fiscal_lsuw, p = 4, exogenous = us_fiscal_ex)

specifies the model and displays the notification of the lower-triangularity.

4. If you want different zero restrictions you have to specify the matrix argument B of the specify_*() function as in teh example:

   > data(us_fiscal_lsuw)
   > data(us_fiscal_ex) 
   > B = matrix(FALSE,3,3)
   > B = matrix(TRUE,3,3)
   > B[1,3] = B[2,3] = FALSE
   > B
    [,1] [,2]  [,3]
   [1,] TRUE TRUE FALSE
   [2,] TRUE TRUE FALSE
   [3,] TRUE TRUE  TRUE
   > specification  = specify_bsvar_sv$new(us_fiscal_lsuw, B= B, p = 4, exogenous = us_fiscal_ex)

   The elements of the B above that are set to TRUE are estimated and those set to FALSE are restricted to zero.

5. The package does not offer sign restrictions. We are just beginning our work on another package that would do that as well as the narrative restrictions!
6. We are just beginning our work on another package that can do things like in the Structural Scenario Analysis paper.

In short, the bsvars package already includes many of the functionality you asked about. More documentation and examples are coming soon. More functionality is upcoming as well, but please be patient. I believe in less than two years, all these things will be covered by my packages :)

I wish you successful implementations, and I am looking forward to hearing about your projects and outputs!

Cheers, Tomasz @donotdespair

[msh855]

i really think this bit would create a lot of confusion:

# run the burn-in
burn_in        = estimate(specification, 10)

# estimate the model
posterior      = estimate(burn_in, 10, thin = 2)

I think better change the names between burn-in and estimations.

[donotdespair]

Do you mean changing the names of the estimate method or the LHS object names, that is, estimations instead of posterior? This is a relatively easy bit to change in the doc! 😸

[msh855]

Yes, I mean the estimate method. It is used to do two different things and I find it confusing. In python, the way this would have worked would be something like this:

 burn_in   = estimate(specification, 10)
posterior  = burn_in(10, thin = 2)

[donotdespair]

Hey @msh855 thanks for this clarification. Yeah, I get that this can be confusing at first. The idea here is that the models used here are non-linear and the MCMC convergence might be slow for some data. the construction of the estimate method is designed to give the possibility of continuing running the MCMC if the last round does not provide convergent outcome.

For instance, lets run a 1000 draws for the burn in and 10000 for the estimation:

# run the burn-in
burn_in        = estimate(specification, 1000)

# estimate the model
posterior      = estimate(burn_in, 10000, thin = 2)

The estimate method is designed in such a way that if a posterior output is provided as the first argument burn_in or posterior then the algo takes the last MCMC draw and simply continues running the MCMC. If the model specification is provided there specification then it takes the starting values from this object and begins the MCMC estimation.

Now, let's suppose that one takes the output in posterior and figures out that the MCMC has not converged. Then, we can simply run the following line:

posterior_next      = estimate(posterior, 10000, thin = 2)

To continue the estimation. This MCMC, saved in posterior_next, has more chances of converging.

I designed it this way based on my experience of working with multivariate dynamic Bayesian models. Sometimes, it's just so surprising how long the convergence takes, despite efficient algos. The bigger the model, the larger the number of parameters, the slower the convergence.

I will not change the name of the routine estimate but I will work on materials explaining it. Say, a vignette.

How does this sound to you?

Cheers, T