Commit e029d466 authored by MichelJuillard's avatar MichelJuillard
Browse files

adding discretionary policy for linear quadratic models, thanks to code provided by Junior

parent 7e771cee
function info = discretionary_policy(var_list)
% Copyright (C) 2007-2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ oo_ M_
options_.discretionary_policy = 1;
oldoptions = options_;
options_.order = 1;
info = stoch_simul(var_list);
if options_.noprint == 0
disp_steady_state(M_,oo_)
for i=M_.orig_endo_nbr:M_.endo_nbr
if strmatch('mult_',M_.endo_names(i,:))
disp(sprintf('%s \t\t %g',M_.endo_names(i,:), ...
oo_.dr.ys(i)));
end
end
end
%oo_ = evaluate_planner_objective(oo_.dr,M_,oo_,options_);
options_ = oldoptions;
\ No newline at end of file
function [dr,ys,info]=discretionary_policy_1(oo_,Instruments)
% Copyright (C) 2007-2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_
persistent Hold
options_ = set_default_option(options_,'noprint',0);
options_ = set_default_option(options_,'qz_criterium',1.000001);
options_ = set_default_option(options_,'solve_maxit',3000);
options_ = set_default_option(options_,'discretion_tol',1e-7);
% safeguard against issues like running ramsey policy first and then running discretion
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;
else
M_.orig_model = M_;
end
beta = options_.planner_discount;
exo_nbr = M_.exo_nbr;
if isfield(M_,'orig_model')
orig_model = M_.orig_model;
endo_nbr = orig_model.endo_nbr;
endo_names = orig_model.endo_names;
lead_lag_incidence = orig_model.lead_lag_incidence;
MaxLead = orig_model.maximum_lead;
MaxLag = orig_model.maximum_lag;
else
endo_names = M_.endo_names;
endo_nbr = M_.endo_nbr;
MaxLag=M_.maximum_lag;
MaxLead=M_.maximum_lead;
lead_lag_incidence = M_.lead_lag_incidence;
end
[U,Uy,W] = feval([M_.fname,'_objective_static'],zeros(endo_nbr,1),[], M_.params);
if any(any(Uy~=0))
error(['discretionary_policy: the objective function must have zero ' ...
'first order derivatives'])
end
W=reshape(W,endo_nbr,endo_nbr);
klen = MaxLag + MaxLead + 1;
iyv=lead_lag_incidence';
% Find the jacobian
z = repmat(zeros(endo_nbr,1),1,klen);
z = z(nonzeros(iyv)) ;
it_ = MaxLag + 1 ;
if exo_nbr == 0
oo_.exo_steady_state = [] ;
end
[junk,jacobia_] = feval([M_.fname '_dynamic'],z, [oo_.exo_simul ...
oo_.exo_det_simul], M_.params, zeros(endo_nbr,1), it_);
if any(junk~=0)
error(['discretionary_policy: the model must be written in deviation ' ...
'form and not have constant terms'])
end
eq_nbr= size(jacobia_,1);
instr_nbr=endo_nbr-eq_nbr;
instr_id=nan(instr_nbr,1);
for j=1:instr_nbr
vj=deblank(Instruments(j,:));
vj_id=strmatch(vj,endo_names,'exact');
if ~isempty(vj_id)
instr_id(j)=vj_id;
else
error([mfilename,':: instrument ',vj,' not found'])
end
end
Indices={'lag','0','lead'};
iter=1;
for j=1:numel(Indices)
eval(['A',Indices{j},'=zeros(eq_nbr,endo_nbr);'])
if strcmp(Indices{j},'0')||(strcmp(Indices{j},'lag') && MaxLag)||(strcmp(Indices{j},'lead') && MaxLead)
[junk,row,col]=find(lead_lag_incidence(iter,:));
eval(['A',Indices{j},'(:,row)=jacobia_(:,col);'])
iter=iter+1;
end
end
B=jacobia_(:,nnz(iyv)+1:end);
%%% MAIN ENGINE %%%
qz_criterium = options_.qz_criterium;
solve_maxit = options_.solve_maxit;
discretion_tol = options_.discretion_tol;
if ~isempty(Hold)
[H,G,info]=discretionary_policy_engine(Alag,A0,Alead,B,W,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,Hold);
else
[H,G,info]=discretionary_policy_engine(Alag,A0,Alead,B,W,instr_id,beta,solve_maxit,discretion_tol,qz_criterium);
end
if info
dr=[];
return
else
Hold=H;
% Hold=[]; use this line if persistent command is not used.
end
% update the following elements
LLI=lead_lag_incidence;
LLI(MaxLag,:)=any(H);
LLI=LLI';
tmp=find(LLI);
LLI(tmp)=1:numel(tmp);
M_.lead_lag_incidence = LLI';
% set the state
dr=oo_.dr;
dr.ys =zeros(endo_nbr,1);
dr=set_state_space(dr,M_);
order_var=dr.order_var;
T=H(order_var,order_var);
dr.ghu=G(order_var,:);
Selection=any(T);
dr.ghx=T(:,Selection);
ys=NondistortionarySteadyState(M_);
dr.ys=ys; % <--- dr.ys =zeros(NewEndo_nbr,1);
function ys=NondistortionarySteadyState(M_)
if exist([M_.fname,'_steadystate.m'],'file')
eval(['ys=',M_.fname,'_steadystate.m;'])
else
ys=zeros(M_.endo_nbr,1);
end
function [H,G,retcode]=discretionary_policy_engine(AAlag,AA0,AAlead,BB,bigw,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,H00,verbose)
% Solves the discretionary problem for a model of the form:
% AAlag*yy_{t-1}+AA0*yy_t+AAlead*yy_{t+1}+BB*e=0, with W the weight on the
% variables in vector y_t and instr_id is the location of the instruments
% in the yy_t vector.
% We use the Dennis algorithm and so we need to re-write the model in the
% form A0*y_t=A1*y_{t-1}+A2*y_{t+1}+A3*x_t+A4*x_{t+1}+A5*e_t, with W the
% weight on the y_t vector and Q the weight on the x_t vector of
% instruments.
% Copyright (C) 2007-2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<12
verbose=0;
if nargin<11
H00=[];
if nargin<10
qz_criterium=1.000001;
if nargin<9
discretion_tol=sqrt(eps);
if nargin<8
solve_maxit=3000;
if nargin<7
beta=.99;
if nargin<6
error([mfilename,':: Insufficient number of input arguments'])
elseif nargin>12
error([mfilename,':: Number of input arguments cannot exceed 12'])
end
end
end
end
end
end
end
[A0,A1,A2,A3,A4,A5,W,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id);
% aux is a logical index of the instruments which appear with lags in the
% model. Their location in the state vector is instr_id(aux)
% endo_augm_id is index (not logical) of locations of the augmented vector
% of non-instrumental variables
AuxiliaryVariables_nbr=sum(aux);
H0=zeros(endo_nbr+AuxiliaryVariables_nbr);
if ~isempty(H00)
H0(1:endo_nbr,1:endo_nbr)=H00;clear H00
end
H10=H0(endo_augm_id,endo_augm_id);
F10=H0(instr_id,endo_augm_id);
yy_nbr=numel(endo_augm_id);
[H1,F1,retcode]=Iterate(H10,F10);
switch retcode
case 2 % nan
retcode=3;
retcode(2)=10000;
if verbose
disp([mfilename,':: NAN elements in the solution'])
end
case 1% maxiter
if verbose
disp([mfilename,':: Maximum Number of Iterations reached'])
end
case 0
BadEig=max(abs(eig(H1)))>qz_criterium;
if BadEig
retcode=3;
retcode(2)=100*max(abs(eig(H1)));
if verbose
disp([mfilename,':: Some eigenvalues greater than qz_criterium, Model potentially unstable'])
end
end
end
if retcode(1)
H=[];
G=[];
else
F2=-(Q+A3DPD*A3)\(A3DPD*A5);
H2=Dinv*(A5+A3*F2);
H=zeros(endo_nbr+AuxiliaryVariables_nbr);
G=zeros(endo_nbr+AuxiliaryVariables_nbr,exo_nbr);
H(endo_augm_id,endo_augm_id)=H1;
H(instr_id,endo_augm_id)=F1;
G(endo_augm_id,:)=H2;
G(instr_id,:)=F2;
% Account for auxilliary variables
H(:,instr_id(aux))=H(:,end-(AuxiliaryVariables_nbr-1:-1:0));
H=H(1:endo_nbr,1:endo_nbr);
G=G(1:endo_nbr,:);
end
function [H1,F1,retcode]=Iterate(H10,F10)
iter=0;
H1=H10;
F1=F10;
while 1
iter=iter+1;
P=SylvesterDoubling(W+beta*F1'*Q*F1,beta*H1',H1,discretion_tol,solve_maxit);
if any(any(isnan(P)))
P=SylvesterHessenbergSchur(W+beta*F1'*Q*F1,beta*H1',H1);
if any(any(isnan(P)))
retcode=2;
return
end
end
D=A0-A2*H1-A4*F1;
Dinv=D\eye(yy_nbr);
A3DPD=A3'*Dinv'*P*Dinv;
F1=-(Q+A3DPD*A3)\(A3DPD*A1);
H1=Dinv*(A1+A3*F1);
[rcode,NQ]=CheckConvergence([H1;F1]-[H10;F10],iter,solve_maxit,discretion_tol);
if rcode
break
else
if verbose
disp(NQ)
end
end
H10=H1;
F10=F1;
end
retcode=rcode-1;
end
end
function [rcode,NQ]=CheckConvergence(Q,iter,MaxIter,crit)
NQ=max(max(abs(Q)));% norm(Q); seems too costly
if isnan(NQ)
rcode=3;
elseif iter>MaxIter;
rcode=2;
elseif NQ<crit
rcode=1;
else
rcode=0;
end
end
function [A00,A11,A22,A33,A44,A55,WW,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id)
[eq_nbr,endo_nbr]=size(AAlag);
exo_nbr=size(BB,2);
y=setdiff(1:endo_nbr,instr_id);
instr_nbr=numel(instr_id);
A0=AA0(:,y);
A1=-AAlag(:,y);
A2=-AAlead(:,y);
A3=-AA0(:,instr_id);
A4=-AAlead(:,instr_id);
A5=-BB;
W=bigw(y,y);
Q=bigw(instr_id,instr_id);
% Adjust for possible lags in instruments by creating auxiliary equations
A6=-AAlag(:,instr_id);
aux=any(A6);
AuxiliaryVariables_nbr=sum(aux);
ny=eq_nbr;
m=eq_nbr+AuxiliaryVariables_nbr;
A00=zeros(m);A00(1:ny,1:ny)=A0;A00(ny+1:end,ny+1:end)=eye(AuxiliaryVariables_nbr);
A11=zeros(m);A11(1:ny,1:ny)=A1;A11(1:ny,ny+1:end)=A6(:,aux);
A22=zeros(m);A22(1:ny,1:ny)=A2;
A33=zeros(m,instr_nbr);A33(1:ny,1:end)=A3;A33(ny+1:end,aux)=eye(AuxiliaryVariables_nbr);
A44=zeros(m,instr_nbr);A44(1:ny,1:end)=A4;
A55=zeros(m,exo_nbr);A55(1:ny,1:end)=A5;
WW=zeros(m);WW(1:ny,1:ny)=W;
endo_augm_id=setdiff(1:endo_nbr+AuxiliaryVariables_nbr,instr_id);
end
function v= SylvesterDoubling (d,g,h,tol,maxit)
% DOUBLES Solves a Sylvester equation using doubling
%
% [v,info] = doubles (g,d,h,tol,maxit) uses a doubling algorithm
% to solve the Sylvester equation v = d + g v h
v = d;
for i =1:maxit,
vadd = g*v*h;
v = v+vadd;
if norm (vadd,1) <= (tol*norm(v,1))
break;
end
g = g*g;
h = h*h;
end
end
function v = SylvesterHessenbergSchur(d,g,h)
%
% DSYLHS Solves a discrete time sylvester equation using the
% Hessenberg-Schur algorithm
%
% v = DSYLHS(g,d,h) computes the matrix v that satisfies the
% discrete time sylvester equation
%
% v = d + g'vh
if size(g,1) >= size(h,1)
[u,gbarp] = hess(g');
[t,hbar] = schur(h);
[vbar] = sylvest_private(gbarp,u'*d*t,hbar,1e-15);
v = u*vbar*t';
else
[u,gbar] = schur(g);
[t,hbarp] = hess(h');
[vbar] = sylvest_private(hbarp,t'*d'*u,gbar,1e-15);
v = u*vbar'*t';
end
end
function v = sylvest_private(g,d,h,tol)
%
% SYLVEST Solves a Sylvester equation
%
% solves the Sylvester equation
% v = d + g v h
% for v where both g and h must be upper block triangular.
% The output info is zero on a successful return.
% The input tol indicates when an element of g or h should be considered
% zero.
[m,n] = size(d);
v = zeros(m,n);
w = eye(m);
i = 1;
temp = [];
%First handle the i = 1 case outside the loop
if i< n,
if abs(h(i+1,i)) < tol,
v(:,i)= (w - g*h(i,i))\d(:,i);
i = i+1;
else
A = [w-g*h(i,i) (-g*h(i+1,i));...
-g*h(i,i+1) w-g*h(i+1,i+1)];
C = [d(:,i); d(:,i+1)];
X = A\C;
v(:,i) = X(1:m,:);
v(:,i+1) = X(m+1:2*m, :);
i = i+2;
end
end
%Handle the rest of the matrix with the possible exception of i=n
while i<n,
b= i-1;
temp = [temp g*v(:,size(temp,2)+1:b)]; %#ok<AGROW>
if abs(h(i+1,i)) < tol,
v(:,i) = (w - g*h(i,i))\(d(:,i) + temp*h(1:b,i));
i = i+1;
else
A = [w - g*h(i,i) (-g*h(i+1,i)); ...
-g*h(i,i+1) w - g*h(i+1,i+1)];
C = [d(:,i) + temp*h(1:b,i); ...
d(:,i+1) + temp*h(1:b,i+1)];
X = A\C;
v(:,i) = X(1:m,:);
v(:,i+1) = X(m+1:2*m, :);
i = i+2;
end
end
%Handle the i = n case if i=n was not in a 2-2 block
if i==n,
b = i-1;
temp = [temp g*v(:,size(temp,2)+1:b)];
v(:,i) = (w-g*h(i,i))\(d(:,i) + temp*h(1:b,i));
end
end
\ No newline at end of file
......@@ -249,6 +249,9 @@ options_.initial_date.freq = 1;
options_.initial_date.period = 1;
options_.initial_date.subperiod = 0;
% discretionary policy
options_.discretionary_policy = 0;
% SWZ SBVAR
options_.ms.freq = 1;
options_.ms.initial_subperiod = 1;
......
......@@ -62,6 +62,12 @@ check_model;
if PI_PCL_solver
[oo_.dr, info] = PCL_resol(oo_.steady_state,0);
elseif options_.discretionary_policy
if ~options_.linear
error(['discretionary_policy solves only linear_quadratic ' ...
'problems']);
end
[oo_.dr,ys,info] = discretionary_policy_1(oo_,options_.instruments);
else
[oo_.dr, info] = resol(oo_.steady_state,0);
end
......
......@@ -216,6 +216,53 @@ RamseyPolicyStatement::writeOutput(ostream &output, const string &basename) cons
output << "ramsey_policy(var_list_);\n";
}
DiscretionaryPolicyStatement::DiscretionaryPolicyStatement(const SymbolList &symbol_list_arg,
const OptionsList &options_list_arg) :
symbol_list(symbol_list_arg),
options_list(options_list_arg)
{
}
void
DiscretionaryPolicyStatement::checkPass(ModFileStructure &mod_file_struct)
{
mod_file_struct.discretionary_policy_present = true;
/* Fill in option_order of mod_file_struct
Since discretionary policy needs one further order of derivation (for example, for 1st order
approximation, it needs 2nd derivatives), we add 1 to the order declared by user */
OptionsList::num_options_t::const_iterator it = options_list.num_options.find("order");
if (it != options_list.num_options.end())
{
int order = atoi(it->second.c_str());
if (order > 1)
{
cerr << "ERROR: discretionary_policy: order > 1 is not yet implemented" << endl;
exit(EXIT_FAILURE);
}
mod_file_struct.order_option = max(mod_file_struct.order_option, order + 1);
}
// Fill in mod_file_struct.partial_information
it = options_list.num_options.find("partial_information");
if (it != options_list.num_options.end() && it->second == "1")
mod_file_struct.partial_information = true;
// Option k_order_solver (implicit when order >= 3)
it = options_list.num_options.find("k_order_solver");
if ((it != options_list.num_options.end() && it->second == "1")
|| mod_file_struct.order_option >= 3)
mod_file_struct.k_order_solver = true;
}
void
DiscretionaryPolicyStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
symbol_list.writeOutput("var_list_", output);
output << "discretionary_policy(var_list_);\n";
}
EstimationStatement::EstimationStatement(const SymbolList &symbol_list_arg,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg) :
......
......@@ -103,6 +103,18 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const;
};
class DiscretionaryPolicyStatement : public Statement
{
private:
const SymbolList symbol_list;
const OptionsList options_list;
public:
DiscretionaryPolicyStatement(const SymbolList &symbol_list_arg,
const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class RplotStatement : public Statement