%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code is to draw impulse responses not only against standard shocks but also expectation shock %
% @Ippei Fujiwara and Heedon Kang(December 28, 2007)                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%warning off MATLAB:dividebyzero
warning off all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition of Parameters                                                                           %
% setting (all you need to specify are just 4 parameters below)                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Determine which shock is expectation shock
% 2. Period for expectation is to be realized. Period must be larger than 0
% 3. Length of period to simulate
% 4. Impulse responses in level difference(impulse_index=1) or percentage deviation(impulse_index=2)
expshock='z';
period=4;
simperiod=20;
impulse_index=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start of the main program                                                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;                                    % To close all the previous figures that are unnecessary
[info,jacobimat]=stoch_simul_es(var_list_);   % To obtain Jacobian matrix
close all;                                    % To close all the unnecessary figures from simulation
%clc;                                         % To clears the command window and homes the cursor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note that matrix iy_ notes that where a variable appears in the model                              %                              
% When there is only one lag and one lead, the matrix has 3 rows                                     % 
% The first rows describe which variables appear with a lag,                                         %
% the second row which variable appear in the model at the current period,                           %
% the third one, which variable appears with a lead                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numx=size(iy_,2);                             % Number of variables at the second row
nums=size(Sigma_e_);                          % (1*2) matrix : [Number of shocks, Number of shocks]                         
numz=numx-nums;                               % Number of variables not related with shocks
numeq=size(jacobimat,1);                      % Number of equations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix difference equation must take a form like A*Z(+1)+B*Z+C*Z(-1)=0                             %
% namely max lead as well as max lag must be one.                                                    %
%                                                                                                    %
% Rational Expectation Model Solution Method                                                         %
%  1. alpha0*E_t(Z(t+1)-Z*) + alpha1*(Z(t)-Z*) + alpha2*(Z(t-1)-Z*)                                  %
%     + beta0*(S(t+1)-S*) + beta1*(S(t)-S*)                         = 0                              %
%  2. S(t) = S* + P*(S(t-1)-S*) + C*epsilon(t)                                                       %
%                                                                                                    %
% What we would like to obtain here are alpha0, alpha1, alpha2, beta0, beta_1, P and C               %
% In Dynare, jacobimat is organized as A*X(-1)+B*X+C*X(+1)+D*shocks, where X=(Z,S)                   %
% in which zero column is eliminated                                                                 %
% In Dynare, jacobian matrix is calculated only with non-zero element at iy_ matrix                  %
% Thus, we need to add zero vectors and to produce new jacobian matrices                             %
% related to lead, current, and lag variables  matching with iy_ matrix                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(iy_,1)~=3
    error('too many lags or too many leads')
end

ABC=jacobimat(:,1:(end-nums));                % Divide Jacobian Matrix into two parts such as ABC and DDD
DDD=jacobimat(:,(end-nums+1):end);             
shockrows=find(sum(DDD,2)~=0);                % Where shocks are located in the Jacobian
varrows=[1:numeq];                            % (1*N) vector, N is the number of equations
varrows=setdiff(varrows, shockrows);          % SETDIFF(A,B) when A and B are vectors returns the values
                                              % in A that are not in B.
zero_A=find(iy_(1,:)==0);                     % Find cells with 0 from the first row of iy_ matrix
zero_B=find(iy_(2,:)==0);                     % Find cells with 0 from the second row of iy_ matrix
zero_C=find(iy_(3,:)==0);                     % Find cells with 0 from the third row of iy_ matrix

ncolumnA=numx-size(zero_A,2);                 % Number of cell with non-zero element in A matrix  
ncolumnB=numx-size(zero_B,2);                 % Number of cell with non-zero element in B matrix  
ncolumnC=numx-size(zero_C,2);                 % Number of cell with non-zero element in C matrix  

AAA=jacobimat(:,[1:ncolumnA]);                % Take a part of Jacobian matrix related to A, B, and C matrix
BBB=jacobimat(:,[ncolumnA+1:ncolumnA+ncolumnB]);
CCC=jacobimat(:,[ncolumnA+ncolumnB+1:ncolumnA+ncolumnB+ncolumnC]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This part is fixed from Fujiwara(2006), since there was a bug regarding if-procedure.              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:size(zero_A,2);
    AAA=[AAA(:,1:zero_A(i)-1) zeros(numeq,1) AAA(:,zero_A(i):end)];
end
clear i;                                      % In order to use 'i' at the next step, we need to dump it out

for i=1:size(zero_B,2);
    BBB=[BBB(:,1:zero_B(i)-1) zeros(numeq,1) BBB(:,zero_B(i):end)];
end
clear i;

for i=1:size(zero_C,2);
    CCC=[CCC(:,1:zero_C(i)-1) zeros(numeq,1) CCC(:,zero_C(i):end)];
end
clear i;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Divide AAA, BBB, and CCC matrix into two parts:                                                    %
% one with endogenous variables and the other with shocks                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AAA1=AAA([varrows],:);
AAA2=AAA([shockrows],:);
BBB1=BBB([varrows],:);
BBB2=BBB([shockrows],:);
CCC1=CCC([varrows],:);
CCC2=CCC([shockrows],:);

shockcolumns=find(sum(BBB2,1)~=0);
varcolumns=[1:numx];
varcolumns=setdiff(varcolumns, shockcolumns);

AAA11=AAA1(:,[varcolumns]);
AAA12=AAA1(:,[shockcolumns]);
AAA21=AAA2(:,[varcolumns]);
AAA22=AAA2(:,[shockcolumns]);

BBB11=BBB1(:,[varcolumns]);
BBB12=BBB1(:,[shockcolumns]);
BBB21=BBB2(:,[varcolumns]);
BBB22=BBB2(:,[shockcolumns]);

CCC11=CCC1(:,[varcolumns]);
CCC12=CCC1(:,[shockcolumns]);

%check
if sum(sum(CCC2))~=0
    error('C2 must be zero matrix')
end

if sum(sum(AAA12))~=0
    error('A12 must be zero matrix')
end

if sum(sum(AAA21))~=0
    error('A21 must be zero matrix')
end

if sum(sum(BBB21))~=0
    error('B21 must be zero matrix')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Since the order at the Jacobian matrix is different from one at the policy function,               %
% we need to reorder variables to be consistent with ghx and ghu                                     %
% Note : It is necessary to check why two orders are different at Dynare!!!!                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
order_var=dr_.order_var;
nnn=1;
for i=1:size(order_var,1);
    fvar=find(varcolumns==order_var(i));
    if fvar>0
        fffvar(nnn)=fvar;
        nnn=nnn+1;
    end
end
clear nnn i;

nnn=1;
for i=1:size(order_var,1);
    fshock=find(shockcolumns==order_var(i));
    if fshock>0
        fffshock(nnn)=fshock;
        nnn=nnn+1;
    end
end
clear nnn i;

CCC11=CCC11(:,[fffvar]);
BBB11=BBB11(:,[fffvar]);
AAA11=AAA11(:,[fffvar]);
CCC12=CCC12(:,[fffshock]);
BBB12=BBB12(:,[fffshock]);

AAA22=AAA22(:,[fffshock]);
BBB22=BBB22(:,[fffshock]);

% Compute alpha_0, alpha_1, alpha_2, beta_0, and beta_1
alpha_0=CCC11;
alpha_1=BBB11;
alpha_2=AAA11;
beta_0=CCC12;
beta_1=BBB12;

% Construct P matrix, this may be redundant
rhos=(-1)*sum(AAA22,1)./sum(BBB22,1);
PPP=zeros(nums);
for i=1:nums;
    PPP(i,i)=rhos(i);
end
clear i;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We now compute matrices A and B in                                                                 %
% Z=Zss+A*Z(-1)+B*S+D1*e                                                                             %
% S=Sss+P*S(-1)    +D2*e                                                                             %
% In Dynare                                                                                          %
% X=Xss+GHX*Xstate(-1)+GHU*e                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GHX=dr_.ghx;
GHU=dr_.ghu;

nstatic=dr_.nstatic;          % In Fujiwara(2006), the static variable is r 
npred=dr_.npred;              
nfwrd=dr_.nfwrd;              

if nstatic+npred+nfwrd~=numx
    error('number of variables is not consistent')     
end

GHX=[zeros(numx,nstatic) GHX zeros(numx,nfwrd)];

for i=1:size(varcolumns,2);
    ffvar=find(order_var==varcolumns(i));
    ffffvar(i)=ffvar;
end
clear i;
ffffvar=sort(ffffvar);

for i=1:size(shockcolumns,2);
    ffshock=find(order_var==shockcolumns(i));
    ffffshock(i)=ffshock;
end
clear i;
ffffshock=sort(ffffshock);

A=GHX([ffffvar],[ffffvar]);
C12=GHX([ffffvar],[ffffshock]);
P=GHX([ffffshock],[ffffshock]);
zeromatrix=GHX([ffffshock],[ffffvar]);
D1=GHU([ffffvar],:);
D2=GHU([ffffshock],:);

B=C12*inv(P);

%check
%P=PPP
if max(max(abs(P-PPP)))>1e-8
    error('matrices are not consistent')
end

%zeromatrix must consist of zeros
if max(max(abs(zeromatrix)))>1e-8
    error('elements of zeromatrix must be zero')
end

%C12-B*P=0
if max(max(abs(C12-B*P)))>1e-8
    error('matrices are not consistent')
end

%D1-B*D2=0
if max(max(abs(D1-B*D2)))>1e-8
    error('matrices are not consistent')
end

%alpha_0*A^2+alpha_1*A+alpha_2=0
if max(max(abs(alpha_0*A^2+alpha_1*A+alpha_2)))>1e-8
    error('matrices are not consistent')
end

%(beta_0+alpha_0*B)*P+beta_1+alpha_1*B+alpha_0*A*B=0
if max(max(abs((beta_0+alpha_0*B)*P+beta_1+alpha_1*B+alpha_0*A*B)))>1e-8
    error('matrices are not consistent')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Expectation shock process accounting                                                               %
% 1. Find out which shocks the expectation shocks are                                                %
% 2. Create new beta_0_new, beta_1_new, and P_new                                                    %
% 3. Create new B_new                                                                                %
%    ---> FindandcheckB.m                                                                            %
%       ---> findal.m                                                                                %
% 4. Make new D2_new                                                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Figure out which shocks the expectation shocks are %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lengthexp=size(expshock,2);
lengthlgy=size(lgy_,2);

expmat=zeros(lengthexp,lengthlgy);
for i=1:lengthexp;
    findletter=find(lgy_(:,i)==expshock(i));
    nn=size(findletter,1);
    for ii=1:nn;
        expmat(i,ii)=findletter(ii);
    end        
end
clear i nn ii;

sumexpmat=sum(expmat);
flagcolumns=find(sumexpmat(1,:)~=0); % Among the element of sumexpmat, find indices of nonzero elements
expmat=expmat(:,[flagcolumns]);

for i=1:size(expmat,2);
    findrows=find(expmat==expmat(1,i));
    if size(findrows,1)==size(expmat,1)
        shocknumber=expmat(1,i);
    end
end

findshock1=find(order_var==shocknumber);
target=find(ffffshock==findshock1);  %this is the place of expectational shock in P

alphas = [alpha_0, alpha_1, alpha_2];
betas = [beta_0, beta_1];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Create betas_new and P_new %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% zerobeta=zeros(size(beta_0,1),period);
zerobeta=zeros(rows(beta_0),period);
beta_0_new1=beta_0(:,[1:target]);
beta_0_new2=beta_0(:,[target+1:end]);
beta_0_new=[beta_0_new1 zerobeta beta_0_new2];

beta_1_new1=beta_1(:,[1:target]);
beta_1_new2=beta_1(:,[target+1:end]);
beta_1_new=[beta_1_new1 zerobeta beta_1_new2];

betas_new=[beta_0_new beta_1_new];

% create P_new matrix
P_new1=P([1:target],[1:target]);
P_new2=P([target+1:end],[target+1:end]);
P_new=zeros(size(beta_0,2)+period);
ematrix=zeros(period);
ematrix([2:end],[1:period-1])=eye(period-1);
P_new([1:target],[1:target])=P_new1;
P_new([target+1:target+1+period-1],[target+1:target+1+period-1])=ematrix;
P_new(target,target+period)=1;
P_new([target+period+1:end],[target+period+1:end])=P_new2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Create B_new               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[B_new] = FindandcheckB(A, P_new, alphas, betas_new);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. We need to make D2_new     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D2_new1=D2([1:target],[1:target]);
D2_new2=D2([target+1:end],[target+1:end]);
D2_new=zeros(size(beta_0,2)+period);
ssmatrix=eye(period)*D2(target,target);
D2_new([1:target],[1:target])=D2_new1;
D2_new([target+1:target+1+period-1],[target+1:target+1+period-1])=ssmatrix;
D2_new([target+period+1:end],[target+period+1:end])=D2_new2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Drawing impulse responses                                                                          %
% According to differencestatespace.tex (since this is not % deviation model),                       %
% Z=A*Z(-1)+B*S  (<--Z=Zss+A*Z(-1)+B*S+D1*e)                                                         %
% S=P*S(-1)+D2*e (<--S=Sss+P*S(-1)    +D2*e)                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
order_var_title=lgy_([order_var],:);
var_title=order_var_title([ffffvar],:);
shock_title=order_var_title([ffffshock],:);
numberfig=size(A,1)+1;
numberfigures=floor(numberfig/9)+1;

% steady state values for endogenous variables
ZSS1=ys_([order_var],:);
ZSS=ZSS1([ffffvar],:);

% case without any reverse in the future
e=zeros(size(P_new,1),simperiod);
e(target+1,1)=1;

S(:,1) = D2_new*e(:,1);
Z(:,1) = B_new*S(:,1);
for ii = 2:simperiod;
    S(:,ii)=P_new*S(:,ii-1)+D2_new*e(:,ii);
    Z(:,ii)=A*Z(:,ii-1)+B_new*S(:,ii);
end

% case with reverse in the future
e2=zeros(size(P_new,1),simperiod);
e2(target+1,1)=1;
e2(target,1+period)=-1;

S2(:,1) = D2_new*e2(:,1);
Z2(:,1) = B_new*S2(:,1);
for ii = 2:simperiod;
    S2(:,ii)=P_new*S2(:,ii-1)+D2_new*e2(:,ii);
    Z2(:,ii)=A*Z2(:,ii-1)+B_new*S2(:,ii);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% There is an option that we can choose.                                                             %
% Impulse responses in level difference(impulse_index=1) or percentage deviation(impulse_index=2)    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zeroline=zeros(1,simperiod);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Case with level difference  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if impulse_index==1;

    if numberfigures == 1
        figure(11)
        subplot(3,3,1)
        hold on
        plot(S(target,:),'b-', 'LineWidth', 2);
        plot(S2(target,:),'r-', 'LineWidth', 0.5);
        plot([zeroline(:)],'g:');
        title([shock_title(target,:)]);
        axis tight
        for i=1:size(A,1)
            subplot(3,3,i+1)
            hold on
            plot([Z(i,:)],'b-', 'LineWidth', 2);
            plot([Z2(i,:)],'r-', 'LineWidth', 0.5);
            plot([zeroline(:)],'g:');
            title([var_title(i,:)]);
            axis tight
        end
        suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
     else
        for ii=1:numberfigures
            if ii == 1
                figure(10+ii)
                subplot(3,3,1)
                hold on
                plot(S(target,:),'b-', 'LineWidth', 2);
                plot(S2(target,:),'r-', 'LineWidth', 0.5);
                plot([zeroline(:)],'g:');
                title([shock_title(target,:)]);
                axis tight
                for iii=1:8
                    subplot(3,3,iii+1)
                    hold on
                    plot([Z(iii,:)],'b-', 'LineWidth', 2);
                    plot([Z2(iii,:)],'r-', 'LineWidth', 0.5);
                    plot([zeroline(:)],'g:');
                    title([var_title(iii,:)]);
                    axis tight
                end
                suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
                clear iii
            elseif ii < numberfigures
                figure(10+ii)
                for iii=1:9
                    subplot(3,3,iii)
                    hold on
                    plot([Z(iii+8+(ii-2)*9,:)],'b-', 'LineWidth', 2);
                    plot([Z2(iii+8+(ii-2)*9,:)],'r-', 'LineWidth', 0.5);
                    plot([zeroline(:)],'g:');
                    title([var_title(iii+8+(ii-2)*9,:)]);
                    axis tight
                end
                suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
                clear iii
            else
                figure(10+ii) 
                for iii=1:(size(A,1)-8-(ii-2)*9)    
                    subplot(3,3,iii)
                    hold on
                    plot([Z(iii+8+(ii-2)*9,:)],'b-','LineWidth', 2);
                    plot([Z2(iii+8+(ii-2)*9,:)],'r-','LineWidth', 0.5);
                    plot([zeroline(:)],'g:');
                    title([var_title(iii+8+(ii-2)*9,:)]);
                    axis tight
                end
                suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
                clear iii      
            end
        end
    end
    clear i ii;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Case with percentage difference  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    for zz=1:size(Z,2);
        Z(:,zz)=Z(:,zz)./ZSS;
    end    
    for zz2=1:size(Z,2);
        Z2(:,zz2)=Z2(:,zz2)./ZSS;
    end
    
    if numberfigures == 1
        figure(11)
        subplot(3,3,1)
        hold on
        plot(S(target,:),'b-', 'LineWidth', 2);
        plot(S2(target,:),'r-', 'LineWidth', 0.5);
        plot([zeroline(:)],'g:');
        title([shock_title(target,:)]);
        axis tight
        for i=1:size(A,1)
            subplot(3,3,i+1)
            hold on
            plot([Z(i,:)],'b-', 'LineWidth', 2);
            plot([Z2(i,:)],'r-', 'LineWidth', 0.5);
            plot([zeroline(:)],'g:');
            title([var_title(i,:)]);
            axis tight
        end
        suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
     else
        for ii=1:numberfigures
            if ii == 1
                figure(10+ii)
                subplot(3,3,1)
                hold on
                plot(S(target,:),'b-', 'LineWidth', 2);
                plot(S2(target,:),'r-', 'LineWidth', 0.5);
                plot([zeroline(:)],'g:');
                title([shock_title(target,:)]);
                axis tight
                for iii=1:8
                    subplot(3,3,iii+1)
                    hold on
                    plot([Z(iii,:)],'b-', 'LineWidth', 2);
                    plot([Z2(iii,:)],'r-', 'LineWidth', 0.5);
                    plot([zeroline(:)],'g:');
                    title([var_title(iii,:)]);
                    axis tight
                end
                suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
                clear iii
            elseif ii < numberfigures
                figure(10+ii)
                for iii=1:9
                    subplot(3,3,iii)
                    hold on
                    plot([Z(iii+8+(ii-2)*9,:)],'b-', 'LineWidth', 2);
                    plot([Z2(iii+8+(ii-2)*9,:)],'r-', 'LineWidth', 0.5);
                    plot([zeroline(:)],'g:');
                    title([var_title(iii+8+(ii-2)*9,:)]);
                    axis tight
                end
                suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
                clear iii
            else
                figure(10+ii) 
                for iii=1:(size(A,1)-8-(ii-2)*9)    
                    subplot(3,3,iii)
                    hold on
                    plot([Z(iii+8+(ii-2)*9,:)],'b-','LineWidth', 2);
                    plot([Z2(iii+8+(ii-2)*9,:)],'r-','LineWidth', 0.5);
                    plot([zeroline(:)],'g:');
                    title([var_title(iii+8+(ii-2)*9,:)]);
                    axis tight
                end
                suptitle(['Expectation  Shock  to  ' shock_title(target,:) '(period  '  num2str(period) ')'])
                clear iii      
            end
        end
    end
    clear i ii;
end