diff --git a/matlab/+hank/+bbeg/check_steady_state.m b/matlab/+hank/+bbeg/check_steady_state.m
new file mode 100644
index 0000000000000000000000000000000000000000..9e5a2903089cc545a90c76d2f859413b3165c764
--- /dev/null
+++ b/matlab/+hank/+bbeg/check_steady_state.m
@@ -0,0 +1,74 @@
+function [G_ss, F_ss] = check_steady_state(M_,sp,om,X,x_bar_f,dH,ind_var_exo, ind_eq_exo)
+   [sp.s, sp.N_a, sp.N_theta, sp.n_a_i, sp.n_theta_i] = set_state_matrix(sp.a, sp.theta);
+   [om.s, om.N_a, om.N_theta, om.n_a_i, sp.n_theta_i] = set_state_matrix(om.a, om.theta);
+   H_ = M_.heterogeneity(1);
+   % Build the values and coefficients of the individual policy functions
+   [x_bar] = compute_x_bar(x_bar_f, sp.s);
+   % Selection matrix for the states in the individual variables
+   iam1 = M_.heterogeneity(1).state_var;
+   if ind_var_exo
+      iam1 = setdiff(iam1,ind_var_exo);
+   end
+   N_x = size(x_bar,1);
+   [p] = compute_p(iam1, N_x);
+   a_bar = p*x_bar;
+   [Mu] = compute_Mu(sp);
+   % Compute Φ and ̃Φ with useful basis matrices
+   Phi = speye(sp.N_theta);
+   Phi_tilde = speye(sp.N_theta);
+   n_a = size(sp.a,2);
+   B_tilde_bar = cell(n_a,1);
+   for i=n_a:-1:1
+      B_tilde_bar{i} = funbas(sp.fbase{i},a_bar(i,:)')';
+      Phi = kron(funbas(sp.fbase{i}, om.a{i})', Phi);
+      Phi_tilde = kron(funbas(sp.fbase{i}, sp.a{i})', Phi_tilde);
+   end
+   % Compute ̃Φᵉ
+   N_sp = size(sp.s,1);
+   hadamard_prod = ones(sp.N_a,N_sp);
+   prod_a_1_to_jm1 = 1;
+   prod_a_jp1_to_n_a = sp.N_a;
+   for j=1:n_a
+      prod_a_jp1_to_n_a = prod_a_jp1_to_n_a/sp.n_a_i(j);
+      hadamard_prod = hadamard_prod .* kron(kron(ones(prod_a_1_to_jm1,1), B_tilde_bar{j}), ones(prod_a_jp1_to_n_a,1));
+      prod_a_1_to_jm1 = prod_a_1_to_jm1*sp.n_a_i(j);
+   end
+   if size(Mu,1) == 1 || size(Mu,2) == 1
+      Phi_tilde_e = kron(hadamard_prod, Mu);
+   else
+      Phi_tilde_e = kron(hadamard_prod, ones(sp.N_theta,1)).*kron(ones(sp.N_a), Mu);
+   end
+   % LU decomposition of Phi_tilde 
+   [L_Phi_tilde,U_Phi_tilde,P_Phi_tilde] = lu(Phi_tilde);
+   % Computing x_bar^#
+   x_bar_dash = x_bar/U_Phi_tilde;
+   x_bar_dash = x_bar_dash/L_Phi_tilde;
+   x_bar_dash = x_bar_dash*P_Phi_tilde;
+   % Computing residuals
+   % Aggregate model
+   Ix = x_bar_dash*Phi*dH;
+   [G_ss, y, x] = G(M_, X, Ix);
+   % Heterogenous model
+   n_a = size(sp.n_a_i,2);
+   N_sp = size(sp.s,1);
+   ind_var_endo = setdiff(1:H_.endo_nbr, ind_var_exo);
+   x_h = zeros(H_.endo_nbr,1);
+   xe_h = zeros(H_.endo_nbr,1);
+   theta_h = zeros(H_.exo_nbr,1);
+   % Compute the relevant subset of Y for individual policy functions iY, which
+   % designates the relevant indices in the vector [x_,x,xe,θ,X_,X,Xe,Θ]. Note
+   % that the latter also includes auxiliary variables
+   F_ss = zeros(H_.endo_nbr-numel(ind_eq_exo),N_sp);
+   for j=1:N_sp
+      if isempty(ind_var_exo) && isempty(ind_eq_exo)
+         theta_h = sp.s(j, n_a+1:end);
+         am1_h = sp.s(j, 1:n_a);
+      else
+         am1_h = sp.s(j, :);
+         x_h(ind_var_exo) = sp.s(j,n_a+1:end);
+      end
+      x_h(ind_var_endo) = x_bar(:, j);
+      xe_h(ind_var_endo) = x_bar_dash*Phi_tilde_e(:, j);
+      [F_ss(:,j)] = F(M_, y, x, theta_h, am1_h, x_h, xe_h, ind_eq_exo);
+   end
+end
\ No newline at end of file
diff --git a/matlab/+hank/+bbeg/het_stoch_simul.m b/matlab/+hank/+bbeg/het_stoch_simul.m
index 66c3ee3226c6f63279483bf3c1a7ca3f75098b86..ba7be25db796e9cccad5feb05893b3523b784436 100644
--- a/matlab/+hank/+bbeg/het_stoch_simul.m
+++ b/matlab/+hank/+bbeg/het_stoch_simul.m
@@ -14,7 +14,6 @@ function [X_hat, x_hat_dash, om_hat, Phi] = het_stoch_simul(T,M_,sp,om,X,x_bar_f
    a_bar = p*x_bar;
    [Mu] = compute_Mu(sp);
    % Compute Φ and ̃Φ with useful basis matrices
-   n_a = size(sp.a,2);
    Phi = speye(sp.N_theta);
    Phi_tilde = speye(sp.N_theta);
    n_a = size(sp.a,2);
@@ -145,7 +144,7 @@ function [X_hat, x_hat_dash, om_hat, Phi] = het_stoch_simul(T,M_,sp,om,X,x_bar_f
    [P] = compute_p(iAm1, N_X);
    X_hat_m1 = zeros(N_X,1);
    E = zeros(T*N_Theta,1);
-   E(1:N_Theta) = 1;
+   E(1:N_Theta) = 100*sqrt(M_.Sigma_e(1));
    iY = iY-N_het;
    % Now iY contains the indices of relevant aggregate variables for individual
    % policy functions in the vector [X_,X,Xe,Θ]. Note that the latter includes
diff --git a/matlab/+hank/+bbeg/private/F.m b/matlab/+hank/+bbeg/private/F.m
new file mode 100644
index 0000000000000000000000000000000000000000..107ed2c8053135c6a6ff3c34c21de585c9b046c9
--- /dev/null
+++ b/matlab/+hank/+bbeg/private/F.m
@@ -0,0 +1,22 @@
+function [F_ss] = F(M_, y, x, theta_h, am1_h, x_h, xe_h, ind_eq_exo)
+   H_ = M_.heterogeneity(1);
+   % --- Build the extended endogenous variable vector at the steady state --- %
+   yh = zeros(3*H_.endo_nbr,1);
+   % - State variables - %
+   iam1 = H_.state_var;
+   yh(iam1) = am1_h;
+   % - Contemporaneous endogenous variables -
+   ix = H_.endo_nbr+1:2*H_.endo_nbr;
+   yh(ix) = x_h;
+   % - Forward endogenous variables - %
+   ixe = H_.endo_nbr+ix;
+   yh(ixe) = xe_h;
+   % - Exogenous variables - %
+   xh = theta_h;
+   % --- Compute the first-order derivatives of F --- %
+   [F_ss,~,~] = feval([M_.fname '.sparse.dynamic_het1_resid'], y, x, M_.params, [], yh, xh, []);
+   if ~isempty(ind_eq_exo)
+      ind_eq = setdiff(1:size(F_ss,1),ind_eq_exo);
+      F_ss = F_ss(ind_eq,:);
+   end
+end
diff --git a/matlab/+hank/+bbeg/private/G.m b/matlab/+hank/+bbeg/private/G.m
new file mode 100644
index 0000000000000000000000000000000000000000..d0ab769ec45271a42b4f0f0be31ba5c15b18e42d
--- /dev/null
+++ b/matlab/+hank/+bbeg/private/G.m
@@ -0,0 +1,36 @@
+function [G_ss, y, x] = G(M_, X, Ix)
+
+   % --- Build the extended endogenous variable vector at the steady state --- %
+   y = zeros(3*M_.endo_nbr,1);
+   % - Variables resulting from an aggregation operator - %
+   iIx_ = [];
+   ix = [];
+   for i=1:size(M_.heterogeneity_aggregates,1)
+      if (M_.heterogeneity_aggregates{i,1} == 'sum')
+         iIx_(end+1) = M_.aux_vars(M_.heterogeneity_aggregates{i,2}).endo_index;
+         ix(end+1) = M_.heterogeneity_aggregates{i,3};
+      end
+   end
+   yagg = Ix(ix);
+   % t-1 variables
+   iX_ = 1:M_.endo_nbr;
+   % Exclude the indices of auxiliary variables stemming from aggregation
+   % operators, e.g. SUM
+   iX_ = setdiff(iX_, iIx_);
+   y(iX_) = X;
+   y(iIx_) = yagg;
+   % t variables
+   iX = M_.endo_nbr+iX_;
+   iIx = M_.endo_nbr+iIx_;
+   y(iX) = X;
+   y(iIx) = yagg;
+   % t+1 variables
+   iXe = M_.endo_nbr+iX;
+   iIxe = M_.endo_nbr+iIx;
+   y(iXe) = X;
+   y(iIxe) = yagg;
+   % - Exogenous variables - %
+   x = zeros(M_.exo_nbr,1);
+   % --- Compute the first-order derivatives of G --- %
+   [G_ss,~,~] = feval([M_.fname '.sparse.dynamic_resid'], y, x, M_.params, [], yagg);
+end
\ No newline at end of file
diff --git a/tests/hank/main.m b/tests/hank/main.m
index 347ca732efa5a3814e813bbc1cedc99b59cdf7d6..e9f014edf7a1e81a030ba26fd64f7147a95d31d3 100644
--- a/tests/hank/main.m
+++ b/tests/hank/main.m
@@ -167,8 +167,8 @@ om.fbase = {fbase_om, fbase_theta};
 T = 300;
 Z = 1.;
 X = [R, W, C, Y, K, Z];
+[G_ss, F_ss] = hank.bbeg.check_steady_state(M_,sp,om,X,x_bar_f,dH,[],[]);
 [X_hat, x_hat_dash, om_hat, Phi] = hank.bbeg.het_stoch_simul(T,M_,sp,om,X,x_bar_f,Lambda,dH,[],[]);
-% [X_hat, x_hat_dash, om_hat, Phi, p, x_bar_dash, Lambda_a, Mu, Phi_a] = hank.bbeg.het_stoch_simul(T,M_,sp,om,X,x_bar_f,Lambda,dH);
 % --------------------------------------------------------------------- %
 % -------------- Responses to aggregate shocks ------------------------ %
 % --------------------------------------------------------------------- %
diff --git a/tests/hank_exo/main.m b/tests/hank_exo/main.m
index 0e519788f8dda59fc461861c18b60d802316f1a8..c4a4dd2f7f27f66463e2f50c226a82630a1b73d5 100644
--- a/tests/hank_exo/main.m
+++ b/tests/hank_exo/main.m
@@ -11,7 +11,7 @@ dynare ks.mod;
 % ------------------------- Steady-State using Rouwenhorst ------------------------------ %
 % --------------------------------------------------------------------------------------- %
 % Coarse grid for predetermined variables: a_sp in BBEG
-a_sp= curved_grid(59,0.,500,2.5);
+a_sp= curved_grid(60,0.,500,2.5);
 fbase_sp = fundef({'spli',a_sp,0,2});
 a_sp = funnode(fbase_sp);
 % Dense grid for predetermined variables: a_ω in BBEG 
@@ -150,6 +150,7 @@ X = [R, W, C, Y, K, Z];
 ie = find(matches(M_.heterogeneity(1).endo_names,'e')==1);
 ind_var_exo = [ie];
 ind_eq_exo = [3];
+[G_ss, F_ss] = hank.bbeg.check_steady_state(M_,sp,om,X,x_bar_f,dH,ind_var_exo, ind_eq_exo);
 [X_hat, x_hat_dash, om_hat, Phi] = hank.bbeg.het_stoch_simul(T,M_,sp,om,X,x_bar_f,Lambda,dH,ind_var_exo, ind_eq_exo);
 
 % --------------------------------------------------------------------- %
@@ -161,7 +162,8 @@ figure;
 plot(1:T, X_hat(iK,:));
 xlabel('$t$','Interpreter','latex');
 ylabel('$\widehat{K}_t$','Interpreter','latex');
-title('IRF');
+title('Aggregate capital response to a unit aggregate productivity shock');
+saveas(gcf,'irf_K.png')
 % Impact response of the individual function for savings to an aggregate TFP
 % shock
 ik = find(matches(M_.heterogeneity(1).endo_names,'k')==1);
@@ -195,4 +197,98 @@ for t=1:T
    dOm_da_dtheta{t}(:,2:end) = dOm_dtheta(:,2:end)-dOm_dtheta(:,1:end-1);
 end
 time_lapse_plot_with_slider_distribution_agg(a_om, dOm_da_dtheta, T);
-time_lapse_plot_with_slider_distribution(a_om, dOm_da_dtheta, T);
\ No newline at end of file
+time_lapse_plot_with_slider_distribution(a_om, dOm_da_dtheta, T);
+% Building up a GIF - Disaggregated Distribution
+fig = figure;
+im = cell(T,1);
+for t=1:T
+   distr_e = sum(dOm_da_dtheta{t}, 2);
+   plot(a_om, dOm_da_dtheta{t}(1,:) ./ distr_e(1), 'DisplayName', 'Lowest income');
+   hold on;
+   plot(a_om, dOm_da_dtheta{t}(4,:) ./ distr_e(4), 'DisplayName', 'Median income');
+   plot(a_om, dOm_da_dtheta{t}(7,:) ./ distr_e(7), 'DisplayName', 'Top income');
+   hold off;
+   title({'Asset distribution';['t = ',num2str(t)]});
+   legend('show');
+   xlim([-1, 200]);
+   ylim([0,0.015]);
+   drawnow
+   frame = getframe(fig);
+   im{t} = frame2im(frame);
+end
+filename = "distr_disaggregated.gif"; % Specify the output file name
+for idx = 1:T
+    [A,map] = rgb2ind(im{idx},256);
+    if idx == 1
+        imwrite(A,map,filename,"gif",LoopCount=Inf, ...
+                DelayTime=0.01)
+    else
+        imwrite(A,map,filename,"gif",WriteMode="append", ...
+                DelayTime=0.01)
+    end
+end
+% Building up a GIF - Aggregated Distribution
+fig = figure;
+im = cell(T,1);
+for t=1:T
+   distr_a = sum(dOm_da_dtheta{t}, 1);
+   plot(a_om, distr_a, 'DisplayName', 'All')
+   title({'Asset distribution';['t = ',num2str(t)]});
+   legend('show');
+   xlim([-1, 200]);
+   ylim([0,0.0075]);
+   drawnow
+   frame = getframe(fig);
+   im{t} = frame2im(frame);
+end
+filename = "distr_aggregated.gif"; % Specify the output file name
+for idx = 1:T
+    [A,map] = rgb2ind(im{idx},256);
+    if idx == 1
+        imwrite(A,map,filename,"gif",LoopCount=Inf, ...
+                DelayTime=0.01)
+    else
+        imwrite(A,map,filename,"gif",WriteMode="append", ...
+                DelayTime=0.01)
+    end
+end
+% Building up a GIF - Policy functions
+fig = figure;
+im = cell(T,1);
+imid = ceil(n_theta/2);
+ic = find(matches(M_.heterogeneity(1).endo_names,'c')==1);
+N_om = size(Phi,2);
+N_theta = n_theta;
+N_a = N_om/N_theta;
+cf_lowest = pol_c{1}(a_om);
+cf_median = pol_c{imid}(a_om);
+cf_top = pol_c{end}(a_om);
+for t=1:T/2
+   c_hat_t = reshape(x_hat_dash(:,ic,t).'*Phi, N_theta, N_a);
+   plot(a_om, cf_lowest'+c_hat_t(1,:), 'DisplayName', 'Lowest income');
+   hold on;
+   plot(a_om, cf_median'+c_hat_t(imid,:), 'DisplayName', 'Median income');
+   plot(a_om, cf_top'+c_hat_t(end,:), 'DisplayName', 'Top income');
+   hold off;
+   title({'Response of the consumption policy function to a TFP shock';['t = ',num2str(t)]});
+   xlabel('$a_{t-1}$','Interpreter','latex')
+   ylabel('$\overline{c} + \hat{c}_{t}$','Interpreter','latex')
+   legend('show');
+   % title('Impact response of the consumption individual policy function to TFP shock');
+   xlim([0, 150]);
+   ylim([0, 10]);
+   drawnow
+   frame = getframe(fig);
+   im{t} = frame2im(frame);
+end
+filename = "policy.gif"; % Specify the output file name
+for idx = 1:T/2
+    [A,map] = rgb2ind(im{idx},256);
+    if idx == 1
+        imwrite(A,map,filename,"gif",LoopCount=Inf, ...
+                DelayTime=0.05)
+    else
+        imwrite(A,map,filename,"gif",WriteMode="append", ...
+                DelayTime=0.05)
+    end
+end
\ No newline at end of file