% Find the stochastic steady state
%% Make sure that the steadystate file takes ys_ as global and assigns this
%% as steady state value
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREAMBLE
% name of model and loading of symbsolve produced information matrix
% clear klmx_d_f
clear all
verbose=1;
warning('off','MATLAB:declareGlobalBeforeUse')

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

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,:);
disp('Give vector of variances')
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
bSig=vech2(Sigma_m);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Msho=[Msh;zeros(matrices.neq-matrices.jumping,nsh)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HERE SHOULD START THE LOOP
zz=1;
maxiter = 500;

ys_old=ys_
ys_seq=ys_;
%%%%%%%%%%%%%%%%%%%
% first round
    [Pi,PHI,AA,BB,numpre,aqz,bqz,Qqz,Zqz,inz22,algo]=firstorder_stuff(namemod,matrices,0);
    [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);
    
    % Unconditional mean
    F4x=[zeros(size(M,2),size(bSig,1));F4];
    V=(eye(size(PHIh,1))-PHIh)\Mh*bSig;%+Mh*bSig;
    S=(eye(size(PHI,1))-(PHI))\(PHI2*V+F4x*bSig);
    ys_=ys_old+(Pi*S+Pi2*V);
%     Iteration
progress(0,'');

while norm(abs((ys_+1e-8)./(ys_old+1e-8))-1)>1e-2 & zz<maxiter
    global ys_
    ys_old=ys_
texto=['Doing iteration ',num2str(zz,'%i'),' out of ',num2str(maxiter,'%i')];
    progress(zz/maxiter,texto);
    [Pi,PHI,AA,BB,numpre,aqz,bqz,Qqz,Zqz,inz22,algo]=firstorder_stuff(namemod,matrices,0);
    [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);
    
    % Unconditional mean
        F4x=[zeros(size(M,2),size(bSig,1));F4];
    V=(eye(size(PHIh,1))-PHIh)\Mh*bSig;%+Mh*bSig;
    S=(eye(size(PHI,1))-(PHI))\(PHI2*V+F4x*bSig);
    ys_=ys_old+(Pi*S+Pi2*V);
    ys_seq(:,zz+1)=ys_;
    zz=zz+1;
end
