stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

corrections to optimal policy code required by making M_, options_, oo_ and it_ local 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1914 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2008-07-02 07:42:50 +00:00
parent 8a94ae8ba1
commit d0b87d6c32
7 changed files with 104 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

188
matlab/dr1.m
View File

@ -66,102 +66,102 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
end end
old_solve_algo = options_.solve_algo; old_solve_algo = options_.solve_algo;
% options_.solve_algo = 1; % options_.solve_algo = 1;
oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,0); oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,0,M_,options_,oo_,it_);
options_.solve_algo = old_solve_algo; options_.solve_algo = old_solve_algo;
[junk,junk,multbar] = ramsey_static(oo_.steady_state); [junk,junk,multbar] = ramsey_static(oo_.steady_state,M_,options_,oo_,it_);
jacobia_=ramsey_dynamic(oo_.steady_state,multbar); [jacobia_,M_] = ramsey_dynamic(oo_.steady_state,multbar,M_,options_,oo_,it_);
klen = M_.maximum_lag + M_.maximum_lead + 1; klen = M_.maximum_lag + M_.maximum_lead + 1;
dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar]; dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar];
if options_.ramsey_policy == 2 % $$$ if options_.ramsey_policy == 2
mask = M_.orig_model.lead_lag_incidence ~= 0; % $$$ mask = M_.orig_model.lead_lag_incidence ~= 0;
incidence_submatrix = M_.lead_lag_incidence(M_.orig_model.maximum_lead+(1:size(mask,1)),1:M_.orig_model.endo_nbr); % $$$ incidence_submatrix = M_.lead_lag_incidence(M_.orig_model.maximum_lead+(1:size(mask,1)),1:M_.orig_model.endo_nbr);
k = nonzeros((incidence_submatrix.*mask)'); % $$$ k = nonzeros((incidence_submatrix.*mask)');
nl = nnz(M_.lead_lag_incidence); % $$$ nl = nnz(M_.lead_lag_incidence);
k = [k; nl+(1:M_.exo_nbr)']; % $$$ k = [k; nl+(1:M_.exo_nbr)'];
kk = reshape(1:(nl+M_.exo_nbr)^2,nl+M_.exo_nbr,nl+M_.exo_nbr); % $$$ kk = reshape(1:(nl+M_.exo_nbr)^2,nl+M_.exo_nbr,nl+M_.exo_nbr);
kk2 = kk(k,k); % $$$ kk2 = kk(k,k);
% $$$
k1 = find(M_.orig_model.lead_lag_incidence'); % $$$ k1 = find(M_.orig_model.lead_lag_incidence');
y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1); % $$$ y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1);
[f,fJ,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_); % $$$ [f,fJ,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_);
% $$$
% looking for dynamic variables that are both in the original model % $$$ % looking for dynamic variables that are both in the original model
% and in the optimal policy model % $$$ % and in the optimal policy model
k1 = k1+nnz(M_.lead_lag_incidence(1:M_.orig_model.maximum_lead,1:M_.orig_model.endo_nbr)); % $$$ k1 = k1+nnz(M_.lead_lag_incidence(1:M_.orig_model.maximum_lead,1:M_.orig_model.endo_nbr));
hessian = sparse([],[],[],size(jacobia_,1),(nl+M_.exo_nbr)^2,nnz(fh)); % $$$ hessian = sparse([],[],[],size(jacobia_,1),(nl+M_.exo_nbr)^2,nnz(fh));
hessian(M_.orig_model.endo_nbr+(1:size(fh,1)),kk2) = fh; % $$$ hessian(M_.orig_model.endo_nbr+(1:size(fh,1)),kk2) = fh;
options_.order = 2; % $$$ options_.order = 2;
elseif options_.ramsey_policy == 3 % $$$ elseif options_.ramsey_policy == 3
maxlag1 = M_.orig_model.maximum_lag; % $$$ maxlag1 = M_.orig_model.maximum_lag;
maxlead1 = M_.orig_model.maximum_lead; % $$$ maxlead1 = M_.orig_model.maximum_lead;
endo_nbr1 = M_.orig_model.endo_nbr; % $$$ endo_nbr1 = M_.orig_model.endo_nbr;
lead_lag_incidence1 = M_.orig_model.lead_lag_incidence; % $$$ lead_lag_incidence1 = M_.orig_model.lead_lag_incidence;
y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1); % $$$ y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1);
k1 = find(M_.orig_model.lead_lag_incidence'); % $$$ k1 = find(M_.orig_model.lead_lag_incidence');
[f,fj,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_); % $$$ [f,fj,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_);
nrj = size(fj,1); % $$$ nrj = size(fj,1);
% $$$
iy = M_.lead_lag_incidence; % $$$ iy = M_.lead_lag_incidence;
kstate = oo_.dr.kstate; % $$$ kstate = oo_.dr.kstate;
inv_order_var = oo_.dr.inv_order_var; % $$$ inv_order_var = oo_.dr.inv_order_var;
offset = 0; % $$$ offset = 0;
i3 = zeros(0,1); % $$$ i3 = zeros(0,1);
i4 = find(kstate(:,2) <= M_.maximum_lag+1); % $$$ i4 = find(kstate(:,2) <= M_.maximum_lag+1);
kstate1 = kstate(i4,:); % $$$ kstate1 = kstate(i4,:);
kk2 = zeros(0,1); % $$$ kk2 = zeros(0,1);
% lagged variables % $$$ % lagged variables
for i=2:M_.maximum_lag + 1 % $$$ for i=2:M_.maximum_lag + 1
i1 = find(kstate(:,2) == i); % $$$ i1 = find(kstate(:,2) == i);
k1 = kstate(i1,:); % $$$ k1 = kstate(i1,:);
i2 = find(oo_.dr.order_var(k1(:,1)) <= M_.orig_model.endo_nbr); % $$$ i2 = find(oo_.dr.order_var(k1(:,1)) <= M_.orig_model.endo_nbr);
i3 = [i3; i2+offset]; % $$$ i3 = [i3; i2+offset];
offset = offset + size(k1,1); % $$$ offset = offset + size(k1,1);
i4 = find(kstate1(:,2) == i); % $$$ i4 = find(kstate1(:,2) == i);
kk2 = [kk2; i4]; % $$$ kk2 = [kk2; i4];
end % $$$ end
i2 = find(oo_.dr.order_var(k1(:,1)) > M_.orig_model.endo_nbr); % $$$ i2 = find(oo_.dr.order_var(k1(:,1)) > M_.orig_model.endo_nbr);
j2 = k1(i2,1); % $$$ j2 = k1(i2,1);
nj2 = length(j2); % $$$ nj2 = length(j2);
k2 = offset+(1:nj2)'; % $$$ k2 = offset+(1:nj2)';
offset = offset + length(i2); % $$$ offset = offset + length(i2);
i3 = [i3; ... % $$$ i3 = [i3; ...
find(M_.orig_model.lead_lag_incidence(M_.orig_model.maximum_lag+1:end,:)')+offset]; % $$$ find(M_.orig_model.lead_lag_incidence(M_.orig_model.maximum_lag+1:end,:)')+offset];
i3 = [i3; (1:M_.exo_nbr)'+length(i3)]; % $$$ i3 = [i3; (1:M_.exo_nbr)'+length(i3)];
ni3 = length(i3); % $$$ ni3 = length(i3);
nrfj = size(fj,1); % $$$ nrfj = size(fj,1);
jacobia_ = zeros(nrfj+length(j2),ni3); % $$$ jacobia_ = zeros(nrfj+length(j2),ni3);
jacobia_(1:nrfj,i3) = fj; % $$$ jacobia_(1:nrfj,i3) = fj;
jacobia_(nrfj+(1:nj2),1:size(oo_.dr.ghx,2)) = oo_.dr.ghx(j2,:); % $$$ jacobia_(nrfj+(1:nj2),1:size(oo_.dr.ghx,2)) = oo_.dr.ghx(j2,:);
jacobia_(nrfj+(1:nj2),k2) = eye(nj2); % $$$ jacobia_(nrfj+(1:nj2),k2) = eye(nj2);
kk1 = reshape(1:ni3^2,ni3,ni3); % $$$ kk1 = reshape(1:ni3^2,ni3,ni3);
hessian = zeros(nrfj+length(j2),ni3^2); % $$$ hessian = zeros(nrfj+length(j2),ni3^2);
hessian(1:nrfj,kk1(i3,i3)) = fh; % $$$ hessian(1:nrfj,kk1(i3,i3)) = fh;
% $$$
k = find(any(M_.lead_lag_incidence(1:M_.maximum_lag, ... % $$$ k = find(any(M_.lead_lag_incidence(1:M_.maximum_lag, ...
M_.orig_model.endo_nbr+1:end))); % $$$ M_.orig_model.endo_nbr+1:end)));
if maxlead1 > maxlag1 % $$$ if maxlead1 > maxlag1
M_.lead_lag_incidence = [ [zeros(maxlead1-maxlag1,endo_nbr1); ... % $$$ M_.lead_lag_incidence = [ [zeros(maxlead1-maxlag1,endo_nbr1); ...
lead_lag_incidence1] ... % $$$ lead_lag_incidence1] ...
[M_.lead_lag_incidence(M_.maximum_lag+(1:maxlead1), ... % $$$ [M_.lead_lag_incidence(M_.maximum_lag+(1:maxlead1), ...
k); zeros(maxlead1,length(k))]]; % $$$ k); zeros(maxlead1,length(k))]];
elseif maxlag1 > maxlead1 % $$$ elseif maxlag1 > maxlead1
M_.lead_lag_incidence = [ [lead_lag_incidence1; zeros(maxlag1- ... % $$$ M_.lead_lag_incidence = [ [lead_lag_incidence1; zeros(maxlag1- ...
maxlead1,endo_nbr1);] ... % $$$ maxlead1,endo_nbr1);] ...
[M_.lead_lag_incidence(M_.maximum_lag+(1:maxlead1), ... % $$$ [M_.lead_lag_incidence(M_.maximum_lag+(1:maxlead1), ...
k); zeros(maxlead1,length(k))]]; % $$$ k); zeros(maxlead1,length(k))]];
else % maxlag1 == maxlead1 % $$$ else % maxlag1 == maxlead1
M_.lead_lag_incidence = [ lead_lag_incidence1 % $$$ M_.lead_lag_incidence = [ lead_lag_incidence1
[M_.lead_lag_incidence(M_.maximum_lag+(1:maxlead1), ... % $$$ [M_.lead_lag_incidence(M_.maximum_lag+(1:maxlead1), ...
k); zeros(maxlead1,length(k))]]; % $$$ k); zeros(maxlead1,length(k))]];
end % $$$ end
M_.maximum_lag = max(maxlead1,maxlag1); % $$$ M_.maximum_lag = max(maxlead1,maxlag1);
M_.maximum_endo_lag = M_.maximum_lag; % $$$ M_.maximum_endo_lag = M_.maximum_lag;
M_.maximum_lead = M_.maximum_lag; % $$$ M_.maximum_lead = M_.maximum_lag;
M_.maximum_endo_lead = M_.maximum_lag; % $$$ M_.maximum_endo_lead = M_.maximum_lag;
% $$$
M_.endo_names = strvcat(M_.orig_model.endo_names, M_.endo_names(endo_nbr1+k,:)); % $$$ M_.endo_names = strvcat(M_.orig_model.endo_names, M_.endo_names(endo_nbr1+k,:));
M_.endo_nbr = endo_nbr1+length(k); % $$$ M_.endo_nbr = endo_nbr1+length(k);
end % $$$ end
else else
klen = M_.maximum_lag + M_.maximum_lead + 1; klen = M_.maximum_lag + M_.maximum_lead + 1;
iyv = M_.lead_lag_incidence'; iyv = M_.lead_lag_incidence';
@ -190,7 +190,7 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
save([M_.fname '_debug.mat'],'jacobia_') save([M_.fname '_debug.mat'],'jacobia_')
end end
dr=set_state_space(dr); dr=set_state_space(dr,M_);
kstate = dr.kstate; kstate = dr.kstate;
kad = dr.kad; kad = dr.kad;
kae = dr.kae; kae = dr.kae;

2
matlab/dynare_estimation.m
View File

@ -127,7 +127,7 @@ end
bayestopt_.penalty = 1e8; % penalty bayestopt_.penalty = 1e8; % penalty
dr = set_state_space([]); dr = set_state_space([],M_);
nstatic = dr.nstatic; nstatic = dr.nstatic;
npred = dr.npred; npred = dr.npred;
nspred = dr.nspred; nspred = dr.nspred;

2
matlab/dynare_estimation_init.m
View File

@ -129,7 +129,7 @@ end
bayestopt_.penalty = 1e8; % penalty bayestopt_.penalty = 1e8; % penalty
dr = set_state_space([]); dr = set_state_space([],M_);
nstatic = dr.nstatic; nstatic = dr.nstatic;
npred = dr.npred; npred = dr.npred;
nspred = dr.nspred; nspred = dr.nspred;

2
matlab/osr1.m
View File

@ -22,7 +22,7 @@ function osr1(i_params,i_var,weights)
exe =zeros(M_.exo_nbr,1); exe =zeros(M_.exo_nbr,1);
oo_.dr = set_state_space(oo_.dr); oo_.dr = set_state_space(oo_.dr,M_);
% check if ys is steady state % check if ys is steady state
fh = str2func([M_.fname '_static']); fh = str2func([M_.fname '_static']);

8
matlab/ramsey_dynamic.m
View File

@ -1,4 +1,4 @@
function J = ramsey_dynamic(ys,lbar) function [J,M_] = ramsey_dynamic(ys,lbar,M_,options_,oo_,it_)
% function J = ramsey_dynamic(ys,lbar) % function J = ramsey_dynamic(ys,lbar)
% ramsey_dynamic sets up the Jacobian of the expanded model for optimal % ramsey_dynamic sets up the Jacobian of the expanded model for optimal
@ -18,8 +18,6 @@ function J = ramsey_dynamic(ys,lbar)
% Gnu Public License. % Gnu Public License.
global M_ options_ it_
% retrieving model parameters % retrieving model parameters
endo_nbr = M_.endo_nbr; endo_nbr = M_.endo_nbr;
i_endo_nbr = 1:endo_nbr; i_endo_nbr = 1:endo_nbr;
@ -52,14 +50,14 @@ function J = ramsey_dynamic(ys,lbar)
y = repmat(ys,1,max_lag+max_lead+1); y = repmat(ys,1,max_lag+max_lead+1);
k = find(i_leadlag'); k = find(i_leadlag');
it_ = 1;
% retrieving derivatives of the objective function % retrieving derivatives of the objective function
[U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr), M_.params); [U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr), M_.params);
Uy = Uy'; Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr); Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% retrieving derivatives of original model % retrieving derivatives of original model
[f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr), M_.params, it_); [f,fJ,fH] = feval([fname '_dynamic'],y(k),[oo_.exo_simul oo_.exo_det_simul], M_.params, it_);
instr_nbr = endo_nbr - size(f,1); instr_nbr = endo_nbr - size(f,1);
mult_nbr = endo_nbr-instr_nbr; mult_nbr = endo_nbr-instr_nbr;

5
matlab/ramsey_static.m
View File

@ -1,4 +1,4 @@
function [resids, rJ,mult] = ramsey_static(x) function [resids, rJ,mult] = ramsey_static(x,M_,options_,oo_,it_)
% function [resids, rJ,mult] = ramsey_static(x) % function [resids, rJ,mult] = ramsey_static(x)
% Computes the static first order conditions for optimal policy % Computes the static first order conditions for optimal policy
@ -17,7 +17,6 @@ function [resids, rJ,mult] = ramsey_static(x)
% part of DYNARE, copyright Dynare Team (2003-2007) % part of DYNARE, copyright Dynare Team (2003-2007)
% Gnu Public License. % Gnu Public License.
global M_ options_ it_
% recovering usefull fields % recovering usefull fields
endo_nbr = M_.endo_nbr; endo_nbr = M_.endo_nbr;
@ -47,7 +46,7 @@ function [resids, rJ,mult] = ramsey_static(x)
k = find(i_lag'); k = find(i_lag');
it_ = 1; it_ = 1;
% [f,fJ,fH] = feval([fname '_dynamic'],y(k),ex); % [f,fJ,fH] = feval([fname '_dynamic'],y(k),ex);
[f,fJ] = feval([fname '_dynamic'],y(k),ex, M_.params, it_); [f,fJ] = feval([fname '_dynamic'],y(k),[oo_.exo_simul oo_.exo_det_simul], M_.params, it_);
% indices of Lagrange multipliers % indices of Lagrange multipliers
inst_nbr = endo_nbr - size(f,1); inst_nbr = endo_nbr - size(f,1);
i_mult = [endo_nbr+1:2*endo_nbr-inst_nbr]'; i_mult = [endo_nbr+1:2*endo_nbr-inst_nbr]';

6
matlab/set_state_space.m
View File

@ -1,5 +1,5 @@
function dr=set_state_space(dr) function dr=set_state_space(dr,M_)
% function dr = set_state_space(dr) % function dr = set_state_space(dr,M_)
% finds the state vector for structural state space representation % finds the state vector for structural state space representation
% sets many fields of dr % sets many fields of dr
% %
@ -18,8 +18,6 @@ function dr=set_state_space(dr)
% part of DYNARE, copyright Dynare Team (1996-2007) % part of DYNARE, copyright Dynare Team (1996-2007)
% Gnu Public License. % Gnu Public License.
global M_ oo_ options_ it_
xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1; xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1;
klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1; klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;