%% COMPUTE SECOND ORDER ACCURATE MEAN BASED ON DYNARE OUTPUT
%% IT ALLOWS TO COMPUTE BEVERIDGE-NELSON DECOMPOSITION TO SUBTRACT LONG RUN
%% TO CHECK BW RESULT
%% cond_mom_dyn4(force_comparison,BN,hor);



function cond_mom_dyn4(force_comp,BN,hor);
persistent info_
if isempty(info_)
    fprintf('=================================================================================================\n')
    fprintf('|| cond_mom_dyn4 info:\n|| Input: 1/0 switch: conditional on saved starting value/on own steady state (defoult=0); \n||\t\t  1/0 switch: Beveridge-Nelson decomposition (default no); \n||\t\t Periods (defoult 4000)\n')
    fprintf('|| output: oo_ (also pased to base)\n')
    fprintf('|| NEED first to save a reference model for the conditional mean:\n')
    fprintf('|| save oo_opt.mat oo_ M_ -MAT \n')
    fprintf('=================================================================================================\n')
    info_=1;
end
if nargin<2
    BN=0; %Beveridge Nelson decomposition for unit root
end

% GENERATE MOMENTS OF MODEL JUST RUN
if nargin<3
hor=4000; % should be so large that beta^hor ~ 0
end
variance_dynare;
% get information on model from workbase
oo_=evalin('base','oo_');
M_=evalin('base','M_');
varnames_neworder=M_.endo_names(oo_.dr.order_var,:);
[xx,pos_rams]=find_incell(varnames_neworder,'\<ramsey');
%% MUST HAVE A DISCOUNT FACTOR CALLED nbeta (the policy discount factor)
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,A_LR]=cond_mom_dyn4_STSTGAP_sub1(oo_.dr,M_,hor,BN);

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',sprintf('Steady-state gap correction\n'))
    refss=input('Give name of steady state oo_ structure file (with .mat ext) [empty if no comparison]==> ','s');
    if ~isempty(refss)
        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
    end
end

if nargin>0
    if force_comp==1;
        if exist('oo_opt.mat','file')
            oo_new=load('oo_opt');
            oo_ref=oo_new.oo_;
            M_ref=oo_new.M_;
        end
    end
end
if exist('oo_new','var')
    athis=str2cell(M_.endo_names(oo_.dr.order_var,:));    
    bref=str2cell(M_ref.endo_names(oo_ref.dr.order_var,:));
    [ii,ref_memb_jj]=ismember(str2cell(M_.endo_names(oo_.dr.order_var,:)),str2cell(M_ref.endo_names(oo_ref.dr.order_var,:)));
    [i1,ithis,iref]=intersect(athis,bref);% av(i2)=av(i3)
    ss_ref=oo_ref.steady_state(oo_ref.dr.order_var);
    
    ss_this=oo_.steady_state(oo_.dr.order_var);
    ss_ref(pos_rams,:)=ss_this(pos_rams,:);
    oo_ref_mean=oo_ref.mean(oo_ref.dr.order_var);
%     oo_ref_mean(pos_rams,:)=ss_this(pos_rams,:); % WAS THIS A TYPO? WHY
%     IS THIS HERE? 1 Dec 2017
    oo_ref_var=oo_ref.var(oo_ref.dr.order_var,oo_ref.dr.order_var);
    if exist('oo_ref','var') %% IF THE PREVIOUS ECONOMY WERE AT THE SS, THIS WOULD BE THE HAT VARIABLE AT TIME ZERO.
      try
          ss_gap=ss_this*0;
          ss_gap(ithis)=ss_ref(iref)-ss_this(ithis);
%         ss_gap=ss_ref(ref_memb_jj)-ss_this;
      catch
          keyboard
      end
    else
        cprintf('String','%s','NEED TO SAVE A REFERENCE MODEL AS oo_opt.mat CONTAINING oo_ and M_ ');
        return
        
    end
    %% moments of the economy used as starting point
    yy=zeros(size(A_mx,1),1);
    Y_mean_ref=yy;
    Y_mean_ref(ithis)=oo_ref_mean(iref);
    Y_mean_ref_state=Y_mean_ref(state_pos);
    
    
    
    %     ss_gap_states=ss_gap(state_pos);
    yy=zeros(size(A_mx,1),1);
    %% Set the initial value of variables at their unconditional mean (of ref model) in level
    yy(:,1)=Y_mean_ref+ss_gap;
    initial_v=Y_mean_ref+ss_gap;
    % This would be ok if start from the unconditional moments;
%          var_ref=oo_ref_var(ref_memb_jj,ref_memb_jj);
%          varc=var_ref(state_pos,state_pos);
    % This instead is better if start from a specific point (not averaging)     
    varc=yy(:,1)*yy(:,1)';
    varc=varc(state_pos,state_pos);
    
    results0=zeros(size(yy,1),2);
else
    yy=zeros(size(A_mx,1),1);
    initial_v=yy;
    results0=zeros(size(yy,1),2);
    varc=zeros(length(state_pos),length(state_pos));
    
end
sum_ey=yy(:,1)*0;%
if ~exist('nbeta','var')||isempty(nbeta)
    error('Missing parameter nbeta (discount factor): must define one in the mod file');
end
% progress('init');
yy_series=zeros(size(yy,1),hor);
for jj=1:hor

    yy=A_mx*yy+C_mx*varc(:)+D_mx*M_.Sigma_e(:)+Delta_mx;
    varc=S_mx1*varc*S_mx1'+S_mx2*M_.Sigma_e*S_mx2';
    sum_ey=sum_ey+nbeta^(jj-1)*yy;
    yy_series(:,jj)=yy;
end




if BN==1
    tmp=A_LR*initial_v;
    results=[yy(:,end)+tmp,sum_ey(:,end)];%The long run term is multiplied by beta^{very large number}=0
else
    results=[yy(:,end),sum_ey(:,end)];
end
%%%%%%%%%%%%%%%%%%%%%%%%



yy_series=yy_series(dr_.inv_order_var,:);
yy=results(:,end-1);
sum_ey=results(:,end);

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;
oo_.mean_disc=sum_ey;
oo_.yy_series=yy_series;
assignin('base','oo_',oo_);

