%% COMPUTE SECOND ORDER ACCURATE MEAN BASED ON DYNARE OUTPUT
%% IT USES BEVERIDGE-NELSON DECOMPOSITION TO SUBTRACT LONG RUN
function oo_=cond_mom_dyn4(varargin);
vrb=evalin('base','exist(''verbose'',''var'')');
if ~vrb
    verbose=0;% one if inputs from keyboard
else
    verbose=evalin('base','verbose');
end
if isempty(verbose)
    verbose=0;
end

hor=3000; % should be so large that beta^hor ~ 0
oo_=evalin('base','oo_');

M_=evalin('base','M_');
pos_nbeta=find_incell(M_.param_names,'nbeta');
nbeta=M_.params(pos_nbeta);

dr_=oo_.dr;

npred=M_.nspred;
nstatic=M_.nstatic;
lgy_=M_.endo_names;
A_mx_0 = [zeros(size(lgy_,1),nstatic),dr_.ghx,zeros(size(lgy_,1),size(lgy_,1)-npred-nstatic)];
%% B-N long term expectation (first order)
A_mx=A_mx_0-A_mx_0^hor;
B_mx = dr_.ghu;
C_mx=0.5*dr_.ghxx;
D_mx=0.5*dr_.ghuu;
E_mx=dr_.ghxu;
Delta_mx=0.5*dr_.ghs2;

% y=A*y(-1)+B*u+C*kron(y(-1),y(-1))+D*kron(u,u)+E*kron(y(-1),u)+Delta_mx;
%%%%
% reordering
% rows of the matrices are re-ordered see oo_.dr
%%%%%%%%%%%%
% State transition matrices
S_mx1=dr_.ghx(nstatic+1:nstatic+size(dr_.ghx,2),:);
S_mx1=S_mx1-S_mx1^hor;
S_mx2=dr_.ghu(nstatic+1:nstatic+size(dr_.ghx,2),:);





%%%%%%%%%%%% just implresp to check ordering and results with Dynare's
%%%%%%%%%%%% output
if 0
    
    % time 0
    innov=zeros(size(S_mx2,2),1);
    innov(1)=1;
    state_v=S_mx2*innov;
    var=zeros(size(S_mx1));
    yfo=dr_.ghu*innov;
    for jj=1:hor
        state_v(:,end+1)=S_mx1*state_v(:,end);
        yfo(:,end+1)=dr_.ghx*state_v(:,end-1);
        var=S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
    end
    yfo=yfo(dr_.inv_order_var,:);
    %variance
    varall=A_mx*var*A_mx'+dr_.ghu*M_.Sigma_e*dr_.ghu';
    varall=varall(dr_.inv_order_var,dr_.inv_order_var);
end
%%%%%%%%%% conditional mean of model


% sum_varall=zeros(rows(yy),rows(yy));

parallel_=1;
leng=100;
innov_series=randn(size(M_.Sigma_e,1),leng);
repet=leng;

kk=4;
kk=min(kk,repet);
if parallel_
    
    if isempty(getCurrentWorker)
        pp=gcp('nocreate');
        %     listatt={[namefile,'.mod'], [namefile,'_steadystate.m']}; %THIS LIST IS MODEL SPECIFIC (in part)
        
        if isempty(pp)
            %% configure pool
            c=parcluster('local');
            
            job1 = createJob(c);
            parpool(c,kk)
        else
            c=parcluster('local');
            
            job1 = createJob(c);
            %         addAttachedFiles(pp, listatt);
        end
        %     pctRunOnAll eval(['global doloop ']);
    end
    %% Make sure there is at least one point per processor
    
    %% Submit job (contained in the wrapper "nodes_part"
    %     hbar = parfor_progressbar(N,'Computing...');
    
    spmd
        zwaitbar=1;
        
        
        % generate distribution of this matrix (i.e. rangeA is subdivided into
        % partitions of equal size distributed over the nodes)
        %         D = codistributed(rangeA ,codistributor1d(2));
        DSH= codistributed(innov_series ,codistributor1d(2));
        %         repet_loc=getLocalPart(D);
        innov_series_=getLocalPart(DSH);
        repet_loc=size(innov_series_,2);
        %%%%%%%%%%%%%% COMPUTATION
        
        for ii=1:repet_loc
            
            if labindex==0
                            progress('init');
%                 upd = textprogressbar(n, 'barlength', 20, ...
%                     'updatestep', 10, ...
%                     'startmsg', 'Waiting... ',...
%                     'endmsg', ' Yay!', ...
%                     'showbar', true, ...
%                     'showremtime', true, ...
%                     'showactualnum', true, ...
%                     'barsymbol', '+', ...
%                     'emptybarsymbol', '-');
                
            end
            sz_burn=10;
            pre_shocks=randn(size(innov_series,1),sz_burn);
            
            yy=zeros(size(A_mx,1),hor+3);
            sum_ey=yy(:,1);
            var=zeros(size(S_mx1));
            S=zeros(size(S_mx1,1),1);
            results0=zeros(size(yy,1),2,repet_loc);
            for uu=1:size(pre_shocks,2);
                
                yy(:,uu+1)=A_mx*yy(:,uu)+C_mx*vec(S(:,uu)*S(:,uu)')+D_mx*vec(pre_shocks(:,uu)*pre_shocks(:,uu)')+Delta_mx;
                %                 var=S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
                sum_ey=sum_ey+nbeta^(uu-1)*yy(:,uu+1);
                S(:,uu+1)=S_mx1*S(:,uu)+S_mx2*chol(M_.Sigma_e)*pre_shocks(:,uu);
            end
            
            
            %             S=S_mx2*chol(M_.Sigma_e)*innov_series_(:,ii);
            var=S(:,end)*S(:,end)';
            %             yy(:,1)=D_mx*vec(innov_series_(:,ii)*innov_series_(:,ii)')+Delta_mx;
            yy=yy(:,4:end);
            %             sum_ey=yy(:,1);%
            for jj=1:hor
                
                yy(:,jj+1)=A_mx*yy(:,jj)+C_mx*var(:)+D_mx*M_.Sigma_e(:)+Delta_mx;
                var=S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
                sum_ey=sum_ey+nbeta^(jj-1)*yy(:,jj+1);
                %MK - the following 2 lines are copied from cond_mom_dyn4_1.m
                %                 varall=A_mx*var*A_mx'+dr_.ghu*M_.Sigma_e*dr_.ghu';
                %                 sum_varall=sum_varall+nbeta^(jj-1)*varall;
                
            end
            results0(:,:,ii)=[yy(:,end),sum_ey(:,end)];
            if labindex==0 %swithc this off in nested parallel counters
                                progress(zwaitbar/repet_loc,['']);
                                zwaitbar=zwaitbar+1;
%                 upd(i);%update progress bar
            end
            %         hbar.iterate(1);   % update progress by one iteration
        end
        
        
        op=results0;
        
    end
    
    %% VERY IMPORTANT:  The dimension in the first entry depends on whether have singleton on second dimension or not if yes then 1 if no then 2
    results=cat(3,op{:}); %
    %     close(hbar);   %close progress bar
end

yy=mean(results(:,end-1,:),3);
sum_ey=mean(results(:,end,:),3);

% presample draws;

% for ii=1:repet
%     y=B_mx*chol(M_.Sigma_e)*innov_series(:,ii);
%     var=y*y';
%     yy=D_mx*vec(innov_series(:,ii)*innov_series(:,ii)')+Delta_mx;
%     sum_ey=yy(:,1);%
%     for jj=1:hor
%
%         yy(:,jj+1)=A_mx*yy(:,jj)+C_mx*var(:)+D_mx*M_.Sigma_e(:)+Delta_mx;
%         var=S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
%         sum_ey=sum_ey+nbeta^(jj-1)*yy(:,jj+1);
%         %MK - the following 2 lines are copied from cond_mom_dyn4_1.m
%         varall=A_mx*var*A_mx'+dr_.ghu*M_.Sigma_e*dr_.ghu';
%         sum_varall=sum_varall+nbeta^(jj-1)*varall;
%
%     end
% end

% % alternatively can add the constant later
% % Constant=(eye(size(A_mx,1))-A_mx)\Delta_mx;
% % sum_ey=sum_ey+1/(1-nbeta)*Constant;
% % yy(:,end)=yy(:,end)+Constant;
% variance_dynare;
%MK - the following 2 lines are copied from cond_mom_dyn4_1.m
% varall=varall(dr_.inv_order_var,dr_.inv_order_var);
% sum_varall=sum_varall(dr_.inv_order_var,dr_.inv_order_var);
sum_ey=sum_ey(dr_.inv_order_var,:);
yy=yy(dr_.inv_order_var,:);
if verbose
    disp('NOTICE THAT  UNCONDITIONAL/CONDITIONAL SUMS  IS A FACTOR OF (1-BETA)')
    %     fprintf('variable\t sum of cond.\t ergodic\n')
    %     disp([lgy_,num2str([sum_ey,yy(:,end)])])
    
    mytable(M_.endo_names,{'sum of cond.','ergodic*(1-BETA)'},[sum_ey,yy(:,end)]);
    
end
% for comparison the following is the asymptotic mean
% replace dynare mean with ergodic
% oo_.var=varall;
oo_.mean=yy(:,end);
oo_.mean_disc=sum_ey(:,end);
assignin('base','oo_',oo_);

