From 40733b85bb0848be1e91ead1b469214bc89ae952 Mon Sep 17 00:00:00 2001
From: michel <michel@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Sat, 30 Jun 2007 06:33:17 +0000
Subject: [PATCH] v4: new function for ramsey_policy; correction in
 ramsey_dynamic.m

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1333 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/ramsey_dynamic.m |  6 ++--
 matlab/ramsey_policy.m  | 75 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 78 insertions(+), 3 deletions(-)
 create mode 100644 matlab/ramsey_policy.m

diff --git a/matlab/ramsey_dynamic.m b/matlab/ramsey_dynamic.m
index 26d7377d5..3a960a095 100644
--- a/matlab/ramsey_dynamic.m
+++ b/matlab/ramsey_dynamic.m
@@ -133,7 +133,7 @@ function J = ramsey_dynamic(ys,lbar)
       kc2 = i_leadlag1(i,kc(ic2));
       ir = find(i_leadlag(max_leadlag1+2-j,:) > 0 );
       kr1 = i_leadlag(max_leadlag1+2-j,ir);
-      J(ir,kc2) = J(ir,kc2) + beta^(j-max_lag-1)...
+      J(ir,kc2) = J(ir,kc2) + beta^(j-max_lead-1)...
 	  *reshape(lbarfH(iH(kr1,kc1)),length(kr1),length(kc1));
     end
   end
@@ -142,7 +142,7 @@ function J = ramsey_dynamic(ys,lbar)
     kc = i_leadlag1(max_lag+i,endo_nbr+(1:exo_nbr));
     ir = find(i_leadlag(leadlag_nbr+1-i,:) > 0);
     kr1 = i_leadlag(leadlag_nbr+1-i,ir);
-    J(ir,kc) = beta^(i-max_lag-1)...
+    J(ir,kc) = beta^(i-max_lead-1)...
 	*reshape(lbarfH(iH(kr1,n_dyn+(1:exo_nbr))),length(kr1), ...
 			exo_nbr);
   end
@@ -152,7 +152,7 @@ function J = ramsey_dynamic(ys,lbar)
     kc1 = i_leadlag(leadlag_nbr+1-i,ic1);
     ic2 = find(i_leadlag1(max_lag+i,endo_nbr+exo_nbr+(1:mult_nbr)) > 0);
     kc2 = i_leadlag1(max_lag+i,endo_nbr+exo_nbr+ic2);
-    J(ic1,kc2) = beta^(i-max_lag-1)*fJ(ic2,kc1)';
+    J(ic1,kc2) = beta^(i-max_lead-1)*fJ(ic2,kc1)';
   end
 
   % Jacobian of original equations
diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m
new file mode 100644
index 000000000..ec3a4f737
--- /dev/null
+++ b/matlab/ramsey_policy.m
@@ -0,0 +1,75 @@
+function info = ramsey_policy(var_list)
+  global options_ oo_ M_
+  
+  oldoptions = options_;
+  options_.ramsey_policy = 1;
+  options_.order = 1;
+  info = stoch_simul(var_list);
+  
+  return
+  if info
+    return
+  end
+  
+  options_.ramsey_policy = 2;
+  info = stoch_simul(var_list);
+  dr = oo_.dr;
+  
+  orig_model = M_.orig_model;
+  endo_nbr1 = orig_model.endo_nbr;
+  exo_nbr1 = M_.exo_nbr;
+  
+  npred = dr.npred;
+  nspred = dr.nspred;
+  inv_order_var1 = dr.inv_order_var(1:endo_nbr1);
+  nstatic = dr.nstatic;
+  
+  ghx = dr.ghx;
+  ghx1 = ghx(inv_order_var1,:);
+  [ghx2,ghu2] = transition_matrix(dr);
+  ghu = dr.ghu;
+  ghu1 = ghu(inv_order_var1,:);
+  ghxx = dr.ghxx;
+  ghxx1 = ghxx(inv_order_var1,:);
+  ghxx2 = [ghxx(nstatic+(1:npred),:); zeros(nspred-npred,size(ghxx,2))];
+  ghxu = dr.ghxu;
+  ghxu1 = ghxu(inv_order_var1,:);
+  ghxu2 = [ghxu(nstatic+(1:npred),:); zeros(nspred-npred,size(ghxu,2))];
+  ghuu = dr.ghuu;
+  ghuu1 = ghuu(inv_order_var1,:);
+  ghuu2 = [ghuu(nstatic+(1:npred),:); zeros(nspred-npred,size(ghuu,2))];
+  ghs2 = dr.ghs2;  
+  ghs21 = ghs2(inv_order_var1,:);
+  ghs22 = [ghs2(nstatic+(1:npred),:); zeros(nspred-npred,size(ghs2,2))];
+  
+  fname = [M_.fname '_objective_static'];
+  [ubar,uj,uh] = feval(fname,dr.ys(1:endo_nbr1),zeros(exo_nbr1,1));
+  
+  bet = options_.planner_discount;
+  wbar = ubar/(1-bet);
+  % wx = ux*gx + bet*wx*hx
+  % wx*(I-bet*hx) = ux*gx
+  wx = uj*ghx1/(eye(nspred)-bet*ghx2);
+  % wu = ux*gu + bet*wx*hu
+  wu = uj*ghu1 + bet*wx*ghu2;
+  % wxx = uxx*kron(gx,gx)+ux*ghxx + bet(wxx*kron(hx,hx) +
+  %       wx*hxx)
+  wxx = (uh*kron(ghx1,ghx1) + uj*ghxx1 + bet*wx*ghxx2)/(eye(nspred^2)-bet*kron(ghx2,ghx2));
+  wxu = uh*kron(ghx1,ghu1) + uj*ghxu1 + bet*(wx*ghxu2+wxx*kron(ghx2,ghu2));
+  wuu = uh*kron(ghu1,ghu1) + uj*ghuu1 + bet*(wx*ghuu2+wxx*kron(ghu2,ghu2));
+  % ws2 = ux*ghs2 + bet*(ws2+wx*hs2+wuu*Sigma_u)
+  ws2 = (uj*ghs21 + bet*(wx*ghs22+wuu*M_.Sigma_e(:)))/(1-bet);
+  oo_.welfare.wbar = wbar;
+  oo_.welfare.wx = wx;
+  oo_.welfare.wu = wu;
+  oo_.welfare.wxx = wxx;
+  oo_.welfare.wxu = wxu;
+  oo_.welfare.wuu = wuu;
+  oo_.welfare.ws2 = ws2;
+
+  disp(' ')
+  disp(' ')
+  disp(['Welfare at the deterministic steady state is ' num2str(wbar+0.5* ...
+						  ws2)]);
+  
+  options_ = oldoptions;
\ No newline at end of file