provisions for Gabaix in solver, stoch_simul, likelihood and smoother

Manually resolved conflicts from 4.7 JRC branch

matlab/DsgeSmoother.m

@@ -104,8 +104,11 @@ if ~options_.smoother_redux
     oldoo.restrict_columns = oo_.dr.restrict_columns;
     oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
     oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns;
-    
-    [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+    if options_.gabaix
+        [T,R,SteadyState,info,M_,oo_] = gabaix_dynare_resolve(M_,options_,oo_);
+    else
+        [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+    end
     
     %get location of observed variables and requested smoothed variables in
     %decision rules
@@ -113,7 +116,11 @@ if ~options_.smoother_redux
     
 else
     if ~options_.occbin.smoother.status
-        [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+        if options_.gabaix
+            [T,R,SteadyState,info,M_,oo_] = gabaix_dynare_resolve(M_,options_,oo_,'restrict');
+        else
+            [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+        end
     else
         [T,R,SteadyState,info,M_,oo_,~,~,~, T0, R0] = ...
             occbin.dynare_resolve(M_,options_,oo_,[],'restrict');

matlab/default_option_values.m

@@ -340,6 +340,7 @@ options_.k_order_solver = false; % by default do not use k_order_perturbation bu
 options_.partial_information = false;
 options_.ACES_solver = false;
 options_.conditional_variance_decomposition = [];
+options_.gabaix = false;
 
 % Ramsey policy
 options_.ramsey_policy = false;

matlab/dsge_likelihood.m

@@ -149,7 +149,7 @@ Hess        = [];
 
 % Ensure that xparam1 is a column vector.
 % (Don't do the transformation if xparam1 is empty, otherwise it would become a
-%  0×1 matrix, which create issues with older MATLABs when comparing with [] in
+%  0×1 matrix, which create issues with older MATLABs when comparing with [] in
 %  check_bounds_and_definiteness_estimation)
 if ~isempty(xparam1)
     xparam1 = xparam1(:);
@@ -208,7 +208,11 @@ if DynareOptions.occbin.likelihood.status
     occbin_.info= {DynareOptions, DynareResults, Model, occbin_options, TTx, RRx, CCx,T0,R0};
 else
     % Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
-    [T,R,SteadyState,info,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+    if DynareOptions.gabaix
+        [T,R,SteadyState,info,Model,DynareResults] = gabaix_dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+    else
+        [T,R,SteadyState,info,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+    end
     occbin_.status = false;
 end
 

matlab/gabaix_dynare_resolve.m

+function [A,B,SteadyState,info,M_,oo_] = gabaix_dynare_resolve(Model,DynareOptions,DynareResults, restricted, A, B, SteadyState)
+
+nargs=nargin;
+if nargin<=4
+    [A,B,SteadyState,info,M_,oo_] = dynare_resolve(Model,DynareOptions,DynareResults);
+    if info(1)
+        return
+    end
+    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);
+%     [~,~,hl1_,~] = get_deriv(Mbase_,SteadyState); % only apply behavioral M to leads
+    Cbarmat = hm1_; % lag
+    Bbarmat = h_;
+    Abarmat = hl1_; % lead
+    
+    % 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)];
+    iyv = M_.lead_lag_incidence';
+    iyv = iyv(:);
+    iyr0 = find(iyv) ;
+    it_ = M_.maximum_lag + 1;
+    z = repmat(oo_.dr.ys,1,klen);
+    task=1; % 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);
+%     if info(1)
+%         A=[];
+%         B=[];
+%         return
+%     end
+    
+%     invmat = inv(Bbarmat+Abarmat*A);
+%     A = -invmat*Cbarmat;
+    A = -(Bbarmat+Abarmat*A)\Cbarmat;
+%     max_eig = max(abs(eig(A)));
+%     if isempty(DynareOptions.qz_criterium)
+%         DynareOptions.qz_criterium = 1+1e-6;
+%     end
+%     if max_eig>DynareOptions.qz_criterium
+%         info(1)=3;
+%         info(2)=(max_eig-DynareOptions.qz_criterium)^2;
+%         A=[];
+%         B=[];
+%         return
+%     end
+    
+%     E = [zeros(length(A)) -Cbarmat; eye(length(A)) zeros(length(A))];
+%     D = [Abarmat Bbarmat; zeros(length(A)) eye(length(A))];
+%     [ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);
+
+%     B = -invmat*Jbarmat_;
+    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)) ;
+else
+    A = A(DynareResults.dr.inv_order_var,DynareResults.dr.inv_order_var);
+    B = B(DynareResults.dr.inv_order_var,:);
+    if isempty(restricted)
+        % we allow for this case, although it just provides same A.B in
+        % input
+        nargs=3;
+    else
+        nargs=4;
+    end
+end
+
+switch nargs
+  case 3
+    endo_nbr = Model.endo_nbr;
+    iv = (1:endo_nbr)';
+  case 4
+    iv = DynareResults.dr.restrict_var_list;
+  otherwise
+    error('gabaix_dynare_resolve:: Error in the calling sequence!')
+end
+
+if nargout==1
+    A = A(iv,iv);
+    return
+end
+
+A = A(DynareResults.dr.order_var(iv),DynareResults.dr.order_var(iv));
+B = B(DynareResults.dr.order_var(iv),:);
+
+

matlab/stoch_simul.m

@@ -24,7 +24,7 @@ end
 
 % The TeX option crashes MATLAB R2014a run with "-nodisplay" option
 % (as is done from the testsuite).
-% Since we can’t directly test whether "-nodisplay" has been passed,
+% Since we cant directly test whether "-nodisplay" has been passed,
 % we test for the "TOP_TEST_DIR" environment variable, which is set
 % by the testsuite.
 % Note that it was not tested whether the crash happens with more
@@ -100,6 +100,13 @@ elseif options_.discretionary_policy
         error('discretionary_policy: only linear-quadratic problems can be solved');
     end
     [~,info,M_,oo_] = discretionary_policy_1(options_.instruments,M_,options_,oo_);
+elseif isfield(options_,'gabaix') && options_.gabaix
+    if options_.logged_steady_state %if steady state was previously logged, undo this
+        oo_.dr.ys=exp(oo_.dr.ys);
+        oo_.steady_state=exp(oo_.steady_state);
+        options_.logged_steady_state=0;
+    end
+    [~,~,~,info,M_,oo_] = gabaix_dynare_resolve(M_,options_,oo_);
 else
     if options_.logged_steady_state %if steady state was previously logged, undo this
         oo_.dr.ys=exp(oo_.dr.ys);