diff --git a/src/ModFile.cc b/src/ModFile.cc
index c3e6d798e5f94dc66c770eee1505b3d1abec523e..00e0d2dfad7cc7d85fd2a9fbd9730f19061d40ba 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -984,8 +984,10 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
   // May be later modified by a shocks block
   mOutputFile << "M_.sigma_e_is_diagonal = true;" << endl;
 
-  // Initialize M_.det_shocks
-  mOutputFile << "M_.det_shocks = [];" << endl;
+  // Initialize M_.det_shocks and M_.heteroskedastic_shocks
+  mOutputFile << "M_.det_shocks = [];" << endl
+              << "M_.heteroskedastic_shocks.Qvalue_orig = [];" << endl
+              << "M_.heteroskedastic_shocks.Qscale_orig = [];" << endl;
 
   auto to_matlab_logical = [](bool m) { return m ? "true" : "false"; };
 
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 10eedf400a913ec108f8ad028a8df9af440cab63..bc063e2e4a02b9b7d1700f0f7f1c1c337c1cf578 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -870,7 +870,7 @@ ParsingDriver::add_heteroskedastic_shock(const string &var, const vector<pair<in
 
   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] });
+    v.emplace_back(periods[i].first, periods[i].second, values[i]);
 
   if (scales)
     heteroskedastic_shocks_scales[symb_id] = v;
diff --git a/src/Shocks.cc b/src/Shocks.cc
index 1ab09572dfafae470d36309f55ea3e19703ccb86..08a4a46486bc689cbafe1c959c300c064824404a 100644
--- a/src/Shocks.cc
+++ b/src/Shocks.cc
@@ -733,29 +733,33 @@ HeteroskedasticShocksStatement::HeteroskedasticShocksStatement(bool overwrite_ar
 void
 HeteroskedasticShocksStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
+  // NB: The first initialization of the fields is done in ModFile::writeMOutput()
   if (overwrite)
-    output << "M_.heteroskedastic_shocks.Qhet = [];" << endl;
+    output << "M_.heteroskedastic_shocks.Qvalue_orig = [];" << endl
+           << "M_.heteroskedastic_shocks.Qscale_orig = [];" << endl;
 
-  for (const auto &[var, vec] : values)
+  for (const auto &[symb_id, vec] : values)
     {
-      string varname = symbol_table.getName(var);
+      int tsid = symbol_table.getTypeSpecificID(symb_id);
       for (const auto &[period1, period2, value] : vec)
         {
-          output << "M_.heteroskedastic_shocks.Qhet." << varname << ".time_value = " << period1 << ":" << period2 << ";" << endl
-                 << "M_.heteroskedastic_shocks.Qhet." << varname << ".value = ";
+          output << "M_.heteroskedastic_shocks.Qvalue_orig = [M_.heteroskedastic_shocks.Qvalue_orig; struct('exo_id', "
+                 << tsid+1 << ",'periods',"
+                 << period1 << ":" << period2 << ",'value',";
           value->writeOutput(output);
-          output << ";" << endl;
+          output << ")];" << endl;
         }
     }
-  for (const auto &[var, vec] : scales)
+  for (const auto &[symb_id, vec] : scales)
     {
-      string varname = symbol_table.getName(var);
+      int tsid = symbol_table.getTypeSpecificID(symb_id);
       for (const auto &[period1, period2, scale] : vec)
         {
-          output << "M_.heteroskedastic_shocks.Qhet." << varname << ".time_scale = " << period1 << ":" << period2 << ";" << endl
-                 << "M_.heteroskedastic_shocks.Qhet." << varname << ".scale = ";
+          output << "M_.heteroskedastic_shocks.Qscale_orig = [M_.heteroskedastic_shocks.Qscale_orig; struct('exo_id', "
+                 << tsid+1 << ",'periods',"
+                 << period1 << ":" << period2 << ",'scale',";
           scale->writeOutput(output);
-          output << ";" << endl;
+          output << ")];" << endl;
         }
     }
 }