%% COMPUTE SECOND ORDER ACCURATE MEAN BASED ON DYNARE OUTPUT
%% IT USES BEVERIDGE-NELSON DECOMPOSITION TO SUBTRACT LONG RUN
% Input: nothing
% output: oo_ (also pased to base)
% NEED first to save a reference model for the conditional mean:
% save oo_opt.mat oo_ M_ -MAT



function oo_=cond_mom_dyn4(varargin)

parallel_=1; %e.g. for debuggin set it to 0
kk=4; % number of nodes
repet=1000; %number of initial conditions
hor=1000; % should be so large that beta^hor ~ 0

oo_=evalin('base','oo_');
M_=evalin('base','M_');
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

if ~isempty(verbose) && verbose
    cprintf('String','%s','Steady-state gap correction\n')
    refss=input('Give name of steady state oo_ structure file (with .mat ext) [empty if no comparison]==> ','s');
    if ~exist(refss,'file')
        cprintf('String','%s',['File ',refss,' does not exist']);
        return
    else
        oo_new=load(refss);
        oo_ref=oo_new.oo_;
        M_ref=oo_new.M_;
        
    end
else
    if exist('oo_opt.mat','file')
        oo_new=load('oo_opt');
        oo_ref=oo_new.oo_;
        M_ref=oo_new.M_;
    end
end

[ii,ref_memb_jj]=ismember(str2cell(M_.endo_names(oo_.dr.inv_order_var,:)),str2cell(M_ref.endo_names(oo_ref.dr.inv_order_var,:)));
ss_ref=oo_ref.steady_state(oo_ref.dr.inv_order_var);
ss_this=oo_.steady_state(oo_.dr.inv_order_var);
if exist('oo_ref','var')
    ss_gap=ss_this-ss_ref(ref_memb_jj);
else
    cprintf('String','%s','NEED TO SAVE A REFERENCE MODEL AS oo_opt.mat CONTAINING oo_ and M_ ');
    return
    
end





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

dr_=oo_.dr;

[A_mx,B_mx,C_mx,D_mx,E_mx,Delta_mx,S_mx1,S_mx2,state_pos]=cond_mom_dyn4_STSTGAP_sub1(oo_.dr,M_,hor);
[A_mx_ref,B_mx_ref,C_mx_ref,D_mx_ref,E_mx_ref,Delta_mx_ref,S_mx_ref1,S_mx_ref2,state_pos_ref]=cond_mom_dyn4_STSTGAP_sub1(oo_ref.dr,M_ref,hor);
ss_gap_states=ss_gap(state_pos);





%%%%%%%%%%%% 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));
%% CORE


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


kk=min(kk,repet);
progb=1;
if parallel_
    pp=gcp('nocreate');
    if isempty(getCurrentWorker)
        
        %     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);
            %
        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...');
    %    addAttachedFiles(pp, {'cond_mom_dyn4_STSTGAP_sub1.m','cond_mom_dyn4_STSTGAP_sub2.m'});
%     Itrans*[1:size(M_ref.endo_names(state_pos_ref,:),1)]'
        Itrans=eye(size(A_mx_ref,1));
        Itrans=Itrans(:,state_pos_ref);
        Itrans=Itrans(ref_memb_jj,:);
        Itrans=Itrans(state_pos,:);
    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
        %% common part of parallel and not parallel (but must be explicitly written for transparency
        %%%%%%%%%%%%%%%%%%%%%%%%

    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_ref=zeros(size(A_mx_ref,1),hor+sz_burn);
        sum_ey_ref=yy_ref(:,1);
        var_ref=zeros(size(S_mx_ref1));
        S_ref=zeros(size(S_mx_ref1,1),1);
        

        for uu=1:size(pre_shocks,2)
            
            yy_ref(:,uu+1)=A_mx_ref*yy_ref(:,uu)+C_mx_ref*vec(S_ref(:,uu)*S_ref(:,uu)')+D_mx_ref*vec(pre_shocks(:,uu)*pre_shocks(:,uu)')+Delta_mx_ref;
            %                 var=S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
            sum_ey_ref=sum_ey_ref+nbeta^(uu-1)*yy_ref(:,uu+1);
            S_ref(:,uu+1)=S_mx_ref1*S_ref(:,uu)+S_mx_ref2*chol(M_ref.Sigma_e)*pre_shocks(:,uu);
        end
        
        
        %             S=S_mx2*chol(M_.Sigma_e)*innov_series_(:,ii);
        var=Itrans*S_ref(:,end)*S_ref(:,end)'*Itrans';
        var_ss=ss_gap_states*ss_gap_states';
        %             var=(var);
        %             yy(:,1)=D_mx*vec(innov_series_(:,ii)*innov_series_(:,ii)')+Delta_mx;
        yy=yy_ref(ref_memb_jj,sz_burn:end);
        %             yy(:,1)=yy(:,1);
        results0=zeros(size(yy,1),2,repet_loc);
        
        sum_ey=yy(:,1)*0;%
        jj=1;
        yy_into_new_env= A_mx*yy(:,jj)+C_mx*var(:)+D_mx*M_.Sigma_e(:)+Delta_mx;
% % %%%% Note that yy=(Y_t-Ybar); After shocks then must be
% yy=(Y_t-Ybar)=(Y_t-Ybar_ref)-(Ybar-Ybar_ref); and note that
% ss_gap=Ybar-Ybar_ref;

        yy(:,jj+1)=yy_into_new_env-ss_gap ;
%         var=-var_ss+S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
        var=yy(state_pos,jj+1)*yy(state_pos,jj+1)'; %This is not in expectation but is initial value of states
        sum_ey=sum_ey+nbeta^(jj-1)*yy(:,jj+1);
        for jj=2: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==progb %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
else
    
    %% common part of parallel and not parallel (but must be explicitly written for transparency
    %%%%%%%%%%%%%%%%%%%%%%%%

    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_ref=zeros(size(A_mx_ref,1),hor+sz_burn);
        sum_ey_ref=yy_ref(:,1);
        var_ref=zeros(size(S_mx_ref1));
        S_ref=zeros(size(S_mx_ref1,1),1);

        for uu=1:size(pre_shocks,2)
            
            yy_ref(:,uu+1)=A_mx_ref*yy_ref(:,uu)+C_mx_ref*vec(S_ref(:,uu)*S_ref(:,uu)')+D_mx_ref*vec(pre_shocks(:,uu)*pre_shocks(:,uu)')+Delta_mx_ref;
            %                 var=S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
            sum_ey_ref=sum_ey_ref+nbeta^(uu-1)*yy_ref(:,uu+1);
            S_ref(:,uu+1)=S_mx_ref1*S_ref(:,uu)+S_mx_ref2*chol(M_ref.Sigma_e)*pre_shocks(:,uu);
        end
        
        
        %             S=S_mx2*chol(M_.Sigma_e)*innov_series_(:,ii);
        var=Itrans*S_ref(:,end)*S_ref(:,end)'*Itrans';
        var_ss=ss_gap_states*ss_gap_states';
        %             var=(var);
        %             yy(:,1)=D_mx*vec(innov_series_(:,ii)*innov_series_(:,ii)')+Delta_mx;
        yy=yy_ref(ref_memb_jj,sz_burn:end);
        %             yy(:,1)=yy(:,1);
        results0=zeros(size(yy,1),2,repet_loc);
        
        sum_ey=yy(:,1)*0;%
        jj=1;
        yy_into_new_env= A_mx*yy(:,jj)+C_mx*var(:)+D_mx*M_.Sigma_e(:)+Delta_mx;
% % %%%% Note that yy=(Y_t-Ybar); After shocks then must be
% yy=(Y_t-Ybar)=(Y_t-Ybar_ref)-(Ybar-Ybar_ref); and note that
% ss_gap=Ybar-Ybar_ref;

        yy(:,jj+1)=yy_into_new_env-ss_gap ;
%         var=-var_ss+S_mx1*var*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
        var=yy(state_pos,jj+1)*yy(state_pos,jj+1)'; %This is not in expectation but is initial value of states
        sum_ey=sum_ey+nbeta^(jj-1)*yy(:,jj+1);
        for jj=2: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==progb %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
if parallel_==1
    results=cat(3,op{:}); %
else
    results=op;
end
%     close(hbar);   %close progress bar

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_);

