if ~exist('verbose','var')
verbose=1;% one if inputs from keyboard
% else
%     verbose=0;
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');
else 
    opty='';
end
if ~isempty(regexp(opty,'.'))
    eval(['load ',opty]);
    opty.mean=opty.mean(oo_.dr.order_var,:);
    for jj=1:size(opty.lgy_,1);
        eval([opty.lgy_(jj,:),'=',num2str(opty.mean(jj,:)),';']);
    end
%     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


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

    varall=dr_.ghx*var*dr_.ghx'+dr_.ghu*M_.Sigma_e*dr_.ghu';
    varall=varall(dr_.inv_order_var,dr_.inv_order_var);
% 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;
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
if verbose
savopty=input('Give name to save ref. model (empty if don''t save)>> ', 's');
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']);
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