Merge branch 'master' into 'master'

occbin fixes and two small new features

See merge request Dynare/dynare!1883

matlab/+occbin/kalman_update_algo_1.m

@@ -194,7 +194,9 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
         else
             opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
         end
-%         opts_simul.init_regime=regimes_(1);
+        if not(options_.occbin.filter.use_relaxation)
+            opts_simul.init_regime=regimes_(1);
+        end
         if M_.occbin.constraint_nbr==1
             myregimestart = [regimes_.regimestart];
         else

matlab/+occbin/kalman_update_algo_3.m

@@ -128,7 +128,7 @@ if M_.occbin.constraint_nbr==1
 else
     myregime = [regimes_.regime1 regimes_.regime2];
 end
-regime_hist = {regimes0};
+regime_hist = {regimes0(1)};
 if M_.occbin.constraint_nbr==1
     regime_end = regimes0(1).regimestart(end);
 end
@@ -180,7 +180,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
             CC(:,2)=ss.C(my_order_var,1);
         end
         newguess=0;
-        regime_hist(niter) = {regimes_};
+        regime_hist(niter) = {regimes_(1)};
         if M_.occbin.constraint_nbr==1
             regime_end(niter) = regimes_(1).regimestart(end);
         end
@@ -199,7 +199,9 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
             opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
         end
         %         end
-        %         opts_simul.init_regime=regimes_(1); %% why don't we use this ???
+        if not(options_.occbin.filter.use_relaxation)
+            opts_simul.init_regime=regimes_(1); 
+        end
         if M_.occbin.constraint_nbr==1
             myregimestart = [regimes_.regimestart];
         else
@@ -212,7 +214,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
         regimes_ = out.regime_history;
         if niter>1
             for kiter=1:niter-1
-                is_periodic(kiter) = isequal(regime_hist{kiter}, regimes_);
+                is_periodic(kiter) = isequal(regime_hist{kiter}, regimes_(1));
             end
             is_periodic =  any(is_periodic);
             if is_periodic
@@ -288,7 +290,13 @@ if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations &&
         CC(:,2) = CCx(:,end);
         [a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
         opts_simul.SHOCKS(1,:) = etahat(:,2)';
-        opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
+        if occbin_options.opts_algo.restrict_state_space
+            tmp=zeros(M_.endo_nbr,1);
+            tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
+            opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
+        else
+            opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
+        end
         if M_.occbin.constraint_nbr==1
             myregimestart = [regimes_.regimestart];
         else

matlab/+occbin/set_default_options.m

@@ -160,6 +160,7 @@ if ismember(flag,{'shock_decomp','all'})
 end
 
 if ismember(flag,{'simul','all'})
+    options_occbin_.simul.algo_truncation = 1;
     options_occbin_.simul.debug = false;
     options_occbin_.simul.curb_retrench=false;
     options_occbin_.simul.endo_init=zeros(M_.endo_nbr,1);
@@ -174,6 +175,7 @@ if ismember(flag,{'simul','all'})
     options_occbin_.simul.check_ahead_periods=200;
     options_occbin_.simul.periodic_solution=true;
     options_occbin_.simul.piecewise_only = false;
+    options_occbin_.simul.reset_regime_in_new_period = false;
     options_occbin_.simul.restrict_state_space=false;
     options_occbin_.simul.SHOCKS=zeros(1,M_.exo_nbr);
     options_occbin_.simul.waitbar=true;

matlab/+occbin/solve_one_constraint.m

@@ -163,6 +163,15 @@ for shock_period = 1:n_shocks_periods
         
         % get the hypothesized piece wise linear solution
         if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:))
+            if iter==1 && opts_simul_.reset_regime_in_new_period
+                binding_indicator=false(size(binding_indicator));
+                binding_indicator_history{iter}=binding_indicator;
+                % analyse violvec and isolate contiguous periods in the other regime.
+                [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
+                regime_history(shock_period).regime = regime;
+                regime_history(shock_period).regimestart = regime_start;
+            end
+            
             [zdatalinear_, SS_out.T(:,:,shock_period), SS_out.R(:,:,shock_period), SS_out.C(:,shock_period), SS, update_flag]=occbin.mkdatap_anticipated_dyn(nperiods_0,DM,...
                 regime_start(end)-1,binding_indicator,...
                 data.exo_pos,data.shocks_sequence(shock_period,:),endo_init,update_flag);
@@ -243,22 +252,26 @@ for shock_period = 1:n_shocks_periods
         
     end
     
-    if regime_change_this_iteration ==1 && max_iter>1
-        disp_verbose(['occbin solver:: period ' int2str(shock_period) '::'],opts_simul_.debug)
-        if is_periodic
-            disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
-            if periodic_solution
-                disp_verbose(['Max error:' num2str(merr) '.'],opts_simul_.debug)
+    if regime_change_this_iteration ==1 
+        if max_iter>opts_simul_.algo_truncation
+            disp_verbose(['occbin solver:: period ' int2str(shock_period) '::'],opts_simul_.debug)
+            if is_periodic
+                disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
+                if periodic_solution
+                    disp_verbose(['Max error:' num2str(merr) '.'],opts_simul_.debug)
+                else
+                    if opts_simul_.waitbar; dyn_waitbar_close(hh); end
+                    error_flag = 310;
+                    return
+                end
             else
+                disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
                 if opts_simul_.waitbar; dyn_waitbar_close(hh); end
-                error_flag = 310;
+                error_flag = 311;
                 return
             end
         else
-            disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
-            if opts_simul_.waitbar; dyn_waitbar_close(hh); end
-            error_flag = 311;
-            return
+            binding_indicator = binding_indicator_history{end};
         end
     end
     if any(error_code_period)

matlab/+occbin/solve_two_constraints.m

@@ -70,11 +70,11 @@ if solve_DM %recompute solution matrices
     [DM.Cbarmat ,DM.Bbarmat, DM.Abarmat, DM.Jbarmat] = occbin.get_deriv(M00_,data.ys);
     
     M10_ = M00_;
-    M10_.params(strmatch(M_.params(M_.occbin.pswitch(1)),M00_.param_names,'exact'))= 1;
+    M10_.params(M_.occbin.pswitch(1))= 1;
     [DM.Cbarmat10, DM.Bbarmat10, DM.Abarmat10, DM.Jbarmat10, DM.Dbarmat10] = occbin.get_deriv(M10_,data.ys);
     
     M01_ = M00_;
-    M01_.params(strmatch(M_.params(M_.occbin.pswitch(2)),M00_.param_names,'exact'))= 1;
+    M01_.params(M_.occbin.pswitch(2))= 1;
     [DM.Cbarmat01, DM.Bbarmat01, DM.Abarmat01, DM.Jbarmat01, DM.Dbarmat01] = occbin.get_deriv(M01_,data.ys);
     
     M11_ = M00_;
@@ -177,6 +177,17 @@ for shock_period = 1:n_shocks_periods
         regime_history(shock_period).regime2 = regime_2;
         regime_history(shock_period).regimestart2 = regime_start_2;
         if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:)) % first period or shock happening
+            if iter==1 && opts_simul_.reset_regime_in_new_period
+                binding_indicator=false(size(binding_indicator));
+                binding_indicator_history{iter}=binding_indicator;
+                % analyse violvec and isolate contiguous periods in the other regime.
+                [regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
+                regime_history(shock_period).regime1 = regime_1;
+                regime_history(shock_period).regimestart1 = regime_start_1;
+                [regime_2, regime_start_2, error_code_period(2)]=occbin.map_regime(binding_indicator(:,2),opts_simul_.debug);
+                regime_history(shock_period).regime2 = regime_2;
+                regime_history(shock_period).regimestart2 = regime_start_2;
+            end
             Tmax=max([regime_start_1(end) regime_start_2(end)])-1;
             [zdatalinear_, SS_out.T(:,:,shock_period), SS_out.R(:,:,shock_period), SS_out.C(:,shock_period), SS, update_flag]=occbin.mkdatap_anticipated_2constraints_dyn(nperiods_0,...
                 DM,Tmax,...
@@ -271,22 +282,28 @@ for shock_period = 1:n_shocks_periods
             binding_indicator_history{iter}=binding_indicator;
         end
     end
-    if regime_change_this_iteration && max_iter>1
-        disp_verbose(['occbin solver: period ' int2str(shock_period) ':'],opts_simul_.debug)
-        if is_periodic
-            disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
-            if periodic_solution
-                disp_verbose(['Max error:' num2str(min_err) '.'],opts_simul_.debug)
+    if regime_change_this_iteration
+        if max_iter>opts_simul_.algo_truncation
+            disp_verbose(['occbin solver: period ' int2str(shock_period) ':'],opts_simul_.debug)
+            if is_periodic
+                disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
+                if periodic_solution
+                    disp_verbose(['Max error:' num2str(min_err) '.'],opts_simul_.debug)
+                else
+                    error_flag = 310;
+                    if opts_simul_.waitbar; dyn_waitbar_close(hh); end
+                    return;
+                end
             else
-                error_flag = 310;
+                disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
+                error_flag = 311;
                 if opts_simul_.waitbar; dyn_waitbar_close(hh); end
                 return;
             end
         else
-            disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
-            error_flag = 311;
-            if opts_simul_.waitbar; dyn_waitbar_close(hh); end
-            return;
+            % if max_iter <= truncation, we force indicator to equal the
+            % last guess
+            binding_indicator = binding_indicator_history{end};
         end
     end
     if any(error_code_period)

matlab/+occbin/write_regimes_to_xls.m

@@ -35,14 +35,14 @@ end
 xls_filename = options_.occbin.write_regimes.filename;
 
 if isfield(regime_history,'regime')
-    Header = {'time', 'regime sequence', 'starting period of regime'};
+    Header = {'time', 'regime_sequence', 'starting_period_of_regime'};
     for tp=1:length(T)
         xlsmat{tp,1}=T(tp);
         xlsmat{tp,2}=int2str(regime_history(tp).regime);
         xlsmat{tp,3}=int2str(regime_history(tp).regimestart);
     end
 else
-    Header = {'time', 'regime sequence 1', 'starting period of regime 1', 'regime sequence 2', 'starting period of regime 2'};
+    Header = {'time', 'regime_sequence_1', 'starting_period_of_regime_1', 'regime_sequence_2', 'starting_period_of_regime_2'};
     for tp=1:length(T)
         xlsmat{tp,1}=T(tp);
         xlsmat{tp,2}=int2str(regime_history(tp).regime1);

matlab/DsgeSmoother.m

@@ -366,6 +366,7 @@ else
     if options_.occbin.smoother.status
         % reconstruct occbin smoother
         if length_varargin>0
+            % sequence of regimes is provided in input
             isoccbin=1;
         else
             isoccbin=0;
@@ -431,6 +432,8 @@ else
                 bbb(oo_.dr.inv_order_var,k) = out.zpiece(1,:);
             end
         end
+        % do not overwrite accurate computations using reduced st. space
+        bbb(oo_.dr.restrict_var_list,:)=ahat;
         ahat0=ahat;
         ahat=bbb;
         if ~isempty(P)

matlab/missing_DiffuseKalmanSmootherH3_Z.m

@@ -212,6 +212,13 @@ else
         R0=R;
         
     end
+    if ~isinf(occbin_options.first_period_occbin_update)
+        % initialize state matrices (otherwise they are set to 0 for
+        % t<first_period_occbin_update!)
+        TTT=repmat(T0,1,1,smpl+1);
+        RRR=repmat(R0,1,1,smpl+1);
+        CCC=repmat(zeros(length(T0),1),1,smpl+1);
+    end
     
 end