diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index 173f6947ef01712c074580ddab3e137720ee5633..417e7ea9eb330212d1300c2dd8ee70364bb97975 100644
--- a/matlab/+occbin/set_default_options.m
+++ b/matlab/+occbin/set_default_options.m
@@ -185,6 +185,7 @@ if ismember(flag,{'simul','all'})
options_occbin_.simul.curb_retrench=false;
options_occbin_.simul.endo_init=zeros(M_.endo_nbr,1);
options_occbin_.simul.full_output=true;
+ options_occbin_.simul.gabaix = false;
options_occbin_.simul.init_regime=[];
options_occbin_.simul.init_binding_indicator=false(0);
options_occbin_.simul.exo_pos=1:M_.exo_nbr;
diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index 2f7834df9d589ea14b7c7f497f30411d432eaad2..947510c285b5008cd1ab1940f531966546a4c37c 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -71,6 +71,13 @@ if solve_DM
DM.n_vars = M_base.endo_nbr;
DM.n_exo = M_base.exo_nbr;
+ if opts_simul_.gabaix
+ % set behavioral matrices
+ for k=1:length(M_base.gabaix_weights)
+ mbeha0 = get_param_by_name_locally(['GABAIX_' M_base.gabaix_weights{k}],M_base);
+ M_base = set_param_value_locally(M_base,M_base.gabaix_weights{k},mbeha0);
+ end
+ end
% get the matrices holding the first derivatives for the model
% each regime is treated separately
[DM.Cbarmat, DM.Bbarmat, DM.Abarmat, DM.Jbarmat] = occbin.get_deriv(M_base,data.ys);
@@ -78,6 +85,14 @@ if solve_DM
Mstar_= M_base;
Mstar_.params(M_.occbin.pswitch(1))= 1;
[DM.Cstarbarmat, DM.Bstarbarmat, DM.Astarbarmat, DM.Jstarbarmat, DM.Dstarbarmat] = occbin.get_deriv(Mstar_,data.ys);
+ if opts_simul_.gabaix
+ % re-set behavioral matrices to compute constant term!
+ for k=1:length(M_base.gabaix_weights)
+ M_base = set_param_value_locally(M_base,M_base.gabaix_weights{k},1);
+ Mstar_ = set_param_value_locally(Mstar_,M_base.gabaix_weights{k},1);
+ end
+ [~,~,~,~,DM.Dstarbarmat] = occbin.get_deriv(Mstar_,data.ys);
+ end
[DM.decrulea,DM.decruleb]=occbin.get_pq(dr_base);
@@ -322,7 +337,7 @@ for shock_period = 1:n_shocks_periods
error_flag = 313;
else
disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
- error_flag = 311;
+ error_flag = 311;
end
if opts_simul_.waitbar; dyn_waitbar_close(hh); end
return
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index eab9c269b6614321f7e9cbd40a36567957834c52..7f7e509179c8ba9edf61b13eb4a873d9cbff801d 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -67,20 +67,49 @@ M00_ = M_;
data.ys = dr.ys;
if solve_DM %recompute solution matrices
+ if opts_simul_.gabaix
+ % set behavioral matrices
+ for k=1:length(M00_.gabaix_weights)
+ mbeha0 = get_param_by_name_locally(['GABAIX_' M00_.gabaix_weights{k}],M00_);
+ M00_ = set_param_value_locally(M00_,M00_.gabaix_weights{k},mbeha0);
+ end
+ end
[DM.Cbarmat ,DM.Bbarmat, DM.Abarmat, DM.Jbarmat] = occbin.get_deriv(M00_,data.ys);
M10_ = M00_;
M10_.params(M_.occbin.pswitch(1))= 1;
[DM.Cbarmat10, DM.Bbarmat10, DM.Abarmat10, DM.Jbarmat10, DM.Dbarmat10] = occbin.get_deriv(M10_,data.ys);
-
+ if opts_simul_.gabaix
+ % compute the constant term with the original non Gabaix params
+ for k=1:length(M00_.gabaix_weights)
+ M10_ = set_param_value_locally(M10_,M00_.gabaix_weights{k},1);
+ end
+ [~,~,~,~,DM.Dbarmat10] = get_deriv(M10_,data.ys);
+ end
+
M01_ = M00_;
M01_.params(M_.occbin.pswitch(2))= 1;
[DM.Cbarmat01, DM.Bbarmat01, DM.Abarmat01, DM.Jbarmat01, DM.Dbarmat01] = occbin.get_deriv(M01_,data.ys);
+ if opts_simul_.gabaix
+ % compute the constant term with the original non Gabaix params
+ for k=1:length(M00_.gabaix_weights)
+ M01_ = set_param_value_locally(M01_,M00_.gabaix_weights{k},1);
+ end
+ [~,~,~,~,DM.Dbarmat01] = get_deriv(M01_,data.ys);
+ end
M11_ = M00_;
M11_.params(M_.occbin.pswitch(1))= 1;
M11_.params(M_.occbin.pswitch(2))= 1;
[DM.Cbarmat11, DM.Bbarmat11, DM.Abarmat11, DM.Jbarmat11, DM.Dbarmat11] = occbin.get_deriv(M11_,data.ys);
+ if opts_simul_.gabaix
+ % compute the constant term with the original non Gabaix params
+ for k=1:length(M00_.gabaix_weights)
+ M11_ = set_param_value_locally(M11_,M00_.gabaix_weights{k},1);
+ M00_ = set_param_value_locally(M00_,M00_.gabaix_weights{k},1); % also re-set the base model
+ end
+ [~,~,~,~,DM.Dbarmat11] = get_deriv(M11_,data.ys);
+ end
[DM.decrulea,DM.decruleb]=occbin.get_pq(dr);
update_flag=true;
diff --git a/matlab/+occbin/solver.m b/matlab/+occbin/solver.m
index 3dc465ba3cfeb102ed811a6e78cd4fade92d275e..28f7466f489ae81fb99b7d80060a24c3195ea559 100644
--- a/matlab/+occbin/solver.m
+++ b/matlab/+occbin/solver.m
@@ -19,7 +19,7 @@ function [oo_, out, ss] = solver(M_,oo_,options_)
% - R: [n_vars by n_exo by n_shock_period] array of shock response matrices
% - C: [n_vars by n_shock_period] array of constants
-% Copyright © 2021 Dynare Team
+% Copyright � 2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -64,10 +64,32 @@ if solve_dr
print_info(error_flag, options_.noprint, options_)
return;
end
+ dr0 = dr;
+ if options_.gabaix
+ % set behvioral form
+ options_.occbin.simul.gabaix=true;
+ [A,B,~,error_flag] = gabaix_dynare_resolve(M_,options_,oo_);
+ if ~error_flag
+ endo_nbr = M_.endo_nbr;
+ nstatic = M_.nstatic;
+ nspred = M_.nspred;
+ ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
+ dr.ghu=B ;
+ dr.ghx=A(:,ic) ;
+ end
+ else
+ options_.occbin.simul.gabaix=false;
+ end
oo_.dr = dr;
sto_dr=dr;
sto_M=M_;
+ out.error_flag=error_flag;
+ if error_flag
+ print_info(error_flag, options_.noprint, options_)
+ return;
+ end
else
+ dr0=oo_.dr;
oo_.dr=sto_dr;
end
@@ -105,3 +127,7 @@ end
out.exo_pos = options_.occbin.simul.exo_pos;
oo_.occbin.simul=out;
+
+if options_.gabaix
+ oo_.dr = dr0;
+end
diff --git a/matlab/check.m b/matlab/check.m
index ec9ec65136bc93e71309c656321b868b8df5eb70..a9dc5e670bdb7993c89f7fd1d3a37c4c5e1590f4 100644
--- a/matlab/check.m
+++ b/matlab/check.m
@@ -43,6 +43,13 @@ end
oo.dr=set_state_space(oo.dr,M,options);
[dr,info,M,~] = resol(1,M,options,oo);
+if options.gabaix
+ % set behvioral form
+ [~,~,~,info,M,oo1] = gabaix_dynare_resolve(M,options,oo);
+ dr.eigval=oo1.dr.eigval;
+ dr.edim=oo1.dr.edim;
+ dr.full_rank=oo1.dr.full_rank;
+end
if info(1) ~= 0 && info(1) ~= 3 && info(1) ~= 4
print_info(info, 0, options);
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 5ca97e08610c15937ec936774f94dce2ca9aaf32..d63adf65448a439e55618447e17358b4831b135a 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -348,6 +348,7 @@ options_.ACES_solver = false;
options_.conditional_variance_decomposition = [];
options_.gabaix = false;
options_.gabaix_check_long_run = true;
+options_.gabaix_long_run = true;
% Ramsey policy
options_.ramsey_policy = false;
diff --git a/matlab/gabaix_dynare_resolve.m b/matlab/gabaix_dynare_resolve.m
index f1bf41a899a32ddc9d3817fa36ef1f136b753881..81b926e768fb014262a3a43282932673e06d4f97 100644
--- a/matlab/gabaix_dynare_resolve.m
+++ b/matlab/gabaix_dynare_resolve.m
@@ -9,19 +9,22 @@ if nargin<=4
A = A(DynareResults.dr.inv_order_var,DynareResults.dr.inv_order_var);
B = B(DynareResults.dr.inv_order_var,:);
- [hm1_,h_,~,Jbarmat_] = get_deriv(Model,SteadyState);
Mbase_=M_;
for k=1:length(Model.gabaix_weights)
mbeha0 = get_param_by_name_locally(['GABAIX_' Model.gabaix_weights{k}],Mbase_);
Mbase_ = set_param_value_locally(Mbase_,Model.gabaix_weights{k},mbeha0);
end
- [hm1_,h_,hl1_,Jbarmat_] = get_deriv(Mbase_,SteadyState);
+ [hm1_,h_,hl1_,Jbarmat_] = occbin.get_deriv(Mbase_,SteadyState);
% [~,~,hl1_,~] = get_deriv(Mbase_,SteadyState); % only apply behavioral M to leads
Cbarmat = hm1_; % lag
Bbarmat = h_;
Abarmat = hl1_; % lead
- if options_.gabaix_check_long_run
+ endo_nbr = Model.endo_nbr;
+ nstatic = Model.nstatic;
+ nspred = Model.nspred;
+ ic = [ nstatic+(1:nspred) endo_nbr+(1:size(oo_.dr.ghx,2)-nspred) ]';
+ if DynareOptions.gabaix_check_long_run || DynareOptions.gabaix_long_run
% check BK of behavioral model, to make sure it solves in forward looking mode
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)];
@@ -30,16 +33,30 @@ if nargin<=4
iyr0 = find(iyv) ;
it_ = M_.maximum_lag + 1;
z = repmat(oo_.dr.ys,1,klen);
- task=1; % only check eigvals
+ task=not(DynareOptions.gabaix_long_run); % only check eigvals
[~,jacobia_] = feval([M_.fname '.dynamic'],z(iyr0),exo_simul, ...
Mbase_.params, oo_.dr.ys, it_);
[dr,info] = dyn_first_order_solver(jacobia_,Mbase_,oo_.dr,DynareOptions,task);
+ oo_.dr.eigval = dr.eigval;
+ oo_.dr.edim = dr.edim;
+ if ~DynareOptions.gabaix_long_run
+ oo_.dr.full_rank = dr.full_rank;
+ else
+ oo_.dr.full_rank = info(1)==0;
+ end
if info(1)
A=[];
B=[];
return
else
- A = -(Bbarmat+Abarmat*A)\Cbarmat;
+ if DynareOptions.gabaix_long_run
+ oo_.dr.ghu=dr.ghu;
+ oo_.dr.ghx=dr.ghx;
+ B(DynareResults.dr.order_var,:) =oo_.dr.ghu;
+ A(DynareResults.dr.order_var,DynareResults.dr.order_var(ic))= oo_.dr.ghx ;
+ else
+ A = -(Bbarmat+Abarmat*A)\Cbarmat;
+ end
end
else
@@ -56,13 +73,11 @@ if nargin<=4
return
end
end
- B = -(Bbarmat+Abarmat*A)\Jbarmat_;
- endo_nbr = Model.endo_nbr;
- nstatic = Model.nstatic;
- nspred = Model.nspred;
- ic = [ nstatic+(1:nspred) endo_nbr+(1:size(oo_.dr.ghx,2)-nspred) ]';
- oo_.dr.ghu=B(DynareResults.dr.order_var,:) ;
- oo_.dr.ghx=A(DynareResults.dr.order_var,DynareResults.dr.order_var(ic)) ;
+ if ~DynareOptions.gabaix_long_run
+ B = -(Bbarmat+Abarmat*A)\Jbarmat_;
+ oo_.dr.ghu=B(DynareResults.dr.order_var,:) ;
+ oo_.dr.ghx=A(DynareResults.dr.order_var,DynareResults.dr.order_var(ic)) ;
+ end
else
A = A(DynareResults.dr.inv_order_var,DynareResults.dr.inv_order_var);
B = B(DynareResults.dr.inv_order_var,:);
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index 4a3efcfe8cfc52e80820a66b086d0151db2c23c3..6a0ff9fe1bc73fe5366688a6651ad743bc540c9b 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -242,9 +242,17 @@ if fload==0
%try stoch_simul([]);
try
if ~ isempty(options_.endogenous_prior_restrictions.moment)
+ if options_.gabaix
+ [Tt,Rr,SteadyState,info,M_,oo_] = gabaix_dynare_resolve(M_,options_,oo_);
+ else
[Tt,Rr,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+ end
+ else
+ if options_.gabaix
+ [Tt,Rr,SteadyState,info,M_,oo_] = gabaix_dynare_resolve(M_,options_,oo_,'restrict');
else
[Tt,Rr,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+ end
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T','var')