diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 66177ad8b..fb70d98ae 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -8886,6 +8886,15 @@ Dynare also has the ability to estimate Bayesian VARs:
     See ``bvar-a-la-sims.pdf``, which comes with Dynare distribution,
     for more information on this command.
 
+.. command:: bvar_irf ;
+
+    |br| Computes the impulse responses of an estimated BVAR model, using
+    Minnesota priors.
+
+    See ``bvar-a-la-sims.pdf``, which comes with Dynare distribution,
+    for more information on this command.
+
+
 Estimation based on moments
 ===========================
 
diff --git a/license.txt b/license.txt
index be78e6d1a..72cb108f8 100644
--- a/license.txt
+++ b/license.txt
@@ -88,8 +88,8 @@ License: public-domain-aim
 
 Files: matlab/optimization/bfgsi1.m matlab/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m
  matlab/optimization/numgrad3.m matlab/optimization/numgrad3_.m matlab/optimization/numgrad5.m
- matlab/optimization/numgrad5_.m matlab/optimization/csminwel1.m matlab/bvar_density.m
- matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m matlab/partial_information/qzswitch.m
+ matlab/optimization/numgrad5_.m matlab/optimization/csminwel1.m matlab/+bvar/density.m
+ matlab/+bvar/toolbox.m matlab/partial_information/PI_gensys.m matlab/partial_information/qzswitch.m
  matlab/partial_information/qzdiv.m
 Copyright: 1993-2009 Christopher Sims
            2006-2023 Dynare Team
diff --git a/matlab/bvar_density.m b/matlab/+bvar/density.m
similarity index 94%
rename from matlab/bvar_density.m
rename to matlab/+bvar/density.m
index dd469a79b..ba69bf9bb 100644
--- a/matlab/bvar_density.m
+++ b/matlab/+bvar/density.m
@@ -1,5 +1,5 @@
-function bvar_density(maxnlags)
-% function bvar_density(maxnlags)
+function density(maxnlags)
+% function density(maxnlags)
 % computes the density of a bayesian var
 %
 % INPUTS
@@ -12,7 +12,7 @@ function bvar_density(maxnlags)
 %    none
 
 % Copyright © 2003-2007 Christopher Sims
-% Copyright © 2007-2017 Dynare Team
+% Copyright © 2007-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -34,7 +34,7 @@ global oo_
 oo_.bvar.log_marginal_data_density=NaN(maxnlags,1);
 
 for nlags = 1:maxnlags
-    [ny, nx, posterior, prior] = bvar_toolbox(nlags);
+    [ny, nx, posterior, prior] = bvar.toolbox(nlags);
     oo_.bvar.posterior{nlags}=posterior;
     oo_.bvar.prior{nlags}=prior;
 
diff --git a/matlab/bvar_forecast.m b/matlab/+bvar/forecast.m
similarity index 80%
rename from matlab/bvar_forecast.m
rename to matlab/+bvar/forecast.m
index 1ceb1a799..2ca24f4dd 100644
--- a/matlab/bvar_forecast.m
+++ b/matlab/+bvar/forecast.m
@@ -1,5 +1,5 @@
-function bvar_forecast(nlags)
-% function bvar_forecast(nlags)
+function forecast(nlags)
+% function forecast(nlags)
 % builds forecats for a bvar model
 %
 % INPUTS
@@ -11,7 +11,7 @@ function bvar_forecast(nlags)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2007-2017 Dynare Team
+% Copyright © 2007-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,9 +31,9 @@ function bvar_forecast(nlags)
 global options_ oo_ M_
 
 if options_.forecast == 0
-    error('bvar_forecast: you must specify "forecast" option')
+    error('bvar.forecast: you must specify "forecast" option')
 end
-[ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
+[ny, nx, posterior, prior, forecast_data] = bvar.toolbox(nlags);
 
 sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic);
 sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);
@@ -116,7 +116,7 @@ sims_with_shocks_up_conf = sims_with_shocks_sort(:, :, sort_idx(2));
 
 OutputDirectoryName = CheckPath('graphs',M_.dname);
 
-dyn_graph=dynare_graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});
+dyn_graph=bvar.graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});
 
 for i = 1:ny
     dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
@@ -166,20 +166,20 @@ for i = 1:length(options_.varobs)
     name = options_.varobs{i};
 
     sims = squeeze(sims_with_shocks(:,i,:));
-    eval(['oo_.bvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']);
-    eval(['oo_.bvar.forecast.with_shocks.Median.' name ' = median(sims, 2);']);
-    eval(['oo_.bvar.forecast.with_shocks.Var.' name ' = var(sims, 0, 2);']);
-    eval(['oo_.bvar.forecast.with_shocks.HPDsup.' name ' = sims_with_shocks_up_conf(:,i);']);
-    eval(['oo_.bvar.forecast.with_shocks.HPDinf.' name ' = sims_with_shocks_down_conf(:,i);']);
+    oo_.bvar.forecast.with_shocks.Mean.(name) = mean(sims, 2);
+    oo_.bvar.forecast.with_shocks.Median.(name) = median(sims, 2);
+    oo_.bvar.forecast.with_shocks.Var.(name) = var(sims, 0, 2);
+    oo_.bvar.forecast.with_shocks.HPDsup.(name) = sims_with_shocks_up_conf(:,i);
+    oo_.bvar.forecast.with_shocks.HPDinf.(name) = sims_with_shocks_down_conf(:,i);
 
     sims = squeeze(sims_no_shock(:,i,:));
-    eval(['oo_.bvar.forecast.no_shock.Mean.' name ' = sims_no_shock_mean(:,i);']);
-    eval(['oo_.bvar.forecast.no_shock.Median.' name ' = sims_no_shock_median(:,i);']);
-    eval(['oo_.bvar.forecast.no_shock.Var.' name ' = var(sims, 0, 2);']);
-    eval(['oo_.bvar.forecast.no_shock.HPDsup.' name ' = sims_no_shock_up_conf(:,i);']);
-    eval(['oo_.bvar.forecast.no_shock.HPDinf.' name ' = sims_no_shock_down_conf(:,i);']);
+    oo_.bvar.forecast.no_shock.Mean.(name) = sims_no_shock_mean(:,i);
+    oo_.bvar.forecast.no_shock.Median.(name) = sims_no_shock_median(:,i);
+    oo_.bvar.forecast.no_shock.Var.(name) = var(sims, 0, 2);
+    oo_.bvar.forecast.no_shock.HPDsup.(name) = sims_no_shock_up_conf(:,i);
+    oo_.bvar.forecast.no_shock.HPDinf.(name) = sims_no_shock_down_conf(:,i);
 
-    if exist('rmse')
-        eval(['oo_.bvar.forecast.rmse.' name ' = rmse(i);']);
+    if exist('rmse','var')
+        oo_.bvar.forecast.rmse.(name) = rmse(i);
     end
 end
diff --git a/matlab/dynare_graph_init.m b/matlab/+bvar/graph_init.m
similarity index 100%
rename from matlab/dynare_graph_init.m
rename to matlab/+bvar/graph_init.m
diff --git a/matlab/bvar_irf.m b/matlab/+bvar/irf.m
similarity index 84%
rename from matlab/bvar_irf.m
rename to matlab/+bvar/irf.m
index ec8d8e0f7..05ec9b31e 100644
--- a/matlab/bvar_irf.m
+++ b/matlab/+bvar/irf.m
@@ -1,5 +1,6 @@
-function bvar_irf(nlags,identification)
-% builds IRFs for a bvar model
+function irf(nlags,identification)
+% function irf(nlags,identification)
+% builds IRFs for a BVAR model
 %
 % INPUTS
 %    nlags            [integer]     number of lags for the bvar
@@ -11,7 +12,7 @@ function bvar_irf(nlags,identification)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2007-2017 Dynare Team
+% Copyright © 2007-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -34,7 +35,7 @@ if nargin==1
     identification = 'Cholesky';
 end
 
-[ny, nx, posterior, prior] = bvar_toolbox(nlags);
+[ny, nx, posterior, prior] = bvar.toolbox(nlags);
 
 S_inv_upper_chol = chol(inv(posterior.S));
 
@@ -77,6 +78,8 @@ for draw=1:options_.bvar_replic
         StructuralMat = Sigma_lower_chol;
     elseif strcmpi(identification,'SquareRoot')
         StructuralMat = sqrtm(Sigma);
+    else
+        error('Unknown identification option. Valid choices are ''Cholesky'' and ''SquareRoot''.')
     end
 
     % Build the IRFs...
@@ -137,10 +140,10 @@ for i=1:ny
     shock_name = options_.varobs{i};
     for j=1:ny
         variable_name = options_.varobs{j};
-        eval(['oo_.bvar.irf.Mean.' variable_name '.' shock_name ' = posterior_mean_irfs(' int2str(j) ',' int2str(i) ',:);'])
-        eval(['oo_.bvar.irf.Median.' variable_name '.' shock_name ' = posterior_median_irfs(' int2str(j) ',' int2str(i) ',:);'])
-        eval(['oo_.bvar.irf.Var.' variable_name '.' shock_name ' = posterior_variance_irfs(' int2str(j) ',' int2str(i) ',:);'])
-        eval(['oo_.bvar.irf.Upper_bound.' variable_name '.' shock_name ' = posterior_up_conf_irfs(' int2str(j) ',' int2str(i) ',:);'])
-        eval(['oo_.bvar.irf.Lower_bound.' variable_name '.' shock_name ' = posterior_down_conf_irfs(' int2str(j) ',' int2str(i) ',:);'])
+        oo_.bvar.irf.Mean.(variable_name).(shock_name) = posterior_mean_irfs(j,i,:);
+        oo_.bvar.irf.Median.(variable_name).(shock_name) = posterior_median_irfs(j,i,:);
+        oo_.bvar.irf.Var.(variable_name).(shock_name) = posterior_variance_irfs(j,i,:);
+        oo_.bvar.irf.Upper_bound.(variable_name).(shock_name) = posterior_up_conf_irfs(j,i,:);
+        oo_.bvar.irf.Lower_bound.(variable_name).(shock_name) = posterior_down_conf_irfs(j,i,:);
     end
 end
\ No newline at end of file
diff --git a/matlab/bvar_toolbox.m b/matlab/+bvar/toolbox.m
similarity index 100%
rename from matlab/bvar_toolbox.m
rename to matlab/+bvar/toolbox.m
diff --git a/preprocessor b/preprocessor
index 978789d02..6fde1b059 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 978789d02a4a27c479094512e791576414c54a73
+Subproject commit 6fde1b05981ba3351ae8662121111bb89bcf6882