Skip to content

GitLab

Commit 927da0e4 authored by michel's avatar michel
Browse files

v4: added code for Ramsey policy 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1124 ac1d8469-bf42-47a9-8791-bf33cf982152
parent d0351756
Hide whitespace changes
Inline Side-by-side
matlab/dr1.m
View file @ 927da0e4
......@@ -50,97 +50,50 @@ if M_.exo_nbr == 0
oo_.exo_steady_state = [] ;
end
tempex = oo_.exo_simul;
it_ = M_.maximum_lag + 1;
% expanding system for Optimal Linear Regulator
if options_.olr
z = repmat(zeros(M_.endo_nbr,1),1,klen);
if isfield(M_,'orig_model')
orig_model = M_.orig_model;
M_.endo_nbr = orig_model.endo_nbr;
M_.endo_names = orig_model.endo_names;
M_.lead_lag_incidence = orig_model.lead_lag_incidence;
M_.maximum_lead = orig_model.maximum_lead;
M_.maximum_endo_lead = orig_model.maximum_endo_lead;
M_.maximum_lag = orig_model.maximum_lag;
M_.maximum_endo_lag = orig_model.maximum_endo_lag;
end
oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,0);
[junk,junk,multbar] = ramsey_static(oo_.steady_state);
jacobia_=ramsey_dynamic(oo_.steady_state,multbar);
klen = M_.maximum_lag + M_.maximum_lead + 1;
dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar];
else
klen = M_.maximum_lag + M_.maximum_lead + 1;
iyv = M_.lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
it_ = M_.maximum_lag + 1 ;
if M_.exo_nbr == 0
oo_.exo_steady_state = [] ;
end
tempex = oo_.exo_simul;
it_ = M_.maximum_lag + 1;
z = repmat(dr.ys,1,klen);
end
z = z(iyr0) ;
if options_.order == 1
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,tempex);
elseif options_.order == 2
z = z(iyr0) ;
if options_.order == 1
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,tempex);
elseif options_.order == 2
[junk,jacobia_,hessian] = feval([M_.fname '_dynamic'],z,...
[oo_.exo_simul ...
oo_.exo_det_simul]);
end
oo_.exo_simul = tempex ;
tempex = [];
% expanding system for Optimal Linear Regulator
if options_.olr
bet = options_.olr_beta;
jacobia1 = [];
n_inst = size(options_.olr_inst,1);
if ~isfield(olr_state_,'done')
olr_state_.done = 1;
olr_state_.old_M_.maximum_endo_lag = M_.maximum_endo_lag;
olr_state_.old_M_.maximum_endo_lead = M_.maximum_endo_lead;
olr_state_.old_M_.endo_nbr = M_.endo_nbr;
olr_state_.old_M_.lead_lag_incidence = M_.lead_lag_incidence;
for i=1:M_.endo_nbr
temp = ['mult_' int2str(i)];
lgoo_.endo_simul = strvcat(lgoo_.endo_simul,temp);
end
M_.endo_nbr = 2*M_.endo_nbr-n_inst;
M_.maximum_endo_lag = max(M_.maximum_endo_lag,M_.maximum_endo_lead);
M_.maximum_endo_lead = M_.maximum_endo_lag;
end
nj = olr_state_.old_M_.endo_nbr-n_inst;
offset_min = M_.maximum_endo_lag - olr_state_.old_M_.maximum_endo_lag;
offset_max = M_.maximum_endo_lead - olr_state_.old_M_.maximum_endo_lead;
newiy = zeros(2*M_.maximum_endo_lag+1,nj+olr_state_.old_M_.endo_nbr);
jacobia_ = jacobia_(1:nj,:);
for i=1:2*M_.maximum_endo_lag+1
if i > offset_min & i <= 2*M_.maximum_endo_lag+1-offset_max
[junk,k1,k2] = find(olr_state_.old_M_.lead_lag_incidence(i-offset_min,:));
if i == M_.maximum_endo_lag+1
jacobia1 = [jacobia1 [jacobia_(:,k2); 2*options_.olr_w]];
else
jacobia1 = [jacobia1 [jacobia_(:,k2); ...
zeros(olr_state_.old_M_.endo_nbr,length(k1))]];
end
newiy(i,k1) = ones(1,length(k1));
end
i1 = 2*M_.maximum_endo_lag+2-i;
if i1 <= 2*M_.maximum_endo_lag+1-offset_max & i1 > offset_min
[junk,k1,k2] = find(olr_state_.old_M_.lead_lag_incidence(i1-offset_min,:));
k3 = find(any(jacobia_(:,k2),2));
x = zeros(olr_state_.old_M_.endo_nbr,length(k3));
x(k1,:) = bet^(-i1+M_.maximum_endo_lag+1)*jacobia_(k3,k2)';
jacobia1 = [jacobia1 [zeros(nj,length(k3)); x]];
newiy(i,k3+olr_state_.old_M_.endo_nbr) = ones(1,length(k3));
end
end
jacobia1 = [jacobia1 [jacobia_(:,end-M_.exo_nbr+1:end); ...
zeros(olr_state_.old_M_.endo_nbr, M_.exo_nbr)]];
newiy = newiy';
newiy = find(newiy(:));
M_.lead_lag_incidence = zeros(M_.endo_nbr*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),1);
M_.lead_lag_incidence(newiy) = [1:length(newiy)]';
M_.lead_lag_incidence =reshape(M_.lead_lag_incidence,M_.endo_nbr,M_.maximum_endo_lag+M_.maximum_endo_lead+1)';
jacobia_ = jacobia1;
clear jacobia1
% computes steady state
resid = feval([M_.fname '_steady'],zeros(olr_state_.old_M_.endo_nbr,1));
if resid'*resid < 1e-12
dr.ys =[dr.ys; zeros(nj,1)];
else
AA = zeros(M_.endo_nbr,M_.endo_nbr);
for i=1:M_.maximum_endo_lag+M_.maximum_endo_lead+1
[junk,k1,k2] = find(M_.lead_lag_incidence(i,:));
AA(:,k1) = AA(:,k1)+jacobia_(:,k2);
end
dr.ys = -AA\[resid; zeros(nj,1)];
oo_.exo_det_simul]);
end
oo_.exo_simul = tempex ;
tempex = [];
end
% end of code section for Optimal Linear Regulator
klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
dr=set_state_space(dr);
kstate = dr.kstate;
kad = dr.kad;
......
matlab/dynare_m.exe
View file @ 927da0e4
No preview for this file type
matlab/ramsey_dynamic.m 0 → 100644
View file @ 927da0e4
function J = ramsey_dynamic(ys,lbar)
% ramsey_dynamic sets up the Jacobian of the expanded model for optimal
% policies. It modifies several fields of M_
% INPUT:
% ys steady state of original endogenous variables
% lbar steady state of Lagrange multipliers
% OUPTUT:
% J jaocobian of expanded model
%
global M_ options_ it_
% retrieving model parameters
endo_nbr = M_.endo_nbr;
i_endo_nbr = 1:endo_nbr;
endo_names = M_.endo_names;
% exo_nbr = M_.exo_nbr+M_.exo_det_nbr;
% exo_names = vertcat(M_.exo_names,M_.exo_det_names);
exo_nbr = M_.exo_nbr;
exo_names = M_.exo_names;
i_leadlag = M_.lead_lag_incidence;
max_lead = M_.maximum_lead;
max_endo_lead = M_.maximum_endo_lead;
max_lag = M_.maximum_lag;
max_endo_lag = M_.maximum_endo_lag;
leadlag_nbr = max_lead+max_lag+1;
fname = M_.fname;
instr_names = options_.olr_inst;
instr_nbr = size(options_.olr_inst,1);
mult_nbr = endo_nbr-instr_nbr;
% discount factor
beta = options_.planner_discount;
% storing original values
orig_model.endo_nbr = endo_nbr;
orig_model.endo_names = endo_names;
orig_model.lead_lag_incidence = i_leadlag;
orig_model.maximum_lead = max_lead;
orig_model.maximum_endo_lead = max_endo_lead;
orig_model.maximum_lag = max_lag;
orig_model.maximum_endo_lag = max_endo_lag;
y = repmat(ys,1,max_lag+max_lead+1);
k = find(i_leadlag');
it_ = 1;
% retrieving derivatives of the objective function
[U,Uy,Uyy] = feval([fname '_objective'],ys,zeros(1,exo_nbr));
% retrieving derivatives of original model
[f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr));
% parameters for expanded model
endo_nbr1 = 2*endo_nbr-instr_nbr+exo_nbr;
max_lead1 = max_lead + max_lag;
max_lag1 = max_lead1;
max_leadlag1 = max_lead1;
% adding new variables names
endo_names1 = endo_names;
% adding shocks to endogenous variables
endo_names1 = strvcat(endo_names1, exo_names);
% adding multipliers names
for i=1:mult_nbr;
endo_names1 = strvcat(endo_names1,['mult_' int2str(i)]);
end
% expanding matrix of lead/lag incidence
%
% multipliers lead/lag incidence
i_mult = [];
for i=1:leadlag_nbr
i_mult = [any(fJ(:,nonzeros(i_leadlag(i,:))) ~= 0,2)' ; i_mult];
end
% putting it all together:
% original variables, exogenous variables made endogenous, multipliers
%
% number of original dynamic variables
n = nnz(i_leadlag);
% numbering columns of dynamic multipliers to be put in the last columns
% of the new Jacobian
i_leadlag1 = [cumsum(i_leadlag(1:max_lag,:),1); ...
repmat(i_leadlag(max_lag+1,:),leadlag_nbr,1); ...
flipud(cumsum(flipud(i_leadlag(max_lag+2:end,:)),1))];
i_leadlag1 = i_leadlag1';
k = find(i_leadlag1 > 0);
n = length(k);
i_leadlag1(k) = 1:n;
i_leadlag1 = i_leadlag1';
i_mult = i_mult';
k = find(i_mult > 0);
i_mult(k) = n+leadlag_nbr*exo_nbr+(1:length(k));
i_mult = i_mult';
i_leadlag1 = [ i_leadlag1 ...
[zeros(max_lead,exo_nbr);...
reshape(n+(1:leadlag_nbr*exo_nbr),exo_nbr,leadlag_nbr)'; ...
zeros(max_lag,exo_nbr)] ...
[zeros(max_lead,mult_nbr);...
i_mult;...
zeros(max_lag,mult_nbr)]];
i_leadlag1 = i_leadlag1';
k = find(i_leadlag1 > 0);
n = length(k);
i_leadlag1(k) = 1:n;
i_leadlag1 = i_leadlag1';
% building Jacobian of expanded model
jacobian = zeros(endo_nbr+mult_nbr,nnz(i_leadlag1)+exo_nbr);
% derivatives of f.o.c. w.r. to endogenous variables
% to be rearranged further down
lbarfH = lbar'*fH;
% indices of Hessian columns
n1 = nnz(i_leadlag)+exo_nbr;
iH = reshape(1:n1^2,n1,n1);
J = zeros(endo_nbr1,nnz(i_leadlag1)+exo_nbr);
% second order derivatives of objective function
J(1:endo_nbr,i_leadlag1(max_leadlag1+1,1:endo_nbr)) = Uyy;
% loop on lead/lags in expanded model
for i=1:2*max_leadlag1 + 1
% index of variables at the current lag in expanded model
kc = find(i_leadlag1(i,i_endo_nbr) > 0);
t1 = max(1,i-max_leadlag1);
t2 = min(i,max_leadlag1+1);
% loop on lead/lag blocks of relevant 1st order derivatives
for j = t1:t2
% derivatives w.r. endogenous variables
ic = find(i_leadlag(i-j+1,:) > 0 );
kc1 = i_leadlag(i-j+1,ic);
[junk,ic1,ic2] = intersect(ic,kc);
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)...
*reshape(lbarfH(iH(kr1,kc1)),length(kr1),length(kc1));
end
end
% derivatives w.r. aux. variables for lead/lag exogenous shocks
for i=1:leadlag_nbr
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)...
*reshape(lbarfH(iH(kr1,endo_nbr+(1:exo_nbr))),length(kr1), ...
exo_nbr);
end
% derivatives w.r. Lagrange multipliers
for i=1:leadlag_nbr
ic1 = find(i_leadlag(leadlag_nbr+1-i,:) > 0);
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)';
end
% Jacobian of original equations
%
% w.r. endogenous variables
ir = endo_nbr+(1:endo_nbr-instr_nbr);
for i=1:leadlag_nbr
ic1 = find(i_leadlag(i,:) > 0);
kc1 = i_leadlag(i,ic1);
ic2 = find(i_leadlag1(max_lag+i,:) > 0);
kc2 = i_leadlag1(max_lag+i,ic2);
[junk,junk,ic3] = intersect(ic1,ic2);
J(ir,kc2(ic3)) = fJ(:,kc1);
end
% w.r. exogenous variables
J(ir,nnz(i_leadlag1)+(1:exo_nbr)) = fJ(:,nnz(i_leadlag)+(1:exo_nbr));
% auxiliary variable for exogenous shocks
ir = 2*endo_nbr-instr_nbr+(1:exo_nbr);
kc = i_leadlag1(leadlag_nbr,endo_nbr+(1:exo_nbr));
J(ir,kc) = eye(exo_nbr);
J(ir,nnz(i_leadlag1)+(1:exo_nbr)) = -eye(exo_nbr);
% eliminating empty columns
% getting indices of nonzero entries
m = find(i_leadlag1');
n = length(m);
k = 1:size(J,2);
for i=1:n
if sum(abs(J(:,i))) < 1e-8
k(i) = 0;
m(i) = 0;
end
end
J = J(:,nonzeros(k));
i_leadlag1 = zeros(size(i_leadlag1))';
i_leadlag1(nonzeros(m)) = 1:nnz(m);
i_leadlag1 = i_leadlag1';
% setting expanded model parameters
% storing original values
M_.endo_nbr = endo_nbr1;
M_.endo_names = endo_names1;
M_.lead_lag_incidence = i_leadlag1;
M_.maximum_lead = max_lead1;
M_.maximum_endo_lead = max_lead1;
M_.maximum_lag = max_lag1;
M_.maximum_endo_lag = max_lag1;
M_.orig_model = orig_model;
matlab/ramsey_static.m 0 → 100644
View file @ 927da0e4
function [resids, rJ,mult] = ramsey_static(x)
% computes the static first order conditions for optimal policy
global M_ options_ it_
% recovering usefull fields
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
fname = M_.fname;
inst_nbr = M_.inst_nbr;
% i_endo_no_inst = M_.endogenous_variables_without_instruments;
max_lead = M_.maximum_lead;
max_lag = M_.maximum_lag;
beta = options_.planner_discount;
% indices of all endogenous variables
i_endo = [1:endo_nbr]';
% indices of endogenous variable except instruments
% i_inst = M_.instruments;
% indices of Lagrange multipliers
i_mult = [endo_nbr+1:2*endo_nbr-inst_nbr]';
% lead_lag incidence matrix for endogenous variables
i_lag = M_.lead_lag_incidence;
% value and Jacobian of objective function
ex = zeros(1,M_.exo_nbr);
[U,Uy,Uyy] = feval([fname '_objective'],x(i_endo),ex);
% value and Jacobian of dynamic function
y = repmat(x(i_endo),1,max_lag+max_lead+1);
k = find(i_lag');
it_ = 1;
% [f,fJ,fH] = feval([fname '_dynamic'],y(k),ex);
[f,fJ] = feval([fname '_dynamic'],y(k),ex);
% derivatives of Lagrangian with respect to endogenous variables
% res1 = Uy;
A = zeros(endo_nbr,endo_nbr-inst_nbr);
for i=1:max_lag+max_lead+1
% select variables present in the model at a given lag
[junk,k1,k2] = find(i_lag(i,:));
% res1(k1) = res1(k1) + beta^(max_lag-i+1)*fJ(:,k2)'*x(i_mult);
A(k1,:) = A(k1,:) + beta^(max_lag-i+1)*fJ(:,k2)';
end
i_inst = var_index(options_.olr_inst);
k = setdiff(1:size(A,1),i_inst);
mult = -A(k,:)\Uy(k);
resids = [f; Uy(i_inst)+A(i_inst,:)*mult];
rJ = [];
return;
% Jacobian of first order conditions
n = nnz(i_lag)+exo_nbr;
iH = reshape(1:n^2,n,n);
rJ = zeros(2*endo_nbr-inst_nbr,2*endo_nbr-inst_nbr);
rJ(i_endo,i_endo) = Uyy;
for i=1:max_lag+max_lead+1
% select variables present in the model at a given lag
[junk,k1,k2] = find(i_lag(i,:));
k3 = length(k2);
rJ(k1,k1) = rJ(k1,k1) + beta^(max_lag-i+1)*reshape(fH(:,iH(k2,k2))'*x(i_mult),k3,k3);
rJ(k1,i_mult) = rJ(k1,i_mult) + beta^(max_lag-1+1)*fJ(:,k2)';
rJ(i_mult,k1) = rJ(i_mult,k1) + fJ(:,k2);
end
% rJ = 1e-3*rJ;
% rJ(209,210) = rJ(209,210)+1-1e-3;
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment