function []=impulse_responses(parameters,wage_grids_outputs,output_vf_large,dirname,...
            hcutoff,apctile,shock_t,ver)

KTlarge = parameters.KTlarge;
mu_l2_shift = parameters.mu_l2_shift;       
mu_l2 = parameters.mu_l2;         
TT = parameters.TT; 

vernum = ver;  % 1: v1, 2: v2, 3: u1, 4: u2
shock_v1=0;
shock_v2=0;
shock_u1=0;
shock_u2=0;

warning('off','MATLAB:MKDIR:DirectoryExists');
subdirname = strcat(dirname,'/t',num2str(shock_t));
mkdir(subdirname);
warning('on','MATLAB:MKDIR:DirectoryExists');

if vernum == 1
    shock_v1=-0.1;
    fname = '/v1';
elseif vernum == 2
    shock_v2=-0.1;
    fname = '/v2';
elseif vernum == 3
    shock_u1=-0.104;
    fname = '/u1';
elseif vernum == 4
    shock_u2=-0.104;
    fname = '/u2';
end

%% Simulating all responses
for k=1:size(KTlarge,1)
    for u=1:size(mu_l2_shift,1)
        parameters1 = parameters;
        parameters1.KT = KTlarge(k,:);
        parameters1.mu_l2  = mu_l2.*mu_l2_shift(u,:);
        simul_output = simul_all(parameters1,wage_grids_outputs,output_vf_large(k,u));
        simul_output_large(k,u)=simul_output;
    end
end

simul_output_noshock = simul_all_het(parameters,simul_output_large);
lowA = simul_output_noshock.A_simul(:,1)<=prctile(simul_output_noshock.A_simul(:,1),apctile);
higA = simul_output_noshock.A_simul(:,1)>=prctile(simul_output_noshock.A_simul(:,1),100-apctile);

for k=1:size(KTlarge,1)
    for u=1:size(mu_l2_shift,1)
        parameters1 = parameters;
        parameters1.KT = KTlarge(k,:);
        parameters1.mu_l2  = mu_l2.*mu_l2_shift(u,:);
        IR_output = simul_all_IR(parameters1,simul_output_large(k,u),wage_grids_outputs,output_vf_large(k,u),shock_u1,shock_u2,shock_v1,shock_v2,shock_t);
        IR_output_large(k,u)=IR_output;
    end
end

simul_output_IR = simul_all_het(parameters,IR_output_large);

%% Aggregate average responses by age
% Families with kids
C_diff  = (mean(simul_output_IR.C_simul)'-mean(simul_output_noshock.C_simul)')./mean(simul_output_noshock.C_simul)';
H1_diff = (mean(simul_output_IR.H1_simul)'-mean(simul_output_noshock.H1_simul)')./mean(simul_output_noshock.H1_simul)';
H2_diff = (mean(simul_output_IR.H2_simul)'-mean(simul_output_noshock.H2_simul)')./mean(simul_output_noshock.H2_simul)';
E1_diff = (mean((simul_output_IR.H1_simul>hcutoff))'-mean((simul_output_noshock.H1_simul>hcutoff))')./mean((simul_output_noshock.H1_simul>hcutoff))';
E2_diff = (mean((simul_output_IR.H2_simul>hcutoff))'-mean((simul_output_noshock.H2_simul>hcutoff))')./mean((simul_output_noshock.H2_simul>hcutoff))';
L1_diff = (mean(simul_output_IR.L1_simul)'-mean(simul_output_noshock.L1_simul)')./mean(simul_output_noshock.L1_simul)';
L2_diff = (mean(simul_output_IR.L2_simul)'-mean(simul_output_noshock.L2_simul)')./mean(simul_output_noshock.L2_simul)';
T1_diff = (mean(simul_output_IR.T1_simul)'-mean(simul_output_noshock.T1_simul)')./mean(simul_output_noshock.T1_simul)';
T2_diff = (mean(simul_output_IR.T2_simul)'-mean(simul_output_noshock.T2_simul)')./mean(simul_output_noshock.T2_simul)';
w1_diff = (mean(exp(simul_output_IR.lw_1_simul))'-mean(exp(simul_output_noshock.lw_1_simul))')./mean(exp(simul_output_noshock.lw_1_simul))';
w2_diff = (mean(exp(simul_output_IR.lw_2_simul))'-mean(exp(simul_output_noshock.lw_2_simul))')./mean(exp(simul_output_noshock.lw_2_simul))';

H2c_diff = zeros(size(C_diff));
Y2c_diff = zeros(size(C_diff));
for tau = 1:TT
    e2both  = (simul_output_IR.H2_simul(:,tau)>hcutoff).*(simul_output_noshock.H2_simul(:,tau)>hcutoff);
    H2c_diff(tau) = (mean(simul_output_IR.H2_simul(e2both==1,tau))'-mean(simul_output_noshock.H2_simul(e2both==1,tau))')./mean(simul_output_noshock.H2_simul(e2both==1,tau))';
    Y2c_diff(tau) = (mean(simul_output_IR.Y2_simul(e2both==1,tau))'-mean(simul_output_noshock.Y2_simul(e2both==1,tau))')./mean(simul_output_noshock.Y2_simul(e2both==1,tau))';
end 

Y1_diff = (mean(simul_output_IR.Y1_simul)'-mean(simul_output_noshock.Y1_simul)')./mean(simul_output_noshock.Y1_simul)';
Y2_diff = (mean(simul_output_IR.Y2_simul)'-mean(simul_output_noshock.Y2_simul)')./mean(simul_output_noshock.Y2_simul)';
Y_diff  = (mean(simul_output_IR.Y_simul)'-mean(simul_output_noshock.Y_simul)')./mean(simul_output_noshock.Y_simul)';
atY_diff = (mean(simul_output_IR.atY_simul)'-mean(simul_output_noshock.atY_simul)')./mean(simul_output_noshock.atY_simul)';

out1 = [C_diff,H1_diff,H2_diff,E1_diff,E2_diff,L1_diff,L2_diff,T1_diff,T2_diff,w1_diff,w2_diff,H2c_diff,Y1_diff,Y2_diff,Y_diff,atY_diff,Y2c_diff]';
means= [mean(simul_output_noshock.C_simul);...
        mean(simul_output_noshock.H1_simul);...
        mean(simul_output_noshock.H2_simul);...
        mean(simul_output_noshock.L1_simul);...
        mean(simul_output_noshock.L2_simul);...
        mean(simul_output_noshock.T1_simul);...
        mean(simul_output_noshock.T2_simul);...
        mean(simul_output_noshock.L1_simul./simul_output_noshock.H1_simul);...
        mean(simul_output_noshock.T1_simul./simul_output_noshock.H1_simul);...
        mean(simul_output_noshock.L2_simul./simul_output_noshock.H2_simul);...
        mean(simul_output_noshock.T2_simul./simul_output_noshock.H2_simul);...
        mean(simul_output_noshock.Y1_simul);...
        mean(simul_output_noshock.Y2_simul);...
        mean(simul_output_noshock.Y_simul);...
        mean(simul_output_noshock.atY_simul);...
        mean(simul_output_noshock.Y1_simul./simul_output_noshock.Y_simul);...
        mean(simul_output_noshock.Y2_simul./simul_output_noshock.Y_simul);...
        mean(simul_output_noshock.Y1_simul./(simul_output_noshock.Y1_simul+simul_output_noshock.Y2_simul));...
        mean(simul_output_noshock.Y2_simul./(simul_output_noshock.Y1_simul+simul_output_noshock.Y2_simul))];

% No kids
C_diff  = (mean(simul_output_IR.C_simul_nk)'-mean(simul_output_noshock.C_simul_nk)')./mean(simul_output_noshock.C_simul_nk)';
H1_diff = (mean(simul_output_IR.H1_simul_nk)'-mean(simul_output_noshock.H1_simul_nk)')./mean(simul_output_noshock.H1_simul_nk)';
H2_diff = (mean(simul_output_IR.H2_simul_nk)'-mean(simul_output_noshock.H2_simul_nk)')./mean(simul_output_noshock.H2_simul_nk)';
E1_diff = (mean((simul_output_IR.H1_simul_nk>hcutoff))'-mean((simul_output_noshock.H1_simul_nk>hcutoff))')./mean((simul_output_noshock.H1_simul_nk>hcutoff))';
E2_diff = (mean((simul_output_IR.H2_simul_nk>hcutoff))'-mean((simul_output_noshock.H2_simul_nk>hcutoff))')./mean((simul_output_noshock.H2_simul_nk>hcutoff))';
L1_diff = (mean(simul_output_IR.L1_simul_nk)'-mean(simul_output_noshock.L1_simul_nk)')./mean(simul_output_noshock.L1_simul_nk)';
L2_diff = (mean(simul_output_IR.L2_simul_nk)'-mean(simul_output_noshock.L2_simul_nk)')./mean(simul_output_noshock.L2_simul_nk)';
T1_diff = (mean(simul_output_IR.T1_simul_nk)'-mean(simul_output_noshock.T1_simul_nk)')./mean(simul_output_noshock.T1_simul_nk)';
T2_diff = (mean(simul_output_IR.T2_simul_nk)'-mean(simul_output_noshock.T2_simul_nk)')./mean(simul_output_noshock.T2_simul_nk)';
w1_diff = (mean(exp(simul_output_IR.lw_1_simul_nk))'-mean(exp(simul_output_noshock.lw_1_simul_nk))')./mean(exp(simul_output_noshock.lw_1_simul_nk))';
w2_diff = (mean(exp(simul_output_IR.lw_2_simul_nk))'-mean(exp(simul_output_noshock.lw_2_simul_nk))')./mean(exp(simul_output_noshock.lw_2_simul_nk))';

H2c_diff = zeros(size(C_diff));
for tau = 1:TT
    e2both  = (simul_output_IR.H2_simul_nk(:,tau)>hcutoff).*(simul_output_noshock.H2_simul_nk(:,tau)>hcutoff);
    H2c_diff(tau) = (mean(simul_output_IR.H2_simul_nk(e2both==1,tau))'-mean(simul_output_noshock.H2_simul_nk(e2both==1,tau))')./mean(simul_output_noshock.H2_simul_nk(e2both==1,tau))';
end 

out2=[C_diff,H1_diff,H2_diff,E1_diff,E2_diff,L1_diff,L2_diff,T1_diff,T2_diff,w1_diff,w2_diff,H2c_diff]';

% With Kids - only low assets 
C_diff  = (mean(simul_output_IR.C_simul(lowA,:))'-mean(simul_output_noshock.C_simul(lowA,:))')./mean(simul_output_noshock.C_simul(lowA,:))';
H1_diff = (mean(simul_output_IR.H1_simul(lowA,:))'-mean(simul_output_noshock.H1_simul(lowA,:))')./mean(simul_output_noshock.H1_simul(lowA,:))';
H2_diff = (mean(simul_output_IR.H2_simul(lowA,:))'-mean(simul_output_noshock.H2_simul(lowA,:))')./mean(simul_output_noshock.H2_simul(lowA,:))';
E1_diff = (mean((simul_output_IR.H1_simul(lowA,:)>hcutoff))'-mean((simul_output_noshock.H1_simul(lowA,:)>hcutoff))')./mean((simul_output_noshock.H1_simul(lowA,:)>hcutoff))';
E2_diff = (mean((simul_output_IR.H2_simul(lowA,:)>hcutoff))'-mean((simul_output_noshock.H2_simul(lowA,:)>hcutoff))')./mean((simul_output_noshock.H2_simul(lowA,:)>hcutoff))';
L1_diff = (mean(simul_output_IR.L1_simul(lowA,:))'-mean(simul_output_noshock.L1_simul(lowA,:))')./mean(simul_output_noshock.L1_simul(lowA,:))';
L2_diff = (mean(simul_output_IR.L2_simul(lowA,:))'-mean(simul_output_noshock.L2_simul(lowA,:))')./mean(simul_output_noshock.L2_simul(lowA,:))';
T1_diff = (mean(simul_output_IR.T1_simul(lowA,:))'-mean(simul_output_noshock.T1_simul(lowA,:))')./mean(simul_output_noshock.T1_simul(lowA,:))';
T2_diff = (mean(simul_output_IR.T2_simul(lowA,:))'-mean(simul_output_noshock.T2_simul(lowA,:))')./mean(simul_output_noshock.T2_simul(lowA,:))';
w1_diff = (mean(exp(simul_output_IR.lw_1_simul(lowA,:)))'-mean(exp(simul_output_noshock.lw_1_simul(lowA,:)))')./mean(exp(simul_output_noshock.lw_1_simul(lowA,:)))';
w2_diff = (mean(exp(simul_output_IR.lw_2_simul(lowA,:)))'-mean(exp(simul_output_noshock.lw_2_simul(lowA,:)))')./mean(exp(simul_output_noshock.lw_2_simul(lowA,:)))';

H2c_diff = zeros(size(C_diff));
for tau = 1:TT
    e2both  = (simul_output_IR.H2_simul(:,tau)>hcutoff).*(simul_output_noshock.H2_simul(:,tau)>hcutoff).*lowA;
    H2c_diff(tau) = (mean(simul_output_IR.H2_simul(e2both==1,tau))'-mean(simul_output_noshock.H2_simul(e2both==1,tau))')./mean(simul_output_noshock.H2_simul(e2both==1,tau))';
end

out3=[C_diff,H1_diff,H2_diff,E1_diff,E2_diff,L1_diff,L2_diff,T1_diff,T2_diff,w1_diff,w2_diff,H2c_diff]';

% with kids - high assets
C_diff  = (mean(simul_output_IR.C_simul(higA,:))'-mean(simul_output_noshock.C_simul(higA,:))')./mean(simul_output_noshock.C_simul(higA,:))';
H1_diff = (mean(simul_output_IR.H1_simul(higA,:))'-mean(simul_output_noshock.H1_simul(higA,:))')./mean(simul_output_noshock.H1_simul(higA,:))';
H2_diff = (mean(simul_output_IR.H2_simul(higA,:))'-mean(simul_output_noshock.H2_simul(higA,:))')./mean(simul_output_noshock.H2_simul(higA,:))';
E1_diff = (mean((simul_output_IR.H1_simul(higA,:)>hcutoff))'-mean((simul_output_noshock.H1_simul(higA,:)>hcutoff))')./mean((simul_output_noshock.H1_simul(higA,:)>hcutoff))';
E2_diff = (mean((simul_output_IR.H2_simul(higA,:)>hcutoff))'-mean((simul_output_noshock.H2_simul(higA,:)>hcutoff))')./mean((simul_output_noshock.H2_simul(higA,:)>hcutoff))';
L1_diff = (mean(simul_output_IR.L1_simul(higA,:))'-mean(simul_output_noshock.L1_simul(higA,:))')./mean(simul_output_noshock.L1_simul(higA,:))';
L2_diff = (mean(simul_output_IR.L2_simul(higA,:))'-mean(simul_output_noshock.L2_simul(higA,:))')./mean(simul_output_noshock.L2_simul(higA,:))';
T1_diff = (mean(simul_output_IR.T1_simul(higA,:))'-mean(simul_output_noshock.T1_simul(higA,:))')./mean(simul_output_noshock.T1_simul(higA,:))';
T2_diff = (mean(simul_output_IR.T2_simul(higA,:))'-mean(simul_output_noshock.T2_simul(higA,:))')./mean(simul_output_noshock.T2_simul(higA,:))';
w1_diff = (mean(exp(simul_output_IR.lw_1_simul(higA,:)))'-mean(exp(simul_output_noshock.lw_1_simul(higA,:)))')./mean(exp(simul_output_noshock.lw_1_simul(higA,:)))';
w2_diff = (mean(exp(simul_output_IR.lw_2_simul(higA,:)))'-mean(exp(simul_output_noshock.lw_2_simul(higA,:)))')./mean(exp(simul_output_noshock.lw_2_simul(higA,:)))';

H2c_diff = zeros(size(C_diff));
for tau = 1:TT
    e2both  = (simul_output_IR.H2_simul(:,tau)>hcutoff).*(simul_output_noshock.H2_simul(:,tau)>hcutoff).*higA;
    H2c_diff(tau) = (mean(simul_output_IR.H2_simul(e2both==1,tau))'-mean(simul_output_noshock.H2_simul(e2both==1,tau))')./mean(simul_output_noshock.H2_simul(e2both==1,tau))';
end

out4=[C_diff,H1_diff,H2_diff,E1_diff,E2_diff,L1_diff,L2_diff,T1_diff,T2_diff,w1_diff,w2_diff,H2c_diff]';

%% Elasticity table
out = [out1;out2;out3;out4];
save(strcat(subdirname,fname),'out','means')
calc_elasticities(out,means,dirname,shock_t,ver);


end