% name of model and loading of symbsolve produced information matrix
% clear klmx_d_f
clear all
    verbose=1;
warning('off','MATLAB:declareGlobalBeforeUse')
global alpha_tilde sig muf;
 sig=10;
 muf=sig/(sig-1);
 fname=input('give name of mod file ==> ','s');
eval(['symbsolve ',fname])
namemod=M_.fname;%'ds_example_lev'_optimal;

eval(['load  matrices_',namemod]);

%  alpha_tilde=0;
% some general reorganization of the information 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lgy_=(matrices.varyt(end:-1:1,:));
M_=struct('param_names',matrices.parnames);
M_.fname=namemod;
M_.endo_names=lgy_;
eval([namemod,'_steadystate;']);
lgy_=str2cell(lgy_);
% keyboard
if 1
eval([namemod,'_M']);
eval(['parameters_' ,namemod]);
nsh=size(Msh,2);
sortinnov=matrices.sortinnov;
% M=zeros(size(PHI,1),size(Msh,2));
Msh=Msh(:,sortinnov);
shocky=str2cell(matrices.shocky(end:-1:1,:));

% optional generate varcovmatrix using key std
perm_inno=matrices.sortinnov;
 inno_shocks=matrices.inno_shock;

innovations=cell2str(shocky);
%innovations=innovations(perm_inno-2,:);
    sprintf('Give vector of %3i variances',size(shocky,1))
    disp(upper('the list of shocks is'))
    disp(shocky(end:-1:1,:))
%    disp(strvcat(upper(['give vector of ',num2str(nsh-1),' standard errors for this list']),...
%         'default is ''std_<name in the list>'''));
     Sigma=input('===>>>>');

     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     if size(Sigma,1)~=size(Sigma,2)
    Sigma=Sigma(:);
end
if isempty(Sigma)

    Sigma=[];
    for jj=1:nsh
        eval(['Sigma=[Sigma;std_',innovations(end-jj+1,3:end),'];']);
    end
end
if verbose
    disp(upper('This is the list of innovations'));
    disp([cell2str(str2cell(innovations(end:-1:1,:))),num2str(Sigma)])
end
% Sigma=Sigma(end:-1:1);
if size(Sigma,2)==1;
    Sigma_m=diag(Sigma).^2;
else
    Sigma_m=Sigma.^2;
end
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Msho=[Msh;zeros(matrices.neq-matrices.jumping,nsh)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
var4mom={ 'Y','GDP','C','INVS','Q_K','R_K','REX','R','PII','PII_H','PII_X','PII_IM','X','IM'};
lgycel=str2cell(strvcat(lgy_));
posinlgy=[];
for kkinx=1:size(var4mom',1);
    tmp=strmatch(var4mom{kkinx},lgycel,'exact');
   posinlgy=[posinlgy,tmp];
end

end
% tmp=[];
% rho_pii=0;rho_p=0;

uu=find_incell(str2cell(M_.param_names),'\<alpha\_tilde\>');
uusig=find_incell(str2cell(M_.param_names),'\<sig\>');
uumuf=find_incell(str2cell(M_.param_names),'\<muf\>');

ww=find_incell((lgy_),'\<Util\>');
wwh=find_incell((lgy_),'\<Ut\_h\>');
wwf=find_incell((lgy_),'\<Ut\_f\>');
drb_pos=find_incell((lgy_),'\<drb\>');
dC_pos=find_incell((lgy_),'\<dC\>');
gg=[];  
zz=1;
global war indet
% alpha_K=.01;
% alpha_N1=1-alpha_K;
res=[];
tmp=[];
zz=1;
tic

plots=1;
 for jj=1;%;10:100:1010;
%        alpha_tilde=0;6.3297;
%       sig=jj;
%       muf=sig/(sig-1);
%       M_.params(uu)=alpha_tilde;
%       M_.params(uusig)=sig;
%       M_.params(uumuf)=muf;
% tic    
[Pi,PHI,AA,BB,numpre,aqz,bqz,Qqz,Zqz,inz22,algo]=firstorder_stuff(namemod,matrices,plots);
hor=500;
 [varY,covy]=variance(PHI,Pi,Msho,(Sigma_m),hor);
 disp(' variance in varY and covariance in covy')
% toc
% sum(sum(abs(PHI)))
% sum(sum(abs(Pi)))
% global C_ss
% tmp(:,zz)=C(1);
% return

if isfield(matrices,'coefcov0');
[P4,F4,PHI1,PHI2,Pi1,Pi2,bigpsi,Ip,xx3,M,Dn3a,Mh,PHIh]=secorder_stuff(namemod,AA,BB,numpre,PHI,Pi,matrices,aqz,bqz,Qqz,Zqz,inz22,0);
% OUTPUT:
%  all second order matrices include the 1/2 of the expansion;
%  P4=constant all variables ;
%  F4= constant of all states;
%  PHI2= coeff. matrix of states on x-product of states
%  PHI1= autorecursive coefficient first order endogenous states (same as last lower block of PHI)
%  Pi1= mapping from endogenous states to all variables (same as last columns of Pi)
%  Pi2= mapping from x-products of states (vech) to all variables
%  bigpsi= x-product of lagged states and contemporaneous shocks
% bet=1
% Sigma(end)=0;
bet=input('Give value of discount factor beta ');
[CondM, CondV,CondACov,justCV,erg_CondM,erg_V,disc]=...
    cond_mom_ur_iter(bet,Dn3a,M,Mh,Pi,Pi2,PHI2,bigpsi,PHI,PHIh,P4,F4,vech2(diag(Sigma)));



soirf=input('Do second order IRFs? ');
if soirf
    ir_so_pot;
end
%  res(zz,1)=jj;
 % uncomment this if using the following part
%  res(zz,1)=jj;
%  res(zz,2)=CondM(ww);
%  res(zz,3)=CondM(wwh);
%  res(zz,4)=CondM(wwf);
%  zz=zz+1;
end
 end
 toc
 if isfield(matrices,'coefcov0');
if size(res,1)>1
     figure;plot(res(:,1),res(:,2:4));
     legend('Aggregate Welfare','Home Welfare','Foreign Welfare');
     XLabel('Home asset');
%      saveas(gcf,'welfare_asset1.fig')
% print(gcf,'-deps2c','-r800','welfare_asset1.eps')
end
fprintf('variable \t Conditional sum E \t Unconditional E\n')
disp([lgy_,num2str([CondM,erg_CondM(:,end)],10)])
 end