matlab/ep/setup_stochastic_perfect_foresight_model_solver.m

@@ -5,7 +5,7 @@ function pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_)
 %  o options_               [struct]    Dynare's options structure
 %  o oo_                    [struct]    Dynare's results structure
 
-% Copyright © 2013-2023 Dynare Team
+% Copyright © 2013-2025 Dynare Team
 %
 % This file is part of Dynare.
 %

matlab/ep/solve_stochastic_perfect_foresight_model.m

-function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model(endo_simul,exo_simul,pfm,nnodes,order)
-
-% Copyright © 2012-2017 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/>.
-
-flag = 0;
-err = 0;
-stop = 0;
-
-params = pfm.params;
-steady_state = pfm.steady_state;
-ny = pfm.ny;
-periods = pfm.periods;
-dynamic_model = pfm.dynamic_model;
-lead_lag_incidence = pfm.lead_lag_incidence;
-lead_lag_incidence_t = transpose(lead_lag_incidence);
-nyp = pfm.nyp;
-nyf = pfm.nyf;
-i_cols_1 = pfm.i_cols_1;
-i_cols_j = pfm.i_cols_j;
-i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
-
-maxit = pfm.maxit_;
-tolerance = pfm.tolerance;
-verbose = pfm.verbose;
-
-number_of_shocks = size(exo_simul,2);
-
-[nodes,weights] = gauss_hermite_weights_and_nodes(nnodes);
-
-if number_of_shocks>1
-    nodes = repmat(nodes,1,number_of_shocks)*chol(pfm.Sigma);
-    % to be fixed for Sigma ~= I
-    for i=number_of_shocks:-1:1
-        rr(i) = {nodes(:,i)};
-        ww(i) = {weights};
-    end
-    nodes = cartesian_product_of_sets(rr{:});
-    weights = prod(cartesian_product_of_sets(ww{:}),2);
-    nnodes = nnodes^number_of_shocks;
-else
-    nodes = nodes*sqrt(pfm.Sigma);
-end
-
-if verbose
-    disp (' -----------------------------------------------------');
-    disp ('MODEL SIMULATION :');
-    fprintf('\n');
-end
-
-z = endo_simul(lead_lag_incidence_t(:)>0);
-[~, jacobian] = dynamic_model(z, exo_simul, params,steady_state, 2);
-
-% Each column of Y represents a different world
-% The upper right cells are unused
-% The first row block is ny x 1
-% The second row block is ny x nnodes
-% The third row block is ny x nnodes^2
-% and so on until size ny x nnodes^order
-world_nbr = nnodes^order;
-Y = repmat(endo_simul(:),1,world_nbr);
-
-% The columns of A map the elements of Y such that
-% each block of Y with ny rows are unfolded column wise
-dimension = ny*(sum(nnodes.^(0:order-1),2)+(periods-order)*world_nbr);
-if order == 0
-    i_upd_r = (1:ny*periods);
-    i_upd_y = i_upd_r + ny;
-else
-    i_upd_r = zeros(dimension,1);
-    i_upd_y = i_upd_r;
-    i_upd_r(1:ny) = (1:ny);
-    i_upd_y(1:ny) = ny+(1:ny);
-    i1 = ny+1;
-    i2 = 2*ny;
-    n1 = ny+1;
-    n2 = 2*ny;
-    for i=2:periods
-        for j=1:nnodes^min(i-1,order)
-            i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods;
-            i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2);
-            i1 = i2+1;
-            i2 = i2+ny;
-        end
-        n1 = n2+1;
-        n2 = n2+ny;
-    end
-end
-if rows(lead_lag_incidence)>2
-    icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
-           find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
-else
-    if nyf
-        icA = [find(lead_lag_incidence(2,:))+world_nbr*ny find(lead_lag_incidence(3,:))+2*world_nbr*ny ]';
-    else
-        icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ]';
-    end
-end
-h1 = clock;
-for iter = 1:maxit
-    A1 = sparse([],[],[],ny*(sum(nnodes.^(0:order-1),2)+1),dimension,(order+1)*world_nbr*nnz(jacobian));
-    res = zeros(ny,periods,world_nbr);
-    i_rows = 1:ny;
-    i_cols = find(lead_lag_incidence');
-    i_cols_p = i_cols(1:nyp);
-    i_cols_s = i_cols(nyp+(1:ny));
-    i_cols_f = i_cols(nyp+ny+(1:nyf));
-    i_cols_Ap = i_cols_p;
-    i_cols_As = i_cols_s;
-    i_cols_Af = i_cols_f - ny;
-    for i = 1:order+1
-        i_w_p = 1;
-        for j = 1:nnodes^(i-1)
-            innovation = exo_simul;
-            if i > 1
-                innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
-            end
-            if i <= order
-                for k=1:nnodes
-                    y = [Y(i_cols_p,i_w_p);
-                         Y(i_cols_s,j);
-                         Y(i_cols_f,(j-1)*nnodes+k)];
-                    [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
-                    if i == 1
-                        % in first period we don't keep track of
-                        % predetermined variables
-                        i_cols_A = [i_cols_As - ny; i_cols_Af];
-                        A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
-                    else
-                        i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
-                        A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
-                    end
-                    res(:,i,j) = res(:,i,j)+weights(k)*d1;
-                    i_cols_Af = i_cols_Af + ny;
-                end
-            else
-                y = [Y(i_cols_p,i_w_p);
-                     Y(i_cols_s,j);
-                     Y(i_cols_f,j)];
-                [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
-                if i == 1
-                    % in first period we don't keep track of
-                    % predetermined variables
-                    i_cols_A = [i_cols_As - ny; i_cols_Af];
-                    A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
-                else
-                    i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
-                    A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
-                end
-                res(:,i,j) = d1;
-                i_cols_Af = i_cols_Af + ny;
-            end
-            i_rows = i_rows + ny;
-            if mod(j,nnodes) == 0
-                i_w_p = i_w_p + 1;
-            end
-            if i > 1
-                if mod(j,nnodes) == 0
-                    i_cols_Ap = i_cols_Ap + ny;
-                end
-                i_cols_As = i_cols_As + ny;
-            end
-        end
-        i_cols_p = i_cols_p + ny;
-        i_cols_s = i_cols_s + ny;
-        i_cols_f = i_cols_f + ny;
-    end
-    nzA = cell(periods,world_nbr);
-    for j=1:world_nbr
-        i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
-        offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
-        offset_r = (j-1)*ny;
-        for i=order+2:periods
-            [d1,jacobian] = dynamic_model(Y(i_rows_y,j), ...
-                                          exo_simul,params, ...
-                                          steady_state,i+1);
-            if i == periods
-                [ir,ic,v] = find(jacobian(:,i_cols_T));
-            else
-                [ir,ic,v] = find(jacobian(:,i_cols_j));
-            end
-            nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]';
-            res(:,i,j) = d1;
-            i_rows_y = i_rows_y + ny;
-            offset_c = offset_c + world_nbr*ny;
-            offset_r = offset_r + world_nbr*ny;
-        end
-    end
-    err = max(abs(res(i_upd_r)));
-    if err < tolerance
-        stop = 1;
-        if verbose
-            fprintf('\n') ;
-            disp([' Total time of simulation        :' num2str(etime(clock,h1))]) ;
-            fprintf('\n') ;
-            disp(' Convergency obtained.') ;
-            fprintf('\n') ;
-        end
-        flag = 0;% Convergency obtained.
-        endo_simul = reshape(Y(:,1),ny,periods+2);
-        break
-    end
-    A2 = [nzA{:}]';
-    A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)];
-    dy = -A\res(i_upd_r);
-    Y(i_upd_y) =   Y(i_upd_y) + dy;
-end
-
-if ~stop
-    if verbose
-        fprintf('\n') ;
-        disp(['     Total time of simulation        :' num2str(etime(clock,h1))]) ;
-        fprintf('\n') ;
-        disp('WARNING : maximum number of iterations is reached (modify options_.simul.maxit).') ;
-        fprintf('\n') ;
-    end
-    flag = 1;% more iterations are needed.
-    endo_simul = 1;
-end
-if verbose
-    disp ('-----------------------------------------------------') ;
-end

matlab/ep/solve_stochastic_perfect_foresight_model_0.m

+function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm)
+
+% Copyright © 2025 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/>.
+
+update_pfm_struct = false;
+update_options_struct = false;
+
+if nargout>4
+    update_pfm_struct = true;
+end
+
+if nargout>5
+    update_options_struct = true;
+end
+
+dynamic_model = pfm.dynamic_model;
+
+periods = pfm.periods;
+
+order = pfm.stochastic_order;
+
+lead_lag_incidence = pfm.lead_lag_incidence;
+lead_lag_incidence_t = transpose(lead_lag_incidence);
+ny = pfm.ny;
+nyp = pfm.nyp;
+nyf = pfm.nyf;
+i_cols_1 = pfm.i_cols_1;
+i_cols_j = pfm.i_cols_j;
+i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
+
+nodes = pfm.nodes;
+weights = pfm.weights;
+nnodes = pfm.nnodes;
+
+if update_pfm_struct
+
+    % make sure that there is a node equal to zero
+    % and permute nodes and weights to have zero first
+    k = find(sum(abs(nodes),2) < 1e-12);
+    if ~isempty(k)
+        nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
+        weights = [weights(k); weights(1:k-1); weights(k+1:end)];
+    else
+        error('there is no nodes equal to zero')
+    end
+
+    pfm.nodes = nodes;
+    pfm.weights = weights;
+
+    if pfm.hybrid_order>0
+        if pfm.hybrid_order==2
+            pfm.h_correction = 0.5*pfm.dr.ghs2(dr.inv_order_var);
+        elseif pfm.hybrid_order>2
+            pfm.h_correction = pfm.dr.g_0(pfm.dr.inv_order_var);
+        else
+            pfm.h_correction = 0;
+        end
+    else
+        pfm.h_correction = 0;
+    end
+
+    z = endo_simul(lead_lag_incidence_t(:)>0);
+    [~, jacobian] = dynamic_model(z, exo_simul, pfm.params, pfm.steady_state, 2);
+
+    world_nbr = nnodes^order;
+
+    % The columns of A map the elements of Y such that
+    % each block of Y with ny rows are unfolded column wise
+    dimension = ny*(sum(nnodes.^(0:order-1),2)+(periods-order)*world_nbr);
+
+    i_upd_r = zeros(dimension,1);
+    i_upd_y = i_upd_r;
+    i_upd_r(1:ny) = (1:ny);
+    i_upd_y(1:ny) = ny+(1:ny);
+    i1 = ny+1;
+    i2 = 2*ny;
+    n1 = ny+1;
+    n2 = 2*ny;
+    for i=2:periods
+        for j=1:nnodes^min(i-1,order)
+            i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods;
+            i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2);
+            i1 = i2+1;
+            i2 = i2+ny;
+        end
+        n1 = n2+1;
+        n2 = n2+ny;
+    end
+
+    if rows(lead_lag_incidence)>2
+        icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
+               find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
+    else
+        if nyf
+            icA = [find(lead_lag_incidence(2,:))+world_nbr*ny find(lead_lag_incidence(3,:))+2*world_nbr*ny ]';
+        else
+            icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ]';
+        end
+    end
+
+    i_rows = 1:ny;
+    i_cols = find(lead_lag_incidence');
+    i_cols_p = i_cols(1:nyp);
+    i_cols_s = i_cols(nyp+(1:ny));
+    i_cols_f = i_cols(nyp+ny+(1:nyf));
+
+    pfm.i_cols_Ap = i_cols_p;
+    pfm.i_cols_As = i_cols_s;
+    pfm.i_cols_Af = i_cols_f - ny;
+    pfm.i_hc = i_cols_f - 2*ny;
+    pfm.i_cols_p = i_cols_p;
+    pfm.i_cols_s = i_cols_s;
+    pfm.i_cols_f = i_cols_f;
+    pfm.i_rows = i_rows;
+
+    pfm.A1 = sparse([],[],[],ny*(sum(nnodes.^(0:order-1),2)+1),dimension,(order+1)*world_nbr*nnz(jacobian));
+    pfm.res = zeros(ny,periods,world_nbr);
+
+    pfm.order = order;
+    pfm.world_nbr = world_nbr;
+
+    pfm.i_cols_1 = i_cols_1;
+    pfm.i_cols_h = i_cols_j;
+    pfm.icA = icA;
+    pfm.i_cols_T = i_cols_T;
+    pfm.i_upd_r = i_upd_r;
+    pfm.i_upd_y = i_upd_y;
+
+    pfm.dimension = dimension;
+
+end
+
+pfm.Y = repmat(endo_simul(:),1,pfm.world_nbr);
+
+if isempty(y)
+    y = repmat(pfm.steady_state, pfm.dimension/pfm.ny, 1);
+end
+
+if update_options_struct
+    % Set algorithm
+    options_.solve_algo = options_.ep.solve_algo;
+    options_.simul.maxit = options_.ep.maxit;
+    [lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), pfm.params);
+    pfm.eq_index = M_.dynamic_mcp_equations_reordering;
+    if options_.ep.solve_algo == 10
+        options_.lmmcp.lb = repmat(lb, pfm.dimension/pfm.ny, 1);
+        options_.lmmcp.ub = repmat(ub, pfm.dimension/pfm.ny, 1);
+    elseif options_.ep.solve_algo == 11
+        options_.mcppath.lb = repmat(lb, pfm.dimension/pfm.ny, 1);
+        options_.mcppath.ub = repmat(ub, pfm.dimension/pfm.ny, 1);
+    end
+end
+
+[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_0, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
+
+endo_simul(:,2) = y(1:ny);

matlab/ep/solve_stochastic_perfect_foresight_model_1.m

-function [errorflag, endo_simul, errorcode, y] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, Options, pfm, order, varargin)
+function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm)
 
-% Copyright © 2012-2022 Dynare Team
+% Copyright © 2012-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -17,87 +17,74 @@ function [errorflag, endo_simul, errorcode, y] = solve_stochastic_perfect_foresi
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if nargin < 6
-    homotopy_parameter = 1;
-else
-    homotopy_parameter = varargin{1};
+update_pfm_struct = false;
+update_options_struct = false;
+
+if nargout>4
+    update_pfm_struct = true;
 end
 
-flag = 0;
-err = 0;
-
-EpOptions = Options.ep;
-
-params = pfm.params;
-steady_state = pfm.steady_state;
-ny = pfm.ny;
-periods = pfm.periods;
-dynamic_model = pfm.dynamic_model;
-lead_lag_incidence = pfm.lead_lag_incidence;
-nyp = pfm.nyp;
-nyf = pfm.nyf;
-i_cols_1 = pfm.i_cols_1;
-i_cols_A1 = pfm.i_cols_A1;
-i_cols_j = pfm.i_cols_j;
-i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
-hybrid_order = pfm.hybrid_order;
-dr = pfm.dr;
-nodes = pfm.nodes;
-weights = pfm.weights;
-nnodes = pfm.nnodes;
-
-maxit = pfm.maxit_;
-tolerance = pfm.tolerance;
-verbose = pfm.verbose;
-
-number_of_shocks = size(exo_simul,2);
-
-% make sure that there is a node equal to zero
-% and permute nodes and weights to have zero first
-k = find(sum(abs(nodes),2) < 1e-12);
-if ~isempty(k)
-    nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
-    weights = [weights(k); weights(1:k-1); weights(k+1:end)];
-else
-    error('there is no nodes equal to zero')
+if nargout>5
+    update_options_struct = true;
 end
-if hybrid_order > 0
-    if hybrid_order == 2
-        h_correction = 0.5*dr.ghs2(dr.inv_order_var);
+
+
+if update_pfm_struct
+
+    periods = pfm.periods;
+
+    order = pfm.stochastic_order;
+
+    ny = pfm.ny;
+    lead_lag_incidence = pfm.lead_lag_incidence;
+    i_cols_1 = pfm.i_cols_1;
+    i_cols_j = pfm.i_cols_j;
+    i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
+
+    nodes = pfm.nodes;
+    weights = pfm.weights;
+    nnodes = pfm.nnodes;
+
+    % Make sure that there is a node equal to zero and permute nodes and weights to have zero first
+    k = find(sum(abs(nodes),2) < 1e-12);
+    if ~isempty(k)
+        nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
+        weights = [weights(k); weights(1:k-1); weights(k+1:end)];
+    else
+        error('there is no nodes equal to zero')
     end
-else
-    h_correction = 0;
-end
 
-if verbose
-    disp (' -----------------------------------------------------');
-    disp ('MODEL SIMULATION :');
-    fprintf('\n');
-end
+    pfm.nodes = nodes;
+    pfm.weights = weights;
 
-% Each column of Y represents a different world
-% The upper right cells are unused
-% The first row block is ny x 1
-% The second row block is ny x nnodes
-% The third row block is ny x nnodes^2
-% and so on until size ny x nnodes^order
-world_nbr = pfm.world_nbr;
-Y = endo_simul(:,2:end-1);
-Y = repmat(Y,1,world_nbr);
-pfm.y0 = endo_simul(:,1);
+    if pfm.hybrid_order > 0
+        if pfm.hybrid_order == 2
+            pfm.h_correction = 0.5*pfm.dr.ghs2(pfm.dr.inv_order_var);
+        elseif pfm.hybrid_order>2
+            pfm.h_correction = pfm.dr.g_0(pfm.dr.inv_order_var);
+        else
+            pfm.h_correction = 0;
+        end
+    else
+        pfm.h_correction = 0;
+    end
+
+    % Each column of Y represents a different world
+    % The upper right cells are unused
+    % The first row block is ny x 1
+    % The second row block is ny x nnodes
+    % The third row block is ny x nnodes^2
+    % and so on until size ny x nnodes^order
+    world_nbr = pfm.world_nbr;
 
-% The columns of A map the elements of Y such that
-% each block of Y with ny rows are unfolded column wise
-% number of blocks
-block_nbr = pfm.block_nbr;
-% dimension of the problem
-dimension = ny*block_nbr;
-pfm.dimension = dimension;
-if order == 0
-    i_upd_r = (1:ny*periods)';
-    i_upd_y = i_upd_r + ny;
-else
-    i_upd_r = zeros(dimension,1);
+    % The columns of A map the elements of Y such that
+    % each block of Y with ny rows are unfolded column wise
+    % number of blocks
+    block_nbr = pfm.block_nbr;
+    % dimension of the problem
+    dimension = ny*block_nbr;
+    pfm.dimension = dimension;
+    i_upd_r = zeros(dimension, 1);
     i_upd_y = i_upd_r;
     i_upd_r(1:ny) = (1:ny);
     i_upd_y(1:ny) = ny+(1:ny);
@@ -116,32 +103,40 @@ else
         n1 = n2+1;
         n2 = n2+ny;
     end
+    icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
+           find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
+
+    pfm.i_rows = 1:ny;
+    pfm.i_cols = find(lead_lag_incidence');
+    pfm.i_cols_1 = i_cols_1;
+    pfm.i_cols_h = i_cols_j;
+    pfm.icA = icA;
+    pfm.i_cols_T = i_cols_T;
+    pfm.i_upd_r = i_upd_r;
+    pfm.i_upd_y = i_upd_y;
+
+end
+
+if isempty(y)
+    y = repmat(pfm.steady_state, pfm.block_nbr, 1);
+end
+
+if update_options_struct
+    % Set algorithm
+    options_.solve_algo = options_.ep.solve_algo;
+    options_.simul.maxit = options_.ep.maxit;
+    [lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), pfm.params);
+    pfm.eq_index = M_.dynamic_mcp_equations_reordering;
+    if options_.ep.solve_algo == 10
+        options_.lmmcp.lb = repmat(lb, pfm.block_nbr, 1);
+        options_.lmmcp.ub = repmat(ub, pfm.block_nbr, 1);
+    elseif options_.ep.solve_algo == 11
+        options_.mcppath.lb = repmat(lb, pfm.block_nbr, 1);
+        options_.mcppath.ub = repmat(ub, pfm.block_nbr, 1);
+    end
 end
-icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
-       find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
-h1 = clock;
-pfm.order = order;
-pfm.world_nbr = world_nbr;
-pfm.nodes = nodes;
-pfm.nnodes = nnodes;
-pfm.weights = weights;
-pfm.h_correction = h_correction;
-pfm.i_rows = 1:ny;
-i_cols = find(lead_lag_incidence');
-pfm.i_cols = i_cols;
-pfm.nyp = nyp;
-pfm.nyf = nyf;
-pfm.hybrid_order = hybrid_order;
-pfm.i_cols_1 = i_cols_1;
-pfm.i_cols_h = i_cols_j;
-pfm.icA = icA;
-pfm.i_cols_T = i_cols_T;
-pfm.i_upd_r = i_upd_r;
-pfm.i_upd_y = i_upd_y;
-
-Options.steady.maxit = 100;
-y = repmat(steady_state,block_nbr,1);
-Options.solve_algo = Options.ep.solve_algo;
-Options.steady.maxit = Options.ep.maxit;
-[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_2, y, Options.simul.maxit, Options.dynatol.f, Options.dynatol.x, Options, exo_simul, pfm);
-endo_simul(:,2) = y(1:ny);
\ No newline at end of file
+
+pfm.y0 = endo_simul(:,1);
+
+[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_1, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
+endo_simul(:,2) = y(1:pfm.ny);

matlab/estimation/dynare_estimation_1.m

@@ -237,6 +237,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
     optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
     for optim_iter = 1:length(optimizer_vec)
         current_optimizer = optimizer_vec{optim_iter};
+        if isnumeric(current_optimizer)
+            if current_optimizer==5
+                if options_.analytic_derivation
+                    old_analytic_derivation = options_.analytic_derivation;
+                    options_.analytic_derivation=-1; %force analytic outer product gradient hessian for each iteration
+                end
+            end
+        end
         [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
         fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
 
@@ -244,6 +252,9 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
             if current_optimizer==5
                 newratflag = new_rat_hess_info.newratflag;
                 new_rat_hess_info = new_rat_hess_info.new_rat_hess_info;
+                if options_.analytic_derivation
+                    options_.analytic_derivation = old_analytic_derivation;
+                end
             elseif current_optimizer==6 %save scaling factor
                 save([M_.dname filesep 'Output' filesep M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
                 options_.mh_jscale = Scale;
@@ -509,7 +520,11 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
                 if error_flag
                     error('%s: I cannot compute the posterior IRFs!',dispString)
                 end
-                oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
+                if options_.occbin.likelihood.status
+                    fprintf('%s: the bayesian_irf option is not compatible with the use of OccBin.',dispString)
+                else
+                    oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
+                end                
             end
             if options_.moments_varendo
                 if error_flag
@@ -538,7 +553,11 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
                         end
                     end
                 end
-                oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
+                if options_.occbin.likelihood.status
+                    fprintf('%s: the moments_varendo option is not compatible with the use of OccBin.',dispString)
+                else
+                    oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
+                end
             end
             if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
                 if error_flag

matlab/estimation/dynare_estimation_init.m

@@ -117,16 +117,21 @@ if isequal(options_.diffuse_filter,1) || (options_.kalman_algo>2)
     end
 end
 
+if strcmp('slice',options_.posterior_sampler_options.posterior_sampling_method)
+    if options_.prior_trunc==0
+        fprintf('\ndynare_estimation_init: slice requires bounded support. Setting options_.prior_trunc=1e-10.\n')
+        options_.prior_trunc=1e-10;
+    end
+end
+
+
 options_=select_qz_criterium_value(options_);
 
 % Set options related to filtered variables.
-if  isequal(options_.filtered_vars,0) && ~isempty(options_.filter_step_ahead)
+if  isequal(options_.filtered_vars,0) && ~isempty(options_.filter_step_ahead) %require filtered var because filter_step_ahead was set
     options_.filtered_vars = 1;
 end
-if ~isequal(options_.filtered_vars,0) && isempty(options_.filter_step_ahead)
-    options_.filter_step_ahead = 1;
-end
-if ~isequal(options_.filtered_vars,0) && isequal(options_.filter_step_ahead,0)
+if ~isequal(options_.filtered_vars,0) && (isempty(options_.filter_step_ahead) || isequal(options_.filter_step_ahead,0)) %require filter_step_ahead because filtered_vars was requested
     options_.filter_step_ahead = 1;
 end
 if ~isequal(options_.filter_step_ahead,0)
@@ -533,12 +538,17 @@ if (options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_fi
             error('IVF-filter: an observable is mapped to a zero variance shock.')
         end
     end
+    if ~isequal(M_.H,0)
+        error('IVF-filter: Measurement erros are not allowed with the inversion filter.')
+    end
 end
 
 if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
     if ~isempty(options_.nk)
-        fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead. Disabling the option.\n')
+        fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead and filtered_variables. Disabling the option.\n')
         options_.nk=[];
+        options_.filter_step_ahead=[];
+        options_.filtered_vars=0;
     end
     if options_.filter_covariance
         fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')

matlab/estimation/get_posterior_parameters.m

@@ -7,7 +7,9 @@ function xparam = get_posterior_parameters(type,M_,estim_params_,oo_,options_,fi
 %   o type              [char]     = 'mode' or 'mean'.
 %   o M_:               [structure] Dynare structure describing the model.
 %   o estim_params_:    [structure] Dynare structure describing the estimated parameters.
-%   o field_1           [char]     optional field like 'mle_'.
+%   o oo_:              [structure] Dynare results structure
+%   o options_:         [structure] Dynare options structure
+%   o field_1           [char]      optional field like 'mle_'.
 %
 % OUTPUTS
 %   o xparam     vector of estimated parameters
@@ -15,7 +17,7 @@ function xparam = get_posterior_parameters(type,M_,estim_params_,oo_,options_,fi
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 2006-2018 Dynare Team
+% Copyright © 2006-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -76,8 +78,6 @@ for i=1:ncn
     m = m+1;
 end
 
-FirstDeep = m;
-
 for i=1:np
     name1 = M_.param_names{estim_params_.param_vals(i,1)};
     xparam(m) = oo_.([field1 type]).parameters.(name1);

matlab/estimation/initial_estimation_checks.m

@@ -86,6 +86,15 @@ if options_.occbin.likelihood.status || options_.occbin.smoother.status
     if options_.fast_kalman_filter
         error('initial_estimation_checks:: Occbin is incompatible with the fast Kalman filter.')        
     end
+    if options_.bayesian_irf
+        error('initial_estimation_checks:: Occbin is incompatible with the bayesian_irf option.')                
+    end
+    if options_.moments_varendo
+        error('initial_estimation_checks:: Occbin is incompatible with the moments_varendo option.')                
+    end
+    if options_.forecast
+        error('initial_estimation_checks:: Occbin is incompatible with the forecast option.')
+    end
 end
 
 if (options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter) || (options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
@@ -161,6 +170,10 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) && options_.mh_nblck<1
     error('initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.')
 end
 
+if options_.mh_drop<0 || options_.mh_drop>=1
+    error('initial_estimation_checks:: mh_drop must be in [0,1).')    
+end
+
 % check and display warnings if steady-state solves static model (except if diffuse_filter == 1) and if steady-state changes estimated parameters
 [oo_.steady_state] = check_steady_state_changes_parameters(M_,estim_params_,oo_,options_, [options_.diffuse_filter==0 options_.diffuse_filter==0] );
 

matlab/estimation/marginal_density.m

@@ -3,10 +3,11 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
 % Computes the marginal density
 %
 % INPUTS
+%   M_               [structure]    Dynare model structure
 %   options_         [structure]    Dynare options structure
 %   estim_params_    [structure]    Dynare estimation parameter structure
-%   M_               [structure]    Dynare model structure
 %   oo_              [structure]    Dynare results structure
+%   bayestopt_       [structure]    describing the priors
 %   outputFolderName [string]       name of folder with results
 %
 % OUTPUTS
@@ -16,7 +17,7 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2005-2023 Dynare Team
+% Copyright © 2005-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -64,7 +65,7 @@ if ~isfield(oo_,'posterior_mode') || (options_.mh_replic && isequal(options_.pos
 end
 
 % save the posterior mean and the inverse of the covariance matrix
-% (usefull if the user wants to perform some computations using
+% (useful if the user wants to perform some computations using
 % the posterior mean instead of the posterior mode ==> ).
 parameter_names = bayestopt_.name;
 save([M_.dname filesep outputFolderName filesep M_.fname '_mean.mat'],'xparam1','hh','parameter_names','SIGMA');

matlab/estimation/online/online_auxiliary_filter.m

@@ -31,6 +31,7 @@ function online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_,
 
 % Set seed for randn().
 options_ = set_dynare_seed_local_options(options_,'default');
+options_.verbosity=0; %particularly suppress warning messages during k_order_pert within the loop
 pruning = options_.particle.pruning;
 variance_update = true;
 online_opt = options_.posterior_sampler_options.current_options;
@@ -41,7 +42,7 @@ SimulationFolder = CheckPath('online', M_.dname);
 bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
 
 % initialization of state particles
-[~, M_, options_, oo_, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+[~, ~, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
 
 order = options_.order;
 mf0 = ReducedForm.mf0;
@@ -78,10 +79,10 @@ xparam = zeros(number_of_parameters,number_of_particles);
 Prior = dprior(bayestopt_, options_.prior_trunc);
 for i=1:number_of_particles
     info = 12042009;
-    while info
+    while info(1)
         candidate = Prior.draw();
-        [info, M_, options_, oo_] = solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
-        if ~info
+        [info] = solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+        if ~info(1)
             xparam(:,i) = candidate(:);
         end
     end
@@ -117,8 +118,8 @@ for t=1:sample_size
     tau_tilde = zeros(1,number_of_particles);
     for i=1:number_of_particles
         % model resolution
-        [info, M_, options_, oo_, ReducedForm] = ...
-            solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+        [info, M_, ReducedForm] = ...
+            solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
         if ~info(1)
             steadystate = ReducedForm.steadystate;
             state_variables_steady_state = ReducedForm.state_variables_steady_state;
@@ -211,8 +212,8 @@ for t=1:sample_size
             candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
             if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
                 % model resolution for new parameters particles
-                [info, M_, options_, oo_, ReducedForm] = ...
-                    solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_) ;
+                [info, M_, ReducedForm] = ...
+                    solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
                 if ~info(1)
                     xparam(:,i) = candidate ;
                     steadystate = ReducedForm.steadystate;
@@ -359,12 +360,7 @@ for t=1:sample_size
         end 
     end
     save(sprintf('%s%sparameters_particles-%u.mat', SimulationFolder, filesep(), t), 'xparam');
-
-    str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)');
-    for l=1:size(xparam,1)
-        str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t));
-    end
-    disp(str)
+    dyntable(options_,'', {'Parameter'; 'Lower Bound (95%)';'Mean';'Upper Bound (95%)'}, bayestopt_.name, [lb95_xparam(:,t), mean_xparam(:,t), ub95_xparam(:,t)], cellofchararraymaxlength(bayestopt_.name)+2, 10, 6)
     disp('')
 
 end

matlab/estimation/prior_posterior_statistics.m

@@ -141,9 +141,9 @@ if options_.filter_covariance
     filter_covariance=1;
 end
 
-smoothed_state_uncertainty=0;
+smoothed_state_uncertainty=false;
 if options_.smoothed_state_uncertainty
-    smoothed_state_uncertainty=1;
+    smoothed_state_uncertainty=true;
 end
 
 % Store the variable mandatory for local/remote parallel computing.
@@ -327,6 +327,9 @@ if options_.smoother
         varlist,'UpdatedVariables',DirectoryName, ...
         '_update',dispString);
     if smoothed_state_uncertainty
+        if isfield(oo_,'Smoother') && isfield(oo_.Smoother,'State_uncertainty')
+            oo_.Smoother=rmfield(oo_.Smoother,'State_uncertainty'); %needs to be removed as classical smoother field has a different format
+        end
         oo_=pm3(M_,options_,oo_,endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',...
             varlist,M_.endo_names_tex,M_.endo_names,...
             varlist,'StateUncertainty',DirectoryName,'_state_uncert',dispString);

matlab/estimation/prior_posterior_statistics_core.m

@@ -76,7 +76,7 @@ naK=myinputs.naK;
 horizon=myinputs.horizon;
 iendo=myinputs.iendo;
 IdObs=myinputs.IdObs; %index of observables
-if horizon
+if horizon && ~options_.occbin.smoother.status
     i_last_obs=myinputs.i_last_obs;
     MAX_nforc1=myinputs.MAX_nforc1;
     MAX_nforc2=myinputs.MAX_nforc2;
@@ -168,7 +168,7 @@ if run_smoother
     stock_smoothed_constant=NaN(endo_nbr,gend,MAX_n_smoothed_constant);
     stock_smoothed_trend=NaN(endo_nbr,gend,MAX_n_smoothed_trend);
     stock_trend_coeff = zeros(endo_nbr,MAX_n_trend_coeff);
-    if horizon
+    if horizon && ~options_.occbin.smoother.status
         stock_forcst_mean= NaN(endo_nbr,horizon,MAX_nforc1);
         stock_forcst_point = NaN(endo_nbr,horizon,MAX_nforc2);
         if ~isequal(M_.H,0)
@@ -230,11 +230,28 @@ for b=fpar:B
             opts_local.occbin.simul.waitbar=0;
             opts_local.occbin.smoother.waitbar = 0;
             opts_local.occbin.smoother.linear_smoother=false; % speed-up
-            [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
-                occbin.DSGE_smoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
-            if oo_.occbin.smoother.error_flag(1)
-                message=get_error_message(oo_.occbin.smoother.error_flag,opts_local);
-                fprintf('\nprior_posterior_statistics: One of the draws failed with the error:\n%s\n',message)
+            if options_.occbin.smoother.inversion_filter
+                dataset_.data=Y';
+                [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, alphahat, etahat] = ...
+                    occbin.IVF_posterior(deep,dataset_,[],options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+                if info(1)
+                    message=get_error_message(info,opts_local);
+                    fprintf('\nprior_posterior_statistics: IVF smoother failed for one of the draws:\n%s\n',message)
+                else
+                    alphatilde = alphahat*nan;
+                    SteadyState=oo_.dr.ys;
+                    trend_coeff = zeros(length(options_.varobs_id),1);
+                    trend_addition=zeros(options_.number_of_observed_variables,gend);                    
+                end
+                %epsilonhat not available as no measurement error allowed
+            else
+                opts_local.verbosity=0;
+                [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
+                    occbin.DSGE_smoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
+                if oo_.occbin.smoother.error_flag(1)
+                    message=get_error_message(oo_.occbin.smoother.error_flag,opts_local);
+                    fprintf('\nprior_posterior_statistics: One of the draws failed with the error:\n%s\n',message)
+                end
             end
         else
             [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
@@ -311,7 +328,7 @@ for b=fpar:B
                 end
             end
         end
-        if horizon
+        if horizon && ~options_.occbin.smoother.status
             yyyy = alphahat(iendo,i_last_obs);
             yf = simulate_posterior_forecasts(yyyy,dr,horizon,false,M_.Sigma_e,1);
             if options_.prefilter
@@ -546,7 +563,7 @@ dyn_waitbar_close(h);
 
 
 function yf=simulate_posterior_forecasts(y0,dr,horizon,stochastic_indicator,Sigma_e,n)
-% function yf=forcst2(y0,horizon,dr,n)
+% function yf=simulate_posterior_forecasts(y0,horizon,dr,n)
 %
 % computes forecasts based on first order model solution, given shocks
 % drawn from the shock distribution, but not including measurement error

matlab/kalman/DsgeSmoother.m

-function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,d] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
-% [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,d] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
+function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,d,info] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
+% [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,d,info] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
 % Estimation of the smoothed variables and innovations.
 %
 % INPUTS
@@ -35,7 +35,14 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
 %                                   about the smoothed state (decision-rule order)
 %   o oo_           [structure] storing the results
 %   o bayestopt_    [structure] describing the priors
-%
+%   o alphahat0     [double]  (m*1) matrix, smoothed endogenous variables
+%                              (a_{0}) for initial period from PKF
+%   o state_uncertainty0 [double] (K,K) matrix storing the uncertainty about 
+%                                   the smoothed state for the initial
+%                                   period from the PKF
+%   o d             [integer]   number of diffuse periods
+%   o info          [1 by 4 double]   error code and penalty
+
 % Notes:
 %   m:  number of endogenous variables (M_.endo_nbr)
 %   T:  number of Time periods (options_.nobs)
@@ -74,9 +81,9 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-alphahat        = [];
-etahat  = [];
-epsilonhat      = [];
+alphahat      = [];
+etahat        = [];
+epsilonhat    = [];
 ahat          = [];
 SteadyState   = [];
 trend_coeff   = [];
@@ -86,8 +93,11 @@ R             = [];
 P             = [];
 PK            = [];
 decomp        = [];
-vobs            = length(options_.varobs);
+vobs          = length(options_.varobs);
 smpl          = size(Y,2);
+alphahat0     = [];
+state_uncertainty0 =[];
+d             = 0;
 
 if ~isempty(xparam1) %not calibrated model
     M_ = set_all_parameters(xparam1,estim_params_,M_);

matlab/model_diagnostics.m

@@ -160,6 +160,7 @@ for b=1:nb
                                     'evaluate', 'static', 'block_decomposed', ['block=' ...
                                 int2str(b)]);
         end
+        n_vars_jacob=size(jacob,2);
     else
         if options_.block
             T = NaN(M_.block_structure_stat.tmp_nbr, 1);

matlab/model_info.m

@@ -279,7 +279,6 @@ function print_line(names,var_index,lead_lag,M_)
         end
     else
         aux_index=find([M_.aux_vars(:).endo_index]==var_index);
-        aux_type=M_.aux_vars(aux_index).type;
         if lead_lag==0
             str = subst_auxvar(var_index, [], M_);
         else

matlab/nonlinear-filters/auxiliary_initialization.m

-function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
-
-% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
-
-% Copyright © 2011-2022 Dynare Team
-%
-% This file is part of Dynare (particles module).
-%
-% 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 particles module 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/>.
-
-persistent init_flag mf1 number_of_particles
-persistent number_of_observed_variables number_of_structural_innovations
-
-% Set default
-if isempty(start)
-    start = 1;
-end
-
-% Set flag for prunning
-%pruning = ParticleOptions.pruning;
-
-% Get steady state and mean.
-%steadystate = ReducedForm.steadystate;
-constant = ReducedForm.constant;
-ss = ReducedForm.ys;
-state_variables_steady_state = ReducedForm.state_variables_steady_state;
-
-% Set persistent variables.
-if isempty(init_flag)
-    mf1 = ReducedForm.mf1;
-    number_of_observed_variables = length(mf1);
-    number_of_structural_innovations = length(ReducedForm.Q);
-    number_of_particles = ParticleOptions.number_of_particles;
-    init_flag = 1;
-end
-
-order = options_.order;
-
-% Set local state space model (first order approximation).
-ghx  = ReducedForm.ghx;
-ghu  = ReducedForm.ghu;
-% Set local state space model (second order approximation).
-ghxx = ReducedForm.ghxx;
-ghuu = ReducedForm.ghuu;
-ghxu = ReducedForm.ghxu;
-ghs2 = ReducedForm.ghs2;
-if (order == 3)
-   ghxxx = ReducedForm.ghxxx;
-   ghuuu = ReducedForm.ghuuu;
-   ghxxu = ReducedForm.ghxxu;
-   ghxuu = ReducedForm.ghxuu;
-   ghxss = ReducedForm.ghxss;
-   ghuss = ReducedForm.ghuss;
-end
-
-% Get covariance matrices
-Q = ReducedForm.Q;
-H = ReducedForm.H;
-if isempty(H)
-    H = 0;
-end
-
-% Get initial condition for the state vector.
-StateVectorMean = ReducedForm.StateVectorMean;
-StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
-state_variance_rank = size(StateVectorVarianceSquareRoot,2);
-%Q_lower_triangular_cholesky = chol(Q)';
-%if pruning
-%    StateVectorMean_ = StateVectorMean;
-%    StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
-%end
-
-% Set seed for randn().
-options_=set_dynare_seed_local_options(options_,'default');
-
-% Initialization of the likelihood.
-const_lik = log(2*pi)*number_of_observed_variables;
-
-% Initialization of the weights across particles.
-weights = ones(1,number_of_particles)/number_of_particles ;
-StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
-%if pruning
-%    StateVectors_ = StateVectors;
-%end
-yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
-%if pruning
-%    yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
-%    [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
-%else
-if (order == 2)
-   tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
-elseif (order == 3)
-   tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ss, options_.threads.local_state_space_iteration_3, false);
-else
-   error('Orders > 3 not allowed');
-end
-%end
-PredictedObservedMean = weights*(tmp(mf1,:)');
-PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
-dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
-PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + H;
-wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
-tau_tilde = weights.*wtilde ;
-tau_tilde = tau_tilde/sum(tau_tilde);
-initial_distribution = resample(StateVectors',tau_tilde',ParticleOptions)' ;