
dim=hor*tsh;
num_pall=1:1:dim;


%Set the shock to save
% disp('List of shocks')
% disp([num2str([1:nsh]'),matrices.shocky])
% not_to_kill=input('Choose the shock you want to keep: ')
% Mshall1=Mshall;
% Mshall1(1:nsh,1:nsh)=zeros(nsh,nsh);
% Mshall2=Mshall1;
% Mshall1(not_to_kill,not_to_kill)=1;
% CSigma=0;
% while size(CSigma,1)~=nsh && size(CSigma,2)~=nsh
%     fprintf('Need %i dimensional vector or matrix\n',nsh)
%     CSigma=input('Give the SQRT(covariance matrix) of the innovations to be simulated ==> ');
% end
% disp('----------------------------------------------------------')
% if any(size(CSigma)==1)
%     CSigma=diag(CSigma);
% end
% Sigma_mx=CSigma^2;

%choose shock to change


%set up a counter
cnt=0;

%set up empty super matrix for results
Y=NaN(neq,S2,S1);
Y2=NaN(neq,S2,S1);
YRE=Y2;
YRE2=Y2;

%set up gain
gain=0.01;
% fraction of periods before shock
inno=.80;
s=0;
SS=S1*S2;





h=waitbar(0,'Please Wait')
set(h,'Position',[496 330.4 288 100])
set(h,'Color',[.5 0 .5])


%Starting the simulation, first the external loop
for jj=1:S1
     
shock_seq2=NaN(nsh,S2);
ee=randn(nsh,S2);
shock_seq2=CSigma*ee;
CSigma_new=CSigma;
CSigma_new(new_s,new_s)=std_new;
shock_seq3=CSigma_new*ee;
    

    TR0=XXX(:,:,ext);
    G_SUB=(A1*TR0+A2);
    ep=CSigma*randn(nsh,1)/100;
    Y(:,1,jj)=G_SUB\Mshall*ep;%yRE(:,end)+randn(size(yRE,1),1)/1000;
    Y2(:,1,jj)=Y(:,1,jj);
    YRE(:,1,jj)=G1\Mshall*ep;
    YRE2(:,1,jj)=YRE(:,1,jj);
    %here I re-initialize each time the starting conditions
    %RET+randn(size(RET))/100;
    %R=eye(length(truenonzero));
    R0=R(:,:,ext);
    prf1=0;
    prf2=0;
    
    
    Pall1(:,:,1)=TR0;
    R1(:,:,1)=R0;
    Pall2(:,:,1)=TR0;
    R2(:,:,1)=R0;
    %star here the inner loop
    for t=2:S2
        
       if t<S2*inno
       epsilon=shock_seq2(:,t);
       y_in=Y(:,t-1,jj);
       [y_out,Ptmp,R1(:,:,t),prf1(t)]=learning_BIS(std_P,mean_P,y_in,R1(:,:,t-1),Pall1(:,:,t-1),Mshall,epsilon,A1,A2,A3,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
       Pall1(:,:,t)=Ptmp;
       Pall2(:,:,t)=Ptmp;
       R1(:,:,t)=R1(:,:,t);
       R2(:,:,t)=R1(:,:,t);
       Y(:,t,jj)=y_out;
       Y2(:,t,jj)=y_out;
       stuck=t;
       YRE(:,t,jj)=RET*YRE(:,t-1,jj)+G1\Mshall*epsilon;       
       YRE2(:,t,jj)=YRE(:,t,jj);  
       prf2(t)=prf1(t);
    if t>200;
    mean_P=mean(Y(:,:,t-200:t,jj),2);%mean(Pall,3);
%     std_P=std(reshape(Pall(njump,truenonzero,:),length(njump)*length(truenonzero),t)')';
    else
    mean_P=Y(:,t,jj);%RET;  
    end
      else
      epsilon=shock_seq3(:,t);
      epsilon2=epsilon;
% epsilon=zeros(nsh,1);
% epsilon2=epsilon;
      if t==ceil(S2*inno)
      epsilon2(new_s,:)=epsilon(new_s)+3*CSigma_new(new_s,new_s);
      end
     if t>200;
    mean_P1=mean(Y(:,t-200:t,jj),2);%mean(Pall,3);
%     std_P=std(reshape(Pall(njump,truenonzero,:),length(njump)*length(truenonzero),t)')';
    else
    mean_P1=Y(:,t,jj);%RET;  
     end
         if t>200;
    mean_P2=mean(Y2(:,t-200:t,jj),2);%mean(Pall,3);
%     std_P=std(reshape(Pall(njump,truenonzero,:),length(njump)*length(truenonzero),t)')';
    else
    mean_P2=Y2(:,t,jj);%RET;  
    end
      
      %case 1, only one shock left
      y_in=Y(:,t-1,jj);
      [y_out,Ptmp,R1(:,:,t),prf1(t)]=learning_BIS(std_P,mean_P1,y_in,R1(:,:,t-1),Pall1(:,:,t-1),Mshall,epsilon2,A1,A2,A3,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
      Pall1(:,:,t)=Ptmp;
      
      Y(:,t,jj)=y_out;
      
      
      y_in=Y2(:,t-1,jj);
      [y_out,Ptmp,R2(:,:,t),prf2(t)]=learning_BIS(std_P,mean_P2,y_in,R2(:,:,t-1),Pall2(:,:,t-1),Mshall,epsilon,A1,A2,A3,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
      Pall2(:,:,t)=Ptmp;
     Y2(:,t,jj)=y_out;
     YRE(:,t,jj)=RET*YRE(:,t-1,jj)+G1\Mshall*epsilon2;   
     YRE2(:,t,jj)=RET*YRE2(:,t-1,jj)+G1\Mshall*epsilon;   
         if t>200;
    mean_P=mean(Y(:,t-200:t,jj),2);%mean(Pall,3);
%     std_P=std(reshape(Pall(njump,truenonzero,:),length(njump)*length(truenonzero),t)')';
    else
    mean_P=Y(:,t,jj);%RET;  
    end
        
      s=s+1;
       end
    waitbar(s/SS,h,sprintf('Done %4.2f\n Use of Projection facility %4.2f\n for draw %1.1f',s/SS*100,sum(prf1)/S2*100,jj))    
            
    end
    prfJ(jj)=sum(prf1);
end
close(h);

block=ceil(S2*inno);
%get the average responses in both regimes
% for jj=1:S2
%     M_Y(:,jj)=mean(Y(:,jj,:),3);
%     M_Y2(:,jj)=mean(Y2(:,jj,:),3);
%     M_YRE(:,jj)=mean(YRE(:,jj,:),3);
%     M_YRE2(:,jj)=mean(YRE2(:,jj,:),3);
%     
% end

% IRF=M_Y-M_Y2;
% IRFRE=M_YRE-M_YRE2;

%get the irf under each regime
TOT_IRF=NaN(neq,S2,S1);
TOT_IRFRE=NaN(neq,S2,S1);
for jj=1:S1
 TOT_IRF(:,:,jj)=Y(:,:,jj)-Y2(:,:,jj);
 TOT_IRFRE(:,:,jj)=YRE(:,:,jj)-YRE2(:,:,jj);
end

%get the top and loer 2.5 % of the irf in any time t

%set the position
UP=round(97.5*S1/100);
DW=round(2.5*S1/100);

for jj=1:neq
    for tt=1:S2
        A=sort(TOT_IRF(jj,tt,:));
        UP_25(jj,tt)=A(:,:,UP);
        LOW_25(jj,tt)=A(:,:,DW);
        ARE=sort(TOT_IRFRE(jj,tt,:));
        UPRE_25(jj,tt)=ARE(:,:,UP);
        LOWRE_25(jj,tt)=ARE(:,:,DW);
   end
end

for jj=1:neq
    
    V_Y(jj,1)=mean(var(Y(jj,:,:),0,2));
    
    V_Y2(jj,1)=mean(var(Y2(jj,:,:),0,2));
    
    V_RE(jj,1)=mean(var(YRE(jj,:,:),0,2));
    V_RE2(jj,1)=mean(var(YRE2(jj,:,:),0,2));
end



