improving Ramsey policy

matlab/add_auxiliary_variables_to_steadystate.m

@@ -24,8 +24,10 @@ k = size(ys,1)+1;
 aux_lead_nbr = 0;
 for i=1:n
     if aux_vars(i).type == 1
+        % lagged variables 
         ys1(k) = ys(aux_vars(i).orig_index);
     elseif aux_vars(i).type == 0
+        % lead variables
         aux_lead_nbr = aux_lead_nbr + 1;
     end
     k = k+1;

matlab/dr1.m

@@ -92,14 +92,20 @@ if options_.ramsey_policy
     if options_.steadystate_flag
         k_inst = [];
         instruments = options_.instruments;
-        for i = 1:size(instruments,1)
+        inst_nbr = size(options_.instruments);
+        for i = 1:inst_nbr
             k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
                                        M_.endo_names,'exact')];
         end
         ys = oo_.steady_state;
-        [inst_val,info1] = dynare_solve('dyn_ramsey_static_', ...
-                                        oo_.steady_state(k_inst),0, ...
-                                        M_,options_,oo_,it_);
+        if inst_nbr == 1
+           nl_func = @(x) dyn_ramsey_static_(x,0,M_,options_,oo_,it_);
+           inst_val = fzero(nl_func,oo_.steady_state(k_inst));
+        else
+            [inst_val,info1] = dynare_solve('dyn_ramsey_static_', ...
+                                            oo_.steady_state(k_inst),0, ...
+                                            M_,options_,oo_,it_);
+        end
         M_.params = evalin('base','M_.params;');
         ys(k_inst) = inst_val;
         [x,check] = feval([M_.fname '_steadystate'],...
@@ -120,6 +126,7 @@ if options_.ramsey_policy
         end
         oo_.steady_state = x;
         [junk,junk,multbar] = dyn_ramsey_static_(oo_.steady_state(k_inst),M_,options_,oo_,it_);
+        oo_.steady_state = [x(1:M_.orig_endo_nbr); multbar];
     else
         xx = oo_.steady_state([1:M_.orig_endo_nbr (M_.orig_endo_nbr+M_.orig_eq_nbr+1):end]);
         [xx,info1] = dynare_solve('dyn_ramsey_static_', ...

matlab/dyn_ramsey_static_.m

@@ -96,11 +96,15 @@ y = f(1:orig_endo_nbr);
 mult = -A\y;
 
 resids1 = y+A*mult;
-[q,r,e] = qr([A y]');
+if inst_nbr == 1
+    r1 = sqrt(resids1'*resids1);
+else
+    [q,r,e] = qr([A y]');
+    r1 = r(end,(orig_endo_nbr-inst_nbr+1:end))';
+end
 if options_.steadystate_flag
-    resids = [r(end,(orig_endo_nbr-inst_nbr+1:end))'];
-    resids = resids1'*resids1;
+    resids = r1;
 else
-    resids = [f(i_mult); r(end,(orig_endo_nbr-inst_nbr+1:end))'];
+    resids = [f(i_mult); r1];
 end
 rJ = [];

matlab/evaluate_planner_objective.m

@@ -106,8 +106,8 @@ else
     Wss = (Uy*gss+beta*(Wuu*M.Sigma_e(:)+Wy*Gss))/(1-beta);
 end
 if options.ramsey_policy
-    yhat = [oo.endo_simul; ...
-            zeros(M.endo_nbr-size(oo.endo_simul,1),1)];
+    yhat = zeros(M.endo_nbr,1);
+    yhat(1:M.orig_endo_nbr) = oo.steady_state(1:M.orig_endo_nbr);
 else
     yhat = oo.endo_simul;
 end
@@ -120,11 +120,15 @@ mexErrCheck('A_times_B_kronecker_C', err);
 mexErrCheck('A_times_B_kronecker_C', err);
 [err,Wyuyhatu] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
-planner_objective_value = Wbar+Wy*yhat+Wu*u+Wyuyhatu ...
+planner_objective_value(1) = Wbar+Wy*yhat+Wu*u+Wyuyhatu ...
     + 0.5*(Wyyyhatyhat + Wuuuu+Wss);
+planner_objective_value(2) = Wbar + 0.5*Wss;
 if ~options.noprint
     disp(' ')
-    disp(['Approximated value of planner objective function: ' ...
-          num2str(planner_objective_value)])
+    disp('Approximated value of planner objective function')
+    disp(['    - with initial Lagrange multipliers set to 0: ' ...
+          num2str(planner_objective_value(1)) ])
+    disp(['    - with initial Lagrange multipliers set to steady state: ' ...
+          num2str(planner_objective_value(2)) ])
     disp(' ')
 end
\ No newline at end of file

switch_dw @ b9f9033d

-Subproject commit 654e1a33387bb2cbf4d48af24ef86fd6f65cc6e4
+Subproject commit b9f9033d6fcb94f5ab3148dec73eed9ed4195d93