diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 180d3e67e6beb71b2adfb9fb215b16194e9f1d8a..ad61be54dd5c7a04a4e93be819400e342c72211d 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -114,14 +114,14 @@ class ParsingDriver;
 %token <string> QUOTED_STRING
 %token QZ_CRITERIUM QZ_ZERO_THRESHOLD DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE PIPE_E PIPE_X PIPE_P
 %token RELATIVE_IRF REPLIC SIMUL_REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY TARGETS
-%token SHOCKS SHOCK_DECOMPOSITION SHOCK_GROUPS USE_SHOCK_GROUPS SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED ENDOGENOUS_TERMINAL_PERIOD
+%token SHOCKS HETEROSKEDASTIC_SHOCKS SHOCK_DECOMPOSITION SHOCK_GROUPS USE_SHOCK_GROUPS SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED ENDOGENOUS_TERMINAL_PERIOD
 %token SMOOTHER SMOOTHER2HISTVAL SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO SOLVER_PERIODS ROBUST_LIN_SOLVE
 %token STDERR STEADY STOCH_SIMUL SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME REALTIME_SHOCK_DECOMPOSITION CONDITIONAL UNCONDITIONAL
 %token TEX RAMSEY_MODEL RAMSEY_POLICY RAMSEY_CONSTRAINTS PLANNER_DISCOUNT PLANNER_DISCOUNT_LATEX_NAME
 %token DISCRETIONARY_POLICY DISCRETIONARY_TOL EVALUATE_PLANNER_OBJECTIVE
 %token <string> TEX_NAME TRUE
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED
-%token VALUES VAR VAREXO VAREXO_DET VARIABLE VAROBS VAREXOBS PREDETERMINED_VARIABLES VAR_EXPECTATION VAR_EXPECTATION_MODEL PLOT_SHOCK_DECOMPOSITION MODEL_LOCAL_VARIABLE
+%token VALUES SCALES VAR VAREXO VAREXO_DET VARIABLE VAROBS VAREXOBS PREDETERMINED_VARIABLES VAR_EXPECTATION VAR_EXPECTATION_MODEL PLOT_SHOCK_DECOMPOSITION MODEL_LOCAL_VARIABLE
 %token WRITE_LATEX_DYNAMIC_MODEL WRITE_LATEX_STATIC_MODEL WRITE_LATEX_ORIGINAL_MODEL WRITE_LATEX_STEADY_STATE_MODEL
 %token XLS_SHEET XLS_RANGE LMMCP BANDPASS_FILTER COLORMAP VAR_MODEL PAC_MODEL QOQ YOY AOA PAC_EXPECTATION TREND_COMPONENT_MODEL
 %left EQUAL_EQUAL EXCLAMATION_EQUAL
@@ -166,6 +166,7 @@ class ParsingDriver;
 %token PARAMETER_CONVERGENCE_CRITERION NUMBER_OF_LARGE_PERTURBATIONS NUMBER_OF_SMALL_PERTURBATIONS
 %token NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION MAX_NUMBER_OF_STAGES
 %token RANDOM_FUNCTION_CONVERGENCE_CRITERION RANDOM_PARAMETER_CONVERGENCE_CRITERION NO_INIT_ESTIMATION_CHECK_FIRST_OBS
+%token HETEROSKEDASTIC_FILTER
 /* Method of Moments */
 %token METHOD_OF_MOMENTS MOM_METHOD
 %token BARTLETT_KERNEL_LAG WEIGHTING_MATRIX WEIGHTING_MATRIX_SCALING_FACTOR ANALYTIC_STANDARD_ERRORS ANALYTIC_JACOBIAN PENALIZED_ESTIMATOR VERBOSE 
@@ -225,6 +226,7 @@ statement : parameters
           | init_param
           | shocks
           | mshocks
+          | heteroskedastic_shocks
           | sigma_e
           | steady
           | check
@@ -1042,6 +1044,24 @@ det_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
                  { driver.add_det_shock($2, $5, $8, false); }
                ;
 
+heteroskedastic_shocks : HETEROSKEDASTIC_SHOCKS ';' heteroskedastic_shock_list END ';'
+                         { driver.end_heteroskedastic_shocks(false); }
+                       | HETEROSKEDASTIC_SHOCKS '(' OVERWRITE ')' ';' heteroskedastic_shock_list END ';'
+                         { driver.end_heteroskedastic_shocks(true); }
+                       | HETEROSKEDASTIC_SHOCKS '(' OVERWRITE ')' ';'  END ';'
+                         { driver.end_heteroskedastic_shocks(true); }
+                       ;
+
+heteroskedastic_shock_list : heteroskedastic_shock_list heteroskedastic_shock_elem
+                           | heteroskedastic_shock_elem
+                           ;
+
+heteroskedastic_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
+                             { driver.add_heteroskedastic_shock($2, $5, $8, false); }
+                           | VAR symbol ';' PERIODS period_list ';' SCALES value_list ';'
+                             { driver.add_heteroskedastic_shock($2, $5, $8, true); }
+                           ;
+
 svar_identification : SVAR_IDENTIFICATION {driver.begin_svar_identification();} ';' svar_identification_list END ';'
                       { driver.end_svar_identification(); }
                     ;
@@ -2071,6 +2091,7 @@ estimation_options : o_datafile
                    | o_occbin_likelihood
                    | o_occbin_smoother
                    | o_no_init_estimation_check_first_obs
+                   | o_heteroskedastic_filter
                    ;
 
 name_value_pair : QUOTED_STRING COMMA QUOTED_STRING
@@ -3778,6 +3799,7 @@ o_use_shock_groups : USE_SHOCK_GROUPS { driver.option_str("plot_shock_decomp.use
 o_colormap : COLORMAP EQUAL symbol { driver.option_num("plot_shock_decomp.colormap",$3); };
 o_icd_colormap : COLORMAP EQUAL symbol { driver.option_num("initial_condition_decomp.colormap",$3); };
 o_no_init_estimation_check_first_obs : NO_INIT_ESTIMATION_CHECK_FIRST_OBS { driver.option_num("no_init_estimation_check_first_obs", "true"); };
+o_heteroskedastic_filter : HETEROSKEDASTIC_FILTER { driver.option_num("heteroskedastic_filter", "true"); };
 
 // Some options to "method_of_moments"
 o_bartlett_kernel_lag : BARTLETT_KERNEL_LAG EQUAL INT_NUMBER { driver.option_num("mom.bartlett_kernel_lag", $3); };
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 9f7c97244b63475dc4a8c3f79cfc5ace8ea5caa6..e87f72262e4f79daab82a6fb1346c20dd8830dc9 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -202,6 +202,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>histval {BEGIN DYNARE_BLOCK; return token::HISTVAL;}
 <INITIAL>filter_initial_state {BEGIN DYNARE_BLOCK; return token::FILTER_INITIAL_STATE;}
 <INITIAL>shocks {BEGIN DYNARE_BLOCK; return token::SHOCKS;}
+<INITIAL>heteroskedastic_shocks {BEGIN DYNARE_BLOCK; return token::HETEROSKEDASTIC_SHOCKS;}
 <INITIAL>shock_groups {BEGIN DYNARE_BLOCK; return token::SHOCK_GROUPS;}
 <INITIAL>init2shocks {BEGIN DYNARE_BLOCK; return token::INIT2SHOCKS;}
 <INITIAL>mshocks {BEGIN DYNARE_BLOCK; return token::MSHOCKS;}
@@ -714,6 +715,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>zero_moments_tolerance {return token::ZERO_MOMENTS_TOLERANCE;}
 <DYNARE_STATEMENT>max_nrows {return token::MAX_NROWS;}
 <DYNARE_STATEMENT>with_epilogue {return token::WITH_EPILOGUE;}
+<DYNARE_STATEMENT>heteroskedastic_filter {return token::HETEROSKEDASTIC_FILTER;}
 
 <DYNARE_STATEMENT>\$[^$]*\$ {
   strtok(yytext+1, "$");
@@ -727,6 +729,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>values {return token::VALUES;}
 <DYNARE_BLOCK>corr {return token::CORR;}
 <DYNARE_BLOCK>periods {return token::PERIODS;}
+<DYNARE_BLOCK>scales {return token::SCALES;}
 <DYNARE_BLOCK>cutoff {return token::CUTOFF;}
 <DYNARE_BLOCK>mfs	{return token::MFS;}
 <DYNARE_BLOCK>balanced_growth_test_tol {return token::BALANCED_GROWTH_TEST_TOL;}
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 05a0c6ec56d25bdaf18b58c49745d04660e2b595..4b0d95200d51d8bde69ca7272653e0b629b07f04 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -809,6 +809,15 @@ ParsingDriver::end_mshocks(bool overwrite)
   det_shocks.clear();
 }
 
+void
+ParsingDriver::end_heteroskedastic_shocks(bool overwrite)
+{
+  mod_file->addStatement(make_unique<HeteroskedasticShocksStatement>(overwrite, heteroskedastic_shocks_values,
+                                                                     heteroskedastic_shocks_scales, mod_file->symbol_table));
+  heteroskedastic_shocks_values.clear();
+  heteroskedastic_shocks_scales.clear();
+}
+
 void
 ParsingDriver::add_det_shock(const string &var, const vector<pair<int, int>> &periods, const vector<expr_t> &values, bool conditional_forecast)
 {
@@ -839,6 +848,30 @@ ParsingDriver::add_det_shock(const string &var, const vector<pair<int, int>> &pe
   det_shocks[symb_id] = v;
 }
 
+void
+ParsingDriver::add_heteroskedastic_shock(const string &var, const vector<pair<int, int>> &periods, const vector<expr_t> &values, bool scales)
+{
+  check_symbol_is_exogenous(var);
+
+  int symb_id = mod_file->symbol_table.getID(var);
+
+  if ((!scales && heteroskedastic_shocks_values.find(symb_id) != heteroskedastic_shocks_values.end())
+      || (scales && heteroskedastic_shocks_scales.find(symb_id) != heteroskedastic_shocks_scales.end()))
+    error("heteroskedastic_shocks: variable " + var + " declared twice");
+
+  if (periods.size() != values.size())
+    error("heteroskedastic_shocks: variable " + var + ": number of periods is different from number of shock values");
+
+  vector<tuple<int, int, expr_t>> v;
+  for (size_t i = 0; i < periods.size(); i++)
+    v.push_back({ periods[i].first, periods[i].second, values[i] });
+
+  if (scales)
+    heteroskedastic_shocks_scales[symb_id] = v;
+  else
+    heteroskedastic_shocks_values[symb_id] = v;
+}
+
 void
 ParsingDriver::add_stderr_shock(const string &var, expr_t value)
 {
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index bbb3d8e5691cc49002e641637a714c67671380fd..4402222ef830e8eff7641e897faf9a76fd5fd687 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -159,6 +159,8 @@ private:
   ShocksStatement::covar_and_corr_shocks_t covar_shocks;
   //! Temporary storage for correlations of shocks
   ShocksStatement::covar_and_corr_shocks_t corr_shocks;
+  //! Temporary storage for values and scales of heteroskedastic_shocks
+  HeteroskedasticShocksStatement::heteroskedastic_shocks_t heteroskedastic_shocks_values, heteroskedastic_shocks_scales;
   //! Temporary storage for Sigma_e rows
   SigmaeStatement::row_t sigmae_row;
   //! Temporary storage for Sigma_e matrix
@@ -434,8 +436,12 @@ public:
   void end_shocks(bool overwrite);
   //! Writes a mshocks statement
   void end_mshocks(bool overwrite);
+  //! Writes a heteroskedastic_shocks statement
+  void end_heteroskedastic_shocks(bool overwrite);
   //! Adds a deterministic shock or a path element inside a conditional_forecast_paths block
   void add_det_shock(const string &var, const vector<pair<int, int>> &periods, const vector<expr_t> &values, bool conditional_forecast);
+  //! Adds a heteroskedastic shock (either values or scales)
+  void add_heteroskedastic_shock(const string &var, const vector<pair<int, int>> &periods, const vector<expr_t> &values, bool scales);
   //! Adds a std error shock
   void add_stderr_shock(const string &var, expr_t value);
   //! Adds a variance shock
diff --git a/src/Shocks.cc b/src/Shocks.cc
index d10f012826d14bcc2ab5dfcf9336aa14b62e0a2f..dd9a770427c76364313c5f5fe3e23bdd3c378045 100644
--- a/src/Shocks.cc
+++ b/src/Shocks.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -721,3 +721,89 @@ Init2shocksStatement::writeJsonOutput(ostream &output) const
     }
   output << "]}";
 }
+
+HeteroskedasticShocksStatement::HeteroskedasticShocksStatement(bool overwrite_arg,
+                                                               const heteroskedastic_shocks_t &values_arg,
+                                                               const heteroskedastic_shocks_t &scales_arg,
+                                                               const SymbolTable &symbol_table_arg)
+  : overwrite{overwrite_arg}, values{values_arg}, scales{scales_arg}, symbol_table{symbol_table_arg}
+{
+}
+
+void
+HeteroskedasticShocksStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  if (overwrite)
+    output << "M_.heteroskedastic_shocks.Qhet = [];" << endl;
+
+  for (const auto &[var, vec] : values)
+    {
+      string varname = symbol_table.getName(var);
+      for (const auto &[period1, period2, value] : vec)
+        {
+          output << "M_.heteroskedastic_shocks.Qhet." << varname << ".time_value = " << period1 << ":" << period2 << ";" << endl
+                 << "M_.heteroskedastic_shocks.Qhet." << varname << ".value = ";
+          value->writeOutput(output);
+          output << ";" << endl;
+        }
+    }
+  for (const auto &[var, vec] : scales)
+    {
+      string varname = symbol_table.getName(var);
+      for (const auto &[period1, period2, scale] : vec)
+        {
+          output << "M_.heteroskedastic_shocks.Qhet." << varname << ".time_scale = " << period1 << ":" << period2 << ";" << endl
+                 << "M_.heteroskedastic_shocks.Qhet." << varname << ".scale = ";
+          scale->writeOutput(output);
+          output << ";" << endl;
+        }
+    }
+}
+
+void
+HeteroskedasticShocksStatement::writeJsonOutput(ostream &output) const
+{
+  output << R"({"statementName": "heteroskedastic_shocks")"
+         << R"(, "overwrite": )" << (overwrite ? "true" : "false")
+         << R"(, "shocks_values": [)";
+  for (auto it = values.begin(); it != values.end(); ++it)
+    {
+      if (it != values.begin())
+        output << ", ";
+      output << R"({"var": ")" << symbol_table.getName(it->first) << R"(", )"
+             << R"("values": [)";
+      for (auto it1 = it->second.begin(); it1 != it->second.end(); ++it1)
+        {
+          if (it1 != it->second.begin())
+            output << ", ";
+          auto [period1, period2, value] = *it1;
+          output << R"({"period1": )" << period1 << ", "
+                 << R"("period2": )" << period2 << ", "
+                 << R"("value": ")";
+          value->writeJsonOutput(output, {}, {});
+          output << R"("})";
+        }
+      output << "]}";
+    }
+  output << R"(], "shocks_scales": [)";
+  for (auto it = scales.begin(); it != scales.end(); ++it)
+    {
+      if (it != scales.begin())
+        output << ", ";
+      output << R"({"var": ")" << symbol_table.getName(it->first) << R"(", )"
+             << R"("scales": [)";
+      for (auto it1 = it->second.begin(); it1 != it->second.end(); ++it1)
+        {
+          if (it1 != it->second.begin())
+            output << ", ";
+          auto [period1, period2, value] = *it1;
+          output << R"({"period1": )" << period1 << ", "
+                 << R"("period2": )" << period2 << ", "
+                 << R"("value": ")";
+          value->writeJsonOutput(output, {}, {});
+          output << R"("})";
+        }
+      output << "]}";
+    }
+  output << "]}";
+}
diff --git a/src/Shocks.hh b/src/Shocks.hh
index 8d17a7b730d897acc583172136c6cdf13b01c0aa..63bc6165fdfecd1b329f95da18b1a9c93fe3781d 100644
--- a/src/Shocks.hh
+++ b/src/Shocks.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -181,4 +181,21 @@ public:
   void writeJsonOutput(ostream &output) const override;
 };
 
+class HeteroskedasticShocksStatement : public Statement
+{
+public:
+  // Maps exo symb_id to list of tuples (period1, period2, value/scale)
+  using heteroskedastic_shocks_t = map<int, vector<tuple<int, int, expr_t>>>;
+private:
+  const bool overwrite;
+  const heteroskedastic_shocks_t values, scales;
+  const SymbolTable &symbol_table;
+public:
+  HeteroskedasticShocksStatement(bool overwrite_arg, const heteroskedastic_shocks_t &values_arg,
+                                 const heteroskedastic_shocks_t &scales_arg,
+                                 const SymbolTable &symbol_table_arg);
+  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
+  void writeJsonOutput(ostream &output) const override;
+};
+
 #endif