diff --git a/matlab/+hank/check_steady_state_input.m b/matlab/+hank/check_steady_state_input.m
new file mode 100644
index 0000000000000000000000000000000000000000..c8b384f8a1460b397dc1b95880e949204459ee52
--- /dev/null
+++ b/matlab/+hank/check_steady_state_input.m
@@ -0,0 +1,347 @@
+% Copyright © 2025 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+% Computes the stochastic simulations of heterogeneus-agent models
+%
+% INPUTS
+% - M_       [structure] Matlab's structure describing the model
+% - options_ [structure] Matlab's structure describing the current options
+% - ss       [structure] Matlab's structure with the model steady-state
+%                       information. It contains the following fields:
+%    - pol [structure] Matlab's structure containing the policy functions
+%                      discretization
+%       - pol.grids [structure]: Matlab's structure containing the nodes of the
+%                                state grids as column vectors.
+%       - pol.values [structure]: Matlab's structure containing the policy
+%                                 function as matrices. Row indices and column
+%                                 indices follow the lexicographic order
+%                                 specified in pol.shocks and pol.states
+%       - pol.shocks [array]: a list containing the order of shock/innovation
+%                             variables used for the rows of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement (resp. varexo(heterogeneity=)
+%                             statement) for discretizes AR(1) processes (resp.
+%                             discretized gaussian i.i.d innovations).
+%       - pol.states [array]: a list containing the order of state variables
+%                             used for the columns of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement.
+%    - shocks [structure]: Matlab's stucture describing the discretization of
+%                          individual shocks or innovations:
+%       - shocks.grids [structure]: Matlab's structure containing the nodes of
+%                                   the shock grids
+%       - shocks.Pi [structure]: Matlab's structure containing the Markov
+%                                matrices if the shock processes are discretized
+%                                AR(1) processes. The field should be absent
+%                                otherwise.
+%       - shocks.w [structure]: Matlab's structure containing the Gauss-Hermite
+%                               weights if the i.i.d gaussian innovation
+%                               processes are discretized
+%    - d [structure]: Matlab's structure describing the steady-state
+%                     distribution
+%       - d.grids [structure]: structure containing the states grids as column
+%                              vectors. If one of the states grid is not
+%                              specified, it is assumed to be identical to
+%                              its policy-function value stored in pol.grids
+%       - d.hist [matrix]: the histogram of the distribution as a matrix with
+%                          shocks/innovation indices as rows and state indices
+%                          as columns. Row and column indices follow the
+%                          lexicographic order specified in d.shocks and
+%                          d.states
+%       - d.shocks [array]: a list containing the order of shock/innovation
+%                           variables used for the rows of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.shocks.
+%       - d.states [array]: a list containing the order of state variables used
+%                           for the columns of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.states.
+%    - agg [structure]: Matlab's tructure containing the steady-state values of
+%                       aggregate variables
+% OUTPUTS
+function [out_ss, sizes] = check_steady_state_input(M_, options_, ss)
+   %% Checks
+   % Retrieve variables information from M_
+   state_symbs = M_.heterogeneity(1).endo_names(M_.heterogeneity(1).state_var);
+   pol_symbs= M_.heterogeneity(1).endo_names;
+   inn_symbs = M_.heterogeneity(1).exo_names;
+   agg_symbs = M_.endo_names(1:M_.orig_endo_nbr);
+   % Initialize output variables
+   sizes = struct;
+   out_ss = struct;
+   if ~isstruct(ss)
+      error('Misspecified steady-state input `ss`: the steady-state input `ss` is not a structure.')
+   end
+   %% Shocks
+   % Check that the field `ss.shocks` exists
+   check_isfield('shocks', ss, 'ss.shocks', ' Models without idiosyncratic shocks but with ex-post heterogeneity over individual state variables cannot be solved yet.')
+   shocks = ss.shocks;
+   % Check that `ss.shocks` is a structure
+   check_isstruct(shocks, 'ss.shocks');
+   % Check that the field `ss.shocks.grids` exists
+   check_isfield('grids', ss.shocks, 'ss.shocks.grids');
+   % Check that `ss.shocks.grids` is a structure
+   check_isstruct(shocks.grids, 'ss.shocks.grids');
+   shock_symbs = fieldnames(shocks.grids);
+   % Check the types of `ss.shocks.grids` values
+   check_fun(shocks.grids, 'ss.shocks.grids', shock_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(shock_symbs)
+      s = shock_symbs{i};
+      if isrow(ss.shocks.grids.(s))
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s)';
+      else
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s);
+      end
+   end
+   % Check shock discretization
+   flag_ar1 = isfield(shocks, 'Pi');
+   flag_gh = isfield(shocks, 'w');
+   if ~flag_ar1 && ~flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both non-specified, while exactly one of them should be set.');
+   end
+   % Check that either individual AR(1) processes or i.i.d gaussion innovation
+   % processes are discretized
+   if flag_ar1 && flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both specified, while only one of them should be.');
+   end
+   % Case of discretized AR(1) exogenous processes
+   if flag_ar1
+      % Check that the grids specification for AR(1) shocks and the
+      % var(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, pol_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).endo_names');
+      % Check that shocks.Pi is a structure
+      check_isstruct(shocks.Pi, 'ss.shocks.Pi');
+      % Check that the grids specification for individual AR(1) shocks, the
+      % information in the M_ structure and the Markov transition matrices
+      % specification are compatible and display a warning if redundancies are
+      % detected
+      check_missingredundant(shocks.Pi, 'ss.shocks.Pi', shock_symbs, options_.hank.nowarningredundant);
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the type of shocks.Pi field values
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) ismatrix(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real matrices');
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      % Check that the matrices provided in shocks.Pi are Markov matrices:
+      %     - all elements are non-negative
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(x >= 0, 'all'), 'contain negative elements')
+      %     - the sum of each row equals 1.
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(abs(sum(x,2)-1) < options_.hank.tol_markov), 'have row sums different from 1.')
+      % Remove AR(1) shock processes from state and policy variables
+      state_symbs = setdiff(state_symbs, shock_symbs);
+      pol_symbs = setdiff(pol_symbs, shock_symbs);
+      % Copy relevant shock-related elements in out_ss 
+      out_ss.shocks.Pi = ss.shocks.Pi;
+   end
+   % Case of discretized i.i.d gaussian innovations
+   if flag_gh
+      % Check that shocks.w is a structure
+      check_isstruct(shocks.w, 'shocks.w');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the varexo(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, inn_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).exo_names');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the Gauss-Hermite weights specification are compatible and display a
+      % warning if redundancies are detected
+      check_missingredundant(shocks.w, 'ss.shocks.w', shock_symbs, options_.hank.nowarningredundant);
+      % Check the type of `ss.shocks.w` elements
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) isvector(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real vectors');
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids 
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) numel(x), 'have incompatible sizes with those of `ss.shocks.grids`', grid_nb);
+      % Check that the arrays provided in shocks.w have the properties of
+      % Gauss-Hermite weights:
+      %     - all elements are non-negative
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(x >= 0.), 'contain negative elements');
+      %     - the sum of Gauss-Hermite weights is 1
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(abs(sum(x)-1) < options_.hank.tol_markov), 'have sums different from 1.');
+      % Make `ss.shocks.w` values column vectors in the output steady-state
+      % structure
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         if isrow(ss.shocks.w.(s))
+            out_ss.shocks.w.(s) = ss.shocks.w.(s)';
+         else
+            out_ss.shocks.w.(s) = ss.shocks.w.(s);
+         end
+      end
+   end
+   %% Policy functions
+   % Check that the `ss.pol` field exists
+   check_isfield('pol', ss, 'ss.pol');
+   pol = ss.pol;
+   % Check that `ss.pol` is a structure
+   check_isstruct(ss.pol, 'ss.pol');
+   % Check that the `ss.pol`.grids field exists
+   check_isfield('grids', ss.pol, 'ss.pol.grids');
+   % Check that `ss.pol.grids` is a structure
+   check_isstruct(ss.pol.grids, 'ss.pol.grids');
+   % Check that the state grids specification and the var(heterogeneity=)
+   % statement are compatible
+   check_missingredundant(pol.grids, 'ss.pol.grids', state_symbs, options_.hank.nowarningredundant);
+   % Check that `ss.pol.grids` values are dense real vectors
+   check_fun(pol.grids, 'ss.pol.grids', state_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Store:
+   %    - the number of states
+   sizes.n_a = numel(state_symbs);
+   %    - the size of the tensor product of the states grids 
+   grid_nb = cellfun(@(s) numel(pol.grids.(s)), state_symbs);
+   sizes.pol.N_a = prod(grid_nb);
+   %    - the size of each state grid
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      sizes.pol.states.(s) = numel(pol.grids.(s));
+   end
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      if isrow(ss.pol.grids.(s))
+         out_ss.pol.grids.(s) = ss.pol.grids.(s)';
+      else
+         out_ss.pol.grids.(s) = ss.pol.grids.(s);
+      end
+   end
+   % Check that the field `ss.pol.values` exists 
+   check_isfield('values', pol, 'ss.pol.values');
+   % Check that `ss.pol.values` is a struct
+   check_isstruct(pol.values, 'ss.pol.values');
+   % Check the missing and redundant variables in `ss.pol.values`
+   check_missingredundant(pol.values, 'ss.pol.values', pol_symbs, options_.hank.nowarningredundant);
+   % Check that `ss.pol.values` values are dense real matrices
+   check_fun(pol.values, 'ss.pol.values', pol_symbs, @(x) isnumeric(x) && isreal(x) && ismatrix(x) && ~issparse(x), 'are not dense real matrices'); 
+   % Check the internal size compatibility of `ss.pol.values`
+   sizes.n_pol = numel(pol_symbs);
+   check_fun(pol.values, 'ss.pol.values', pol_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the sizes of `ss.shocks.grids` elements', sizes.N_e);
+   check_fun(pol.values, 'ss.pol.values', pol_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the sizes of `ss.pol.grids` elements', sizes.pol.N_a);
+   % Copy `ss.pol.values` in `out_ss`
+   out_ss.pol.values = ss.pol.values;
+   % Check the permutation of state variables for policy functions
+   [out_ss.pol.states] = check_permutation(pol, 'states', 'ss.pol', state_symbs);
+   % Check the permutation of shock variables for policy functions
+   [out_ss.pol.shocks] = check_permutation(pol, 'shocks', 'ss.pol', shock_symbs);
+   %% Distribution
+   % Check that the field `ss.d` exists 
+   check_isfield('d', ss, 'ss.d');
+   d = ss.d;
+   % Check that the field `ss.d.hist` exists 
+   check_isfield('hist', ss.d, 'ss.d.hist');
+   % Check the type of `ss.d.hist`
+   if ~(ismatrix(d.hist) && isnumeric(d.hist) && isreal(d.hist) && ~issparse(d.hist))
+      error('Misspecified steady-state input `ss`: `ss.d.hist` is not a dense real matrix.');
+   end
+   % Check the consistency of `ss.d.grids` 
+   if ~isfield(d, 'grids')
+      if ~options_.hank.nowarningdgrids
+         warning(['In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states ', cell2mat(state_symbs), '. The policy grids in `ss.pol.grids` shall be used.']);
+      end
+      % Copy the relevant sizes from the pol field
+      sizes.d = sizes.pol;
+      out_ss.d.grids = out_ss.pol.grids;
+   else
+      % Check that `ss.d.grids` is a struct
+      check_isstruct(d.grids, 'ss.d.grids');
+      % Check redundant variables in `ss.d.grids`
+      d_grids_symbs = fieldnames(d.grids);
+      if isempty(d_grids_symbs)
+         if ~options_.hank.nowarningredundant
+            warning(['In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states ', cell2mat(state_symbs), '. The policy grids in `ss.pol.grids` shall be used.']);
+         end
+         % Copy the relevant sizes from the pol field
+         sizes.d = sizes.pol;
+         out_ss.d.grids = out_ss.pol.grids;
+      else
+         d_grids_symbs_in_states = ismember(d_grids_symbs, state_symbs);
+         if ~all(d_grids_symbs_in_states)
+            if ~options_.hank.nowarningredundant
+               warning(['In the steady-state input `ss.d.states`, the following specification for the states grids in the distribution structure are not useful: ', cell2mat(d_grids_symbs(~d_grids_symbs_in_states))]);
+            end
+            d_grids_symbs = d_grids_symbs(d_grids_symbs_in_states);
+         end
+         % Check the types of `ss.d.grids` elements
+         check_fun(pol.grids, 'ss.d.grids', d_grids_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+         % Store the size of the states grids for the distribution
+         for i=1:numel(d_grids_symbs)
+            s = d_grids_symbs{i};
+            sizes.d.states.(s) = numel(d.grids.(s));
+            out_ss.d.grids.(s) = d.grids.(s);
+         end
+         states_out_of_d = setdiff(state_symbs, d_grids_symbs);
+         if ~isempty(states_out_of_d)
+            if ~options_.hank.nowarningdgrids
+               warning(['hank.bbeg.het_stoch_simul: in the steady-state structure, no distribution-specific grid set for states ', cell2mat(states_out_of_d), '. The policy grids specified in pol.grids shall be used.']);
+            end
+            % Store the sizes of the unspecified states grids from the pol field
+            for i=1:numel(states_out_of_d)
+               s = states_out_of_d{i};
+               sizes.d.states.(s) = sizes.pol.states.(s);
+               out_ss.d.grids.(s) = pol.grids.(s);
+            end
+         end
+      end
+   end
+   % Check the internal size compatibility of the distribution histogram
+   if size(d.hist,1) ~= sizes.N_e
+      error('Misspecified steady-state input `ss`: the number of rows of the histogram matrix `ss.d.hist` is not consistent with the sizes of shocks grids `ss.shocks.grids`');
+   end
+   sizes.d.N_a = prod(structfun(@(x) x, sizes.d.states));
+   if size(d.hist,2) ~= sizes.d.N_a
+      error('Misspecified steady-state input `ss`: the number of columns of the histogram matrix `ss.d.hist` is not consistent with the sizes of states grids `ss.d.grids`/`ss.pol.grids`');
+   end
+   % Copy `ss.d.hist` in `out_ss` 
+   out_ss.d.hist = ss.d.hist;
+   % Check the permutation of state variables in the distribution
+   [out_ss.d.states] = check_permutation(d, 'states', 'ss.d', state_symbs);
+   % Check the permutation of shock variables in the distribution
+   [out_ss.d.shocks] = check_permutation(d, 'shocks', 'ss.d', shock_symbs);
+   %% Aggregate variables
+   % Check that the field `ss.agg` exists
+   check_isfield('agg', ss, 'ss.agg');
+   agg = ss.agg;
+   % Check that `ss.agg` is a structure
+   check_isstruct(agg, 'ss.agg');
+   % Check that the aggregate variables specification and the var statement are
+   % compatible
+   check_missingredundant(agg, 'ss.agg', agg_symbs, options_.hank.nowarningredundant);
+   % Store the number of aggregate variables
+   sizes.agg = numel(fieldnames(agg));
+   % Check the types of `ss.agg` values
+   check_fun(agg, 'ss.agg', agg_symbs, @(x) isreal(x) && isscalar(x), 'are not real scalars');
+   % Copy `ss.agg` into `out_ss`
+   out_ss.agg = agg;
+end
diff --git a/matlab/+hank/private/check_consistency.m b/matlab/+hank/private/check_consistency.m
new file mode 100644
index 0000000000000000000000000000000000000000..4ba1d8e788f3f6ee656a7c142df3f0981747af83
--- /dev/null
+++ b/matlab/+hank/private/check_consistency.m
@@ -0,0 +1,6 @@
+function [] = check_consistency(f1, f2, f1_name, f2_name)
+   f1_in_f2 = ismember(f1,f2);
+   if ~all(f1_in_f2)
+      error('Misspecified steady-state input `ss`: the following variables are mentioned in `%s`, but are missing in `%s`: %s', f1_name, f2_name, cell2mat(f1(~f1_in_f2)));
+   end
+end
\ No newline at end of file
diff --git a/matlab/+hank/private/check_fun.m b/matlab/+hank/private/check_fun.m
new file mode 100644
index 0000000000000000000000000000000000000000..1174c6b71b496e9706892c7d4c6874f32b682d40
--- /dev/null
+++ b/matlab/+hank/private/check_fun.m
@@ -0,0 +1,9 @@
+function [] = check_fun(f, f_name, symbs, fun, text, eq)
+   elt_check = cellfun(@(s) fun(f.(s)), symbs);
+   if nargin == 6
+      elt_check = elt_check == eq;
+   end
+   if ~all(elt_check)
+      error('Misspecified steady-state input `ss`: the following fields in `%s` %s: %s', f_name, text, cell2mat(symbs(~elt_check)));
+   end
+end
\ No newline at end of file
diff --git a/matlab/+hank/private/check_isfield.m b/matlab/+hank/private/check_isfield.m
new file mode 100644
index 0000000000000000000000000000000000000000..4f4f537d3a933857e92de5f099c3fc2e8df723ef
--- /dev/null
+++ b/matlab/+hank/private/check_isfield.m
@@ -0,0 +1,11 @@
+function [] = check_isfield(f, s, f_name, details)
+   if nargin < 4
+      details = '';
+   end
+   if nargin < 3
+      f_name = f;
+   end
+   if ~isfield(s, f)
+      error('Misspecified steady-state input `ss`: the `%s` field is missing.%s', f_name, details);
+   end
+end
\ No newline at end of file
diff --git a/matlab/+hank/private/check_isstruct.m b/matlab/+hank/private/check_isstruct.m
new file mode 100644
index 0000000000000000000000000000000000000000..292d45fd41720fc082a235eb00e8265f2ce3ab2a
--- /dev/null
+++ b/matlab/+hank/private/check_isstruct.m
@@ -0,0 +1,5 @@
+function [] = check_isstruct(f, f_name)
+   if ~isstruct(f)
+      error('Misspecified steady-state input `ss`: the `%s` field should be a structure.', f_name);
+   end
+end
\ No newline at end of file
diff --git a/matlab/+hank/private/check_missingredundant.m b/matlab/+hank/private/check_missingredundant.m
new file mode 100644
index 0000000000000000000000000000000000000000..c076b7a90dc0ba7e53bbb879cfe578bea68ba205
--- /dev/null
+++ b/matlab/+hank/private/check_missingredundant.m
@@ -0,0 +1,13 @@
+function [] = check_missingredundant(f, f_name, symbs, nowarningredundant)
+   fields = fieldnames(f);
+   symbs_in_fields = ismember(symbs, fields);
+   if ~all(symbs_in_fields)
+      error('Misspecified steady-state input `ss`. The following fields are missing in `%s`: %s', f_name, cell2mat(symbs(~symbs_in_fields)));
+   end
+   if ~nowarningredundant
+      fields_in_symbs = ismember(fields, symbs);
+      if ~all(fields_in_symbs)
+         warning('Steady-state input `ss`. The following fields are redundant in `%s`: %s', f_name, cell2mat(fields(~fields_in_symbs)));
+      end
+   end
+end
\ No newline at end of file
diff --git a/matlab/+hank/private/check_permutation.m b/matlab/+hank/private/check_permutation.m
new file mode 100644
index 0000000000000000000000000000000000000000..2a7931d46aff01e850256c84ba15c729719e1479
--- /dev/null
+++ b/matlab/+hank/private/check_permutation.m
@@ -0,0 +1,15 @@
+function [out_order] = check_permutation(s, f_name, s_name, symbs)
+   if ~isfield(s, f_name)
+      out_order = symbs;
+   else
+      f = s.(f_name);
+      if ~isstring(f) && ~iscellstr(f)
+         error('Misspecified steady-state input `ss`: the `%s.%s` field should be a cell array of character vectors or a string array.', s_name, f_name);
+      end
+      err_var = setdiff(f, symbs);
+      if ~isempty(err_var)
+         error('Misspecified steady-state input `ss`: the set of variables of the `%s.%s` field is not consistent with the information in M_ and the other fields in the steady-state structure `ss`. Problematic variables: %s', s_name, f_name, cell2mat(err_var));
+      end
+      out_order = f;
+   end
+end
\ No newline at end of file
diff --git a/meson.build b/meson.build
index 351922ff619aab9d42e549e2521eb12e2df6ed1d..b838bc90f6f3da8b2a3304a3c90eb225d2f6f0c9 100644
--- a/meson.build
+++ b/meson.build
@@ -1896,6 +1896,8 @@ mod_and_m_tests = [
                 'solver-test-functions/variablydimensioned.m',
                 'solver-test-functions/watson.m',
                 'solver-test-functions/wood.m',] },
+  { 'test' : [ 'hank/ks_ar1_ss.m',
+               'hank/ks_iid_ss.m' ] },
   { 'test' : [ 'cyclereduction.m' ] },
   { 'test' : [ 'logarithmicreduction.m' ] },
   { 'test' : [ 'riccatiupdate.m' ] },
diff --git a/tests/hank/expect_error.m b/tests/hank/expect_error.m
new file mode 100644
index 0000000000000000000000000000000000000000..6c065eba2c639cb95442344d98006e06540d9181
--- /dev/null
+++ b/tests/hank/expect_error.m
@@ -0,0 +1,16 @@
+function [incTestFailed] = expect_error(fn, description, testFailed, show_message)
+    if nargin < 3
+        show_message = false;
+    end
+    try
+        fn();
+        incTestFailed = testFailed+1;
+        fprintf('❌ Expected error for: %s\n', description);
+    catch ME
+        fprintf('✔ Correctly threw error for: %s\n', description);
+        if show_message
+            fprintf('   ↪ Error message: %s\n', ME.message);
+        end
+        incTestFailed = testFailed;
+    end
+end
\ No newline at end of file
diff --git a/tests/hank/ks.mod b/tests/hank/ks.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0164957d3abe683752a28c2b0d024c89b751788b
--- /dev/null
+++ b/tests/hank/ks.mod
@@ -0,0 +1,66 @@
+heterogeneity_dimension households;
+
+var(heterogeneity=households) 
+    c  // Consumption
+    a  // Assets
+    e  // Idiosyncratic (log-)efficiency
+    Va // Derivative of the value function w.r.t assets
+;
+
+varexo(heterogeneity=households) eps_e; // Idiosyncratic efficiency shock
+
+var
+    r // Rate of return on capital net of depreciation
+    w // Wage rate
+    Y // Aggregate output
+    K // Aggregate capital
+    Z // Aggregate productivity
+;
+varexo eps_Z; // Aggregate productivity shock
+
+parameters
+    L     // Labor
+    alpha // Share of capital in production fuction
+    beta  // Subjective discount rate of houselholds
+    delta // Capital depreciation rate
+    eis   // Elasticity of intertemporal substitution
+    rho_e // Earning shock persistence
+    sig_e // Earning shock innovation std err
+    rho_Z // Aggregate TFP shock persistence
+    sig_Z // Aggregate TFP shock innovation std err
+    Z_ss  // Aggregate TFP shock average value
+;
+
+model(heterogeneity=households);
+    c^(-1/eis)-beta*Va(+1)=0 ⟂ a>=0;
+    (1+r)*a(-1)+w*e-c-a; 
+    Va = (1+r)*c^(-1/eis);
+    log(e) - rho_e*log(e(-1)) - eps_e;
+end;
+
+model;
+    Z * K(-1)^alpha * L^(1 - alpha) - Y;
+    alpha * Z * (K(-1) / L)^(alpha - 1) - delta - r;
+    (1 - alpha) * Z * (K(-1) / L)^alpha - w;
+    K - SUM(a);
+    log(Z) - rho_Z*log(Z(-1)) - (1-rho_Z)*log(Z_ss) - eps_Z;
+end;
+
+alpha = 0.11;
+beta = 0.9819527880123727;
+eis = 1;
+delta = 0.025;
+L = 1;
+rho_e = 0.966;
+sig_e = 0.5*sqrt(1-rho_e^2);
+rho_Z = 0.8;
+sig_Z = 0.014;
+Z_ss = 0.8816;
+
+shocks(heterogeneity=households);
+    var eps_e; stderr sig_e;
+end;
+
+shocks;
+    var eps_Z; stderr sig_Z;
+end;
\ No newline at end of file
diff --git a/tests/hank/ks_ar1_ss.m b/tests/hank/ks_ar1_ss.m
new file mode 100644
index 0000000000000000000000000000000000000000..7a753cfd270df9d5d6e04151f92aae9985c0b33c
--- /dev/null
+++ b/tests/hank/ks_ar1_ss.m
@@ -0,0 +1,402 @@
+% source_dir = getenv('source_root');
+% addpath([source_dir filesep 'matlab']);
+
+dynare_config;
+
+testFailed = 0;
+
+skipline()
+disp('*** TESTING: hank.check_steady_state_input.m - Discretized AR(1) case ***');
+
+dynare ks.mod;
+
+load 'ks_ar1_ss.mat';
+
+verbose = true;
+
+%% === OUTPUTS ===
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, base_struct);
+    if ~(sizes.n_e == numel(fieldnames(base_struct.shocks.grids)) && ...
+         sizes.n_a == numel(fieldnames(base_struct.pol.grids)) && ...
+         sizes.N_e == numel(base_struct.shocks.grids.e) && ...
+         sizes.n_pol == numel(fieldnames(base_struct.pol.values)) && ...
+         sizes.agg == numel(fieldnames(base_struct.agg)) && ...
+         sizes.shocks.e == numel(base_struct.shocks.grids.e) && ...
+         sizes.pol.N_a == numel(base_struct.pol.grids.a) && ...
+         sizes.pol.states.a == sizes.pol.N_a && ...
+         sizes.d.N_a == sizes.pol.N_a && sizes.d.states.a == sizes.pol.states.a)
+      testFailed = testFailed+1;
+      dprintf('The `sizes` output is not correct!');
+    end
+    if ~(iscolumn(out_ss.shocks.grids.e) && iscolumn(out_ss.pol.grids.a) && ...
+         iscolumn(out_ss.d.grids.a) && ...
+         size(out_ss.shocks.grids.e,1) == sizes.shocks.e && ...
+         size(out_ss.pol.grids.a,1) == sizes.pol.states.a && ...
+         size(out_ss.d.grids.a,1) == sizes.d.states.a && ...
+         numel(fieldnames(out_ss.agg)) == sizes.agg)
+      testFailed = testFailed+1;
+      dprintf('The `out_ss` output is not correct!');
+    end
+catch ME
+    testFailed = testFailed+2;
+    dprintf('Outputs `sizes` or `out_ss` are not correct!')
+end
+
+%% === NON-STRUCT FIELDS ===
+
+% ss is not a struct
+ss = 123;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss is not a struct', testFailed, verbose);
+
+% ss.shocks is not a struct
+ss = base_struct;
+ss.shocks = 1;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks is not a struct', testFailed, verbose);
+
+% ss.shocks.grids is not a struct
+ss = base_struct;
+ss.shocks.grids = 42;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks.grids is not a struct', testFailed, verbose);
+
+% ss.shocks.Pi is not a struct
+ss = base_struct;
+ss.shocks.Pi = "not_a_struct";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks.Pi is not a struct', testFailed, verbose);
+
+% ss.shocks.w is not a struct
+ss = base_struct;
+ss.shocks = rmfield(ss.shocks, 'Pi');
+ss.shocks.w = "not_a_struct";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks.w is not a struct', testFailed, verbose);
+
+% ss.pol is not a struct
+ss = base_struct;
+ss.pol = pi;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.pol is not a struct', testFailed, verbose);
+
+% ss.pol.grids is not a struct
+ss = base_struct;
+ss.pol.grids = "grid";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.pol.grids is not a struct', testFailed, verbose);
+
+% ss.pol.values is not a struct
+ss = base_struct;
+ss.pol.values = "values";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.pol.values is not a struct', testFailed, verbose);
+
+% ss.d is not a struct
+ss = base_struct;
+ss.d = 1;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.d is not a struct', testFailed, verbose);
+
+% ss.d.grids is not a struct
+ss = base_struct;
+ss.d.grids = "grid";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.d.grids is not a struct', testFailed, verbose);
+
+% ss.agg is not a struct
+ss = base_struct;
+ss.agg = 0;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.agg is not a struct', testFailed, verbose);
+
+%% === SHOCKS TESTS ===
+
+% Missing `shocks` field
+ss = base_struct;
+ss = rmfield(ss, 'shocks');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `shocks` field', testFailed, verbose);
+
+% Missing `shocks.grids`
+ss = base_struct;
+ss.shocks = rmfield(ss.shocks, 'grids');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `shocks.grids`', testFailed, verbose);
+
+% Wrong type for `shocks.grids.e`
+ss = base_struct;
+ss.shocks.grids.e = ones(3,3);  % invalid: cell array
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'shocks.grids.e are not dense real vectors', testFailed, verbose);
+
+% Missing both `shocks.Pi` and `shocks.w`
+ss = base_struct;
+ss.shocks = rmfield(ss.shocks, 'Pi');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'No discretization method (Pi or w) provided', testFailed, verbose);
+
+% Both `shocks.Pi` and `shocks.w` provided (mutual exclusivity)
+ss = base_struct;
+ss.shocks.w = struct('e', ones(size(ss.shocks.Pi.e,1),1));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Both Pi and w should not coexist', testFailed, verbose);
+
+% AR(1) Tests
+% Shock grid not in var(heterogeneity=) list
+ss = base_struct;
+ss.shocks.grids.fake_shock = ss.shocks.grids.e;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Shock grid symbol not in heterogeneity declaration', testFailed, verbose);
+
+% Missing Markov matrix for a shock
+ss = base_struct;
+ss.shocks.Pi = rmfield(ss.shocks.Pi, 'e');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing Pi entry for a declared shock', testFailed, verbose);
+
+% Wrong type for `shocks.grids.Pi.e`
+ss = base_struct;
+ss.shocks.Pi.e = speye(numel(ss.shocks.grids.e));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'shocks.Pi.e are not dense real matrices', testFailed, verbose);
+
+% Markov matrix wrong row number
+ss = base_struct;
+ss.shocks.Pi.e = rand(3, 7);  % Should be square matching grid
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Pi row count mismatch with grid', testFailed, verbose);
+
+% Markov matrix wrong column number
+ss = base_struct;
+ss.shocks.Pi.e = rand(7, 3);  % Should be square matching grid
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Pi column count mismatch with grid', testFailed, verbose);
+
+% Markov matrix has negative elements
+ss = base_struct;
+ss.shocks.Pi.e(1,1) = -0.1;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Negative elements in Pi', testFailed, verbose);
+
+% Markov matrix does not sum to 1
+ss = base_struct;
+ss.shocks.Pi.e = ss.shocks.Pi.e * 2;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Rows of Pi not summing to 1', testFailed, verbose);
+
+
+%% === POLICY TESTS ===
+
+% Missing `pol` field
+ss = base_struct;
+ss = rmfield(ss, 'pol');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `pol` field', testFailed, verbose);
+
+% Missing `pol.grids`
+ss = base_struct;
+ss.pol = rmfield(ss.pol, 'grids');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `pol.grids`', testFailed, verbose);
+
+% Missing one state grid
+ss = base_struct;
+ss.pol.grids = rmfield(ss.pol.grids, 'a');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing state grid for a', testFailed, verbose);
+
+% Wrong type for `pol.grids.a`
+ss = base_struct;
+ss.pol.grids.a = ones(3, 3);  % should be a vector
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.grids.a should be a vector', testFailed, verbose);
+
+% Missing `pol.values`
+ss = base_struct;
+ss.pol = rmfield(ss.pol, 'values');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `pol.values`', testFailed, verbose);
+
+% Missing policy value for declared variable
+ss = base_struct;
+ss.pol.values = rmfield(ss.pol.values, 'a');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing policy value for a', testFailed, verbose);
+
+% Wrong type for `pol.values.a`
+ss = base_struct;
+ss.pol.values.a = rand(2, 2, 2); 
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.values.a should be a matrix', testFailed, verbose);
+
+% Incompatible policy matrix row size (shocks × states)
+ss = base_struct;
+ss.pol.values.a = rand(99, size(ss.pol.values.a, 2));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of rows in policy matrix', testFailed, verbose);
+
+% Incompatible policy matrix column size
+ss = base_struct;
+ss.pol.values.a = rand(size(ss.pol.values.a, 1), 99);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of columns in policy matrix', testFailed, verbose);
+
+%% === DISTRIBUTION TESTS ===
+
+% Missing `d` field
+ss = base_struct;
+ss = rmfield(ss, 'd');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing distribution structure `d`', testFailed, verbose);
+
+% Wrong type for `ss.d.grids.a`
+ss = base_struct;
+ss.d.grids.a = rand(3, 3);  % not a vector
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.grids.a should be a vector', testFailed, verbose);
+
+% Missing `d.hist`
+ss = base_struct;
+ss.d = rmfield(ss.d, 'hist');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing d.hist', testFailed, verbose);
+
+% Wrong type for `ss.d.hist`
+ss = base_struct;
+ss.d.hist = rand(2,2,2);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.hist should be a matrix', testFailed, verbose);
+
+% Wrong row size in `d.hist`
+ss = base_struct;
+ss.d.hist = rand(99, size(ss.d.hist,2));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of rows in d.hist', testFailed, verbose);
+
+% Wrong column size in `d.hist`
+ss = base_struct;
+ss.d.hist = rand(size(ss.d.hist,1), 99);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of columns in d.hist', testFailed, verbose);
+
+%% === AGGREGATES TESTS ===
+
+% Missing `agg` field
+ss = base_struct;
+ss = rmfield(ss, 'agg');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `agg` field', testFailed, verbose);
+
+% Remove one required aggregate variable
+ss = base_struct;
+ss.agg = rmfield(ss.agg, 'r');  % remove interest rate
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing aggregate variable r', testFailed, verbose);
+
+% Wrong type for `ss.agg.r`
+ss = base_struct;
+ss.agg.r = 'non-numeric';
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'agg.r should be numeric', testFailed, verbose);
+
+%% === PERMUTATION TESTS ===
+
+% pol.shocks is not a string or cellstr
+ss = base_struct;
+ss.pol.shocks = 123;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.shocks is not a string array or cellstr', testFailed, verbose);
+
+% pol.shocks has wrong names
+ss = base_struct;
+ss.pol.shocks = {'wrong_name'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.shocks not matching exo_names', testFailed, verbose);
+
+% pol.states is not a string or cellstr
+ss = base_struct;
+ss.pol.states = 42;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.states is not a string array or cellstr', testFailed, verbose);
+
+% pol.states has wrong names
+ss = base_struct;
+ss.pol.states = {'bad_state'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.states not matching state_var', testFailed, verbose);
+
+% d.shocks is not a string or cellstr
+ss = base_struct;
+ss.d.shocks = 3.14;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.shocks is not a string array or cellstr', testFailed, verbose);
+
+% d.shocks inconsistent with pol.shocks
+ss = base_struct;
+ss.d.shocks = {'ghost_shock'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.shocks not matching pol.shocks', testFailed, verbose);
+
+% d.states is not a string or cellstr
+ss = base_struct;
+ss.d.states = true;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.states is not a string array or cellstr', testFailed, verbose);
+
+% d.states inconsistent with pol.states
+ss = base_struct;
+ss.d.states = {'phantom_state'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.states not matching pol.states', testFailed, verbose);
+
+%% === REDUNDANCY TESTS ===
+
+options_.hank.nowarningredundant = false;
+
+% shocks.Pi contains redundant entry (not in shocks.grids)
+ss = base_struct;
+ss.shocks.Pi.extra = eye(length(ss.shocks.grids.e));  % not listed in shocks.grids
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant shocks.Pi entry accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant shocks.Pi entry: %s\n', ME.message);
+end
+
+% pol.values contains redundant field (not in model's pol_symbs)
+ss = base_struct;
+ss.pol.values.redundant_var = ones(size(ss.pol.values.a));
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant pol.values field accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant pol.values field: %s\n', ME.message);
+end
+
+% d.grids contains redundant field (not a declared state)
+ss = base_struct;
+ss.d.grids.not_a_state = [0;1;2];  % not in M_.heterogeneity.state_var
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant d.grids field accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant d.grids field: %s\n', ME.message);
+end
+
+% agg contains extra field (not in M_.endo_names)
+ss = base_struct;
+ss.agg.extra_agg = 999;
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant agg field accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant agg field: %s\n', ME.message);
+end
+
+quit(testFailed > 0)
diff --git a/tests/hank/ks_ar1_ss.mat b/tests/hank/ks_ar1_ss.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9062b46082afdf5415596d6c11a4ff2657329601
--- /dev/null
+++ b/tests/hank/ks_ar1_ss.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d1cd897108ad2a9db76fca6ec42889ad8755094d3c7b714e62c3083a973601c7
+size 666774
diff --git a/tests/hank/ks_iid_ss.m b/tests/hank/ks_iid_ss.m
new file mode 100644
index 0000000000000000000000000000000000000000..732630163cdfc4e302944d19efd3a0d6705d9b02
--- /dev/null
+++ b/tests/hank/ks_iid_ss.m
@@ -0,0 +1,392 @@
+source_dir = getenv('source_root');
+addpath([source_dir filesep 'matlab']);
+
+dynare_config;
+
+testFailed = 0;
+
+skipline()
+disp('*** TESTING: hank.check_steady_state_input.m - Discretized i.i.d case ***');
+
+dynare ks.mod;
+
+load 'ks_iid_ss.mat';
+
+verbose = true;
+
+%% === OUTPUTS ===
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, base_struct);
+    if ~(sizes.n_e == numel(fieldnames(base_struct.shocks.grids)) && ...
+         sizes.n_a == numel(fieldnames(base_struct.pol.grids)) && ...
+         sizes.N_e == numel(base_struct.shocks.grids.eps_e) && ...
+         sizes.n_pol == numel(fieldnames(base_struct.pol.values)) && ...
+         sizes.agg == numel(fieldnames(base_struct.agg)) && ...
+         sizes.shocks.eps_e == numel(base_struct.shocks.grids.eps_e) && ...
+         sizes.pol.N_a == N_a && sizes.pol.states.a == numel(base_struct.pol.grids.a) && ...
+         sizes.pol.states.e == numel(base_struct.pol.grids.e) && ...
+         sizes.d.N_a == sizes.pol.N_a && sizes.d.states.a == sizes.pol.states.a)
+      testFailed = testFailed+1;
+      dprintf('The `sizes` output is not correct!');
+    end
+    if ~(iscolumn(out_ss.shocks.grids.eps_e) && iscolumn(out_ss.pol.grids.a) && ...
+         iscolumn(out_ss.pol.grids.e) && iscolumn(out_ss.d.grids.a) && ...
+         iscolumn(out_ss.d.grids.e) && ...
+         size(out_ss.shocks.grids.eps_e,1) == sizes.shocks.eps_e && ...
+         size(out_ss.pol.grids.a,1) == sizes.pol.states.a && ...
+         size(out_ss.pol.grids.e,1) == sizes.pol.states.e && ...
+         size(out_ss.d.grids.a,1) == sizes.d.states.a && ...
+         size(out_ss.d.grids.e,1) == sizes.d.states.e && ...
+         numel(fieldnames(out_ss.agg)) == sizes.agg)
+      testFailed = testFailed+1;
+      dprintf('The `out_ss` output is not correct!');
+    end
+catch ME
+    testFailed = testFailed+2;
+    dprintf('Outputs `sizes` or `out_ss` are not correct!')
+end
+
+%% === NON-STRUCT FIELDS ===
+
+% ss is not a struct
+ss = 123;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss is not a struct', testFailed, verbose);
+
+% ss.shocks is not a struct
+ss = base_struct;
+ss.shocks = 1;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks is not a struct', testFailed, verbose);
+
+% ss.shocks.grids is not a struct
+ss = base_struct;
+ss.shocks.grids = 42;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks.grids is not a struct', testFailed, verbose);
+
+% ss.shocks.w is not a struct
+ss = base_struct;
+ss.shocks.w = "not_a_struct";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.shocks.w is not a struct', testFailed, verbose);
+
+% ss.pol is not a struct
+ss = base_struct;
+ss.pol = pi;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.pol is not a struct', testFailed, verbose);
+
+% ss.pol.grids is not a struct
+ss = base_struct;
+ss.pol.grids = "grid";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.pol.grids is not a struct', testFailed, verbose);
+
+% ss.pol.values is not a struct
+ss = base_struct;
+ss.pol.values = "values";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.pol.values is not a struct', testFailed, verbose);
+
+% ss.d is not a struct
+ss = base_struct;
+ss.d = 1;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.d is not a struct', testFailed, verbose);
+
+% ss.d.grids is not a struct
+ss = base_struct;
+ss.d.grids = "grid";
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.d.grids is not a struct', testFailed, verbose);
+
+% ss.agg is not a struct
+ss = base_struct;
+ss.agg = 0;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'ss.agg is not a struct', testFailed, verbose);
+
+%% === SHOCKS TESTS ===
+
+% Missing `shocks` field
+ss = base_struct;
+ss = rmfield(ss, 'shocks');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `shocks` field', testFailed, verbose);
+
+% Missing `shocks.grids`
+ss = base_struct;
+ss.shocks = rmfield(ss.shocks, 'grids');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `shocks.grids`', testFailed, verbose);
+
+% Wrong type for `shocks.grids.eps_e`
+ss = base_struct;
+ss.shocks.grids.eps_e = ones(3,3);  % invalid: cell array
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'shocks.grids.eps_e are not dense real vectors', testFailed, verbose);
+
+% Missing both `shocks.Pi` and `shocks.w`
+ss = base_struct;
+ss.shocks = rmfield(ss.shocks, 'w');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'No discretization method (Pi or w) provided', testFailed, verbose);
+
+% Both `shocks.Pi` and `shocks.w` provided (mutual exclusivity)
+ss = base_struct;
+ss.shocks.Pi = struct('e', ones(size(ss.pol.grids.e,1)));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Both Pi and w should not coexist', testFailed, verbose);
+
+% AR(1) Tests
+% Shock grid not in var(heterogeneity=) list
+ss = base_struct;
+ss.shocks.grids.fake_shock = ss.shocks.grids.eps_e;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Shock grid symbol not in heterogeneity declaration', testFailed, verbose);
+
+% Missing Markov matrix for a shock
+ss = base_struct;
+ss.shocks.w = rmfield(ss.shocks.w, 'eps_e');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing w entry for a declared shock', testFailed, verbose);
+
+% Wrong type for `shocks.grids.w.e`
+ss = base_struct;
+ss.shocks.w.eps_e = speye(numel(ss.shocks.grids.eps_e));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'shocks.w.eps_e are not dense real vectors', testFailed, verbose);
+
+% Gauss-Hermite weights count
+ss = base_struct;
+ss.shocks.w.eps_e = rand(sizes.N_e, 1);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    ' row count mismatch with grid', testFailed, verbose);
+
+% Gauss-Hermite weights have negative elements
+ss = base_struct;
+ss.shocks.w.eps_e(1) = -0.1;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Negative elements in Pi', testFailed, verbose);
+
+% Gauss-Hermite weights do not sum to 1
+ss = base_struct;
+ss.shocks.w.eps_e = ss.shocks.w.eps_e * 2;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Rows of Pi not summing to 1', testFailed, verbose);
+
+
+%% === POLICY TESTS ===
+
+% Missing `pol` field
+ss = base_struct;
+ss = rmfield(ss, 'pol');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `pol` field', testFailed, verbose);
+
+% Missing `pol.grids`
+ss = base_struct;
+ss.pol = rmfield(ss.pol, 'grids');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `pol.grids`', testFailed, verbose);
+
+% Missing one state grid
+ss = base_struct;
+ss.pol.grids = rmfield(ss.pol.grids, 'a');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing state grid for a', testFailed, verbose);
+
+% Wrong type for `pol.grids.a`
+ss = base_struct;
+ss.pol.grids.a = ones(3, 3);  % should be a vector
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.grids.a should be a vector', testFailed, verbose);
+
+% Missing `pol.values`
+ss = base_struct;
+ss.pol = rmfield(ss.pol, 'values');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `pol.values`', testFailed, verbose);
+
+% Missing policy value for declared variable
+ss = base_struct;
+ss.pol.values = rmfield(ss.pol.values, 'a');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing policy value for a', testFailed, verbose);
+
+% Wrong type for `pol.values.a`
+ss = base_struct;
+ss.pol.values.a = rand(2, 2, 2); 
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.values.a should be a matrix', testFailed, verbose);
+
+% Incompatible policy matrix row size (shocks × states)
+ss = base_struct;
+ss.pol.values.a = rand(99, size(ss.pol.values.a, 2));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of rows in policy matrix', testFailed, verbose);
+
+% Incompatible policy matrix column size
+ss = base_struct;
+ss.pol.values.a = rand(size(ss.pol.values.a, 1), 99);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of columns in policy matrix', testFailed, verbose);
+
+%% === DISTRIBUTION TESTS ===
+
+% Missing `d` field
+ss = base_struct;
+ss = rmfield(ss, 'd');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing distribution structure `d`', testFailed, verbose);
+
+% Wrong type for `ss.d.grids.a`
+ss = base_struct;
+ss.d.grids.a = rand(3, 3);  % not a vector
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.grids.a should be a vector', testFailed, verbose);
+
+% Missing `d.hist`
+ss = base_struct;
+ss.d = rmfield(ss.d, 'hist');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing d.hist', testFailed, verbose);
+
+% Wrong type for `ss.d.hist`
+ss = base_struct;
+ss.d.hist = rand(2,2,2);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.hist should be a matrix', testFailed, verbose);
+
+% Wrong row size in `d.hist`
+ss = base_struct;
+ss.d.hist = rand(99, size(ss.d.hist,2));
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of rows in d.hist', testFailed, verbose);
+
+% Wrong column size in `d.hist`
+ss = base_struct;
+ss.d.hist = rand(size(ss.d.hist,1), 99);
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Wrong number of columns in d.hist', testFailed, verbose);
+
+%% === AGGREGATES TESTS ===
+
+% Missing `agg` field
+ss = base_struct;
+ss = rmfield(ss, 'agg');
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing `agg` field', testFailed, verbose);
+
+% Remove one required aggregate variable
+ss = base_struct;
+ss.agg = rmfield(ss.agg, 'r');  % remove interest rate
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'Missing aggregate variable r', testFailed, verbose);
+
+% Wrong type for `ss.agg.r`
+ss = base_struct;
+ss.agg.r = 'non-numeric';
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'agg.r should be numeric', testFailed, verbose);
+
+%% === PERMUTATION TESTS ===
+
+% pol.shocks is not a string or cellstr
+ss = base_struct;
+ss.pol.shocks = 123;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.shocks is not a string array or cellstr', testFailed, verbose);
+
+% pol.shocks has wrong names
+ss = base_struct;
+ss.pol.shocks = {'wrong_name'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.shocks not matching exo_names', testFailed, verbose);
+
+% pol.states is not a string or cellstr
+ss = base_struct;
+ss.pol.states = 42;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.states is not a string array or cellstr', testFailed, verbose);
+
+% pol.states has wrong names
+ss = base_struct;
+ss.pol.states = {'bad_state'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'pol.states not matching state_var', testFailed, verbose);
+
+% d.shocks is not a string or cellstr
+ss = base_struct;
+ss.d.shocks = 3.14;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.shocks is not a string array or cellstr', testFailed, verbose);
+
+% d.shocks inconsistent with pol.shocks
+ss = base_struct;
+ss.d.shocks = {'ghost_shock'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.shocks not matching pol.shocks', testFailed, verbose);
+
+% d.states is not a string or cellstr
+ss = base_struct;
+ss.d.states = true;
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.states is not a string array or cellstr', testFailed, verbose);
+
+% d.states inconsistent with pol.states
+ss = base_struct;
+ss.d.states = {'phantom_state'};
+testFailed = expect_error(@() hank.check_steady_state_input(M_, options_, ss), ...
+    'd.states not matching pol.states', testFailed, verbose);
+
+%% === REDUNDANCY TESTS ===
+
+options_.hank.nowarningredundant = false;
+
+% shocks.w contains redundant entry (not in shocks.grids)
+ss = base_struct;
+ss.shocks.xtra = ones(length(ss.shocks.grids.eps_e));  % not listed in shocks.grids
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant shocks.Pi entry accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant shocks.Pi entry: %s\n', ME.message);
+end
+
+% pol.values contains redundant field (not in model's pol_symbs)
+ss = base_struct;
+ss.pol.values.redundant_var = ones(size(ss.pol.values.a));
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant pol.values field accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant pol.values field: %s\n', ME.message);
+end
+
+% d.grids contains redundant field (not a declared state)
+ss = base_struct;
+ss.d.grids.not_a_state = [0;1;2];  % not in M_.heterogeneity.state_var
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant d.grids field accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant d.grids field: %s\n', ME.message);
+end
+
+% agg contains extra field (not in M_.endo_names)
+ss = base_struct;
+ss.agg.extra_agg = 999;
+try
+    [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+    disp('✔ Redundant agg field accepted (warning expected)');
+catch ME
+    testFailed = testFailed+1;
+    fprintf('❌ Unexpected error from redundant agg field: %s\n', ME.message);
+end
+
+quit(testFailed > 0)
\ No newline at end of file
diff --git a/tests/hank/ks_iid_ss.mat b/tests/hank/ks_iid_ss.mat
new file mode 100644
index 0000000000000000000000000000000000000000..947828e358deb79ccc8167ab594c5d0646c49a3f
--- /dev/null
+++ b/tests/hank/ks_iid_ss.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:10e29a3e6a0e791e77d5c95abac4a5eddc7ae9c58aefd19989196be3d4fb8abe
+size 666774