Commit 77cdcce6 authored by michel's avatar michel
Browse files

v4: adding Ramsey policy to parser;

changes to Ramsey policy Matlab code

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1133 ac1d8469-bf42-47a9-8791-bf33cf982152
parent e3bb7606
......@@ -54,7 +54,8 @@ while ~done
grad = (feval(FUN,x*ones(1,nv)+tvec,varargin{:})-f0*ones(1,nv))/delta;
end
else % use analytic gradient
grad=feval(gradfun,x,varargin{:});
% grad=feval(gradfun,x,varargin{:});
[f0,grad] = feval(FUN,x,varargin{:});
end
if isreal(grad)
if rcond(grad)<1e-12
......
......@@ -51,7 +51,7 @@ if M_.exo_nbr == 0
end
% expanding system for Optimal Linear Regulator
if options_.olr
if options_.ramsey_policy
if isfield(M_,'orig_model')
orig_model = M_.orig_model;
M_.endo_nbr = orig_model.endo_nbr;
......@@ -62,7 +62,10 @@ if options_.olr
M_.maximum_lag = orig_model.maximum_lag;
M_.maximum_endo_lag = orig_model.maximum_endo_lag;
end
old_solve_algo = options_.solve_algo;
options_.solve_algo = 1;
oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,0);
options_.solve_algo = old_solve_algo;
[junk,junk,multbar] = ramsey_static(oo_.steady_state);
jacobia_=ramsey_dynamic(oo_.steady_state,multbar);
klen = M_.maximum_lag + M_.maximum_lead + 1;
......
No preview for this file type
......@@ -75,8 +75,10 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
error(sprintf('Solve block = %d check = %d\n',i,info));
end
end
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
fvec = feval(func,x,varargin{:});
if max(abs(fvec)) > tolf
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
end
elseif options_.solve_algo == 3
[x,info] = csolve(func,x,'grad_ss',1e-6,500,varargin{:});
end
......
......@@ -54,6 +54,13 @@ function global_initialization()
% Solution
options_.order = 2;
options_.dr_algo = 0;
options_.solve_algo = 2;
% Ramsey policy
options_.planner_discount = 1.0;
% estimation
options_.load_mh_file = 0;
% Misc
options_.conf_sig = 0.9;
......
......@@ -229,7 +229,7 @@ elseif options_.load_mh_file == 1
save([MhDirectoryName '/' M_.fname '_mh_history'],'record');
disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
disp(' ')
elseif options_.load_mh_file == -1
elseif options_.mh_recover == 1
%% The previous metropolis-hastings crashed before the end! I try to
%% recover the saved draws...
disp('MH: Recover mode!')
......
function dr=mult_elimination(void)
% function mult_elimination()
% replaces Lagrange multipliers in Ramsey policy by lagged value of state
% and shock variables
% INPUT
% none
% OUTPUT
% dr: a structure with the new decision rule
%
global M_ options_ oo_
dr = oo_.dr;
nstatic = dr.nstatic;
npred = dr.npred;
order_var = dr.order_var;
nstates = M_.endo_names(order_var(nstatic+(1:npred)),:);
il = strmatch('mult_',nstates);
nil = setdiff(1:dr.npred,il);
m_nbr = length(il);
nm_nbr = length(nil);
AA1 = dr.ghx(:,nil);
AA2 = dr.ghx(:,il);
A1 = dr.ghx(nstatic+(1:npred),nil);
A2 = dr.ghx(nstatic+(1:npred),il);
B = dr.ghu(nstatic+(1:npred),:);
A11 = A1(nil,:);
A21 = A1(il,:);
A12 = A2(nil,:);
A22 = A2(il,:)
[Q1,R1,E1] = qr(A2);
n1 = sum(abs(diag(R1)) > 1e-8);
Q1_12 = Q1(1:nm_nbr,n1+1:end);
Q1_22 = Q1(nm_nbr+1:end,n1+1:end);
[Q2,R2,E2] = qr(Q1_22');
n2 = sum(abs(diag(R2)) > 1e-8);
R2_1 = inv(R2(1:n2,1:n2));
M1(order_var,:) = AA1 - AA2*E2*[R2_1*Q2(:,1:n2)'*Q1_12'; zeros(m_nbr-n2,m_nbr)];
M2(order_var,:) = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*A1; zeros(m_nbr-n2,nil)];
M3(order_var,:) = dr.ghu;
M4(order_var,:) = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*B; zeros(m_nbr-n2,size(B,2))];
endo_nbr = M_.orig_model.endo_nbr;
exo_nbr = M_.exo_nbr;
lead_lag_incidence = M_.lead_lag_incidence(:,1:endo_nbr+exo_nbr);
lead_lag_incidence1 = lead_lag_incidence(1,:) > 0;
maximum_lag = M_.maximum_lag;
for i=1:maximum_lag-1
lead_lag_incidence1 = [lead_lag_incidence1; lead_lag_incidence(i,:)| ...
lead_lag_incidence(i+1,:)];
end
lead_lag_incidence1 = [lead_lag_incidence1; ...
lead_lag_incidence(M_.maximum_lag,:) > 0];
k = find(lead_lag_incidence1');
lead_lag_incidence1 = zeros(size(lead_lag_incidence1'));
lead_lag_incidence1(k) = 1:length(k);
lead_lag_incidence1 = lead_lag_incidence1';
kstate = zeros(0,2);
for i=maximum_lag:-1:1
k = find(lead_lag_incidence(i,:));
kstate = [kstate; [k repmat(i+1,length(k),1)]];
end
dr.M1 = M1;
dr.M2 = M2;
dr.M3 = M3;
dr.M4 = M4;
\ No newline at end of file
......@@ -7,6 +7,7 @@ function J = ramsey_dynamic(ys,lbar)
% OUPTUT:
% J jaocobian of expanded model
%
global M_ options_ it_
% retrieving model parameters
......@@ -24,9 +25,8 @@ function J = ramsey_dynamic(ys,lbar)
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;
% instr_names = options_.olr_inst;
% instr_nbr = size(options_.olr_inst,1);
% discount factor
beta = options_.planner_discount;
......@@ -44,9 +44,14 @@ function J = ramsey_dynamic(ys,lbar)
k = find(i_leadlag');
it_ = 1;
% retrieving derivatives of the objective function
[U,Uy,Uyy] = feval([fname '_objective'],ys,zeros(1,exo_nbr));
[U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr));
Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% retrieving derivatives of original model
[f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr));
instr_nbr = endo_nbr - size(f,1);
mult_nbr = endo_nbr-instr_nbr;
% parameters for expanded model
endo_nbr1 = 2*endo_nbr-instr_nbr+exo_nbr;
......
......@@ -7,8 +7,8 @@ function [resids, rJ,mult] = ramsey_static(x)
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;
% 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;
......@@ -16,22 +16,25 @@ function [resids, rJ,mult] = ramsey_static(x)
% 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]';
% i_inst = M_.instruments;
% 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);
[U,Uy,Uyy] = feval([fname '_objective_static'],x(i_endo),ex);
Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% 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);
% indices of Lagrange multipliers
inst_nbr = endo_nbr - size(f,1);
i_mult = [endo_nbr+1:2*endo_nbr-inst_nbr]';
% derivatives of Lagrangian with respect to endogenous variables
% res1 = Uy;
......@@ -43,10 +46,13 @@ function [resids, rJ,mult] = ramsey_static(x)
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];
% i_inst = var_index(options_.olr_inst);
% k = setdiff(1:size(A,1),i_inst);
% mult = -A(k,:)\Uy(k);
mult = -A\Uy;
% resids = [f; Uy(i_inst)+A(i_inst,:)*mult];
resids1 = Uy+A*mult;
resids = [f; sqrt(resids1'*resids1/endo_nbr)];
rJ = [];
return;
......
......@@ -12,7 +12,6 @@ global it_
%unfinished
jacobian_flag = 0;
options_ = set_default_option(options_,'olr',0);
options_ = set_default_option(options_,'jacobian_flag',1);
info = 0;
......@@ -32,7 +31,7 @@ end
dr.ys = ys;
fh = str2func([M_.fname '_static']);
if options_.linear == 0
if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state]))) > options_.dynatol & options_.olr == 0
if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state]))) > options_.dynatol & options_.ramsey_policy == 0
if options_.steadystate_flag
[dr.ys,check1] = feval([M_.fname '_steadystate'],dr.ys,...
[oo_.exo_steady_state; oo_.exo_det_steady_state]);
......
......@@ -81,6 +81,27 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename) const
output << "stoch_simul(var_list_);\n";
}
RamseyPolicyStatement::RamseyPolicyStatement(const TmpSymbolTable &tmp_symbol_table_arg,
const OptionsList &options_list_arg) :
tmp_symbol_table(tmp_symbol_table_arg),
options_list(options_list_arg)
{
}
void
RamseyPolicyStatement::checkPass(ModFileStructure &mod_file_struct)
{
mod_file_struct.stoch_simul_or_similar_present = true;
}
void
RamseyPolicyStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
tmp_symbol_table.writeOutput("var_list_", output);
output << "options_.ramsey_policy=1;\nstoch_simul(var_list_);\n";
}
EstimationStatement::EstimationStatement(const TmpSymbolTable &tmp_symbol_table_arg,
const OptionsList &options_list_arg) :
tmp_symbol_table(tmp_symbol_table_arg),
......
......@@ -63,7 +63,7 @@ typedef pair<int, Type> ExpObj;
%token QZ_CRITERIUM
%token RELATIVE_IRF REPLIC RPLOT
%token SHOCKS SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER SOLVE_ALGO STDERR STEADY STOCH_SIMUL
%token TEX
%token TEX RAMSEY_POLICY PLANNER_DISCOUNT
%token <string_val> TEX_NAME
%token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL
%token VALUES VAR VAREXO VAREXO_DET VAROBS
......@@ -129,7 +129,8 @@ typedef pair<int, Type> ExpObj;
| olr
| olr_inst
| model_comparison
| planner_objective
| planner_objective
| ramsey_policy
;
......@@ -1039,8 +1040,6 @@ typedef pair<int, Type> ExpObj;
| o_noprint
;
planner_objective : PLANNER_OBJECTIVE { driver.begin_planner_objective(); } hand_side { driver.end_planner_objective($3); } ';'
filename_list : filename {driver.add_mc_filename($1);}
| filename_list COMMA filename {driver.add_mc_filename($3);}
| filename '(' value ')' {driver.add_mc_filename($1, $3);}
......@@ -1058,6 +1057,28 @@ typedef pair<int, Type> ExpObj;
| '.' { $$ = new string("."); }
;
planner_objective : PLANNER_OBJECTIVE { driver.begin_planner_objective(); } hand_side { driver.end_planner_objective($3); } ';'
;
ramsey_policy : RAMSEY_POLICY ';'
{driver.ramsey_policy();}
| RAMSEY_POLICY '(' ramsey_policy_options_list ')' ';'
{driver.ramsey_policy();}
| RAMSEY_POLICY tmp_var_list ';'
{driver.ramsey_policy();}
| RAMSEY_POLICY '(' ramsey_policy_options_list ')' tmp_var_list ';'
{driver.ramsey_policy();}
;
ramsey_policy_options_list :
ramsey_policy_options_list COMMA ramsey_policy_options
| ramsey_policy_options
;
ramsey_policy_options : stoch_simul_options
| o_planner_discount
;
o_dr_algo: DR_ALGO EQUAL INT_NUMBER {driver.option_num("dr_algo", $3);};
o_solve_algo: SOLVE_ALGO EQUAL INT_NUMBER {driver.option_num("solve_algo", $3);};
o_simul_algo: SIMUL_ALGO EQUAL INT_NUMBER {driver.option_num("simul_algo", $3);};
......@@ -1129,6 +1150,7 @@ typedef pair<int, Type> ExpObj;
o_constant : CONSTANT {driver.option_num("noconstant", "0");}
o_noconstant : NOCONSTANT {driver.option_num("noconstant", "1");}
o_mh_recover : MH_RECOVER {driver.option_num("load_mh_file", "-1");}
o_planner_discount : PLANNER_DISCOUNT EQUAL FLOAT_NUMBER {driver.option_num("planner_discount",$3);}
range : NAME ':' NAME
{
......@@ -1137,6 +1159,7 @@ typedef pair<int, Type> ExpObj;
delete $3;
$$ = $1;
}
vec_int_elem : INT_NUMBER
| INT_NUMBER ':' INT_NUMBER
{ $1->append(":"); $1->append(*$3); delete $3; $$ = $1; }
......
......@@ -79,6 +79,7 @@ int sigma_e = 0;
<INITIAL>Sigma_e {BEGIN DYNARE_STATEMENT; sigma_e = 1; return token::SIGMA_E;}
<INITIAL>calib {BEGIN DYNARE_STATEMENT; return token::CALIB;}
<INITIAL>planner_objective {BEGIN DYNARE_STATEMENT; return token::PLANNER_OBJECTIVE;}
<INITIAL>ramsey_policy {BEGIN DYNARE_STATEMENT; return token::RAMSEY_POLICY;}
/* End of a Dynare statement */
<DYNARE_STATEMENT>; {
......@@ -199,6 +200,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT>xls_sheet {return token::XLS_SHEET;}
<DYNARE_STATEMENT>xls_range {return token::XLS_RANGE;}
<DYNARE_STATEMENT>mh_recover {return token::MH_RECOVER;}
<DYNARE_STATEMENT>planner_discount {return token::PLANNER_DISCOUNT;}
<DYNARE_STATEMENT>[\.] {return yy::parser::token_type (yytext[0]);}
......
......@@ -926,6 +926,14 @@ ParsingDriver::end_planner_objective(NodeID expr)
mod_file->addStatement(new PlannerObjectiveStatement(model_tree));
}
void
ParsingDriver::ramsey_policy()
{
mod_file->addStatement(new RamseyPolicyStatement(*tmp_symbol_table, options_list));
tmp_symbol_table->clear();
options_list.clear();
}
NodeID
ParsingDriver::add_model_equal(NodeID arg1, NodeID arg2)
{
......
......@@ -49,6 +49,18 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const;
};
class RamseyPolicyStatement : public Statement
{
private:
const TmpSymbolTable tmp_symbol_table;
const OptionsList options_list;
public:
RamseyPolicyStatement(const TmpSymbolTable &tmp_symbol_table_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
{
private:
......
......@@ -289,6 +289,8 @@ public:
void begin_planner_objective();
//! End a planner objective statement
void end_planner_objective(NodeID expr);
//! ramsey policy statement
void ramsey_policy();
//! Writes token "arg1=arg2" to model tree
NodeID add_model_equal(NodeID arg1, NodeID arg2);
//! Writes token "arg=0" to model tree
......
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