diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 95d932248fe8cb231a21026171de81146c0fac47..8de89f5f7e6c0b4c00bbae1bcf58fedbe061add9 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -659,7 +659,6 @@ options_.particleswarm = particleswarm;
 
 % prior analysis
 options_.prior_mc = 20000;
-options_.prior_analysis_endo_var_list = {};
 
 % did model undergo block decomposition + minimum feedback set computation ?
 options_.block = false;
diff --git a/matlab/estimation/check_prior_analysis_data.m b/matlab/estimation/check_prior_analysis_data.m
index 569f324d22e6cf840889c56b7dc328c8b0b86130..ee400509870971f196387a15b016014c31d056f5 100644
--- a/matlab/estimation/check_prior_analysis_data.m
+++ b/matlab/estimation/check_prior_analysis_data.m
@@ -47,7 +47,6 @@ if ~exist([ M_.dname '/prior/draws'],'dir')
 end
 
 prior_draws_info = dir([ M_.dname '/prior/draws/prior_draws*.mat']);
-date_of_the_last_prior_draw_file = prior_draws_info(end).datenum;
 
 %% Get informations about _posterior_draws files.
 if isempty(prior_draws_info)
@@ -57,6 +56,7 @@ if isempty(prior_draws_info)
     end
     return
 else
+    date_of_the_last_prior_draw_file = prior_draws_info(end).datenum;
     date_of_the_prior_definition = get_date_of_a_file([ M_.dname '/prior/definition.mat']);
     if date_of_the_prior_definition>date_of_the_last_prior_draw_file
         info = 2;
diff --git a/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m
index c0ca04a89f01b2f2c5c8f359224856e5b6084334..ef55dad23f8fddd4f01bfc64b92c1f2cf066ef1c 100644
--- a/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m
+++ b/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m
@@ -61,19 +61,8 @@ else
 end
 
 % Set varlist (vartan)
-if ~posterior
-    if isfield(options_,'varlist')
-        temp = options_.varlist;
-    end
-    options_.varlist = options_.prior_analysis_endo_var_list;
-end
 endo_names=options_.varlist;
 [ivar,vartan ] = get_variables_list(options_, M_);
-if ~posterior
-    if exist('temp','var')
-        options_.varlist = temp;
-    end
-end
 nvar = length(ivar);
 
 % Set the size of the auto-correlation function to zero.
@@ -123,6 +112,9 @@ for file = 1:NumberOfDrawsFiles
         temp=load([M_.dname filesep folder_name filesep M_.fname '_' type '_draws' num2str(file) ]);
     else
         temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        if columns(temp.pdraws)==3 && isstruct(temp.pdraws{1,3})
+            temp.pdraws(:,2)=[];
+        end
     end
     isdrsaved = columns(temp.pdraws)-1;
     NumberOfDraws = rows(temp.pdraws);
diff --git a/matlab/estimation/dsge_simulated_theoretical_covariance.m b/matlab/estimation/dsge_simulated_theoretical_covariance.m
index cb976ff6844320ba345b7c706eaea28d838ce7da..267f0d0e4d1ffebbc095736c703ea3aeb35e2265 100644
--- a/matlab/estimation/dsge_simulated_theoretical_covariance.m
+++ b/matlab/estimation/dsge_simulated_theoretical_covariance.m
@@ -59,19 +59,8 @@ else
 end
 
 % Set varlist (vartan)
-if ~posterior
-    if isfield(options_,'varlist')
-        temp = options_.varlist;
-    end
-    options_.varlist = options_.prior_analysis_endo_var_list;
-end
 endo_names=options_.varlist;
 [ivar,vartan] = get_variables_list(options_,M_);
-if ~posterior
-    if exist('temp','var')
-        options_.varlist = temp;
-    end
-end
 nvar = length(ivar);
 
 if options_.pruning
@@ -117,6 +106,9 @@ for file = 1:NumberOfDrawsFiles
         temp=load([M_.dname filesep folder_name filesep M_.fname '_' type '_draws' num2str(file) ]);
     else
         temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        if columns(temp.pdraws)==3 && isstruct(temp.pdraws{1,3})
+            temp.pdraws(:,2)=[];
+        end
     end
     NumberOfDraws = rows(temp.pdraws);
     isdrsaved = columns(temp.pdraws)-1;
diff --git a/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m b/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m
index 209efad68ad577437ea68592a87ebbd1d28581ae..be017702ff941d73e7af247ff5434fe1477bdfdb 100644
--- a/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m
+++ b/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m
@@ -61,19 +61,8 @@ else
 end
 
 % Set varlist (vartan)
-if ~posterior
-    if isfield(options_,'varlist')
-        temp = options_.varlist;
-    end
-    options_.varlist = options_.prior_analysis_endo_var_list;
-end
 [ivar,vartan,options_] = get_variables_list(options_,M_);
 endo_names=options_.varlist;
-if ~posterior
-    if exist('temp','var')
-        options_.varlist = temp;
-    end
-end
 nvar = length(ivar);
 
 % Set the size of the auto-correlation function to zero.
@@ -127,6 +116,9 @@ for file = 1:NumberOfDrawsFiles
         temp=load([M_.dname filesep folder_name filesep M_.fname '_' type '_draws' num2str(file) ]);
     else
         temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        if columns(temp.pdraws)==3 && isstruct(temp.pdraws{1,3})
+            temp.pdraws(:,2)=[];
+        end
     end
     isdrsaved = columns(temp.pdraws)-1;
     NumberOfDraws = rows(temp.pdraws);
diff --git a/matlab/estimation/prior_analysis.m b/matlab/estimation/prior_analysis.m
index a247ba6f1caf908c92456eb89d258deb594b2651..e9065a4516d3e66519bd66240cdd49fba889abe4 100644
--- a/matlab/estimation/prior_analysis.m
+++ b/matlab/estimation/prior_analysis.m
@@ -19,26 +19,26 @@ function oo_ = prior_analysis(type,arg1,arg2,arg3,options_,M_,oo_,estim_params_)
 info = check_prior_analysis_data(type,M_);
 SampleSize = options_.prior_mc;
 switch info
-  case {0,1,2}
-    MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
-    drsize = size_of_the_reduced_form_model(oo_.dr);
-    if drsize*SampleSize>MaxMegaBytes
-        drsave=0;
-    else
-        drsave=1;
-    end
-    load([M_.dname '/prior/definition.mat']);
-    prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_);
-    clear('bayestopt_');
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
-  case {4,5}
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
-  case 6
-    [ivar,vartan] = get_variables_list(options_,M_);
-    nvar = length(ivar);
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
-  otherwise
-    error('prior_analysis:: Check_prior_analysis_data gave a meaningless output!')
+    case {0,1,2}
+        MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
+        drsize = size_of_the_reduced_form_model(oo_.dr);
+        if drsize*SampleSize>MaxMegaBytes
+            drsave=0;
+        else
+            drsave=1;
+        end
+        load([M_.dname '/prior/definition.mat']);
+        prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_);
+        clear('bayestopt_');
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
+    case {4,5}
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
+    case 6
+        [ivar,vartan] = get_variables_list(options_,M_);
+        nvar = length(ivar);
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
+    otherwise
+        error('prior_analysis:: Check_prior_analysis_data gave a meaningless output!')
 end
 
 
@@ -53,7 +53,7 @@ switch type
   case 'variance'
     if nargin==narg1
         [nvar,vartan] = ...
-            dsge_simulated_theoretical_covariance(SampleSize,M_,options_,oo_,'prior');
+            dsge_simulated_theoretical_covariance(SampleSize,arg3,M_,options_,oo_,'prior');
     end
     oo_ = covariance_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                  vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_,options_);
@@ -64,6 +64,17 @@ switch type
     end
     oo_ = variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                              M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+    if ~all(diag(M_.H)==0)
+        if strmatch(arg1,options_.varobs,'exact')
+            if isoctave && octave_ver_less_than('8.4') %Octave bug #60347
+                observable_name_requested_vars=intersect_stable(vartan,options_.varobs);
+            else
+                observable_name_requested_vars=intersect(vartan,options_.varobs,'stable');
+            end
+            oo_ = variance_decomposition_ME_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
+                [M_.exo_names;'ME'],arg2,observable_name_requested_vars,arg1,options_.mh_conf_sig,oo_,options_);
+        end
+    end
   case 'correlation'
     if nargin==narg1
         [nvar,vartan] = ...
@@ -79,9 +90,9 @@ switch type
     oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                                       arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
     if ~all(diag(M_.H)==0)
-        if strmatch(vartan(arg1,:),options_.varobs,'exact')
+        if strmatch(arg1,options_.varobs,'exact')
             oo_ = conditional_variance_decomposition_ME_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
-                                                              arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+                                                              arg3,[M_.exo_names;'ME'],arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
         end
     end
   otherwise
diff --git a/matlab/estimation/prior_sampler.m b/matlab/estimation/prior_sampler.m
index 172190ce716fbc363cfb2c4a0b7b27304d7f2130..3521df61f3befa36cd0d34382244e1a562a10198 100644
--- a/matlab/estimation/prior_sampler.m
+++ b/matlab/estimation/prior_sampler.m
@@ -6,6 +6,8 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 %   M_          [structure]  Model description.
 %   bayestopt_  [structure]  Prior distribution description.
 %   options_    [structure]  Global options of Dynare.
+%   oo_         [structure]     Dynare structure where the results are saved.
+%   estim_params_       [structure] characterizing parameters to be estimated
 %
 % OUTPUTS:
 %   results     [structure]  Various statistics.
@@ -13,7 +15,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2009-2023 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,14 +35,11 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 % Initialization.
 Prior = dprior(bayestopt_, options_.prior_trunc);
 PriorDirectoryName = CheckPath('prior/draws',M_.dname);
-work = ~drsave;
 iteration = 0;
 loop_indx = 0;
-file_indx = [];
 count_bk_indeterminacy = 0;
 count_bk_unstability = 0;
 count_bk_singularity = 0;
-count_static_var_def = 0;
 count_no_steadystate = 0;
 count_steadystate_file_exit = 0;
 count_dll_problem = 0;
diff --git a/matlab/moments/compute_moments_varendo.m b/matlab/moments/compute_moments_varendo.m
index 4b26fb17a6dd4c7890e2ca9d8983f780f168175b..874e51c9e47f9bcf9ce4c9d56cd2d7133f7e7e2c 100644
--- a/matlab/moments/compute_moments_varendo.m
+++ b/matlab/moments/compute_moments_varendo.m
@@ -64,10 +64,7 @@ if strcmpi(type,'posterior')
 elseif strcmpi(type,'prior')
     posterior = 0;
     if nargin==5
-        var_list_ = options_.prior_analysis_endo_var_list;
-        if isempty(var_list_)
-            options_.prior_analysis_var_list = options_.varobs;
-        end
+        var_list_ = options_.varobs;
     end
     if isfield(oo_,'PriorTheoreticalMoments')
         oo_=rmfield(oo_,'PriorTheoreticalMoments');
@@ -103,7 +100,7 @@ if posterior
 else
     for i=1:NumberOfEndogenousVariables
         for j=i:NumberOfEndogenousVariables
-            oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, [], options_, M_, oo_, estim_params_);
+            oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, NumberOfLags, options_, M_, oo_, estim_params_);
         end
     end
 end
diff --git a/tests/moments/fs2000_post_moments.mod b/tests/moments/fs2000_post_moments.mod
index f97968ca137b1663392da45507cecfa9d80c82e0..b3bc124c18f876c50dd254a559c5c59af63a4b23 100644
--- a/tests/moments/fs2000_post_moments.mod
+++ b/tests/moments/fs2000_post_moments.mod
@@ -285,6 +285,12 @@ for var_iter_1=1:nvars
     end
 end
 
+
+options_.prior_mc=100;
+options_.qz_criterium = 1+1e-6;
+oo_ = compute_moments_varendo('prior',options_,M_,oo_,estim_params_,M_.endo_names(1:M_.orig_endo_nbr));
+
+
 /*
  * The following lines were used to generate the data file. If you want to
  * generate another random data file, comment the "estimation" line and uncomment