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

	gg='4';

	else
	gg=which('dynare');
end 
% gg=which('dynare');
dyny=0;
hor=1000;
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
    
A_mx = [zeros(size(lgy_,1),dr_.nstatic),dr_.ghx,zeros(size(lgy_,1),size(lgy_,1)-dr_.npred-dr_.nstatic)];
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(dr_.nstatic+1:dr_.nstatic+size(dr_.ghx,2),:);
S_mx2=dr_.ghu(dr_.nstatic+1:dr_.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(regexp(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]);
        opt.lgy_=opt.lgy_(bb,:);
        opt.eyy=opt.eyy(bb,:);
        opt.varall=opt.varall(bb,bb);
    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
    
    yy=opt.eyy(dr_.order_var,:);
    varall=opt.varall(dr_.order_var,dr_.order_var);
    var=varall(dr_.nstatic+1:dr_.nstatic+size(dr_.ghx,2),dr_.nstatic+1:dr_.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),1);
yfo=yy;

end

state_v=zeros(dr_.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=dr_.ghx*var*dr_.ghx'+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*0;%the first is zero anyway
sum_varall=zeros(rows(yy),rows(yy));

for jj=1:hor
    
    yy(:,end+1)=A_mx*yy(:,end)+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(:,end);
    %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;

%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)])])
end
% for comparison the following is the asymptotic mean
% replace dynare mean with ergodic
oo_.var=varall;
oo_.mean=yy(:,end);

if verbose
savopty=input('Give name to save ref. model (empty if don''t save)>> ', 's');
%savopty='';
else 
    savopty='';
end
if ~isempty(regexp(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;']); 
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