Make SIS (particle filter) compatible with linear models.

matlab/estimation/dynare_estimation_1.m

@@ -12,7 +12,7 @@ function dynare_estimation_1(var_list_,dname)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2003-2024 Dynare Team
+% Copyright © 2003-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -68,35 +68,31 @@ else
 end
 
 % Set particle filter flag.
-if options_.order > 1
-    if options_.particle.status
-        skipline()
-        disp('Estimation using a non linear filter!')
-        skipline()
-        if ~options_.nointeractive && ismember(options_.mode_compute,[1,3,4]) && ~strcmpi(options_.particle.filter_algorithm,'gf')% Known gradient-based optimizers
-            disp('You are using a gradient-based mode-finder. Particle filtering introduces discontinuities in the')
-            disp('objective function w.r.t the parameters. Thus, should use a non-gradient based optimizer.')
-            fprintf('\nPlease choose a mode-finder:\n')
-            fprintf('\t 0 - Continue using gradient-based method (it is most likely that you will no get any sensible result).\n')
-            fprintf('\t 6 - Monte Carlo based algorithm\n')
-            fprintf('\t 7 - Nelder-Mead simplex based optimization routine (Matlab optimization toolbox required)\n')
-            fprintf('\t 8 - Nelder-Mead simplex based optimization routine (Dynare''s implementation)\n')
-            fprintf('\t 9 - CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm\n')
-            choice = [];
-            while isempty(choice)
-                choice = input('Please enter your choice: ');
-                if isnumeric(choice) && isint(choice) && ismember(choice,[0 6 7 8 9])
-                    if choice
-                        options_.mode_compute = choice;
-                    end
-                else
-                    fprintf('\nThis is an invalid choice (you have to choose between 0, 6, 7, 8 and 9).\n')
-                    choice = [];
+if options_.particle.status
+    skipline()
+    disp('Estimation using a non linear filter!')
+    skipline()
+    if ~options_.nointeractive && ismember(options_.mode_compute,[1,3,4]) && ~strcmpi(options_.particle.filter_algorithm,'gf')% Known gradient-based optimizers
+        disp('You are using a gradient-based mode-finder. Particle filtering introduces discontinuities in the')
+        disp('objective function w.r.t the parameters. Thus, should use a non-gradient based optimizer.')
+        fprintf('\nPlease choose a mode-finder:\n')
+        fprintf('\t 0 - Continue using gradient-based method (it is most likely that you will no get any sensible result).\n')
+        fprintf('\t 6 - Monte Carlo based algorithm\n')
+        fprintf('\t 7 - Nelder-Mead simplex based optimization routine (Matlab optimization toolbox required)\n')
+        fprintf('\t 8 - Nelder-Mead simplex based optimization routine (Dynare''s implementation)\n')
+        fprintf('\t 9 - CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm\n')
+        choice = [];
+        while isempty(choice)
+            choice = input('Please enter your choice: ');
+            if isnumeric(choice) && isint(choice) && ismember(choice,[0 6 7 8 9])
+                if choice
+                    options_.mode_compute = choice;
                 end
+            else
+                fprintf('\nThis is an invalid choice (you have to choose between 0, 6, 7, 8 and 9).\n')
+                choice = [];
             end
         end
-    else
-        error('For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=true;')
     end
 end
 
@@ -579,14 +575,14 @@ end
 
 %Run and store classical smoother if needed
 if options_.smoother && ... %Bayesian smoother requested before
-    (any(bayestopt_.pshape > 0) && options_.mh_replic || ... % Bayesian with MCMC run
-    any(bayestopt_.pshape > 0) && options_.load_mh_file) % Bayesian with loaded MCMC
-    % nothing to do
+        (any(bayestopt_.pshape > 0) && options_.mh_replic || ... % Bayesian with MCMC run
+         any(bayestopt_.pshape > 0) && options_.load_mh_file) % Bayesian with loaded MCMC
+                                                              % nothing to do
 elseif options_.partial_information ||...
-    options_.order>1 %no particle smoother
-    % smoothing not yet supported
+        options_.order>1 %no particle smoother
+                         % smoothing not yet supported
 else
-    %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
+    % ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
     oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
 end
 

matlab/estimation/non_linear_dsge_likelihood.m

@@ -120,7 +120,10 @@ state_variables_idx = restrict_variables_idx(mf0);
 number_of_state_variables = length(mf0);
 
 ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
-ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
+ReducedForm.constant = ReducedForm.steadystate;
+if options_.order==2
+    ReducedForm.constant = ReducedForm.constant + .5*dr.ghs2(restrict_variables_idx);
+end
 ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
 ReducedForm.Q = Q;
 ReducedForm.H = H;
@@ -138,17 +141,19 @@ else
     ReducedForm.use_k_order_solver = false;
     ReducedForm.ghx  = dr.ghx(restrict_variables_idx,:);
     ReducedForm.ghu  = dr.ghu(restrict_variables_idx,:);
-    ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
-    ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
-    ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
-    ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
-    if options_.order==3
-        ReducedForm.ghxxx = dr.ghxxx(restrict_variables_idx,:);
-        ReducedForm.ghuuu = dr.ghuuu(restrict_variables_idx,:);
-        ReducedForm.ghxxu = dr.ghxxu(restrict_variables_idx,:);
-        ReducedForm.ghxuu = dr.ghxuu(restrict_variables_idx,:);
-        ReducedForm.ghxss = dr.ghxss(restrict_variables_idx,:);
-        ReducedForm.ghuss = dr.ghuss(restrict_variables_idx,:);
+    if options_.order>1
+        ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
+        ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
+        ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
+        ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
+        if options_.order==3
+            ReducedForm.ghxxx = dr.ghxxx(restrict_variables_idx,:);
+            ReducedForm.ghuuu = dr.ghuuu(restrict_variables_idx,:);
+            ReducedForm.ghxxu = dr.ghxxu(restrict_variables_idx,:);
+            ReducedForm.ghxuu = dr.ghxuu(restrict_variables_idx,:);
+            ReducedForm.ghxss = dr.ghxss(restrict_variables_idx,:);
+            ReducedForm.ghuss = dr.ghuss(restrict_variables_idx,:);
+        end
     end
 end
 

matlab/nonlinear-filters/sequential_importance_particle_filter.m

@@ -2,7 +2,7 @@ function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,P
 
 % Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
 
-% Copyright © 2011-2022 Dynare Team
+% Copyright © 2011-2025 Dynare Team
 %
 % This file is part of Dynare (particles module).
 %
@@ -47,20 +47,21 @@ else
     % 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
-        % Set local state space model (third order approximation).
-        ghxxx = ReducedForm.ghxxx;
-        ghuuu = ReducedForm.ghuuu;
-        ghxxu = ReducedForm.ghxxu;
-        ghxuu = ReducedForm.ghxuu;
-        ghxss = ReducedForm.ghxss;
-        ghuss = ReducedForm.ghuss;
+    if order>1
+        ghxx = ReducedForm.ghxx;
+        ghuu = ReducedForm.ghuu;
+        ghxu = ReducedForm.ghxu;
+        ghs2 = ReducedForm.ghs2;
+        if order == 3
+            % Set local state space model (third order approximation).
+            ghxxx = ReducedForm.ghxxx;
+            ghuuu = ReducedForm.ghuuu;
+            ghxxu = ReducedForm.ghxxu;
+            ghxuu = ReducedForm.ghxuu;
+            ghxss = ReducedForm.ghxss;
+            ghuss = ReducedForm.ghuss;
+        end
     end
 end
 
@@ -126,20 +127,22 @@ for t=1:sample_size
         if ReducedForm.use_k_order_solver
             tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
         else
-            if order == 2
+            switch order
+              case 1
+                tmp = NaN(rows(ghx), columns(yhat));
+                parfor i=1:columns(yhat)
+                    tmp(:,i) = constant + ghx*yhat(:,i) + ghu*epsilon(:,i)
+                end
+              case 2
                 tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
-            elseif order == 3
+              case 3
                 tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
-            else
+              otherwise
                 error('Order > 3: use_k_order_solver should be set to true');
             end
         end
     end
-    %PredictedObservedMean = tmp(mf1,:)*transpose(weights);
     PredictionError = bsxfun(@minus,Y(:,t),tmp(ReducedForm.mf1,:));
-    %dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
-    %PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
-    %PredictedObservedVariance = H;
     if rcond(H) > 1e-16
         lnw = -.5*(const_lik+sum(PredictionError.*(H\PredictionError),1));
     else
@@ -150,13 +153,13 @@ for t=1:sample_size
     wtilde = weights.*exp(lnw-dfac);
     lik(t) = log(sum(wtilde))+dfac;
     weights = wtilde/sum(wtilde);
-    if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
+    if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*number_of_particles) || ParticleOptions.resampling.status.systematic
         if pruning
-            temp = resample([tmp(ReducedForm.mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
+            temp = resample([tmp(ReducedForm.mf0,:)' tmp_(mf0_,:)'], weights', ParticleOptions);
             StateVectors = temp(:,1:number_of_state_variables)';
             StateVectors_ = temp(:,number_of_state_variables+1:end)';
         else
-            StateVectors = resample(tmp(ReducedForm.mf0,:)',weights',ParticleOptions)';
+            StateVectors = resample(tmp(ReducedForm.mf0,:)', weights', ParticleOptions)';
         end
         weights = ones(1,number_of_particles)/number_of_particles;
     elseif ParticleOptions.resampling.status.none