From 336afa9f2ba76f52dba2b9d62572801e24e9a0f8 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 18 Dec 2011 17:33:08 +0100 Subject: [PATCH] fixing bugs related to new functions --- matlab/AIM_first_order_solver.m | 60 +++++++------- matlab/dyn_first_order_solver.m | 1 + matlab/dyn_ramsey_linearized_foc.m | 124 ----------------------------- matlab/osr1.m | 4 +- matlab/stochastic_solvers.m | 83 ++++++++++--------- 5 files changed, 77 insertions(+), 195 deletions(-) delete mode 100644 matlab/dyn_ramsey_linearized_foc.m diff --git a/matlab/AIM_first_order_solver.m b/matlab/AIM_first_order_solver.m index 962da3dee..82e6ddfbf 100644 --- a/matlab/AIM_first_order_solver.m +++ b/matlab/AIM_first_order_solver.m @@ -1,33 +1,31 @@ -function [dr,info]=AIM_first_order_solver(jacobia,M_,dr) - try - [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr); +function [dr,info]=AIM_first_order_solver(jacobia,M,dr,qz_criterium,nd) + + info = 0; + + [dr,aimcode]=dynAIMsolver1(jacobia,M,dr); - % reuse some of the bypassed code and tests that may be needed - - if aimcode ~=1 - info(1) = aimcode; - info(2) = 1.0e+8; - return - end - [A,B] =transition_matrix(dr); - dr.eigval = eig(A); - nba = nd-sdim; - - nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1); - - if nba ~= nyf - temp = sort(abs(dr.eigval)); - if nba > nyf - temp = temp(nd-nba+1:nd-nyf)-1-options_.qz_criterium; - info(1) = 3; - elseif nba < nyf; - temp = temp(nd-nyf+1:nd-nba)-1-options_.qz_criterium; - info(1) = 4; - end - info(2) = temp'*temp; - return - end - catch - disp(lasterror.message) - error('Problem with AIM solver - Try to remove the "aim_solver" option') + if aimcode ~=1 + info(1) = aimcode; + info(2) = 1.0e+8; + return end + [A,B] =transition_matrix(dr); + dr.eigval = eig(A); + nd = size(dr.kstate,1); + nba = nd-sum( abs(dr.eigval) < qz_criterium ); + + nyf = dr.nfwrd+dr.nboth; + + if nba ~= nyf + temp = sort(abs(dr.eigval)); + if nba > nyf + temp = temp(nd-nba+1:nd-nyf)-1-qz_criterium; + info(1) = 3; + elseif nba < nyf; + temp = temp(nd-nyf+1:nd-nba)-1-qz_criterium; + info(1) = 4; + end + info(2) = temp'*temp; + return + end + diff --git a/matlab/dyn_first_order_solver.m b/matlab/dyn_first_order_solver.m index 98a186a50..ef577cb98 100644 --- a/matlab/dyn_first_order_solver.m +++ b/matlab/dyn_first_order_solver.m @@ -130,6 +130,7 @@ function [dr,info] = dyn_first_order_solver(jacobia,b,M_,dr,options,task) k1 = find(kstate(n4:nd,2) == M_.maximum_endo_lag+1); k2 = find(kstate(1:n3,2) == M_.maximum_endo_lag+2); + dr.gx = gx; dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)]; %lead variables actually present in the model diff --git a/matlab/dyn_ramsey_linearized_foc.m b/matlab/dyn_ramsey_linearized_foc.m deleted file mode 100644 index 2a1949cba..000000000 --- a/matlab/dyn_ramsey_linearized_foc.m +++ /dev/null @@ -1,124 +0,0 @@ -function [jacobia_,dr,info,M_,oo_] = dyn_ramsey_linearized_foc(dr,M_,options_,oo_) -% function [jacobia_,dr,info,M_,oo_] = dyn_ramsey_linearized_foc(dr,M_,options_,oo_) -% computes the Jacobian of the linear approximation of the F.O.C of a -% Ramsey problem -% -% INPUTS -% dr [matlab structure] Decision rules for stochastic simulations. -% M_ [matlab structure] Definition of the model. -% options_ [matlab structure] Global options. -% oo_ [matlab structure] Results -% -% OUTPUTS -% dr [matlab structure] Decision rules for stochastic simulations. -% info [integer] info=1: the model doesn't define current variables uniquely -% info=2: problem in mjdgges.dll info(2) contains error code. -% info=3: BK order condition not satisfied info(2) contains "distance" -% absence of stable trajectory. -% info=4: BK order condition not satisfied info(2) contains "distance" -% indeterminacy. -% info=5: BK rank condition not satisfied. -% info=6: The jacobian matrix evaluated at the steady state is complex. -% M_ [matlab structure] -% options_ [matlab structure] -% oo_ [matlab structure] -% -% ALGORITHM -% ... -% -% SPECIAL REQUIREMENTS -% none. -% - -% Copyright (C) 1996-2009 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 . - -info = 0; - -if isfield(M_,'orig_model') - orig_model = M_.orig_model; - M_.endo_nbr = orig_model.endo_nbr; - M_.orig_endo_nbr = orig_model.orig_endo_nbr; - M_.aux_vars = orig_model.aux_vars; - 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 - - if options_.steadystate_flag - k_inst = []; - instruments = options_.instruments; - for i = 1:size(instruments,1) - k_inst = [k_inst; strmatch(options_.instruments(i,:), ... - M_.endo_names,'exact')]; - end - ys = oo_.steady_state; - [inst_val,info1] = dynare_solve('dyn_ramsey_static_', ... - oo_.steady_state(k_inst),0, ... - M_,options_,oo_,it_); - M_.params = evalin('base','M_.params;'); - ys(k_inst) = inst_val; - [x,check] = feval([M_.fname '_steadystate'],... - ys,[oo_.exo_steady_state; ... - oo_.exo_det_steady_state]); - if size(x,1) < M_.endo_nbr - if length(M_.aux_vars) > 0 - x = add_auxiliary_variables_to_steadystate(x,M_.aux_vars,... - M_.fname,... - oo_.exo_steady_state,... - oo_.exo_det_steady_state,... - M_.params); - else - error([M_.fname '_steadystate.m doesn''t match the model']); - end - end - oo_.steady_state = x; - [junk,junk,multbar] = dyn_ramsey_static_(oo_.steady_state(k_inst),M_,options_,oo_,it_); - else - [oo_.steady_state,info1] = dynare_solve('dyn_ramsey_static_', ... - oo_.steady_state,0,M_,options_,oo_,it_); - [junk,junk,multbar] = dyn_ramsey_static_(oo_.steady_state,M_,options_,oo_,it_); - end - - check1 = max(abs(feval([M_.fname '_static'],... - oo_.steady_state,... - [oo_.exo_steady_state; ... - oo_.exo_det_steady_state], M_.params))) > options_.dynatol ; - if check1 - info(1) = 20; - info(2) = check1'*check1; - return - end - - [jacobia_,M_] = dyn_ramsey_dynamic_(oo_.steady_state,multbar,M_,options_,oo_,it_); - klen = M_.maximum_lag + M_.maximum_lead + 1; - dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar]; - oo_.steady_state = dr.ys; - - 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,:), ... - dr.ys(i))); - end - end - end - diff --git a/matlab/osr1.m b/matlab/osr1.m index 478dff673..0f2079cd9 100644 --- a/matlab/osr1.m +++ b/matlab/osr1.m @@ -40,12 +40,12 @@ end exe =zeros(M_.exo_nbr,1); -dr = set_state_space(oo_.dr,M_); +oo_.dr = set_state_space(oo_.dr,M_); np = size(i_params,1); t0 = M_.params(i_params); -inv_order_var = dr.inv_order_var; +inv_order_var = oo_.dr.inv_order_var; H0 = 1e-4*eye(np); crit = 1e-7; diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m index 08f15e388..97fe44fe5 100644 --- a/matlab/stochastic_solvers.m +++ b/matlab/stochastic_solvers.m @@ -66,35 +66,42 @@ if options_.k_order_solver; return; end -if options_.ramsey_policy - % expanding system for Optimal Linear Regulator - [jacobia_,dr,info,M_,oo_] = dyn_ramsey_linearized_foc(dr,M_,options_,oo_); -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 - - it_ = M_.maximum_lag + 1; - z = repmat(dr.ys,1,klen); - z = z(iyr0) ; - if options_.order == 1 - [junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ... - oo_.exo_det_simul], M_.params, dr.ys, it_); - elseif options_.order == 2 - [junk,jacobia_,hessian1] = feval([M_.fname '_dynamic'],z,... - [oo_.exo_simul ... - oo_.exo_det_simul], M_.params, dr.ys, it_); - if options_.use_dll - % In USE_DLL mode, the hessian is in the 3-column sparse representation - hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ... - size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2)); - end +klen = M_.maximum_lag + M_.maximum_lead + 1; +exo_simul = [repmat(oo_.exo_steady_state',klen,1) repmat(oo_.exo_det_steady_state',klen,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 + +it_ = M_.maximum_lag + 1; +z = repmat(dr.ys,1,klen); +if options_.order == 1 + if (options_.bytecode) + [chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + M_.params, dr.ys, 1); + jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; + else + [junk,jacobia_] = feval([M_.fname '_dynamic'],z(iyr0),exo_simul, ... + M_.params, dr.ys, it_); + end; +elseif options_.order == 2 + if (options_.bytecode) + [chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + M_.params, dr.ys, 1); + jacobia_ = [loc_dr.g1 loc_dr.g1_x]; + else + [junk,jacobia_,hessian1] = feval([M_.fname '_dynamic'],z(iyr0),... + exo_simul, ... + M_.params, dr.ys, 3); + end; + if options_.use_dll + % In USE_DLL mode, the hessian is in the 3-column sparse representation + hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ... + size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2)); end end @@ -119,6 +126,7 @@ nstatic = dr.nstatic; nfwrd = dr.nfwrd; npred = dr.npred; nboth = dr.nboth; +nfwrds = nfwrd+nboth; order_var = dr.order_var; nd = size(kstate,1); nz = nnz(M_.lead_lag_incidence); @@ -135,10 +143,9 @@ if M_.maximum_endo_lead == 0 % If required, use AIM solver if not check only if options_.order == 1 [k1,junk,k2] = find(kstate(:,4)); - temp = -b\jacobia_(:,[k2 nz+1:end]); - dr.ghx = temp(:,1:npred); + dr.ghx(:,k1) = -b\jacobia_(:,k2); if M_.exo_nbr - dr.ghu = temp(:,npred+1:end); + dr.ghu = -b\jacobia_(:,nz+1:end); end dr.eigval = eig(transition_matrix(dr)); dr.rank = 0; @@ -154,16 +161,16 @@ if M_.maximum_endo_lead == 0 'backward models']) end elseif M_.maximum_endo_lag == 0 - % purely forward model - dr.ghx = []; - dr.ghu = -b\jacobia_(:,nz+1:end); +% purely forward model +dr.ghx = []; +dr.ghu = -b\jacobia_(:,nz+1:end); elseif options_.risky_steadystate [dr,info] = dyn_risky_steadystate_solver(oo_.steady_state,M_,dr, ... options_,oo_); else % If required, use AIM solver if not check only if (options_.aim_solver == 1) && (task == 0) - [dr,info] = AIM_first_order_solver(jacobia_,M_,dr,sdim); + [dr,info] = AIM_first_order_solver(jacobia_,M_,dr,options_.qz_criterium); else % use original Dynare solver [dr,info] = dyn_first_order_solver(jacobia_,b,M_,dr,options_,task); @@ -187,12 +194,12 @@ else f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var)))); f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var)))); fudet = sparse(jacobia_(:,nz+M_.exo_nbr+1:end)); - M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nyf-nboth)]); + M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*dr.gx zeros(M_.endo_nbr,nfwrds-nboth)]); M2 = M1*f1; dr.ghud = cell(M_.exo_det_length,1); dr.ghud{1} = -M1*fudet; for i = 2:M_.exo_det_length - dr.ghud{i} = -M2*dr.ghud{i-1}(end-nyf+1:end,:); + dr.ghud{i} = -M2*dr.ghud{i-1}(end-nfwrds+1:end,:); end end