From 956b42fdff783e376fe916e79c53df52bd6f49b1 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 26 Apr 2010 15:04:42 +0200 Subject: [PATCH] corrected bugs in shock_decomposition --- matlab/evaluate_planner_objective.m | 40 +++++++++++++++++++++-------- matlab/global_initialization.m | 5 ++++ matlab/shock_decomposition.m | 4 +-- 3 files changed, 37 insertions(+), 12 deletions(-) diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m index 4ecd38ca0..3f7396506 100644 --- a/matlab/evaluate_planner_objective.m +++ b/matlab/evaluate_planner_objective.m @@ -50,27 +50,46 @@ function oo1 = evaluate_planner_objective(dr,M,oo,options) i_org = (1:M.endo_nbr)'; end ipred = find(lead_lag_incidence(M.maximum_lag,:))'; + orig_ifrwrd = find(M.orig_model.lead_lag_incidence(3,:)); + orig_npred = nnz(M.orig_model.lead_lag_incidence(1,:)); order_var = dr.order_var; ys = oo.dr.ys; [U,Uy,Uyy] = feval([M.fname '_objective_static'],ys,zeros(1,exo_nbr), ... M.params); + z = repmat(ys(1:M.orig_model.endo_nbr),1,3); + z = z(find(M.orig_model.lead_lag_incidence')); + [F,Fy,Fyy] = feval([M.fname '_dynamic'],z,zeros(3,exo_nbr), ... + M.params,2); + mu_ss = oo.dr.ys(M.orig_model.endo_nbr+exo_nbr+(1:size(F,1))); %K = reshape(1:endo_nbr^2,endo_nbr,endo_nbr); %K = K(order_var,order_var); %Uyy = Uyy(:,K(:)); beta = options.planner_discount; + Gy = dr.ghx(nstatic+(1:npred),:); + Gu = dr.ghu(nstatic+(1:npred),:); gy(dr.order_var,:) = dr.ghx; gu(dr.order_var,:) = dr.ghu; - Gy = gy(ipred,:); - Gu = gu(ipred,:); gy = gy(i_org,:); gu = gu(i_org,:); + muFy = mu_ss'*Fy; + muFyy = mu_ss'*Fyy; + Zy = [eye(orig_npred,npred); gy; gy(orig_ifrwrd,:)*Gy; zeros(exo_nbr, ... + npred)]; + Zu = [zeros(orig_npred,exo_nbr); gu; gy(orig_ifrwrd,:)*Gu; eye(exo_nbr, ... + exo_nbr)]; + Zyy = kron(Zy,Zy); + Zuu = kron(Zu,Zu); + Zyu = kron(Zy,Zu); Wbar = U/(1-beta); - Wy = Uy*gy/(eye(npred)-beta*Gy); - Wu = Uy*gu + beta*Wy*Gu; - Wyy = Uyy*kron(gy,gy)/(eye(npred^2)-beta*kron(Gy,Gy)); - Wyu = Uyy*kron(gy,gu)+beta*Wyy*kron(Gy,Gu); - Wuu = Uyy*kron(gu,gu)+beta*Wyy*kron(Gu,Gu); + Wy = (Uy*gy+muFy*Zy)/(eye(npred)-beta*Gy); + Wu = Uy*gu + muFy*Zu + beta*Wy*Gu; + % N1 = Uyy*kron(gy,gy) + % N2 = muFyy*Zyy + format long + Wyy = (Uyy*kron(gy,gy)+muFyy*Zyy)/(eye(npred^2)-beta*kron(Gy,Gy)) + Wyu = Uyy*kron(gy,gu) + muFyy*Zyu + beta*Wyy*kron(Gy,Gu); + Wuu = Uyy*kron(gu,gu) + muFyy*Zuu + beta*Wyy*kron(Gu,Gu); Wss = beta*Wuu*vec(M.Sigma_e)/(1-beta); if options.ramsey_policy @@ -79,17 +98,18 @@ function oo1 = evaluate_planner_objective(dr,M,oo,options) else yhat = oo.endo_simul; end - yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred))); + yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred))) u = oo.exo_simul(1,:)'; planner_objective_value = Wbar+Wy*yhat+Wu*u+Wyu*kron(yhat,u) ... + 0.5*(Wyy*kron(yhat,yhat) + Wuu*kron(u,u)+Wss); + disp(Wyy*kron(yhat,yhat)) oo1.planner_objective_value = planner_objective_value; if ~options.noprint disp(' ') if options.ramsey_policy - disp(['Value of planner objective function under Ramsey policy: ' ... - num2str(planner_objective_value)]) + disp(sprintf('Value of planner objective function under Ramsey policy: %16.8f', ... + planner_objective_value)) else disp(['Value of planner objective function: ' ... num2str(planner_objective_value)]) diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index dbee353e3..6d774b556 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -235,6 +235,11 @@ options_.use_dll = 0; % model evaluated using bytecode.dll options_.bytecode = 0; +% dates for historical time series +options_.initial_date.freq = 1; +options_.initial_date.period = 1; +options_.initial_date.subperiod = 0; + % SWZ SBVAR options_.ms.freq = 1; options_.ms.initial_subperiod = 1; diff --git a/matlab/shock_decomposition.m b/matlab/shock_decomposition.m index 7f7ef1dd6..1278afa0b 100644 --- a/matlab/shock_decomposition.m +++ b/matlab/shock_decomposition.m @@ -38,7 +38,7 @@ nshocks = M_.exo_nbr; % indices of endogenous variables if size(varlist,1) == 0 - varlist = M_.endo_names(M_.orig_endo_nbr); + varlist = M_.endo_names(1:M_.orig_endo_nbr,:); end [i_var,nvar] = varlist_indices(varlist,M_.endo_names); @@ -95,4 +95,4 @@ options_.initial_date.freq = 1; options_.initial_date.period = 1; options_.initial_date.sub_period = 0; -graph_decomp(z,M_.exo_names,i_var,options_.initial_date) +graph_decomp(z,M_.exo_names,M_.endo_names,i_var,options_.initial_date)