diff --git a/matlab/+conditional_forecasts/run.m b/matlab/+conditional_forecasts/run.m
index 9737147557f31084ae9952ba9578fdc8bc84864d..1ac117bb54cd1944c53e9b7c855a67fe1c45cce1 100644
--- a/matlab/+conditional_forecasts/run.m
+++ b/matlab/+conditional_forecasts/run.m
@@ -33,7 +33,7 @@ function forecasts=run(M_,options_,oo_,bayestopt_,estim_params_,constrained_path
 % [1] Results are stored in oo_.conditional_forecast.
 % [2] Use the function conditional_forecasts.plot to plot the results.
 
-% Copyright © 2006-2023 Dynare Team
+% Copyright © 2006-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -161,7 +161,7 @@ else
         options_.qz_criterium = 1+1e-6;
     end
     graph_title='Calibration';
-    if ~isfield(oo_.dr,'kstate')
+    if ~isfield(oo_.dr,'ghx')
         error('You need to call stoch_simul before conditional_forecast')
     end
 end
@@ -308,4 +308,4 @@ for i = 1:EndoSize
     forecasts.uncond.ci.(M_.endo_names{oo_.dr.order_var(i)}) = [tmp(t1,:)' ,tmp(t2,:)' ]';
 end
 forecasts.graph.title = graph_title;
-forecasts.graph.fname = M_.fname;
\ No newline at end of file
+forecasts.graph.fname = M_.fname;
diff --git a/matlab/moments/UnivariateSpectralDensity.m b/matlab/moments/UnivariateSpectralDensity.m
index 419bf0152bd20bd1ab472d82062a80b1a2b754ba..d1963e62ac01acb58f514c0057b33ed83ec471e7 100644
--- a/matlab/moments/UnivariateSpectralDensity.m
+++ b/matlab/moments/UnivariateSpectralDensity.m
@@ -19,7 +19,7 @@ function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list)
 
 % Adapted from th_autocovariances.m.
 
-% Copyright © 2006-2023 Dynare Team
+% Copyright © 2006-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -68,29 +68,11 @@ ghx = oo_.dr.ghx;
 ghu = oo_.dr.ghu;
 nspred = M_.nspred;
 nstatic = M_.nstatic;
-kstate = oo_.dr.kstate;
 order = oo_.dr.order_var;
 iv(order) = 1:length(order);
 nx = size(ghx,2);
-ikx = nstatic+1:nstatic+nspred;
-k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
-i0 = find(k0(:,2) == M_.maximum_lag+1);
-i00 = i0;
-AS = ghx(:,i0);
-ghu1 = zeros(nx,M_.exo_nbr);
-ghu1(i0,:) = ghu(ikx,:);
-for i=M_.maximum_lag:-1:2
-    i1 = find(k0(:,2) == i);
-    n1 = size(i1,1);
-    j1 = zeros(n1,1);
-    j2 = j1;
-    for k1 = 1:n1
-        j1(k1) = find(k0(i00,1)==k0(i1(k1),1));
-        j2(k1) = find(k0(i0,1)==k0(i1(k1),1));
-    end
-    AS(:,j1) = AS(:,j1)+ghx(:,i1);
-    i0 = i1;
-end
+ikx = nstatic+(1:nspred);
+ghu1 = ghu(ikx,:);
 
 [A,B] = kalman_transition_matrix(oo_.dr,ikx',1:nx);
 [~, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
diff --git a/matlab/optimal_policy/mult_elimination.m b/matlab/optimal_policy/mult_elimination.m
index 8961c99c5f0181a2021ec78837efea9640772162..6f351e1bfcb1685cd18e3c81d9febfb001c4d02e 100644
--- a/matlab/optimal_policy/mult_elimination.m
+++ b/matlab/optimal_policy/mult_elimination.m
@@ -12,7 +12,7 @@ function dr=mult_elimination(varlist,M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2003-2018 Dynare Team
+% Copyright © 2003-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -88,12 +88,6 @@ lead_lag_incidence1 = zeros(size(lead_lag_incidence1'));
 lead_lag_incidence1(k) = 1:length(k);
 lead_lag_incidence1 = lead_lag_incidence1';
 
-kstate = zeros(0,2);
-for i=maximum_lag:-1:1
-    k = find(lead_lag_incidence(i,:));
-    kstate = [kstate; [k' repmat(i+1,length(k),1)]];
-end
-
 dr.M1 = M1;
 dr.M2 = M2;
 dr.M3 = M3;
diff --git a/matlab/partial_information/dr1_PI.m b/matlab/partial_information/dr1_PI.m
index d64651c56d7bf3638c6917a13de551f2d8ea17ee..b583d61d98190a32f9fd16d66cedd2da64e6ccd9 100644
--- a/matlab/partial_information/dr1_PI.m
+++ b/matlab/partial_information/dr1_PI.m
@@ -40,7 +40,7 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
 %   none.
 %
 
-% Copyright © 1996-2018 Dynare Team
+% Copyright © 1996-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -158,13 +158,12 @@ if options_.debug
 end
 
 dr=set_state_space(dr,M_);
-kstate = dr.kstate;
 nstatic = M_.nstatic;
 nfwrd = M_.nfwrd;
 nspred = M_.nspred;
 nboth = M_.nboth;
 order_var = dr.order_var;
-nd = size(kstate,1);
+nd = M_.nspred+M_.nsfwrd;
 nz = nnz(M_.lead_lag_incidence);
 
 sdyn = M_.endo_nbr - nstatic;
diff --git a/matlab/shock_decomposition/initial_condition_decomposition.m b/matlab/shock_decomposition/initial_condition_decomposition.m
index 1004645c1c9ba2a3bda999e53ea3265ef787eb00..f5606b05044628406bb0dfb5c58ffc847fccf544 100644
--- a/matlab/shock_decomposition/initial_condition_decomposition.m
+++ b/matlab/shock_decomposition/initial_condition_decomposition.m
@@ -21,7 +21,7 @@ function oo_ = initial_condition_decomposition(M_,oo_,options_,varlist,bayestopt
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2017-2018 Dynare Team
+% Copyright © 2017-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -124,8 +124,7 @@ if ~isfield(oo_,'initval_decomposition') || isequal(varlist,0)
 
     maximum_lag = M_.maximum_lag;
 
-    k2 = dr.kstate(find(dr.kstate(:,2) <= maximum_lag+1),[1 2]);
-    i_state = order_var(k2(:,1))+(min(i,maximum_lag)+1-k2(:,2))*M_.endo_nbr;
+    i_state = order_var(M_.nstatic+(1:M_.nspred));
     for i=1:gend
         if i > 1 && i <= maximum_lag+1
             lags = min(i-1,maximum_lag):-1:1;
diff --git a/matlab/shock_decomposition/realtime_shock_decomposition.m b/matlab/shock_decomposition/realtime_shock_decomposition.m
index 597999fdc5ccca7e07254dfc26c86366a7484fba..a60789a54eaf458765dc7b03b4b108dbb75af010 100644
--- a/matlab/shock_decomposition/realtime_shock_decomposition.m
+++ b/matlab/shock_decomposition/realtime_shock_decomposition.m
@@ -22,7 +22,7 @@ function oo_ = realtime_shock_decomposition(M_,oo_,options_,varlist,bayestopt_,e
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2009-2020 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -204,8 +204,7 @@ for j=presample+1:nobs
 
     maximum_lag = M_.maximum_lag;
 
-    k2 = dr.kstate(find(dr.kstate(:,2) <= maximum_lag+1),[1 2]);
-    i_state = order_var(k2(:,1))+(min(i,maximum_lag)+1-k2(:,2))*M_.endo_nbr;
+    i_state = order_var(M_.nstatic+(1:M_.nspred));
     for i=1:gend+forecast_
         if i > 1 && i <= maximum_lag+1
             lags = min(i-1,maximum_lag):-1:1;
diff --git a/matlab/shock_decomposition/shock_decomposition.m b/matlab/shock_decomposition/shock_decomposition.m
index 790629ee96d69ecef3ba7eb6b73bba51faa557f5..3c30d0aa42fa3b96f146c0c17eff217d2a6486f9 100644
--- a/matlab/shock_decomposition/shock_decomposition.m
+++ b/matlab/shock_decomposition/shock_decomposition.m
@@ -23,7 +23,7 @@ function [oo_,M_] = shock_decomposition(M_,oo_,options_,varlist,bayestopt_,estim
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2009-2023 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -112,8 +112,7 @@ z(:,end,:) = Smoothed_Variables_deviation_from_mean;
 
 maximum_lag = M_.maximum_lag;
 
-k2 = dr.kstate(find(dr.kstate(:,2) <= maximum_lag+1),[1 2]);
-i_state = order_var(k2(:,1))+(min(i,maximum_lag)+1-k2(:,2))*M_.endo_nbr;
+i_state = order_var(M_.nstatic+(1:M_.nspred));
 for i=1:gend
     if i > 1 && i <= maximum_lag+1
         lags = min(i-1,maximum_lag):-1:1;
diff --git a/matlab/stochastic_solver/AIM_first_order_solver.m b/matlab/stochastic_solver/AIM_first_order_solver.m
index ce149e1a6f0678c51f3424a78a4d2a39eaeeef7c..a8e48fbd7dcf1c0b9226879ffbe2df9dbd15a6e6 100644
--- a/matlab/stochastic_solver/AIM_first_order_solver.m
+++ b/matlab/stochastic_solver/AIM_first_order_solver.m
@@ -51,7 +51,7 @@ function [dr,info]=AIM_first_order_solver(jacobia,M_,dr,qz_criterium)
 %! @end deftypefn
 %@eod:
 
-% Copyright © 2001-2017 Dynare Team
+% Copyright © 2001-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -80,7 +80,7 @@ end
 A = kalman_transition_matrix(dr,M_.nstatic+(1:M_.nspred), 1:M_.nspred);
 dr.eigval = eig(A);
 disp(dr.eigval)
-nd = size(dr.kstate,1);
+nd = M_.nspred+M_.nsfwrd;
 nba = nd-sum( abs(dr.eigval) < qz_criterium );
 
 nsfwrd = M_.nsfwrd;
diff --git a/matlab/stochastic_solver/disp_dr.m b/matlab/stochastic_solver/disp_dr.m
index e349d3e946e3596cffc4fd944e75acc51b3668b8..338021e6086e38af2b5770f425ad991dc105118f 100644
--- a/matlab/stochastic_solver/disp_dr.m
+++ b/matlab/stochastic_solver/disp_dr.m
@@ -12,7 +12,7 @@ function disp_dr(M_,options_,dr,order,var_list)
 % OUTPUTS
 % none
 
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,8 +36,6 @@ end
 
 nx =size(dr.ghx,2);
 nu =size(dr.ghu,2);
-k = find(dr.kstate(:,2) <= M_.maximum_lag+1);
-klag = dr.kstate(k,[1 2]);
 k1 = dr.order_var;
 
 if isempty(var_list)
@@ -142,7 +140,7 @@ for k=1:nx
     if isfield(dr,'state_var')
         str1 = subst_auxvar(dr.state_var(k),-1,M_);
     else
-        str1 = subst_auxvar(k1(klag(k,1)),klag(k,2)-M_.maximum_lag-2,M_);
+        str1 = subst_auxvar(k1(M_.nstatic+k),-1,M_);
     end
     if options_.loglinear
         str = sprintf(label_format,['log(',str1,')']);
@@ -177,8 +175,8 @@ if order > 1
     for k = 1:nx
         for j = 1:k
             flag = 0;
-            str1 = sprintf('%s,%s',subst_auxvar(k1(klag(k,1)),klag(k,2)-M_.maximum_lag-2,M_), ...
-                           subst_auxvar(k1(klag(j,1)),klag(j,2)-M_.maximum_lag-2,M_));
+            str1 = sprintf('%s,%s',subst_auxvar(k1(M_.nstatic+k),-1,M_), ...
+                           subst_auxvar(k1(M_.nstatic+j),-1,M_));
             str = sprintf(label_format, str1);
             for i=1:nvar
                 if k == j
@@ -219,7 +217,7 @@ if order > 1
     for k = 1:nx
         for j = 1:nu
             flag = 0;
-            str1 = sprintf('%s,%s',subst_auxvar(k1(klag(k,1)),klag(k,2)-M_.maximum_lag-2,M_), M_.exo_names{j});
+            str1 = sprintf('%s,%s',subst_auxvar(k1(M_.nstatic+k),-1,M_), M_.exo_names{j});
             str = sprintf(label_format,str1);
             for i=1:nvar
                 x = dr.ghxu(ivar(i),(k-1)*nu+j);
diff --git a/matlab/stochastic_solver/k_order_pert.m b/matlab/stochastic_solver/k_order_pert.m
index 573c76e166b541b1a47e545ef357fe542b579a79..41ea2df57bcfc7034592cfddf3362d7e6b29c512 100644
--- a/matlab/stochastic_solver/k_order_pert.m
+++ b/matlab/stochastic_solver/k_order_pert.m
@@ -65,11 +65,9 @@ dr.ghx = dyn_derivs.gy;
 dr.ghu = dyn_derivs.gu;
 
 if options_.loglinear
-    k = find(dr.kstate(:,2) <= M_.maximum_endo_lag+1);
-    klag = dr.kstate(k,[1 2]);
     k1 = dr.order_var;
     dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ...
-             repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1);
+             repmat(dr.ys(k1(M_.nstatic+(1:M_.nspred)))',size(dr.ghx,1),1);
     dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
 end
 
diff --git a/matlab/stochastic_solver/resol.m b/matlab/stochastic_solver/resol.m
index 66916534abadb16c07a04b62761db4431ca3f9af..6364aad7845a82b07d8af074ff655aaee8545c38 100644
--- a/matlab/stochastic_solver/resol.m
+++ b/matlab/stochastic_solver/resol.m
@@ -34,7 +34,7 @@ function [dr, info, params] = resol(check_flag, M_, options_, dr_in, endo_steady
 %   info(1)=24    ->    M_.params has been updated in the steadystate routine and has some NaNs.
 %   info(1)=30    ->    Ergodic variance can't be computed.
 
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -52,9 +52,6 @@ function [dr, info, params] = resol(check_flag, M_, options_, dr_in, endo_steady
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 %preserve only the following fields:
-if isfield(dr_in,'kstate')
-    dr.kstate = dr_in.kstate;
-end
 if isfield(dr_in,'inv_order_var')
     dr.inv_order_var = dr_in.inv_order_var;
 end
diff --git a/matlab/stochastic_solver/set_state_space.m b/matlab/stochastic_solver/set_state_space.m
index 5e5de1da4ad5f14ed773b304b98154376f15bb5b..fbb65c1e408d9e16bb7556415b9fda9a65f039ba 100644
--- a/matlab/stochastic_solver/set_state_space.m
+++ b/matlab/stochastic_solver/set_state_space.m
@@ -1,6 +1,8 @@
 function dr=set_state_space(dr,M_)
 % dr=set_state_space(dr,M_)
-% Write the state space representation of the reduced form solution.
+% This function computes the DR ordering and inverse ordering.
+% It used to compute stuff related to the state-space representation of the reduced form, hence
+%  its name.
 
 %@info:
 %! @deftypefn {Function File} {[@var{dr} =} set_state_space (@var{dr},@var{M_})
@@ -33,7 +35,7 @@ function dr=set_state_space(dr,M_)
 %! @end deftypefn
 %@eod:
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -68,48 +70,8 @@ else
     both_var = [];
     stat_var = setdiff([1:endo_nbr]',fwrd_var);
 end
-order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
-inv_order_var(order_var) = (1:endo_nbr);
 
-% building kmask for z state vector in t+1
-if max_lag > 0
-    kmask = [];
-    if max_lead > 0
-        kmask = lead_lag_incidence(max_lag+2,order_var) ;
-    end
-    kmask = [kmask; lead_lag_incidence(1,order_var)] ;
-else
-    if max_lead==0 %%in this case lead_lag_incidence has no entry max_lag+2
-        error('Dynare currently does not allow to solve purely static models in a stochastic context.')
-    end
-    kmask = lead_lag_incidence(max_lag+2,order_var) ;
-end
-
-kmask = kmask';
-kmask = kmask(:);
-i_kmask = find(kmask);
-nd = nnz(kmask);           % size of the state vector
-kmask(i_kmask) = (1:nd);
-% auxiliary equations
-
-% composition of state vector
-% col 1: variable (index in DR-order);  col 2: lead/lag in z(t+1);
-% col 3: A cols for t+1 (D);            col 4: A cols for t (E)
-kstate = [ repmat([1:endo_nbr]',klen-1,1) kron([klen:-1:2]',ones(endo_nbr,1)) ...
-           zeros((klen-1)*endo_nbr,2)];
-kiy = flipud(lead_lag_incidence(:,order_var))';
-kiy = kiy(:);
-if max_lead > 0
-    kstate(1:endo_nbr,3) = kiy(1:endo_nbr)-nnz(lead_lag_incidence(max_lag+1,:));
-    kstate(kstate(:,3) < 0,3) = 0;
-    kstate(endo_nbr+1:end,4) = kiy(2*endo_nbr+1:end);
-else
-    kstate(:,4) = kiy(endo_nbr+1:end);
-end
-kstate = kstate(i_kmask,:);
-
-dr.order_var = order_var;
-dr.inv_order_var = inv_order_var';
-dr.kstate = kstate;
+dr.order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
+dr.inv_order_var(dr.order_var) = 1:endo_nbr;
 
 dr.transition_auxiliary_variables = [];
diff --git a/matlab/stochastic_solver/simult_.m b/matlab/stochastic_solver/simult_.m
index ed7f2c096fabb7dc60fa924ce64a1967c7e4f9e7..94466df25435d27821bd8a5ff2b49dd8adf714b7 100644
--- a/matlab/stochastic_solver/simult_.m
+++ b/matlab/stochastic_solver/simult_.m
@@ -17,7 +17,7 @@ function y_=simult_(M_,options_,y0,dr,ex_,iorder)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -62,8 +62,7 @@ if options_.k_order_solver
                        options_.pruning);
     y_(dr.order_var,:) = y_;
 else
-    k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
-    k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
+    k2 = M_.nstatic+(1:M_.nspred);
     order_var = dr.order_var;
 
     switch iorder
diff --git a/matlab/stochastic_solver/simultxdet.m b/matlab/stochastic_solver/simultxdet.m
index 074cc4b5e60f81b1b0b1089951e0b96885934a60..7585ddfe96f3e46b87f2ed64d9c9e6c58ff6444c 100644
--- a/matlab/stochastic_solver/simultxdet.m
+++ b/matlab/stochastic_solver/simultxdet.m
@@ -22,7 +22,7 @@ function [y_,int_width,int_width_ME]=simultxdet(y0,ex,ex_det, iorder,var_list,M_
 % The condition size(ex,1)+M_.maximum_lag=size(ex_det,1) must be verified
 %  for consistency.
 
-% Copyright © 2008-2018 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -54,12 +54,7 @@ nx = size(dr.ghu,2);
 y_ = zeros(size(y0,1),iter+ykmin);
 y_(:,1:ykmin) = y0;
 k1 = ykmin:-1:1;
-k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin+1),[1 2]);
-k2 = k2(:,1)+(ykmin+1-k2(:,2))*endo_nbr;
-k3 = M_.lead_lag_incidence(1:ykmin,:)';
-k3 = find(k3(:));
-k4 = dr.kstate(find(dr.kstate(:,2) < ykmin+1),[1 2]);
-k4 = k4(:,1)+(ykmin+1-k4(:,2))*endo_nbr;
+k2 = nstatic+(1:nspred);
 
 nvar = length(var_list);
 if nvar == 0
diff --git a/matlab/stochastic_solver/stochastic_solvers.m b/matlab/stochastic_solver/stochastic_solvers.m
index 62fb72bb18e516236dd3d48efb02915a93dba503..678810d56ec6eef4ab27ee27b2273debf0115ee4 100644
--- a/matlab/stochastic_solver/stochastic_solvers.m
+++ b/matlab/stochastic_solver/stochastic_solvers.m
@@ -25,7 +25,7 @@ function [dr, info] = stochastic_solvers(dr, task, M_, options_, exo_steady_stat
 %                                 info=6 -> The jacobian matrix evaluated at the steady state is complex.
 %                                 info=9 -> k_order_pert was unable to compute the solution
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -192,14 +192,13 @@ if ~isempty(nanrow)
     return
 end
 
-kstate = dr.kstate;
 nstatic = M_.nstatic;
 nfwrd = M_.nfwrd;
 nspred = M_.nspred;
 nboth = M_.nboth;
 nsfwrd = M_.nsfwrd;
 order_var = dr.order_var;
-nd = size(kstate,1);
+nd = M_.nspred+M_.nsfwrd;
 nz = nnz(M_.lead_lag_incidence);
 
 sdyn = M_.endo_nbr - nstatic;
@@ -212,13 +211,12 @@ b(:,cols_b) = jacobia_(:,cols_j);
 if M_.maximum_endo_lead == 0
     % backward models: simplified code exist only at order == 1
     if local_order == 1
-        [k1,~,k2] = find(kstate(:,4));
         if M_.maximum_endo_lag
             dr.state_var = find(M_.lead_lag_incidence(1,:));
         else
             dr.state_var = [];
         end
-        dr.ghx(:,k1) = -b\jacobia_(:,k2);
+        dr.ghx = -b\jacobia_(:,1:nspred);
         if M_.exo_nbr
             dr.ghu =  -b\jacobia_(:,nz+1:end);
         end
diff --git a/tests/decision_rules/third_order/comparison_policy_functions_dynare_mathematica.m b/tests/decision_rules/third_order/comparison_policy_functions_dynare_mathematica.m
index 12271db230c13d714a4e5b3940edf8343e9793b3..a5f542390216479de5a7d80e81955cea612e2252 100644
--- a/tests/decision_rules/third_order/comparison_policy_functions_dynare_mathematica.m
+++ b/tests/decision_rules/third_order/comparison_policy_functions_dynare_mathematica.m
@@ -5,9 +5,7 @@ order=options_.order;
 dr=oo_.dr;
 nx =size(dr.ghx,2);
 nu =size(dr.ghu,2);
-k = find(dr.kstate(:,2) <= M_.maximum_lag+1);
-klag = dr.kstate(k,[1 2]);
-state_vars=dr.order_var(klag(:,1));
+state_vars=dr.order_var(M_.nstatic+(1:M_.nspred));
 %order of endogenous variables in FV et al. paper: c, invest, y, h, r, dp, kp, lambda, pi
 varlist_FV={'C';'I';'Y';'H';'r';'D';'K';'lambda';'phi'};
 for ii=1: length(varlist_FV)