%% COMPUTE SECOND ORDER ACCURATE MEAN BASED ON DYNARE OUTPUT
%% IT USES BEVERIDGE-NELSON DECOMPOSITION TO SUBTRACT LONG RUN

if ~exist('verbose','var')
    verbose=0;% one if inputs from keyboard
else
    %    verbose=1;
end


hor=1000;
% if exist('dynare_version')>0
%
% 	gg='4';
%
% 	else
% 	gg=which('dynare');
% end
% % gg=which('dynare');
% dyny=0;

% if ~isempty(regexp(gg,'4'));

dr_=oo_.dr;
% else
%     M_=struct('Sigma_e',Sigma_e_);
% [a,b]=sort(dr_.order_var);
%     dr_.inv_order_var=b;
%     oo_.steady_state=ys_;
%
% end
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^1000;
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^1000;
S_mx2=dr_.ghu(nstatic+1:nstatic+size(dr_.ghx,2),:);

% state_v=S_mx1*state_v(-1)+S_mx2*u;
%%%%%%%%%%%%%%%%%%%%
% Initial conditions
if verbose
    opty=input('Give name of reference model (empty for no referece)>> ','s');
    %opty='europol_bgg_optimal';
else
    opty='';
end
if ~isempty(opty)
    eval(['load ',opty]);
    
    %     for jj=1:size(lgy_,1);
    %         yy(jj,1)=opt.yy(find_incell(opt.lgy_,deblank(lgy_(jj,:))),1);
    %         for kk=1:size(lgy_,1)
    %             varall(jj,kk)=opt.varall(find_incell(opt.lgy_,deblank(lgy_(jj,:))),find_incell(opt.lgy_,deblank(lgy_(kk,:))));
    %         end
    %     end
    if isempty(find_incell(lgy_,'ramsey'))
        [a1,b1]=find_incell(opt.lgy_,'ramsey');
        %         [a2,b2]=find_incell(opt.lgy_,'dummy');
        b2=[];
        bb=setdiff([1:size(opt.lgy_,1)]',[b1;b2]); % take out the lagrange multipliers
        opt.lgy_=opt.lgy_(bb,:); % only endogenous variables;
        opt.eyy=opt.eyy(bb,:);   % only endogenous variables' means
                opt.varall=opt.varall(bb,bb); % only endogenous variables' variances;
    end
    
    %The Ramsey monetary union case - one ramsey more than in cooperative
    if size(lgy_,1)==size(opt.lgy_,1)+1
        [a1,b1]=find_incell(opt.lgy_,'ramsey');
        if ~isempty(b1);
            
            
            opt.lgy_=[opt.lgy_(1:b1(end),:);['ramsey' num2str(size(b1,1)+1) '         '];opt.lgy_(b1(end)+1:end,:)];
            opt.eyy=[opt.eyy(1:b1(end),:);0;opt.eyy(b1(end)+1:end,:)];
            opt.varall=[opt.varall(1:b1(end),:);zeros(1,size(opt.varall,1));opt.varall(b1(end)+1:end,:)];
            opt.varall=[opt.varall(:,1:b1(end)) zeros(size(opt.varall,1),1) opt.varall(:,b1(end)+1:end)];
        end
    end
    posvars=[];
    for jj=1:size(lgy_,1);
        posvars(jj,1)=find_incell(opt.lgy_,regexprep(['\<',lgy_(jj,:),'\>'],'\s*',''));
    end
    opt.eyy=opt.eyy(posvars);
    opt.varall=opt.varall(posvars,posvars);
    %assign unconditional mean to initial values
    yy=opt.eyy(dr_.order_var,:);
    varall=opt.varall(dr_.order_var,dr_.order_var);
    var=varall(nstatic+1:nstatic+size(dr_.ghx,2),nstatic+1:nstatic+size(dr_.ghx,2));
    
    %     lgy_st=lgy_st(bb,:);
    %     eval(['Matrix_lagrange_',fname_]);
    %    yy=eval(cell2str(str2cell(lgy_),0))';
    %    lgy_st=lgy_(oo_.dr.order_var,:);
    %    lgy_st=lgy_st(oo_.dr.nstatic+1:oo_.dr.nstatic+oo_.dr.npred,:);
    %     [a,b]=find_incell(str2cell(lgy_st),'ramsey');
    %     bb=setdiff([1:size(lgy_st,1)]',b);
    %     lgy_st=lgy_st(bb,:);
    
    %     posend=[];
    %     for jj=1:size(lgy_,1);
    %         posend=[posend;find_incell(str2cell(opty.lgy_),['\<',deblank(lgy_(jj,:)),'\>'])];
    %     end
    %     posstate_ref=[];
    %     posstate=[];
    %     for jj=1:size(lgy_st,1);
    %         tmp=find_incell(str2cell(opty.state),['\<',deblank(lgy_st(jj,:)),'\>']);
    %         posstate_ref=[posstate_ref;tmp];
    %         if ~isempty(tmp)
    %         posstate=[posstate;jj];
    %         end
    %     end
    % %     yy=opty.mean(posend);
    %     var=zeros(oo_.dr.npred,oo_.dr.npred);
    %     var(posstate,posstate)=opty.var(posstate_ref,posstate_ref);
else
    
    var=zeros(size(S_mx1));
    yy=zeros(size(dr_.ghx,1),hor+1);
    yfo=yy;
    
end

state_v=zeros(npred,1);

% yy2=zeros(size(S_mx1,1)^2,1);

%%%%%%%%%%%% 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_ey=yy(:,1)*0;%
 sum_varall=zeros(rows(yy),rows(yy));

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=dr_.ghx*var*dr_.ghx'+dr_.ghu*M_.Sigma_e*dr_.ghu';
         sum_varall=sum_varall+nbeta^(jj-1)*varall;
    
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,{'variance','sum of cond.','ergodic*(1-BETA)'},[(diag(oo_.var)),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);

if verbose
    savopty=input('Give name of ref. model to save (empty if don''t save)>> ', 's');
    %savopty='';
else
    savopty='';
end
if ~isempty(savopty)
    %    state=lgy_(oo_.dr.order_var,:);
    % state=state(oo_.dr.nstatic+1:oo_.dr.nstatic+oo_.dr.npred,:);
    %  opty=struct('mean',yy(:,end),'var',var,'lgy_',M_.endo_names,'state',state,'stst',ys_);
    %  eval(['save ',savopty,'.mat', ' opty -MAT']);
    %         verbose=0;
    %         cond_mom_dyn4;
    % FOLLOWING ADDED BY MARCIN
    opt.lgy_=lgy_;
    opt.eyy=yy(:,end);
    %             opt.varall=varall;
    eval(['save ' M_.fname '.mat opt;']);
    % FOLLOWING IN ORDER TO COMPUTE PENALTY -- ONLY NEED TO SAVE
    % STUFF FOR OPTIMAL
    variance_dynare;
    dr_opt.vc=tot_var(oo_.dr.order_var,oo_.dr.order_var);% 
         dr_opt.lgy2_=lgy_(oo_.dr.order_var,:);
    save dr_opt.mat dr_opt -MAT
end
if 0
    UU=(eye(size(A_mx,1))-A_mx);
    vv=(eye(size(S_mx1,1)^2)-kron(S_mx1,S_mx1))\vec(S_mx2*M_.Sigma_e*S_mx2');
    vvv=reshape(vv,size(S_mx1));
    yy_as=UU\(C_mx*vvv(:)+D_mx*M_.Sigma_e(:)+Delta_mx);
    yy_as=yy_as(dr_.inv_order_var);
    [lgy_,num2str([yy_as,sum_ey,yy(:,end),(oo_.mean-oo_.steady_state)])]
    
end