diff --git a/src/auxiliary_particle_filter.m b/src/auxiliary_particle_filter.m
index 6a4e89feb15fbc2fcdd67a3ab3ba509d1585748d..228aeb3faa4743ad6f961ef9229f321cc2528e71 100644
--- a/src/auxiliary_particle_filter.m
+++ b/src/auxiliary_particle_filter.m
@@ -42,6 +42,7 @@ number_of_particles = ParticleOptions.number_of_particles;
 
 if ReducedForm.use_k_order_solver
     dr = ReducedForm.dr;
+    [udr] = folded_to_unfolded_dr(dr, DynareModel, DynareOptions);
 else
     % Set local state space model (first order approximation).
     ghx  = ReducedForm.ghx;
@@ -96,7 +97,7 @@ for t=1:sample_size
         if ReducedForm.use_k_order_solver
             tmp = 0;
             for i=1:size(nodes)
-                tmp = tmp + nodes_weights(i)*local_state_space_iteration_k(yhat, nodes(i,:)'*ones(1,number_of_particles), dr, Model, DynareOptions);
+                tmp = tmp + nodes_weights(i)*local_state_space_iteration_k(yhat, nodes(i,:)'*ones(1,number_of_particles), dr, Model, DynareOptions, udr);
             end
         else
             tmp = 0;
@@ -121,7 +122,7 @@ for t=1:sample_size
         StateVectors_ = tmp_(mf0,:);
     else
         if ReducedForm.use_k_order_solver
-            tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+            tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
         else
             tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
         end
diff --git a/src/conditional_filter_proposal.m b/src/conditional_filter_proposal.m
index 3b6b65538881e884280168033aab3405bde7408d..b816d139e1290d040e449fc655da1ac1b71598cf 100644
--- a/src/conditional_filter_proposal.m
+++ b/src/conditional_filter_proposal.m
@@ -79,7 +79,7 @@ epsilon = Q_lower_triangular_cholesky*nodes';
 yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2));
 
 if ReducedForm.use_k_order_solver
-    tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+    tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions,folded_to_unfolded_dr(dr, Model, DynareOptions));
 else
     tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
 end
diff --git a/src/gaussian_filter_bank.m b/src/gaussian_filter_bank.m
index d7361207bc4d2f9bca5d81028f656ee83a2cd7b8..0731ecaa5459344b02dad2d503703b09c95b2e90 100644
--- a/src/gaussian_filter_bank.m
+++ b/src/gaussian_filter_bank.m
@@ -81,7 +81,7 @@ StateVectors = sigma_points(1:number_of_state_variables,:);
 epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
 yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
 if ReducedForm.use_k_order_solver
-    tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+    tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, folded_to_unfolded_dr(dr, Model, DynareOptions));
 else
     tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
 end
diff --git a/src/gaussian_mixture_filter_bank.m b/src/gaussian_mixture_filter_bank.m
index 74b9dea477585c22f32fb2ff0e4c1cbf37a16cad..963b3eb2f7e85636f9dd4bb1269bf0c531946855 100644
--- a/src/gaussian_mixture_filter_bank.m
+++ b/src/gaussian_mixture_filter_bank.m
@@ -77,7 +77,7 @@ epsilon = bsxfun(@plus, StructuralShocksSqrtP*nodes3(:,number_of_state_variables
 StateVectors = bsxfun(@plus, StateSqrtP*nodes3(:,1:number_of_state_variables)', StateMu);
 yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
 if ReducedForm.use_k_order_solver
-    tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+    tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions,folded_to_unfolded_dr(dr, Model, DynareOptions));
 else
     tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
 end
diff --git a/src/measurement_equations.m b/src/measurement_equations.m
index 87edc75fc129b184287b71ec7932b3893db58dca..45cee63455038a796b12beb7f2b401a2cf884a01 100644
--- a/src/measurement_equations.m
+++ b/src/measurement_equations.m
@@ -32,7 +32,7 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state;
 number_of_structural_innovations = length(ReducedForm.Q);
 yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
 if ReducedForm.use_k_order_solver
-    tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions);
+    tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, folded_to_unfolded_dr(dr, Model, DynareOptions));
     measure = tmp(mf1,:);
 else
     measure = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
diff --git a/src/nonlinear_kalman_filter.m b/src/nonlinear_kalman_filter.m
index 322eb6e6ec00d2ce40143dc21331c481327a832a..594c751afea59596e6d98ab7a4b73be0ce740d07 100644
--- a/src/nonlinear_kalman_filter.m
+++ b/src/nonlinear_kalman_filter.m
@@ -56,6 +56,7 @@ end
 
 if ReducedForm.use_k_order_solver
     dr = ReducedForm.dr;
+    [udr] = folded_to_unfolded_dr(dr, Model, DynareOptions);
 else
     % Set local state space model (first-order approximation).
     ghx  = ReducedForm.ghx;
@@ -116,7 +117,7 @@ for t=1:sample_size
     epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
     yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
     if ReducedForm.use_k_order_solver
-        tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+        tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions,udr);
     else
         tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
     end
diff --git a/src/online_auxiliary_filter.m b/src/online_auxiliary_filter.m
index 91c51f53027817de20b318b480c823d59bb23a2e..a7b79fa2a83218ea610713ade4f77b3a021509c6 100644
--- a/src/online_auxiliary_filter.m
+++ b/src/online_auxiliary_filter.m
@@ -125,6 +125,7 @@ for t=1:sample_size
             % Set local state space model (second-order approximation).
             if ReducedForm.use_k_order_solver
                 dr = ReducedForm.dr;
+                [udr] = folded_to_unfolded_dr(dr, Model, DynareOptions);
             else
                 constant = ReducedForm.constant;
                 % Set local state space model (first-order approximation).
@@ -138,7 +139,7 @@ for t=1:sample_size
             % particle likelihood contribution
             yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state);
             if ReducedForm.use_k_order_solver
-                tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions);
+                tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions, udr);
             else
                 if pruning
                     yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
@@ -194,7 +195,7 @@ for t=1:sample_size
                     % compute particles likelihood contribution
                     yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state);
                     if ReducedForm.use_k_order_solver
-                        tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+                        tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
                     else
                         if pruning
                             yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state);
diff --git a/src/sequential_importance_particle_filter.m b/src/sequential_importance_particle_filter.m
index 975ac12996ba1a3db4cd7fe861cc7ffdbc24310f..365bbc287ef96d06a15d5cdcfa1f50707c9c5445 100644
--- a/src/sequential_importance_particle_filter.m
+++ b/src/sequential_importance_particle_filter.m
@@ -51,6 +51,7 @@ end
 
 if ReducedForm.use_k_order_solver
     dr = ReducedForm.dr;
+    [udr] = folded_to_unfolded_dr(dr, Model, DynareOptions);
 else
     % Set local state space model (first order approximation).
     ghx  = ReducedForm.ghx;
@@ -106,7 +107,7 @@ for t=1:sample_size
         [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
     else
         if ReducedForm.use_k_order_solver
-            tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+            tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
         else
             tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
         end