diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index 26e23b4a1..878b06ec3 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -42,6 +42,13 @@ function planner_objective_value = evaluate_planner_objective(M_,options_,oo_)
 % The approximated conditional expectation of the planner's objective function starting from the non-stochastic steady-state and allowing for future shocks thus verifies
 % W(y,0,1) = Wbar + 0.5*Wss
 
+% In the discretionary case, the model is assumed to be linear and the utility is assumed to be linear-quadratic. This changes 2 aspects of the results delinated above:
+% 1) the second-order derivatives of the policy and transition functions h and g are zero. 
+% 2) the unconditional expectation of states coincides with its steady-state, which entails E(yhat) = 0
+% Therefore, the unconditional welfare can now be approximated as
+% E(W) = (1 - beta)^{-1} ( Ubar + 0.5 ( U_xx h_y^2 E(yhat^2) + U_xx h_u^2 E(u^2) )
+% As for the conditional welfare, the second-order formula above is still valid, but the derivatives of W no longer contain any second-order derivatives of the policy and transition functions h and g.
+
 % INPUTS
 %   M_:        (structure) model description
 %   options_:  (structure) options
@@ -76,26 +83,91 @@ beta = get_optimal_policy_discount_factor(M_.params, M_.param_names);
 
 ys = oo_.dr.ys;
 planner_objective_value = zeros(2,1);
-if options_.order == 1
-    [U] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
-    planner_objective_value(1) = U/(1-beta);
-    planner_objective_value(2) = U/(1-beta);  
-elseif options_.order == 2
+if options_.ramsey_policy
+    if options_.order == 1
+        [U] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
+        planner_objective_value(1) = U/(1-beta);
+        planner_objective_value(2) = U/(1-beta);  
+    elseif options_.order == 2
+        [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
+
+        Gy = dr.ghx(nstatic+(1:nspred),:);
+        Gu = dr.ghu(nstatic+(1:nspred),:);
+        Gyy = dr.ghxx(nstatic+(1:nspred),:);
+        Gyu = dr.ghxu(nstatic+(1:nspred),:);
+        Guu = dr.ghuu(nstatic+(1:nspred),:);
+        Gss = dr.ghs2(nstatic+(1:nspred),:);
+
+        gy(dr.order_var,:) = dr.ghx;
+        gu(dr.order_var,:) = dr.ghu;
+        gyy(dr.order_var,:) = dr.ghxx;
+        gyu(dr.order_var,:) = dr.ghxu;
+        guu(dr.order_var,:) = dr.ghuu;
+        gss(dr.order_var,:) = dr.ghs2;
+
+        Uyy = full(Uyy);
+
+        Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
+        Uyygugu = A_times_B_kronecker_C(Uyy,gu,gu);
+
+        %% Unconditional welfare
+
+        old_noprint = options_.noprint;
+
+        if ~old_noprint
+            options_.noprint = 1;
+        end
+        var_list = M_.endo_names(dr.order_var(nstatic+(1:nspred)));
+        [info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list); %get decision rules and moments
+        if ~old_noprint
+            options_.noprint = 0;
+        end
+
+        oo_.mean(isnan(oo_.mean)) = options_.huge_number;
+        oo_.var(isnan(oo_.var)) = options_.huge_number;
+
+        Ey = oo_.mean;
+        Eyhat = Ey - ys(dr.order_var(nstatic+(1:nspred)));
+
+        var_corr = Eyhat*Eyhat';
+        Eyhatyhat = oo_.var(:) + var_corr(:);
+        Euu = M_.Sigma_e(:);
+
+        EU = U + Uy*gy*Eyhat + 0.5*((Uyygygy + Uy*gyy)*Eyhatyhat + (Uyygugu + Uy*guu)*Euu + Uy*gss);
+        EW = EU/(1-beta);
+
+        %% Conditional welfare starting from the non-stochastic steady-state
+
+        Wbar = U/(1-beta);
+        Wy = Uy*gy/(eye(nspred)-beta*Gy);
+
+        if isempty(options_.qz_criterium)
+            options_.qz_criterium = 1+1e-6;
+        end
+        %solve Lyapunuv equation Wyy=gy'*Uyy*gy+Uy*gyy+beta*Wy*Gyy+beta*Gy'Wyy*Gy
+        Wyy = reshape(lyapunov_symm(sqrt(beta)*Gy',reshape(Uyygygy + Uy*gyy + beta*Wy*Gyy,nspred,nspred),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, options_.debug),1,nspred*nspred);
+        Wyygugu = A_times_B_kronecker_C(Wyy,Gu,Gu);
+        Wuu = Uyygugu + Uy*guu + beta*(Wyygugu + Wy*Guu);
+        Wss = (Uy*gss + beta*(Wy*Gss + Wuu*M_.Sigma_e(:)))/(1-beta);
+        W = Wbar + 0.5*Wss;
+
+        planner_objective_value(1) = EW;
+        planner_objective_value(2) = W;
+    else
+        %Order k code will go here!
+        fprintf('\nevaluate_planner_objective: order>2 not yet supported\n')
+        planner_objective_value(1) = NaN;
+        planner_objective_value(2) = NaN;
+        return
+    end
+elseif options_.discretionary_policy
+    
     [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
 
     Gy = dr.ghx(nstatic+(1:nspred),:);
     Gu = dr.ghu(nstatic+(1:nspred),:);
-    Gyy = dr.ghxx(nstatic+(1:nspred),:);
-    Gyu = dr.ghxu(nstatic+(1:nspred),:);
-    Guu = dr.ghuu(nstatic+(1:nspred),:);
-    Gss = dr.ghs2(nstatic+(1:nspred),:);
-
     gy(dr.order_var,:) = dr.ghx;
     gu(dr.order_var,:) = dr.ghu;
-    gyy(dr.order_var,:) = dr.ghxx;
-    gyu(dr.order_var,:) = dr.ghxu;
-    guu(dr.order_var,:) = dr.ghuu;
-    gss(dr.order_var,:) = dr.ghs2;
 
     Uyy = full(Uyy);
 
@@ -125,9 +197,10 @@ elseif options_.order == 2
     Eyhatyhat = oo_.var(:) + var_corr(:);
     Euu = M_.Sigma_e(:);
 
-    EU = U + Uy*gy*Eyhat + 0.5*((Uyygygy + Uy*gyy)*Eyhatyhat + (Uyygugu + Uy*guu)*Euu + Uy*gss);
+    EU = U + Uy*gy*Eyhat + 0.5*(Uyygygy*Eyhatyhat + Uyygugu*Euu);
     EW = EU/(1-beta);
 
+
     %% Conditional welfare starting from the non-stochastic steady-state
 
     Wbar = U/(1-beta);
@@ -136,27 +209,23 @@ elseif options_.order == 2
     if isempty(options_.qz_criterium)
         options_.qz_criterium = 1+1e-6;
     end
-    %solve Lyapunuv equation Wyy=gy'*Uyy*gy+Uy*gyy+beta*Wy*Gyy+beta*Gy'Wyy*Gy
-    Wyy = reshape(lyapunov_symm(sqrt(beta)*Gy',reshape(Uyygygy + Uy*gyy + beta*Wy*Gyy,nspred,nspred),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, options_.debug),1,nspred*nspred);
+    %solve Lyapunuv equation Wyy=gy'*Uyy*gy+beta*Gy'Wyy*Gy
+    Wyy = reshape(lyapunov_symm(sqrt(beta)*Gy',reshape(Uyygygy,nspred,nspred),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, options_.debug),1,nspred*nspred);
     Wyygugu = A_times_B_kronecker_C(Wyy,Gu,Gu);
-    Wuu = Uyygugu + Uy*guu + beta*(Wyygugu + Wy*Guu);
-    Wss = (Uy*gss + beta*(Wy*Gss + Wuu*M_.Sigma_e(:)))/(1-beta);
+    Wuu = Uyygugu + beta*Wyygugu;
+    Wss = beta*Wuu*M_.Sigma_e(:)/(1-beta);
     W = Wbar + 0.5*Wss;
 
     planner_objective_value(1) = EW;
     planner_objective_value(2) = W;
-else
-    %Order k code will go here!
-    fprintf('\nevaluate_planner_objective: order>2 not yet supported\n')
-    planner_objective_value(1) = NaN;
-    planner_objective_value(2) = NaN;
-    return
 end
 if ~options_.noprint
     if options_.ramsey_policy
         fprintf('\nApproximated value of unconditional welfare:  %10.8f\n', planner_objective_value(1))
         fprintf('\nApproximated value of conditional welfare:  %10.8f\n', planner_objective_value(2))
     elseif options_.discretionary_policy
-        fprintf('\nApproximated value of unconditional welfare with discretionary policy: %10.8f\n\n', EW)
+        fprintf('\nApproximated value of unconditional welfare with discretionary policy:  %10.8f\n', planner_objective_value(1))
+        fprintf('\nApproximated value of conditional welfare with discretionary policy:  %10.8f\n', planner_objective_value(2))
+
     end
 end
diff --git a/tests/discretionary_policy/Gali_discretion.mod b/tests/discretionary_policy/Gali_discretion.mod
index fc596311e..af32cfa54 100644
--- a/tests/discretionary_policy/Gali_discretion.mod
+++ b/tests/discretionary_policy/Gali_discretion.mod
@@ -1,5 +1,5 @@
 /*
- * This file implements the baseline New Keynesian model of Jordi Galí (2008): Monetary Policy, Inflation,
+ * This file implements the baseline New Keynesian model of Jordi Gal� (2008): Monetary Policy, Inflation,
  * and the Business Cycle, Princeton University Press, Chapter 5
  *
  * This implementation was written by Johannes Pfeifer. 
@@ -9,7 +9,7 @@
  */
 
 /*
- * Copyright (C) 2013-15 Johannes Pfeifer
+ * Copyright (C) 2013-21 Johannes Pfeifer
  *
  * This is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -130,7 +130,7 @@ end
 %Compute theoretical objective function
 V=betta/(1-betta)*(var_pi_theoretical+alpha_x*var_y_gap_theoretical); %evaluate at steady state in first period
 
-if isnan(oo_.planner_objective_value) || abs(V-oo_.planner_objective_value)>1e-10
+if any( [ isnan(oo_.planner_objective_value(2)), abs(V-oo_.planner_objective_value(2))>1e-10 ] )
     error('Computed welfare deviates from theoretical welfare')
 end
 end;
@@ -144,6 +144,6 @@ end;
 V=var_pi_theoretical+alpha_x*var_y_gap_theoretical+ betta/(1-betta)*(var_pi_theoretical+alpha_x*var_y_gap_theoretical); %evaluate at steady state in first period
  
 discretionary_policy(instruments=(i),irf=20,discretionary_tol=1e-12,planner_discount=betta) y_gap pi p u;
-if isnan(oo_.planner_objective_value) || abs(V-oo_.planner_objective_value)>1e-10
+if any( [ isnan(oo_.planner_objective_value(1)), abs(V-oo_.planner_objective_value(1))>1e-10 ] )
     error('Computed welfare deviates from theoretical welfare')
 end