
gg=which('dynare');
dyny=0;
if ~isempty(regexp(gg,'dynare_v4'));
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
if ~exist('sum_ey','var')
    M_.old_sigma=M_.Sigma_e;
elseif isfield(M_,'old_sigma');
    M_.Sigma_e=M_.old_sigma;
end
%loop for calculating var-cov
lng=rows(M_.Sigma_e);

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

%opty=input('Give name of reference model (empty for no referece)>> ','s');
opty='';
if ~isempty(regexp(opty,'.'))
    eval(['load ',opty]);
    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
hor=1000;
% 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
res=[];
for loop=1:1:(lng+1);
    var=zeros(size(S_mx1));

    

yy=zeros(size(dr_.ghx,1),1);
if loop==1
sigma_old=M_.Sigma_e;
else
M_.Sigma_e=zeros(lng,lng);
M_.Sigma_e((loop-1),(loop-1))=sigma_old((loop-1),(loop-1));
end
% disp('VAR-COV Matrix');
% disp(diag(M_.Sigma_e));
% const=dr_.newconst*M_.Sigma_e(:);
Delta_mx=0.5*dr_.newconst*M_.Sigma_e(:);
%%%%%%%%%% conditional mean of model 
hor=500;
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';
    varall=dr_.ghx*var*dr_.ghx'+dr_.ghu*M_.Sigma_e*dr_.ghu';
    sum_ey=sum_ey+nbeta^(jj-1)*yy(:,end);
    sum_varall=sum_varall+nbeta^(jj-1)*varall;
%    display(sum_varall(1,1));
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;
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,:);


% % disp('difference w.r. to Dynare output')
%   [lgy_,num2str([sum_ey,yy(:,end)])]
  
  varlist={'Ut_h','Ut_f','Ce','Ce_s','l','l_s','pii','pii_s','gdp_h','gdp_f','tot'};
  gg=[varlist;varlist];
  gg=gg(:);
  
  pos=[];
  for jj=1:size(varlist,2);
      pos=[pos,find_incell(str2cell(M_.endo_names),['\<', deblank(varlist{jj}),'\>'])];
     eval([deblank(varlist{jj}),'_pos=',num2str(pos(end)),';']); 
 end
  
% xutil=find_incell(str2cell(M_.endo_names),'\<Util\>');
% xuth=find_incell(str2cell(M_.endo_names),'\<Ut_h\>');
% xutf=find_incell(str2cell(M_.endo_names),'\<Ut_h\>');
% xce=find_incell(str2cell(M_.endo_names),'\<Ce\>');
% xces=find_incell(str2cell(M_.endo_names),'\<Ce_s\>');
% xl=find_incell(str2cell(M_.endo_names),'\<l\>');
% xls=find_incell(str2cell(M_.endo_names),'\<l_s\>');
% xtot=find_incell(str2cell(M_.endo_names),'\<tot\>');
% xepsch=find_incell(str2cell(M_.endo_names),'\<epsilon_ch\>');
% xepscf=find_incell(str2cell(M_.endo_names),'\<epsilon_cf\>');
% xvsigh=find_incell(str2cell(M_.endo_names),'\<vsig_h\>');
% xvsigf=find_incell(str2cell(M_.endo_names),'\<vsig_f\>');


%save first the expected mean and variances of variable
% util_mean=sum_ey(xutil,:);
% util_var=sum_varall(xutil,xutil);
% c_mean=sum_ey(xce,:);
% c_var=sum_varall(xce,xce);
% l_mean=sum_ey(xl,:);
% l_var=sum_varall(xl,xl);
% tot_mean=sum_ey(xtot,:);
% tot_var=sum_varall(xtot,xtot);
% c_pref_var=sum_varall(xce,xepsch);
% l_pref_var=sum_varall(xl,xepsch);
% c_ls_var=sum_varall(xce,xvsigh);
% l_ls_var=sum_varall(xl,xvsigh);


if loop==1
    
%save some of the results in an excel file
d1 =str2cell(lgy_(vec([pos;pos]),:));
%s = xlswrite('Results.xls', d1, 'Figure 1.1', 'B1')
end
variables=[];
% _{util_mean,util_var,c_mean,c_var,l_mean,l_var,tot_mean,tot_var}
z=1;
variables=[
    sum_ey(Ut_h_pos)
sum_varall(Ut_h_pos,Ut_h_pos)
sum_ey(Ut_f_pos)
sum_varall(Ut_f_pos,Ut_f_pos)
sum_ey(Ce_pos)
sum_varall(Ce_pos,Ce_pos)
sum_ey(Ce_s_pos)
sum_varall(Ce_s_pos,Ce_s_pos)
sum_ey(l_pos)
sum_varall(l_pos,l_pos)
sum_ey(l_s_pos)
sum_varall(l_s_pos,l_s_pos)
sum_ey(pii_pos)
sum_varall(pii_pos,pii_pos)
sum_ey(pii_s_pos)
sum_varall(pii_s_pos,pii_s_pos)
sum_ey(gdp_h_pos)
sum_varall(gdp_h_pos,gdp_h_pos)
sum_ey(gdp_f_pos)
sum_varall(gdp_f_pos,gdp_f_pos)
sum_ey(tot_pos)
sum_varall(tot_pos,tot_pos)]';
% write2excel(variables,loop);
res=[res;(variables)];

end
rel_res=res./kron(ones(size(res,1),1),res(1,:));

tt=tabtex(strvcat('Variable',cell2str(gg)),'Welfare Decomposition',3,(M_.exo_names),rel_res(2:end,:)*100);
table('Relative',strvcat('Variable',M_.exo_names([1,4:end],:)),strvcat(cell2str(gg)),rel_res([2,5:end],:)',10,10,  3)
table('Absolute',strvcat('Variable','total',(M_.exo_names([1,4:end],:))),strvcat(cell2str(gg)),rel_res([1,2,5:end],:)',10,10,  3)



  
  % for comparison the following is the asymptotic mean

%savopty=input('Give name to save ref. model (empty if don''t save)>> ', '');
savopty='';
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
% %afterloop;
% end;

% variance=sum_ey, mean=yy(:,end) ???? 

