diff --git a/matlab/@dprior/dprior.m b/matlab/@dprior/dprior.m
new file mode 100644
index 0000000000000000000000000000000000000000..94c38e9e8b587774d97b729122aa2ba64f0968c6
--- /dev/null
+++ b/matlab/@dprior/dprior.m
@@ -0,0 +1,376 @@
+classdef dprior
+
+    properties
+        p6 = [];                         % Prior first hyperparameter.
+        p7 = [];                         % Prior second hyperparameter.
+        p3 = [];                         % Lower bound of the prior support.
+        p4 = [];                         % Upper bound of the prior support.
+        lb = [];                         % Truncated prior lower bound.
+        ub = [];                         % Truncated prior upper bound.
+        uniform_index = [];              % Index for the uniform priors.
+        gaussian_index = [];             % Index for the gaussian priors.
+        gamma_index = [];                % Index for the gamma priors.
+        beta_index = [];                 % Index for the beta priors.
+        inverse_gamma_1_index = [];      % Index for the inverse gamma type 1 priors.
+        inverse_gamma_2_index = [];      % Index for the inverse gamma type 2 priors.
+        weibull_index = [];              % Index for the weibull priors.
+        uniform_draws = false;
+        gaussian_draws = false;
+        gamma_draws = false;
+        beta_draws = false;
+        inverse_gamma_1_draws = false;
+        inverse_gamma_2_draws = false;
+        weibull_draws = false;
+    end
+
+    methods
+
+        function o = dprior(BayesInfo, PriorTrunc, Uniform)
+        % Class constructor.
+        %
+        % INPUTS
+        % - BayesInfo    [struct]   Informations about the prior distribution, aka bayestopt_.
+        % - PriorTrunc   [double]   scalar, probability mass to be excluded, aka options_.prior_trunc
+        % - Uniform      [logical]  scalar, produce uniform random deviates on the prior support.
+        %
+        % OUTPUTS
+        % - o            [dprior]   scalar, prior object.
+        %
+        % REQUIREMENTS
+        % None.
+            o.p6 = BayesInfo.p6;
+            o.p7 = BayesInfo.p7;
+            o.p3 = BayesInfo.p3;
+            o.p4 = BayesInfo.p4;
+            bounds = prior_bounds(BayesInfo, PriorTrunc);
+            o.lb = bounds.lb;
+            o.ub = bounds.ub;
+            if nargin>2 && Uniform
+                prior_shape = repmat(5, length(o.p6), 1);
+            else
+                prior_shape = BayesInfo.pshape;
+            end
+            o.beta_index = find(prior_shape==1);
+            if ~isempty(o.beta_index)
+                o.beta_draws = true;
+            end
+            o.gamma_index = find(prior_shape==2);
+            if ~isempty(o.gamma_index)
+                o.gamma_draws = true;
+            end
+            o.gaussian_index = find(prior_shape==3);
+            if ~isempty(o.gaussian_index)
+                o.gaussian_draws = true;
+            end
+            o.inverse_gamma_1_index = find(prior_shape==4);
+            if ~isempty(o.inverse_gamma_1_index)
+                o.inverse_gamma_1_draws = true;
+            end
+            o.uniform_index = find(prior_shape==5);
+            if ~isempty(o.uniform_index)
+                o.uniform_draws = true;
+            end
+            o.inverse_gamma_2_index = find(prior_shape==6);
+            if ~isempty(o.inverse_gamma_2_index)
+                o.inverse_gamma_2_draws = true;
+            end
+            o.weibull_index = find(prior_shape==8);
+            if ~isempty(o.weibull_index)
+                o.weibull_draws = true;
+            end
+        end
+
+        function p = draw(o)
+        % Return a random draw from the prior distribution.
+        %
+        % INPUTS
+        % - o    [dprior]
+        %
+        % OUTPUTS
+        % - p    [double]   m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
+        %
+        % REMARKS
+        % None.
+        %
+        % EXAMPLE
+        %
+        % >> Prior = dprior(bayestopt_, options_.prior_trunc);
+        % >> d = Prior.draw()
+            p = NaN(rows(o.lb), 1);
+            if o.uniform_draws
+                p(o.uniform_index) = rand(length(o.uniform_index),1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
+                out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
+                while ~isempty(out_of_bound)
+                    p(o.uniform_index) = rand(length(o.uniform_index), 1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
+                    out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
+                end
+            end
+            if o.gaussian_draws
+                p(o.gaussian_index) = randn(length(o.gaussian_index), 1).*o.p7(o.gaussian_index) + o.p6(o.gaussian_index);
+                out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
+                while ~isempty(out_of_bound)
+                    p(o.gaussian_index(out_of_bound)) = randn(length(o.gaussian_index(out_of_bound)), 1).*o.p7(o.gaussian_index(out_of_bound)) + o.p6(o.gaussian_index(out_of_bound));
+                    out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
+                end
+            end
+            if o.gamma_draws
+                p(o.gamma_index) = gamrnd(o.p6(o.gamma_index), o.p7(o.gamma_index))+o.p3(o.gamma_index);
+                out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
+                while ~isempty(out_of_bound)
+                    p(o.gamma_index(out_of_bound)) = gamrnd(o.p6(o.gamma_index(out_of_bound)), o.p7(o.gamma_index(out_of_bound)))+o.p3(o.gamma_index(out_of_bound));
+                    out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
+                end
+            end
+            if o.beta_draws
+                p(o.beta_index) = (o.p4(o.beta_index)-o.p3(o.beta_index)).*betarnd(o.p6(o.beta_index), o.p7(o.beta_index))+o.p3(o.beta_index);
+                out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
+                while ~isempty(out_of_bound)
+                    p(o.beta_index(out_of_bound)) = (o.p4(o.beta_index(out_of_bound))-o.p3(o.beta_index(out_of_bound))).*betarnd(o.p6(o.beta_index(out_of_bound)), o.p7(o.beta_index(out_of_bound)))+o.p3(o.beta_index(out_of_bound));
+                    out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
+                end
+            end
+            if o.inverse_gamma_1_draws
+                p(o.inverse_gamma_1_index) = ...
+                    sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index)/2, 2./o.p6(o.inverse_gamma_1_index)))+o.p3(o.inverse_gamma_1_index);
+                out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
+                while ~isempty(out_of_bound)
+                    p(o.inverse_gamma_1_index(out_of_bound)) = ...
+                        sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_1_index(out_of_bound))))+o.p3(o.inverse_gamma_1_index(out_of_bound));
+                    out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
+                end
+            end
+            if o.inverse_gamma_2_draws
+                p(o.inverse_gamma_2_index) = ...
+                    1./gamrnd(o.p7(o.inverse_gamma_2_index)/2, 2./o.p6(o.inverse_gamma_2_index))+o.p3(o.inverse_gamma_2_index);
+                out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
+                while ~isempty(out_of_bound)
+                    p(o.inverse_gamma_2_index(out_of_bound)) = ...
+                        1./gamrnd(o.p7(o.inverse_gamma_2_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_2_index(out_of_bound)))+o.p3(o.inverse_gamma_2_index(out_of_bound));
+                    out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
+                end
+            end
+            if o.weibull_draws
+                p(o.weibull_index) = wblrnd(o.p7(o.weibull_index), o.p6(o.weibull_index)) + o.p3(o.weibull_index);
+                out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
+                while ~isempty(out_of_bound)
+                    p(o.weibull_index(out_of_bound)) = wblrnd(o.p7(o.weibull_index(out_of_bound)), o.p6(o.weibull_index(out_of_bound)))+o.p3(o.weibull_index(out_of_bound));
+                    out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
+                end
+            end
+        end
+
+        function P = draws(o, n)
+        % Return n independent random draws from the prior distribution.
+        %
+        % INPUTS
+        % - o    [dprior]
+        %
+        % OUTPUTS
+        % - P    [double]   m×n matrix, random draw from the prior distribution.
+        %
+        % REMARKS
+        % If the Parallel Computing Toolbox is available, the main loop is run in parallel.
+        %
+        % EXAMPLE
+        %
+        % >> Prior = dprior(bayestopt_, options_.prior_trunc);
+        % >> Prior.draws(1e6)
+            P = NaN(rows(o.lb), 1);
+            parfor i=1:n
+                P(:,i) = draw(o);
+            end
+        end
+
+    end % methods
+end % classdef --*-- Unit tests --*--
+
+%@test:1
+%$ % Fill global structures with required fields...
+%$ prior_trunc = 1e-10;
+%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1);    % Prior shape
+%$ p1 = .4*ones(14,1);                          % Prior mean
+%$ p2 = .2*ones(14,1);                          % Prior std.
+%$ p3 = NaN(14,1);
+%$ p4 = NaN(14,1);
+%$ p5 = NaN(14,1);
+%$ p6 = NaN(14,1);
+%$ p7 = NaN(14,1);
+%$
+%$ for i=1:14
+%$    switch p0(i)
+%$      case 1
+%$        % Beta distribution
+%$        p3(i) = 0;
+%$        p4(i) = 1;
+%$        [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
+%$      case 2
+%$        % Gamma distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
+%$      case 3
+%$        % Normal distribution
+%$        p3(i) = -Inf;
+%$        p4(i) = Inf;
+%$        p6(i) = p1(i);
+%$        p7(i) = p2(i);
+%$        p5(i) = p1(i);
+%$      case 4
+%$        % Inverse Gamma (type I) distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
+%$      case 5
+%$        % Uniform distribution
+%$        [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
+%$        p3(i) = p6(i);
+%$        p4(i) = p7(i);
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
+%$      case 6
+%$        % Inverse Gamma (type II) distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
+%$      case 8
+%$        % Weibull distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
+%$      otherwise
+%$        error('This density is not implemented!')
+%$    end
+%$ end
+%$
+%$ BayesInfo.pshape = p0;
+%$ BayesInfo.p1 = p1;
+%$ BayesInfo.p2 = p2;
+%$ BayesInfo.p3 = p3;
+%$ BayesInfo.p4 = p4;
+%$ BayesInfo.p5 = p5;
+%$ BayesInfo.p6 = p6;
+%$ BayesInfo.p7 = p7;
+%$
+%$ ndraws = 1e5;
+%$ m0 = BayesInfo.p1; %zeros(14,1);
+%$ v0 = diag(BayesInfo.p2.^2); %zeros(14);
+%$
+%$ % Call the tested routine
+%$ try
+%$    % Instantiate dprior object
+%$    o = dprior(BayesInfo, prior_trunc, false);
+%$    % Do simulations in a loop and estimate recursively the mean and the variance.
+%$    for i = 1:ndraws
+%$         draw = o.draw();
+%$         m1 = m0 + (draw-m0)/i;
+%$         m2 = m1*m1';
+%$         v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
+%$         m0 = m1;
+%$    end
+%$    t(1) = true;
+%$ catch
+%$     t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$     t(2) = all(abs(m0-BayesInfo.p1)<3e-3);
+%$     t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3));
+%$ end
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ % Fill global structures with required fields...
+%$ prior_trunc = 1e-10;
+%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1);    % Prior shape
+%$ p1 = .4*ones(14,1);                          % Prior mean
+%$ p2 = .2*ones(14,1);                          % Prior std.
+%$ p3 = NaN(14,1);
+%$ p4 = NaN(14,1);
+%$ p5 = NaN(14,1);
+%$ p6 = NaN(14,1);
+%$ p7 = NaN(14,1);
+%$
+%$ for i=1:14
+%$    switch p0(i)
+%$      case 1
+%$        % Beta distribution
+%$        p3(i) = 0;
+%$        p4(i) = 1;
+%$        [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
+%$      case 2
+%$        % Gamma distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
+%$      case 3
+%$        % Normal distribution
+%$        p3(i) = -Inf;
+%$        p4(i) = Inf;
+%$        p6(i) = p1(i);
+%$        p7(i) = p2(i);
+%$        p5(i) = p1(i);
+%$      case 4
+%$        % Inverse Gamma (type I) distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
+%$      case 5
+%$        % Uniform distribution
+%$        [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
+%$        p3(i) = p6(i);
+%$        p4(i) = p7(i);
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
+%$      case 6
+%$        % Inverse Gamma (type II) distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
+%$      case 8
+%$        % Weibull distribution
+%$        p3(i) = 0;
+%$        p4(i) = Inf;
+%$        [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
+%$        p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
+%$      otherwise
+%$        error('This density is not implemented!')
+%$    end
+%$ end
+%$
+%$ BayesInfo.pshape = p0;
+%$ BayesInfo.p1 = p1;
+%$ BayesInfo.p2 = p2;
+%$ BayesInfo.p3 = p3;
+%$ BayesInfo.p4 = p4;
+%$ BayesInfo.p5 = p5;
+%$ BayesInfo.p6 = p6;
+%$ BayesInfo.p7 = p7;
+%$
+%$ ndraws = 1e5;
+%$
+%$ % Call the tested routine
+%$ try
+%$    % Instantiate dprior object.
+%$    o = dprior(BayesInfo, prior_trunc, false);
+%$    X = o.draws(ndraws);
+%$    m = mean(X, 2);
+%$    v = var(X, 0, 2);
+%$    t(1) = true;
+%$ catch
+%$     t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$     t(2) = all(abs(m-BayesInfo.p1)<3e-3);
+%$     t(3) = all(all(abs(v-BayesInfo.p2.^2)<5e-3));
+%$ end
+%$ T = all(t);
+%@eof:2
diff --git a/matlab/GetOneDraw.m b/matlab/GetOneDraw.m
index ffd5665c90154a4b89ce53849bc8a2d4ca9c7e54..8c92a84ba70dac74984995a35d658ca067dd4c54 100644
--- a/matlab/GetOneDraw.m
+++ b/matlab/GetOneDraw.m
@@ -1,5 +1,5 @@
-function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_)
-% function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_)
+function [xparams, logpost] = GetOneDraw(distribution, M_, estim_params_, oo_, options_, bayestopt_)
+
 % draws one parameter vector and its posterior from MCMC or the prior
 %
 % INPUTS
@@ -18,7 +18,7 @@ function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,baye
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2005-2017 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -35,12 +35,13 @@ function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,baye
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-switch type
-  case 'posterior'
-    [xparams, logpost] = metropolis_draw(0);
-  case 'prior'
-    xparams = prior_draw();
+if isprior(distribution)
+    xparams = distribution.draw();
     if nargout>1
-        logpost = evaluate_posterior_kernel(xparams',M_,estim_params_,oo_,options_,bayestopt_);
+        logpost = evaluate_posterior_kernel(xparams, M_, estim_params_, oo_, options_, bayestopt_);
     end
-end
\ No newline at end of file
+elseif ischar(distribution) && strcmpi(distribution, 'posterior')
+    [xparams, logpost] = metropolis_draw(0);
+else
+    error('GetOneDraw:: Wrong inputs.')
+end
diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m
index 2118704707ca69edfa29d4fd0fd4025effc07b4e..9b5122d14b496dcc09f2d621276f67737f1fd6e7 100644
--- a/matlab/PosteriorIRF_core1.m
+++ b/matlab/PosteriorIRF_core1.m
@@ -23,7 +23,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
 % SPECIAL REQUIREMENTS.
 %   None.
 %
-% Copyright © 2006-2019 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -134,12 +134,17 @@ end
 % Parallel 'while' very good!!!
 stock_param=zeros(MAX_nruns,npar);
 stock_irf_dsge=zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
+
+if strcmp(type, 'prior')
+    Prior = dprior(bayestopt_, options_.prior_trunc);
+end
+
 while fpar<B
     fpar = fpar + 1;
     irun = irun+1;
     irun2 = irun2+1;
     if strcmpi(type,'prior')
-        deep = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
+        deep = Prior.draw();
     else
         deep = x(fpar,:);
     end
@@ -294,7 +299,7 @@ end
 % directory on call machine that contain the model).
 
 myoutput.OutputFileName = [OutputFileName_dsge;
-                    OutputFileName_param;
-                    OutputFileName_bvardsge];
+                           OutputFileName_param;
+                           OutputFileName_bvardsge];
 
 myoutput.nosaddle = nosaddle;
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index 89838b5d9ac4ac483326d524b5e2787aaa664166..0dd48f0714729f7039453dc300fb11b6084a5858 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -40,13 +40,14 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
 %    * identification_analysis
 %    * isoctave
 %    * plot_identification
-%    * prior_draw
+%    * dprior.draw
 %    * set_default_option
 %    * set_prior
 %    * skipline
 %    * vnorm
 % =========================================================================
-% Copyright © 2010-2022 Dynare Team
+
+% Copyright © 2010-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -328,18 +329,18 @@ if options_.discretionary_policy || options_.ramsey_policy
         options_ident.analytic_derivation_mode=-1;
     end
 end
-    
+
 options_.analytic_derivation_mode = options_ident.analytic_derivation_mode; %overwrite setting in options_
 
-% initialize persistent variables in prior_draw
+% Instantiate dprior object (Prior)
 if prior_exist
     if any(bayestopt_.pshape > 0)
         if options_ident.prior_range
             %sample uniformly from prior ranges (overwrite prior specification)
-            prior_draw(bayestopt_, options_.prior_trunc, true);
+            Prior = dprior(bayestopt_, options_.prior_trunc, true);
         else
             %sample from prior distributions
-            prior_draw(bayestopt_, options_.prior_trunc, false);
+            Prior = dprior(bayestopt_, options_.prior_trunc, false);
         end
     else
         options_ident.prior_mc = 1; %only one single point
@@ -495,7 +496,7 @@ if iload <=0
             kk=0;
             while kk<50 && info(1)
                 kk=kk+1;
-                params = prior_draw();
+                params = Prior.draw();
                 options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
                 % perform identification analysis
                 [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
@@ -548,7 +549,7 @@ if iload <=0
         if external_sample
             params = pdraws0(iteration+1,:); % loaded draws
         else
-            params = prior_draw(); % new random draw from prior
+            params = Prior.draw(); % new random draw from prior
         end
         options_ident.tittxt = []; % clear title text for graphs and figures
         % run identification analysis
@@ -692,7 +693,7 @@ if iload <=0
             end
 
             % store results for minimal system
-            if ~options_MC.no_identification_minimal 
+            if ~options_MC.no_identification_minimal
                 if ~error_indicator.identification_minimal
                     STO_dMINIMAL(:,:,run_index)                  = ide_minimal.dMINIMAL;
                     IDE_MINIMAL.cond(iteration,1)                = ide_minimal.cond;
@@ -901,7 +902,7 @@ if SampleSize > 1
                 store_nodisplay = options_.nodisplay;
                 options_.nodisplay = 1;
                 % HIGHEST CONDITION NUMBER
-                [~, jmax] = max(IDE_MOMENTS.cond);                
+                [~, jmax] = max(IDE_MOMENTS.cond);
                 tittxt = 'Draw with HIGHEST condition number';
                 fprintf('\nTesting %s.\n',tittxt);
                 if ~iload
diff --git a/matlab/execute_prior_posterior_function.m b/matlab/execute_prior_posterior_function.m
index 00cbe64e7dc10314adaf79987564641ef8f3c49f..3ac13f177edf5348d62b852de5f60e99349f917b 100644
--- a/matlab/execute_prior_posterior_function.m
+++ b/matlab/execute_prior_posterior_function.m
@@ -17,7 +17,7 @@ function oo_=execute_prior_posterior_function(posterior_function_name,M_,options
 % OUTPUTS
 %   oo_          [structure]     Matlab/Octave structure gathering the results (initialized by dynare, see @ref{oo_}).
 
-% Copyright © 2013-2015 Dynare Team
+% Copyright © 2013-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -47,21 +47,21 @@ end
 %Create function handle
 functionhandle=str2func(posterior_function_name);
 
-prior = true;
 n_draws=options_.sampling_draws;
-% Get informations about the _posterior_draws files.
+
 if strcmpi(type,'posterior')
-    %% discard first mh_drop percent of the draws:
+    % Get informations about the _posterior_draws files.
+    % discard first mh_drop percent of the draws:
     CutSample(M_, options_, estim_params_);
-    %% initialize metropolis draws
-    options_.sub_draws=n_draws; %set draws for sampling; changed value is not returned to base workspace
-    [error_flag,~,options_]= metropolis_draw(1,options_,estim_params_,M_);
+    % initialize metropolis draws
+    options_.sub_draws = n_draws; % set draws for sampling; changed value is not returned to base workspace
+    [error_flag, ~, options_] = metropolis_draw(1, options_, estim_params_, M_);
     if error_flag
         error('EXECUTE_POSTERIOR_FUNCTION: The draws could not be initialized')
     end
-    n_draws=options_.sub_draws;
-    prior = false;
+    n_draws = options_.sub_draws;
 elseif strcmpi(type,'prior')
+    % Get informations about the prior distribution.
     if isempty(bayestopt_)
         if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
             [xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
@@ -72,38 +72,36 @@ elseif strcmpi(type,'prior')
     if exist([M_.fname '_prior_restrictions.m'])
         warning('prior_function currently does not support endogenous prior restrictions. They will be ignored. Consider using a prior_function with nobs=1.')
     end
-    prior_draw(bayestopt_, options_.prior_trunc);
+    Prior = dprior(bayestopt_, options_.prior_trunc);
 else
     error('EXECUTE_POSTERIOR_FUNCTION: Unknown type!')
 end
 
-%get draws for later use
-first_draw=GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
-parameter_mat=NaN(n_draws,length(first_draw));
-parameter_mat(1,:)=first_draw;
-for draw_iter=2:n_draws
-    parameter_mat(draw_iter,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
+if strcmpi(type, 'prior')
+    parameter_mat = Prior.draws(n_draws);
+else
+    parameter_mat = NaN(length(bayestopt_.p6), n_draws);
+    for i = 1:n_draws
+        parameter_mat(:,i) = GetOneDraw(type, M_, estim_params_, oo_, options_, bayestopt_);
+    end
 end
 
-% get output size
+% Get output size
 try
-    junk=functionhandle(parameter_mat(1,:),M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info);
+    junk = functionhandle(parameter_mat(:,1), M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info);
 catch err
     fprintf('\nEXECUTE_POSTERIOR_FUNCTION: Execution of prior/posterior function led to an error. Execution cancelled.\n')
     rethrow(err)
 end
 
-%initialize cell with number of columns
-results_cell=cell(n_draws,size(junk,2));
+% Initialize cell with number of columns
+results_cell = cell(n_draws, columns(junk));
 
-%% compute function on draws
-for draw_iter = 1:n_draws
-    M_ = set_all_parameters(parameter_mat(draw_iter,:),estim_params_,M_);
-    [results_cell(draw_iter,:)]=functionhandle(parameter_mat(draw_iter,:),M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info);
+% Evaluate function on each draw
+for i = 1:n_draws
+    M_ = set_all_parameters(parameter_mat(:,i), estim_params_, M_);
+    [results_cell(i,:)] = functionhandle(parameter_mat(:,i), M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info);
 end
 
-if prior
-    oo_.prior_function_results = results_cell;
-else
-    oo_.posterior_function_results = results_cell;
-end
+% Save results under oo_
+oo_.(sprintf('%s_function_results', type)) = results_cell;
diff --git a/matlab/isprior.m b/matlab/isprior.m
new file mode 100644
index 0000000000000000000000000000000000000000..a7a014b14b88b09057f67d91d50180f1f199662e
--- /dev/null
+++ b/matlab/isprior.m
@@ -0,0 +1,22 @@
+function b = isprior(o)
+
+% Return True iff o is a dprior object.
+
+% Copyright © 2023 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+b = isa(o, 'dprior');
diff --git a/matlab/particles b/matlab/particles
index 7926e75aa525b73c91bbb944d19556c43a8fbbe6..bb16df7f2e6ddbc81e23b9661ee100f5e8dc675a 160000
--- a/matlab/particles
+++ b/matlab/particles
@@ -1 +1 @@
-Subproject commit 7926e75aa525b73c91bbb944d19556c43a8fbbe6
+Subproject commit bb16df7f2e6ddbc81e23b9661ee100f5e8dc675a
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
index 7925f482b63098b3aaae4a30b40aa90476a2d938..0722858c4cdfa9a91d5c829a0fb322f1f7b315ea 100644
--- a/matlab/posterior_sampler_core.m
+++ b/matlab/posterior_sampler_core.m
@@ -36,7 +36,7 @@ function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatl
 % See the comments in the posterior_sampler.m funtion.
 
 
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -86,8 +86,6 @@ oo_ = myinputs.oo_;
 if whoiam
     % initialize persistent variables in priordens()
     priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
-    % initialize persistent variables in prior_draw()
-    prior_draw(bayestopt_,options_.prior_trunc);
 end
 
 MetropolisFolder = CheckPath('metropolis',M_.dname);
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 52ba42060acc6b7543f37ed79dd630e557d09f66..60a909d2c6d6993707dcf33bd667aa21c01a2f26 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -38,7 +38,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -113,7 +113,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
     fprintf(fidlog,['  Number of simulations per block: ' int2str(nruns(1)) '\n']);
     fprintf(fidlog,' \n');
     if isempty(d)
-        prior_draw(bayestopt_,options_.prior_trunc);
+        Prior = dprior(bayestopt_, options_.prior_trunc);
     end
     if options_.mh_initialize_from_previous_mcmc.status
         PreviousFolder0 = options_.mh_initialize_from_previous_mcmc.directory;
@@ -132,7 +132,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
             MetropolisFolder0 = fileparts(RecordFile0);
             PreviousFolder0=fileparts(MetropolisFolder0);
             [~, ModelName0]=fileparts(PreviousFolder0);
-        else            
+        else
             %% check for proper filesep char in user defined paths
             PreviousFolder0=strrep(PreviousFolder0,'\',filesep);
             MetropolisFolder0 = [PreviousFolder0 filesep 'metropolis'];
@@ -195,7 +195,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
             trial = 1;
             while validate == 0 && trial <= 10
                 if isempty(d)
-                    candidate = prior_draw();
+                    candidate = Prior.draw();
                 else
                     if isfield(options_,'mh_init_scale')
                         if trial==1
@@ -228,7 +228,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
                     if options_.nointeractive
                         disp(['Estimation::mcmc: I reduce mh_init_scale_factor by 10 percent:'])
                         if isfield(options_,'mh_init_scale')
-                           options_.mh_init_scale = .9*options_.mh_init_scale; 
+                           options_.mh_init_scale = .9*options_.mh_init_scale;
                            fprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.\n',options_.mh_init_scale)
                         else
                             options_.mh_init_scale_factor = .9*options_.mh_init_scale_factor;
@@ -301,7 +301,6 @@ if ~options_.load_mh_file && ~options_.mh_recover
     if options_.mh_initialize_from_previous_mcmc.status
         record.InitialSeeds = record0.LastSeeds;
     else
-        
         for j=1:NumberOfBlocks
             % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
             set_dynare_seed(options_.DynareRandomStreams.seed+j);
@@ -502,7 +501,7 @@ elseif options_.mh_recover
     % How many mh-files are saved in this block?
     ExistingDrawsInLastMCFile=zeros(NumberOfBlocks,1); %initialize: no MCMC draws of current MCMC are in file from last run
     % Check whether last present file is a file included in the last MCMC run
-    
+
     update_record=0;
     for k=1:NumberOfBlocks
         FirstBlock = k;
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index 7d04d28450e4df46fbf949220aa6eeaf5116c63f..64cf7685a7c77bb63390733fbdcffe44cfc577a9 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -30,7 +30,7 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
 % SPECIAL REQUIREMENTS.
 %   None.
 
-% Copyright © 2005-2020 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -109,8 +109,8 @@ if ~strcmpi(type,'prior')
         logpost=myinputs.logpost;
     end
 else
-    x=(prior_draw(bayestopt_,options_.prior_trunc))';
-    options_=select_qz_criterium_value(options_);
+    Prior = dprior(bayestopt_,options_.prior_trunc);
+    options_ = select_qz_criterium_value(options_);
 end
 if whoiam
     Parallel=myinputs.Parallel;
@@ -180,9 +180,9 @@ end
 if naK
     stock_filter_step_ahead =NaN(length(options_.filter_step_ahead),endo_nbr,gend+max(options_.filter_step_ahead),MAX_naK);
 end
-stock_param = NaN(MAX_nruns,size(x,2));
-stock_logpo = NaN(MAX_nruns,1);
-stock_ys = NaN(MAX_nruns,endo_nbr);
+stock_param = NaN(MAX_nruns, length(bayestopt_.p6));
+stock_logpo = NaN(MAX_nruns, 1);
+stock_ys = NaN(MAX_nruns, endo_nbr);
 if filter_covariance
     stock_filter_covariance = zeros(endo_nbr,endo_nbr,gend+1,MAX_filter_covariance);
 end
@@ -196,7 +196,8 @@ for b=fpar:B
         iter=1;
         logpo=[];
         while (isempty(logpo) || isinf(logpo)) && iter<1000
-            [deep, logpo] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
+            deep = Prior.draw();
+            logpo = evaluate_posterior_kernel(deep, M_, estim_params_, oo_, options_, bayestopt_)
             iter=iter+1;
         end
         if iter==1000
@@ -214,7 +215,7 @@ for b=fpar:B
 
     if run_smoother
         [dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_);
-        if ismember(info(1),[3,4]) 
+        if ismember(info(1),[3,4])
             opts_local.qz_criterium = 1 + (opts_local.qz_criterium-1)*10; %loosen tolerance, useful for big models where BK conditions were tested on restricted state space
             [dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_);
         end
diff --git a/matlab/prior_sampler.m b/matlab/prior_sampler.m
index 75b1896fd30c70f17eb76214eb1e6cd0713f0a33..1303cbcf02784a45d506f4f30e83aa26b3787263 100644
--- a/matlab/prior_sampler.m
+++ b/matlab/prior_sampler.m
@@ -13,7 +13,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2009-2018 Dynare Team
+% Copyright © 2009-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,7 +31,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 % Initialization.
-prior_draw(bayestopt_, options_.prior_trunc);
+Prior = dprior(bayestopt_, options_.prior_trunc);
 PriorDirectoryName = CheckPath('prior/draws',M_.dname);
 work = ~drsave;
 iteration = 0;
@@ -94,10 +94,10 @@ while iteration < NumberOfSimulations
         dyn_waitbar(iteration/NumberOfSimulations,hh,'Please wait. Prior sampler...');
     end
     loop_indx = loop_indx+1;
-    params = prior_draw();
-    M_ = set_all_parameters(params,estim_params_,M_);
-    [T,R,~,INFO,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
-    dr=oo_.dr;
+    params = Prior.draw();
+    M_ = set_all_parameters(params, estim_params_, M_);
+    [T, R, ~, INFO, M_, oo_] = dynare_resolve(M_, options_, oo_, 'restrict');
+    dr = oo_.dr;
     if ~INFO(1)
         INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_);
     end
@@ -200,4 +200,4 @@ m1 = m0 + (newobs'-m0)/iter;
 qq = m1*m1';
 s1 = s0 + ( (newobs'*newobs-qq-s0) + (iter-1)*(m0*m0'-qq') )/iter;
 mu = m1;
-sigma = s1;
\ No newline at end of file
+sigma = s1;
diff --git a/matlab/slice_sampler.m b/matlab/slice_sampler.m
index a50c12ab6c4944fd329379ba1b613e21d541d010..73fc6c3479a641dde7a861e86407a1269ec9506e 100644
--- a/matlab/slice_sampler.m
+++ b/matlab/slice_sampler.m
@@ -24,7 +24,7 @@ function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetapri
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -58,6 +58,8 @@ neval = zeros(npar,1);
 % % % fname0=fname;
 fname = [ int2str(sampler_options.curr_block)];
 
+Prior = dprior(varargin{6},varargin{3}.prior_trunc);
+
 it=0;
 while it<npar
     it=it+1;
@@ -80,8 +82,8 @@ while it<npar
         icount=0;
         while ~isfinite(fxold) && icount<1000
             icount=icount+1;
-            theta = prior_draw();
-            if all(theta(:) >= thetaprior(:,1)) && all(theta(:) <= thetaprior(:,2))
+            theta = Prior.draw();
+            if all(theta >= thetaprior(:,1)) && all(theta <= thetaprior(:,2))
                 fxold = -feval(objective_function,theta,varargin{:});
             end
         end