% initialization of time, RE equilibrium and LS learning vector of
% variables

hor=input('Length of simulation ==>');
if isempty(hor)
    hor=3000;
end
randn('seed',1999)
t=1;

yRE=RET*y_ss*0+G1\Mshall*CSigma*randn(nsh,1);

%% training period for initial conditions of learning based on RE

train=100;
for t=2:train;
yRE(:,t)=RET*yRE(:,t-1)+G1\Mshall*CSigma*randn(nsh,1);
end
%% extract PC of states;
states_all=yRE(truenonzero,:);
nfactors=nsh+1;rank(covy(truenonzero,truenonzero))
FAVAR=0;
LAMBDA=1;
FAVAR=input('USE PC to shrink regressors?[1/0] ==> ');% initialization of the PC decomposition for YW learning; if LAMBDA=0, all states are used
if FAVAR
[fact,LAMBDA]=extract(states_all',nfactors-1);
% LAMBDA=1;
LAMBDA=LAMBDA/size(states_all,1);
end
% Use only LAMBDA to reduce the states to PC
%%% use only the last two of training period
% yRE=yRE(:,end-1:end);
y=yRE;
t=1;
%% INITIALIZATION OF COV MATRIX AND TRANSITION MATRIX (CHOOSE RE EQUILIBRIUM OR ELSE)

Pall=RET;%+randn/1000;eye(neq)/3;
R=eye(length(truenonzero));%*mean(diag(Sigma_mx)); covy(truenonzero,truenonzero);
prf=0;% projection facility count;
converg=vec(Pall(njump,truenonzero,t)-PHI2(njump,truenonzero));
% njump=1:4;
% picknojumps=1;
y_record=y;
%  h=waitbar(0,'Please Wait')
%  set(h,'Position',[496 330.4 288 100])
%  set(h,'Color',[1 1 1])
 progress('init'); %initializes the progress bar with title 'Please wait...'
%  progress('init',sprintf('Please Wait\n')); %initializes the progress bar with title TITLE
shock_seq=CSigma*randn(nsh,hor);
% mute_sh=setdiff([1:nsh],inno);
mean_y=y;
y_YW=yRE;
A_YW=Pall(njump,truenonzero);
prf_YW=prf;
converg_YW=converg;
gain=0.5;
states_YW=LAMBDA'*y_YW(truenonzero,:);
Z=zeros(length(njump),size(states_YW,1));
V=zeros(size(states_YW,1),size(states_YW,1));
Z(:,:,1)=gain*y_YW(njump,2:end)*states_YW(:,1:end-1)';
V(:,:,1)=  gain*states_YW(:,2:end)*states_YW(:,2:end)'+gain*(1-gain)*states_YW(:,1:end-1)*states_YW(:,1:end-1)';

std_y=sqrt(vv);

std_P=ones(length(njump)*length(truenonzero),1).*(vec(RET(njump,truenonzero))+(.01)); %initialization of standard error of parameter estimates
mean_P=y_ss(truenonzero);%RET(njump,truenonzero);
for t=2:hor;
        cnv=norm((converg(:,t-1)),2);
        cnv_YW=norm((converg_YW(:,t-1)),2);
%          waitbar(t/hor,h,sprintf('Done %4.2f; Use Proj.Fac.= %4.2f;\n ||estimates-RE||_{2}= %4.2f \nUse Proj.Fac. YW= %4.2f;\n ||estimates YW-RE||_{2}= %4.2f ',t/hor*100,sum(prf(1:t-1))/t*100,cnv,sum(prf_YW(1:t-1))/t*100,cnv_YW));
% msg=sprintf('Done %4.2f; Use Proj.Fac.= %4.2f; ||estimates-RE||_{2}= %4.2f  Use Proj.Fac. YW= %4.2f; ||estimates YW-RE||_{2}= %4.2f \n',t/hor*100,sum(prf(1:t-1))/t*100,cnv,sum(prf_YW(1:t-1))/t*100,cnv_YW);
    progress(t/hor,sprintf('Done %3.1f; Use PF= %3.1f; ||LS-RE||_{2}= %4.2f ',t/hor*100,sum(prf(1:t-1))/t*100,cnv)); %sets the RATE and the TITLE

    gain=1/(1+t^(0.85));
    epsilon=shock_seq(:,t);

    
    
    y_in=y(:,t-1);
    
    
    [y_out,Ptmp,R(:,:,t),prf(t)]=learning_BIS(std_P,mean_P,y_in,R(:,:,t-1),Pall(:,:,t-1),Mshall,epsilon,A1,A2,A3,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
    if t<1000
    [y_YW,A_YW(:,:,t),LAMBDA,Z,V,prf_YW(t)]=learning_BIS_YW_zm(Z,V,std_y,y_YW,A_YW(:,:,t-1),Mshall,t,epsilon,A1,A2,A3,LAMBDA,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
    else
     A_YW(:,:,t)=A_YW(:,:,t-1);   
     prf_YW(t)=0;
    end
    if FAVAR==1;
    [fact,LAMBDA]=extract(y_YW(truenonzero,:)',nfactors-1);
    y_YW(:,t)=y_YW(:,t-1);
    Z(:,:,t)=Z(:,:,t-1);
    V(:,:,t)=V(:,:,t-1);
    end
    
    Pall(:,:,t)=Ptmp;

    
    y(:,t)=y_out;
    
    if t>200;
    mean_P=(mean(y(truenonzero,t-200:t),2));%mean(Pall,3);
    std_P=std(reshape(Pall(njump,truenonzero,:),length(njump)*length(truenonzero),t)')';
    else
      mean_P=y(truenonzero,t);%RET;  
    end
    
    
    
    Trans=(A1*Ptmp+A2)\(-A3);
    converg(:,t)=vec(Trans(njump,truenonzero)-RET(njump,truenonzero));
     Trans_YW=(A1*picknojumps*A_YW(:,:,end)*picknonzero+A2)\(-A3);
    converg(:,t)=vec(Trans(njump,truenonzero)-RET(njump,truenonzero));
    converg_YW(:,t)=vec(Trans_YW(njump,truenonzero)-RET(njump,truenonzero));
    yRE(:,t)=RET*yRE(:,t-1)+G1\Mshall*epsilon;
    y_record(:,t)=y(:,t);
end
% close(h)
% progress('close')
TransF=(A1*Pall(:,:,end)+A2)\(-A3);
TransF_YW=(A1*picknojumps*A_YW(:,:,end)*picknonzero+A2)\(-A3);
G1LS=(A1*Pall(:,:,end)+A2);
% MF=(A1*Pall(:,:,end)+A2)\[Msh;zeros(neq-nsh,nsh)];
% MRE=G1\[Msh;zeros(neq-nsh,nsh)];


%% rearrange variables names

disp('Response to shocks: LS/RE')
mytable(new_varnames(njump,:),new_varnames([1:nsh],:),[TransF(njump,1:nsh)./RET(njump,1:nsh)]);
disp('----------------------------------------------------------------------------------------------------------')

disp('Response to shocks: YW/RE')
mytable(new_varnames(njump,:),new_varnames([1:nsh],:),[TransF_YW(njump,1:nsh)./RET(njump,1:nsh)]);
disp('----------------------------------------------------------------------------------------------------------')


figure(1);
trim=floor(.25*hor);
subplot(2,1,1)
plot(1:size(converg,2)-trim,converg(:,trim+1:end),'>')
title('convergence LS')
subplot(2,1,2)
plot(1:size(converg,2),prf,'o')
title('use of projection facility LS')

figure(2);
subplot(2,1,1)
plot(1:size(converg,2)-trim,y_record(:,trim+1:end),'r--',1:size(converg,2)-trim,yRE(:,trim+1:end),'k-.')
title('LS-learning (--) vs RE (-.)')
subplot(2,1,2)
plot(1:size(converg,2)-trim,(y_record(:,trim+1:end)-repmat(mean(y_record(:,trim+1:end),2),1,size(y_record(:,trim+1:end),2))-yRE(:,trim+1:end)))
title('LS learning - RE ')


% figure(3);
% trim=floor(.25*hor);
% subplot(2,1,1)
% plot(1:size(converg_YW,2)-trim,converg_YW(:,trim+1:end),'>')
% title('convergence YW')
% subplot(2,1,2)
% plot(1:size(converg_YW,2),prf_YW,'o')
% title('use of projection facilityYW ')
% 
% figure(4);
% subplot(2,1,1)
% plot(1:size(converg_YW,2)-trim,y_YW(:,trim+1:end),'r--',1:size(converg_YW,2)-trim,yRE(:,trim+1:end),'k-.')
% title('YW-learning (--) vs RE (-.)')
% subplot(2,1,2)
% plot(1:size(converg_YW,2)-trim,(y_YW(:,trim+1:end)-repmat(mean(y_YW(:,trim+1:end),2),1,size(y_YW(:,trim+1:end),2))-yRE(:,trim+1:end)))
% title('YW learning - RE ')
%%
if do_irf
    y_IRF(:,1)=y_ss*0;
    y_YW_IRF(:,1)=y_ss*0;
yRE_IRF(:,1)=y_ss*0;
epsilon=zeros(nsh,1);

epsilon(inno)=1;
R_IRF=R(:,:,end);
Pall_IRF=Pall(:,:,end);
A_YW_IRF=A_YW(:,:,end);
prf_IRF=0;
prf_YW_IRF=0;
mean_y_IRF=y_YW_IRF*0;
    for t=2:irf_length;
        
        
        
        y_in=y_IRF(:,t-1);
        
        
        [y_out,Ptmp,R_IRF(:,:,t),prf_IRF(t)]=learning_BIS(std_P,mean_P,y_in,R_IRF(:,:,t-1),Pall_IRF(:,:,t-1),Mshall,epsilon,A1,A2,A3,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
        [y_YW_IRF,A_YW_IRF(:,:,t),LAMBDA,Z,V,prf_YW_IRF(t)]=learning_BIS_YW_zm(Z,V,std_y,y_YW_IRF,A_YW_IRF(:,:,t-1),Mshall,t,epsilon,A1,A2,A3,LAMBDA,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);

        Pall_IRF(:,:,t)=Ptmp;
        y_IRF(:,t)=y_out;
        
        
        
        
        yRE_IRF(:,t)=RET*yRE_IRF(:,t-1)+G1\Mshall*epsilon;
        epsilon=epsilon*0;
    end
    
    disp('LIST OF VARIABLES')
    [num2str([1:neq]'),new_varnames]
    nn=input('Give vector of variables # for IRF ==>');
    rows=floor(sqrt(length(nn)));
    cols=ceil(length(nn)/rows);
   
    figure(5);% IFR
    
    subplot(3,1,1)
    plot(1:irf_length,y_IRF(nn,:))
    title('LS-learning  ')
    subplot(3,1,2)
    plot(1:irf_length,yRE_IRF(nn,:))
    title('RE solution ')
    legend(new_varnames(nn,:))
     subplot(3,1,3)
        plot(1:irf_length,y_IRF(nn,:)-yRE_IRF(nn,:))
          legend(new_varnames(nn,:))
    title('LS-RE solutions ')
    
        figure(6);% IFR
    
    subplot(3,1,1)
    plot(1:irf_length,y_YW_IRF(nn,:))
    title('YW-learning  ')
    subplot(3,1,2)
    plot(1:irf_length,yRE_IRF(nn,:))
    title('RE solution ')
    legend(new_varnames(nn,:))
     subplot(3,1,3)
        plot(1:irf_length,y_YW_IRF(nn,:)-yRE_IRF(nn,:),'-r')
    title('YW-RE solutions ')

    % plot(1:size(converg,2)-trim,(y_record(:,trim+1:end)-repmat(mean(y_record(:,trim+1:end),2),1,size(y_record(:,trim+1:end),2))-yRE(:,trim+1:end)))
    % title('LS learning - RE ')
end
clearvars y_YW Z V