diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index 196ae7d5e354c6ba4e7f0faaab749344357f753d..8a4515441b32b250c036f9d071e06d533c77d6b0 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -321,6 +321,14 @@ options_ident = set_default_option(options_ident,'analytic_derivation_mode', opt
     %  1: kronecker products method to compute analytical derivatives as in Iskrev (2010) (only for order=1)
     % -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification_numerical_objective.m (previously thet2tau.m)
     % -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the identification Jacobians are then computed analytically as with 0
+
+if options_.discretionary_policy || options_.ramsey_policy
+    if options_ident.analytic_derivation_mode~=-1
+        fprintf('dynare_identification: discretionary_policy and ramsey_policy require analytic_derivation_mode=-1. Resetting the option.')
+        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
diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m
index 453cb4ef4fde031bdf996e0853ba4bdc915109a2..669ce2fb3a0084268c8266592b9ba0d53a9f5880 100644
--- a/matlab/get_identification_jacobians.m
+++ b/matlab/get_identification_jacobians.m
@@ -159,9 +159,9 @@ if order == 1
     [~, g1 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1);
     %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
     DYNAMIC = [Yss;
-               vec(g1(oo.dr.order_var,:))]; %add steady state and put rows of g1 in DR order
+               vec(g1)]; %add steady state and put rows of g1 in DR order
     dDYNAMIC = [oo.dr.derivs.dYss;
-                reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state
+                reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state
     REDUCEDFORM = [Yss;
                    vec(oo.dr.ghx);
                    dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu))]; %in DR order
@@ -177,11 +177,11 @@ elseif order == 2
     %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
     %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
     DYNAMIC = [Yss;
-               vec(g1(oo.dr.order_var,:));
-               vec(g2(oo.dr.order_var,:))]; %add steady state and put rows of g1 and g2 in DR order
+               vec(g1);
+               vec(g2)]; %add steady state and put rows of g1 and g2 in DR order
     dDYNAMIC = [oo.dr.derivs.dYss;
-                reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3));  %reshape dg1 in DR order
-                reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order
+                reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3));  %reshape dg1 in DR order
+                reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order
     REDUCEDFORM = [Yss;
                    vec(oo.dr.ghx);
                    dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu));
@@ -204,13 +204,13 @@ elseif order == 3
     %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
     %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
     DYNAMIC = [Yss;
-               vec(g1(oo.dr.order_var,:));
-               vec(g2(oo.dr.order_var,:));
-               vec(g3(oo.dr.order_var,:))]; %add steady state and put rows of g1 and g2 in DR order
+               vec(g1);
+               vec(g2);
+               vec(g3)]; %add steady state and put rows of g1 and g2 in DR order
     dDYNAMIC = [oo.dr.derivs.dYss;
-                reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3));  %reshape dg1 in DR order
-                reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3));
-                reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order
+                reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3));  %reshape dg1 in DR order
+                reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3));
+                reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order
     REDUCEDFORM = [Yss;
                    vec(oo.dr.ghx);
                    dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu));
diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m
index 2deb29da74655fa1dd9924b1bc149394d84daddf..f23a58d17c70aeabdca3a96a5559058542485aa6 100644
--- a/matlab/get_perturbation_params_derivs.m
+++ b/matlab/get_perturbation_params_derivs.m
@@ -344,9 +344,17 @@ if analytic_derivation_mode == -1
     %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
     dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options);
     ind_Yss = 1:endo_nbr;
-    ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
+    if options.discretionary_policy || options.ramsey_policy
+        ind_g1 = ind_Yss(end) + (1:M.eq_nbr*yy0ex0_nbr);
+    else
+        ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
+    end
     DERIVS.dYss = dYss_g(ind_Yss, :); %in tensor notation, wrt selected model parameters only
-    DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
+    if options.discretionary_policy || options.ramsey_policy
+        DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[M.eq_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
+    else
+        DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
+    end
     if order > 1
         ind_g2 = ind_g1(end) + (1:endo_nbr*yy0ex0_nbr^2);
         DERIVS.dg2 = reshape(sparse(dYss_g(ind_g2,:)),[endo_nbr, yy0ex0_nbr^2*modparam_nbr]); %blockwise in matrix notation, i.e. [dg2_dp1 dg2_dp2 ...], where dg2_dpj has dimension endo_nbr by yy0ex0_nbr^2
diff --git a/matlab/get_perturbation_params_derivs_numerical_objective.m b/matlab/get_perturbation_params_derivs_numerical_objective.m
index 5193bc323ae10697f7866cb1a9a05ab6724340ac..c888eef48b4f7a6e0e452329dbbefa6b174a674c 100644
--- a/matlab/get_perturbation_params_derivs_numerical_objective.m
+++ b/matlab/get_perturbation_params_derivs_numerical_objective.m
@@ -51,7 +51,7 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
 
 %% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
 M = set_all_parameters(params,estim_params,M);
-[~,info,M,options,oo] = resol(0,M,options,oo);
+[~,info,M,options,oo] = compute_decision_rules(M,options,oo);
 Sigma_e = M.Sigma_e;
 
 if info(1) > 0
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index ff773195cf7810d42f5618469f271a4b3b0f7140..5f2b47660e73354c3306b91669e7a95998704c1d 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -133,7 +133,7 @@ no_identification_minimal     = options_ident.no_identification_minimal;
 no_identification_spectrum    = options_ident.no_identification_spectrum;
 
 %Compute linear approximation and fill dr structure
-[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
+[oo_.dr,info,M_,options_,oo_] = compute_decision_rules(M_,options_,oo_);
 
 if info(1) == 0 %no errors in solution
     % Compute parameter Jacobians for identification analysis
diff --git a/matlab/identification_numerical_objective.m b/matlab/identification_numerical_objective.m
index 66112ca6b3960abe936cb355125a92cb9135d806..f38575af9d352bc89e1ca4893311138d0152f964 100644
--- a/matlab/identification_numerical_objective.m
+++ b/matlab/identification_numerical_objective.m
@@ -76,7 +76,7 @@ else
 end
 
 %% compute Kalman transition matrices and steady state with updated parameters
-[~,info,M,options,oo] = resol(0,M,options,oo);
+[~,info,M,options,oo] = compute_decision_rules(M,options,oo);
 options = rmfield(options,'options_ident');
 pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0);
 
diff --git a/tests/discretionary_policy/dennis_1_estim.mod b/tests/discretionary_policy/dennis_1_estim.mod
index f869a69e661681d0f2e4b473a72ccfd8bed145e0..ea05d899e8a5b6fe084e7bb1f9732ec0e875b2d4 100644
--- a/tests/discretionary_policy/dennis_1_estim.mod
+++ b/tests/discretionary_policy/dennis_1_estim.mod
@@ -36,8 +36,10 @@ estimated_params;
 end;
 
 options_.plot_priors=0;
-estimation(order = 1, datafile = dennis_simul, mh_replic = 2000, mh_nblocks=1,smoother,bayesian_irf,moments_varendo) y i pi pi_c q;
+estimation(order = 1, datafile = dennis_simul, mh_replic = 2000, mh_nblocks=1,smoother,bayesian_irf,moments_varendo, conditional_variance_decomposition=[1,2]) y i pi pi_c q;
 
 if max(abs(oo_.posterior.optimization.mode - [1; 0.3433])) > 0.025
   error('Posterior mode too far from true parameter values');
 end
+
+identification;
\ No newline at end of file
diff --git a/tests/optimal_policy/Ramsey/Ramsey_Example_estimation.mod b/tests/optimal_policy/Ramsey/Ramsey_Example_estimation.mod
index b51233fc1bcc83b141632b7751cdd56dcae9f090..4b7845c17a985bb32512c41f9aaa0c7263bc5ffa 100644
--- a/tests/optimal_policy/Ramsey/Ramsey_Example_estimation.mod
+++ b/tests/optimal_policy/Ramsey/Ramsey_Example_estimation.mod
@@ -222,7 +222,7 @@ end;
         ramsey_model(instruments=(R),planner_discount=beta,planner_discount_latex_name=$\beta$); 
         
         //conduct stochastic simulations of the Ramsey problem
-        stoch_simul(order=1,irf=20,periods=500) pi_ann log_h R_ann log_C Z r_real;
+        stoch_simul(TeX,order=1,irf=20,periods=500) pi_ann log_h R_ann log_C Z r_real;
         evaluate_planner_objective;
         
         @# if Estimation_under_Ramsey==1
@@ -234,6 +234,7 @@ end;
             varobs log_C;
                 
             estimation(datafile=ramsey_simulation,mode_compute=5,mh_nblocks=1,mh_replic=0);
+            identification(parameter_set=posterior_mode);
         @# endif        
     @# endif
 @# endif