use M_.occbin.filter_initial_state to perform smoothers running particles from...

use M_.occbin.filter_initial_state to perform smoothers running particles from smoothed estimate of states in t=0

matlab/estimation/prior_posterior_statistics_core.m

@@ -171,8 +171,8 @@ if run_smoother
     stock_smoothed_trend=NaN(endo_nbr,gend,MAX_n_smoothed_trend);
     stock_trend_coeff = zeros(endo_nbr,MAX_n_trend_coeff);
     if options_.occbin.smoother.status
-%         stock_occbin_regime = struct();
-%         stock_occbin_realtime_regime = struct();
+        %         stock_occbin_regime = struct();
+        %         stock_occbin_realtime_regime = struct();
     end
     if horizon
         stock_forcst_mean= NaN(endo_nbr,horizon,MAX_nforc1);
@@ -200,6 +200,11 @@ if smoothed_state_uncertainty
 end
 
 opts_local = options_;
+state_distrib=false;
+if options_.occbin.smoother.status && not(opts_local.lik_init==2 && opts_local.Harvey_scale_factor==0)
+    state_distrib=true;
+end
+
 for b=fpar:B
     if strcmpi(type,'prior')
         iter=1;
@@ -246,10 +251,45 @@ for b=fpar:B
                 message=get_error_message(oo_.occbin.smoother.error_flag,opts_local);
                 fprintf('\nprior_posterior_statistics: One of the draws failed with the error:\n%s\n',message)
                 continue
-            else
-                stock_occbin_regime(:,irun(5))=oo_.occbin.smoother.regime_history;
-                stock_occbin_realtime_regime(:,irun(5))=oo_.occbin.smoother.realtime_regime_history;
             end
+            if not(opts_local.lik_init==2 && opts_local.Harvey_scale_factor==0)
+                state_uncertainty1 = state_uncertainty0(oo_.dr.inv_order_var,oo_.dr.inv_order_var);
+                alphahat0 = a0T(oo_.dr.inv_order_var);
+                alphahat1 = alphahat0(M_.state_var);
+                state_uncertainty1 = state_uncertainty1(M_.state_var,M_.state_var);
+                [U,X] = svd(0.5*(state_uncertainty1+state_uncertainty1'));
+                % P= U*X*U'; % symmetric matrix!
+                is = find(diag(X)>options_.kalman_tol);
+                StateVectorVarianceSquareRoot = chol(X(is,is))';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+
+                % Get the rank of StateVectorVarianceSquareRoot
+                state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+                U = U(:,is);
+
+                if not(isempty(is))
+                    alphahat01 = U(:,is)*StateVectorVarianceSquareRoot*randn(state_variance_rank,1)+alphahat1;
+                else
+                    alphahat01 = alphahat1;
+                end
+                a0TU = a0T*0;
+                a0TU(M_.state_var) = alphahat01;
+                a0TU = a0TU(oo_.dr.order_var);
+                if isfield(M_.occbin,'filter_initial_state') && ~isempty(M_.occbin.filter_initial_state)
+                    M_local = M_;
+                    M_local.filter_initial_state = M_.occbin.filter_initial_state;
+                    opts_local1 = opts_local;
+                    opts_local1.lik_init = 2;
+                    opts_local1.Harvey_scale_factor = 0;
+                    %                 a0T,state_uncertainty0
+                    for jj=1:length(M_.state_var)
+                        M_local.params(strcmp([M_.endo_names{M_.state_var(jj)} 'init'],M_.param_names)) = alphahat01(jj);
+                    end
+                    [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_,a0T,state_uncertainty0] = ...
+                        occbin.DSGE_smoother(deep,gend,Y,data_index,missing_value,M_local,oo_,opts_local1,bayestopt_,estim_params_);
+                end
+            end
+            stock_occbin_regime(:,irun(5))=oo_.occbin.smoother.regime_history;
+            stock_occbin_realtime_regime(:,irun(5))=oo_.occbin.smoother.realtime_regime_history;
         else
             [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_,a0T,state_uncertainty0] = ...
                 DsgeSmoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
@@ -334,7 +374,7 @@ for b=fpar:B
             if options_.prefilter
                 % add mean
                 yf(:,IdObs) = yf(:,IdObs)+repmat(mean_varobs, ...
-                                                 horizon+maxlag,1);
+                    horizon+maxlag,1);
                 % add trend, taking into account that last point of sample is still included in forecasts and only cut off later
                 yf(:,IdObs) = yf(:,IdObs)+((options_.first_obs-1)+gend+[1-maxlag:horizon]')*trend_coeff'-...
                     repmat(mean(trend_coeff*[options_.first_obs:options_.first_obs+gend-1],2)',length(1-maxlag:horizon),1); %center trend
@@ -358,7 +398,7 @@ for b=fpar:B
             else
                 % add trend, taking into account that last point of sample is still included in forecasts and only cut off later
                 yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat(((options_.first_obs-1)+gend+[1-maxlag:horizon]')* ...
-                                                       trend_coeff',[1,1,1]);
+                    trend_coeff',[1,1,1]);
             end
             if options_.loglinear
                 yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
@@ -561,21 +601,21 @@ end
 myoutput.ifil=ifil;
 if RemoteFlag==1
     myoutput.OutputFileName = [OutputFileName_smooth;
-                        OutputFileName_update;
-                        OutputFileName_inno;
-                        OutputFileName_error;
-                        OutputFileName_filter_step_ahead;
-                        OutputFileName_param;
-                        OutputFileName_forc_mean;
-                        OutputFileName_forc_point;
-                        OutputFileName_forc_point_ME;
-                        OutputFileName_filter_covar;
-                        OutputFileName_trend_coeff;
-                        OutputFileName_init_state;
-                        OutputFileName_smoothed_trend;
-                        OutputFileName_smoothed_constant;
-                        OutputFileName_state_uncert;
-                        OutputFileName_init_state_uncert];
+        OutputFileName_update;
+        OutputFileName_inno;
+        OutputFileName_error;
+        OutputFileName_filter_step_ahead;
+        OutputFileName_param;
+        OutputFileName_forc_mean;
+        OutputFileName_forc_point;
+        OutputFileName_forc_point_ME;
+        OutputFileName_filter_covar;
+        OutputFileName_trend_coeff;
+        OutputFileName_init_state;
+        OutputFileName_smoothed_trend;
+        OutputFileName_smoothed_constant;
+        OutputFileName_state_uncert;
+        OutputFileName_init_state_uncert];
 end
 
 dyn_waitbar_close(h);
@@ -622,7 +662,7 @@ endo_nbr=length(y0);
 yf = zeros(endo_nbr,1+horizon,n);
 yf(:,1,:,:) = repmat(y0,[1,1,n]);
 
-for iter=1:horizon    
+for iter=1:horizon
     if stochastic_indicator
         yf(:,iter+1,:) = dr.ghx*squeeze(yf(k2,iter,:))+B1*squeeze(e(:,:,iter));
     else