diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 0b2c82500186b6856f968822e6ec0c370c38cb97..f5f2f8e9aea9be44b0b7d70b2cd73cc7a77bd62a 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4016,13 +4016,14 @@ DynamicModel::fillTrendComponentModelTable() const
 }
 
 void
-DynamicModel::fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &ECr,
+DynamicModel::fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &A0r,
+                                       map<string, map<tuple<int, int, int>, expr_t>> &A0starr,
                                        ExprNode::subst_table_t &diff_subst_table) const
 {
   for (const auto & it : trend_component_model_table.getEqNums())
     {
       int i = 0;
-      map<tuple<int, int, int>, expr_t> EC;
+      map<tuple<int, int, int>, expr_t> A0, A0star;
       vector<int> target_lhs = trend_component_model_table.getTargetLhs(it.first);
       vector<int> nontarget_eqnums = trend_component_model_table.getNonTargetEqNums(it.first);
       vector<int> undiff_nontarget_lhs = getUndiffLHSForPac(it.first, diff_subst_table);
@@ -4038,8 +4039,9 @@ DynamicModel::fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, exp
       i = 0;
       for (auto eqn : it.second)
         if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
-          equations[eqn]->arg2->fillErrorCorrectionRow(i++, parsed_undiff_nontarget_lhs, target_lhs, EC);
-      ECr[it.first] = EC;
+          equations[eqn]->arg2->fillErrorCorrectionRow(i++, parsed_undiff_nontarget_lhs, target_lhs, A0, A0star);
+      A0r[it.first] = A0;
+      A0starr[it.first] = A0star;
     }
 }
 
@@ -4113,11 +4115,11 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel(StaticModel &static_mode
 void
 DynamicModel::fillTrendComponentmodelTableAREC(ExprNode::subst_table_t &diff_subst_table) const
 {
-  map<string, map<tuple<int, int, int>, expr_t>> ARr, ECr;
+  map<string, map<tuple<int, int, int>, expr_t>> ARr, A0r, A0starr;
   fillAutoregressiveMatrix(ARr, true);
   trend_component_model_table.setAR(ARr);
-  fillErrorComponentMatrix(ECr, diff_subst_table);
-  trend_component_model_table.setEC(ECr);
+  fillErrorComponentMatrix(A0r, A0starr, diff_subst_table);
+  trend_component_model_table.setA0(A0r, A0starr);
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index a4b665c13c6179ace5664faaf36e1dbefa57a105..03b168998d7ed82f0ae24b1792870cbcd4aefa30 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -323,7 +323,7 @@ public:
   void fillAutoregressiveMatrix(map<string, map<tuple<int, int, int>, expr_t>> &ARr, bool is_trend_component_model) const;
 
   //! Fill Error Component Matrix for trend_component_model
-  void fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &ECr, ExprNode::subst_table_t &diff_subst_table) const;
+  void fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &A0r, map<string, map<tuple<int, int, int>, expr_t>> &A0starr, ExprNode::subst_table_t &diff_subst_table) const;
 
   //! Fill the Trend Component Model Table
   void fillTrendComponentModelTable() const;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index c19698abf2cbd0e62ea50c92c4213b07fe884007..502563ee341cac399cbd95ce729c9ecb517bc55d 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -722,7 +722,7 @@ NumConstNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<i
 }
 
 void
-NumConstNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+NumConstNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
 }
 
@@ -2045,7 +2045,7 @@ VariableNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<i
 }
 
 void
-VariableNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+VariableNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
 }
 
@@ -3891,9 +3891,9 @@ UnaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<in
 }
 
 void
-UnaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+UnaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
-  arg->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
+  arg->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
 }
 
 void
@@ -5857,7 +5857,8 @@ BinaryOpNode::fillErrorCorrectionRowHelper(expr_t arg1, expr_t arg2,
                                            int eqn,
                                            const vector<int> &nontarget_lhs,
                                            const vector<int> &target_lhs,
-                                           map<tuple<int, int, int>, expr_t> &EC) const
+                                           map<tuple<int, int, int>, expr_t> &A0,
+                                           map<tuple<int, int, int>, expr_t> &A0star) const
 {
   if (op_code != BinaryOpcode::times)
     return;
@@ -5872,64 +5873,130 @@ BinaryOpNode::fillErrorCorrectionRowHelper(expr_t arg1, expr_t arg2,
   if (tmp.size() != 1)
     return;
 
+  auto *vn = dynamic_cast<VariableNode *>(arg2);
   auto *bopn = dynamic_cast<BinaryOpNode *>(arg2);
-  if (bopn == nullptr || bopn->op_code != BinaryOpcode::minus)
-    return;
-
-  arg2->collectDynamicVariables(SymbolType::endogenous, endogs);
-  if (endogs.size() != 2)
+  if ((bopn == nullptr || bopn->op_code != BinaryOpcode::minus)
+      && vn == nullptr)
     return;
 
-  arg2->collectDynamicVariables(SymbolType::exogenous, endogs);
-  arg2->collectDynamicVariables(SymbolType::parameter, endogs);
-  if (endogs.size() != 2)
+  if (bopn != nullptr)
     {
-      cerr << "ERROR in model; expecting param*endog or param*(endog-endog)" << endl;
-      exit(EXIT_FAILURE);
-    }
+      arg2->collectDynamicVariables(SymbolType::endogenous, endogs);
+      if (endogs.size() != 2)
+        return;
 
-  auto *vn1 = dynamic_cast<VariableNode *>(bopn->arg1);
-  auto *vn2 = dynamic_cast<VariableNode *>(bopn->arg2);
-  if (vn1 == nullptr || vn2 == nullptr)
-    {
-      cerr << "ERROR in model; expecting param*endog or param*(endog-endog)" << endl;
-      exit(EXIT_FAILURE);
-    }
+      arg2->collectDynamicVariables(SymbolType::exogenous, endogs);
+      arg2->collectDynamicVariables(SymbolType::parameter, endogs);
+      if (endogs.size() != 2)
+        {
+          cerr << "ERROR in model; expecting param*endog or param*(endog-endog)" << endl;
+          exit(EXIT_FAILURE);
+        }
 
-  int endog1 = vn1->symb_id;
-  int endog2 = vn2->symb_id;
-  int orig_endog2 = endog2;
+      auto *vn1 = dynamic_cast<VariableNode *>(bopn->arg1);
+      auto *vn2 = dynamic_cast<VariableNode *>(bopn->arg2);
+      if (vn1 == nullptr || vn2 == nullptr)
+        {
+          cerr << "ERROR in model; expecting param*endog or param*(endog-endog)" << endl;
+          exit(EXIT_FAILURE);
+        }
 
-  bool isauxvar1 = datatree.symbol_table.isAuxiliaryVariable(endog1);
-  endog1 = isauxvar1 ?
-    datatree.symbol_table.getOrigSymbIdForDiffAuxVar(endog1) : endog1;
+      int endog1 = vn1->symb_id;
+      int endog2 = vn2->symb_id;
+      int orig_endog1 = endog1;
+      int orig_endog2 = endog2;
+
+      bool isauxvar1 = datatree.symbol_table.isAuxiliaryVariable(endog1);
+      endog1 = isauxvar1 ?
+        datatree.symbol_table.getOrigSymbIdForDiffAuxVar(endog1) : endog1;
+
+      bool isauxvar2 = datatree.symbol_table.isAuxiliaryVariable(endog2);
+      endog2 = isauxvar2 ?
+        datatree.symbol_table.getOrigSymbIdForDiffAuxVar(endog2) : endog2;
+
+      int A0_max_lag = vn1->lag;
+      int A0star_max_lag = vn2->lag;
+      int A0_colidx = -1;
+      int A0star_colidx = -1;
+      if (find(nontarget_lhs.begin(), nontarget_lhs.end(), endog1) != nontarget_lhs.end()
+          && find(target_lhs.begin(), target_lhs.end(), endog2) != target_lhs.end())
+        {
+          A0_colidx = (int) distance(nontarget_lhs.begin(), find(nontarget_lhs.begin(), nontarget_lhs.end(), endog1));
+          int tmp_lag = vn1->lag;
+          if (isauxvar1)
+            tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(orig_endog1);
+          if (tmp_lag < A0_max_lag)
+            A0_max_lag = tmp_lag;
+
+          A0star_colidx = (int) distance(target_lhs.begin(), find(target_lhs.begin(), target_lhs.end(), endog2));
+          tmp_lag = vn2->lag;
+          if (isauxvar2)
+            tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(orig_endog2);
+          if (tmp_lag < A0star_max_lag)
+            A0star_max_lag = tmp_lag;
+        }
+      else
+        return;
 
-  bool isauxvar2 = datatree.symbol_table.isAuxiliaryVariable(endog2);
-  endog2 = isauxvar2 ?
-    datatree.symbol_table.getOrigSymbIdForDiffAuxVar(endog2) : endog2;
+      if (A0.find(make_tuple(eqn, -A0_max_lag, A0_colidx)) != A0.end())
+        {
+          cerr << "BinaryOpNode::fillErrorCorrectionRowHelper: Error filling A0 matrix: "
+               << "lag/symb_id encountered more than once in equtaion" << endl;
+          exit(EXIT_FAILURE);
+        }
 
-  int max_lag = vn2->lag;
-  int colidx = -1;
-  if (find(nontarget_lhs.begin(), nontarget_lhs.end(), endog1) != nontarget_lhs.end()
-      && find(target_lhs.begin(), target_lhs.end(), endog2) != target_lhs.end())
-    {
-      colidx = (int) distance(target_lhs.begin(), find(target_lhs.begin(), target_lhs.end(), endog2));
-      int tmp_lag = vn2->lag;
-      if (isauxvar2)
-        tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(orig_endog2);
-      if (tmp_lag < max_lag)
-        max_lag = tmp_lag;
+      if (A0star.find(make_tuple(eqn, -A0star_max_lag, A0star_colidx)) != A0star.end())
+        {
+          cerr << "BinaryOpNode::fillErrorCorrectionRowHelper: Error filling A0star matrix: "
+               << "lag/symb_id encountered more than once in equtaion" << endl;
+          exit(EXIT_FAILURE);
+        }
+      A0[make_tuple(eqn, -A0_max_lag, A0_colidx)] = arg1;
+      A0star[make_tuple(eqn, -A0star_max_lag, A0star_colidx)] = arg1;
     }
   else
-    return;
-
-  if (EC.find(make_tuple(eqn, -max_lag, colidx)) != EC.end())
     {
-      cerr << "BinaryOpNode::fillErrorCorrectionRowHelper: Error filling EC matrix: "
-           << "lag/symb_id encountered more than once in equtaion" << endl;
-      exit(EXIT_FAILURE);
+      arg2->collectDynamicVariables(SymbolType::endogenous, endogs);
+      if (endogs.size() != 1)
+        return;
+
+      arg2->collectDynamicVariables(SymbolType::exogenous, endogs);
+      arg2->collectDynamicVariables(SymbolType::parameter, endogs);
+      if (endogs.size() != 1)
+        {
+          cerr << "ERROR in model; expecting param*endog or param*(endog-endog)" << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      int endog1 = vn->symb_id;
+      int orig_endog1 = endog1;
+
+      bool isauxvar1 = datatree.symbol_table.isAuxiliaryVariable(endog1);
+      endog1 = isauxvar1 ?
+        datatree.symbol_table.getOrigSymbIdForDiffAuxVar(endog1) : endog1;
+
+      int max_lag = vn->lag;
+      int colidx = -1;
+      if (find(nontarget_lhs.begin(), nontarget_lhs.end(), endog1) != nontarget_lhs.end())
+        {
+          colidx = (int) distance(nontarget_lhs.begin(), find(nontarget_lhs.begin(), nontarget_lhs.end(), endog1));
+          int tmp_lag = vn->lag;
+          if (isauxvar1)
+            tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(orig_endog1);
+          if (tmp_lag < max_lag)
+            max_lag = tmp_lag;
+        }
+      else
+        return;
+
+      if (A0.find(make_tuple(eqn, -max_lag, colidx)) != A0.end())
+        {
+          cerr << "BinaryOpNode::fillErrorCorrectionRowHelper: Error filling A0 matrix: "
+               << "lag/symb_id encountered more than once in equtaion" << endl;
+          exit(EXIT_FAILURE);
+        }
+      A0[make_tuple(eqn, -max_lag, colidx)] = arg1;
     }
-  EC[make_tuple(eqn, -max_lag, colidx)] = arg1;
 }
 
 void
@@ -5960,12 +6027,14 @@ BinaryOpNode::replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table)
 }
 
 void
-BinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+BinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
-  fillErrorCorrectionRowHelper(arg1, arg2, eqn, nontrend_lhs, trend_lhs, EC);
-  fillErrorCorrectionRowHelper(arg2, arg1, eqn, nontrend_lhs, trend_lhs, EC);
-  arg1->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
-  arg2->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
+  int A0size = A0.size();
+  fillErrorCorrectionRowHelper(arg1, arg2, eqn, nontrend_lhs, trend_lhs, A0, A0star);
+  if (A0size == A0.size())
+    fillErrorCorrectionRowHelper(arg2, arg1, eqn, nontrend_lhs, trend_lhs, A0, A0star);
+  arg1->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
+  arg2->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
 }
 
 void
@@ -6974,11 +7043,11 @@ TrinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<
 }
 
 void
-TrinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+TrinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
-  arg1->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
-  arg2->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
-  arg3->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
+  arg1->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
+  arg2->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
+  arg3->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
 }
 
 void
@@ -7651,7 +7720,7 @@ AbstractExternalFunctionNode::fillAutoregressiveRow(int eqn, const vector<int> &
 }
 
 void
-AbstractExternalFunctionNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+AbstractExternalFunctionNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
   cerr << "External functions not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
@@ -9193,7 +9262,7 @@ VarExpectationNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<t
 }
 
 void
-VarExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+VarExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
   cerr << "Var Expectation not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
@@ -9719,7 +9788,7 @@ PacExpectationNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<t
 }
 
 void
-PacExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const
+PacExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const
 {
   cerr << "Pac Expectation not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 6c7f4aff1e93f0396c1d7c56981e88112e356e60..dcbdc144ad90668bcbf7e13e92fd65a83a67d358 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -618,7 +618,8 @@ class ExprNode
 
       //! Fills the EC matrix structure
       virtual void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs,
-                                          map<tuple<int, int, int>, expr_t> &EC) const = 0;
+                                          map<tuple<int, int, int>, expr_t> &A0,
+                                          map<tuple<int, int, int>, expr_t> &A0star) const = 0;
 
       //! Finds equations where a variable is equal to a constant
       virtual void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const = 0;
@@ -742,7 +743,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -836,7 +837,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -958,7 +959,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -1097,8 +1098,8 @@ public:
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRowHelper(expr_t arg1, expr_t arg2,
                                     int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs,
-                                    map<tuple<int, int, int>, expr_t> &AR) const;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+                                    map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -1217,7 +1218,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -1348,7 +1349,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -1567,7 +1568,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
@@ -1672,7 +1673,7 @@ public:
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
   void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+  void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 267b8242ac8e9b08c6cecc90930bd5b4846092e4..526492463c79ef7159d37a77eadc95aa3ee6fa71 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -111,9 +111,11 @@ TrendComponentModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> A
 }
 
 void
-TrendComponentModelTable::setEC(map<string, map<tuple<int, int, int>, expr_t>> EC_arg)
+TrendComponentModelTable::setA0(map<string, map<tuple<int, int, int>, expr_t>> A0_arg,
+                                map<string, map<tuple<int, int, int>, expr_t>> A0star_arg)
 {
-  EC = move(EC_arg);
+  A0 = move(A0_arg);
+  A0star = move(A0star_arg);
 }
 
 map<string, vector<string>>
@@ -249,7 +251,7 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
   if (names.empty())
     return;
 
-  string filename = "+" + basename + "/trend_component_ar_ec.m";
+  string filename = "+" + basename + "/trend_component_ar_a0.m";
   ofstream ar_ec_output;
   ar_ec_output.open(filename, ios::out | ios::binary);
   if (!ar_ec_output.is_open())
@@ -257,8 +259,8 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
       cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  ar_ec_output << "function [ar, ec] = trend_component_ar_ec(model_name, params)" << endl
-            << "%function [ar, ec] = trend_component_ar_ec(model_name, params)" << endl
+  ar_ec_output << "function [AR, A0, A0star] = trend_component_ar_a0(model_name, params)" << endl
+            << "%function [AR, A0, A0star] = trend_component_ar_a0(model_name, params)" << endl
             << "% File automatically generated by the Dynare preprocessor" << endl << endl;
 
   for (const auto &name : names)
@@ -328,32 +330,49 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
 
       ar_ec_output << "if strcmp(model_name, '" << name << "')" << endl
                 << "    % AR" << endl
-                << "    ar = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << getMaxLag(name) << ");" << endl;
+                << "    AR = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << getMaxLag(name) << ");" << endl;
       for (const auto & it : AR.at(name))
         {
           int eqn, lag, lhs_symb_id;
           tie (eqn, lag, lhs_symb_id) = it.first;
           int colidx = (int) distance(nontarget_lhs_vec.begin(), find(nontarget_lhs_vec.begin(), nontarget_lhs_vec.end(), lhs_symb_id));
-          ar_ec_output << "    ar(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
+          ar_ec_output << "    AR(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
           it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
         }
 
-      int ec_lag = 0;
-      for (const auto & it : EC.at(name))
-        if (get<1>(it.first) > ec_lag)
-          ec_lag = get<1>(it.first);
+      int a0_lag = 0;
+      for (const auto & it : A0.at(name))
+        if (get<1>(it.first) > a0_lag)
+          a0_lag = get<1>(it.first);
       ar_ec_output << endl
-                   << "    % EC" << endl
-                   << "    ec = zeros(" << nontarget_lhs_vec.size() << ", " << target_lhs_vec.size() << ", " << ec_lag << ");" << endl;
-      for (const auto & it : EC.at(name))
+                   << "    % A0" << endl
+                   << "    A0 = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << a0_lag << ");" << endl;
+      for (const auto & it : A0.at(name))
         {
           int eqn, lag, colidx;
           tie (eqn, lag, colidx) = it.first;
-          ar_ec_output << "    ec(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
+          ar_ec_output << "    A0(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
           it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
         }
+
+      int a0star_lag = 0;
+      for (const auto & it : A0star.at(name))
+        if (get<1>(it.first) > a0star_lag)
+          a0star_lag = get<1>(it.first);
+      ar_ec_output << endl
+                   << "    % A0star" << endl
+                   << "    A0star = zeros(" << nontarget_lhs_vec.size() << ", " << target_lhs_vec.size() << ", " << a0star_lag << ");" << endl;
+      for (const auto & it : A0star.at(name))
+        {
+          int eqn, lag, colidx;
+          tie (eqn, lag, colidx) = it.first;
+          ar_ec_output << "    A0star(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
+          it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_ec_output << ";" << endl;
+        }
+
       ar_ec_output << "    return" << endl
                 << "end" << endl << endl;
     }
diff --git a/src/SubModel.hh b/src/SubModel.hh
index b3883f5ecc6afe5077f739e440ca657a74b5c0e3..5444cf125221d19c667375a7c8056afc5e5f1838 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -46,7 +46,7 @@ private:
   map<string, vector<expr_t>> lhs_expr_t;
   map<string, vector<int>> target_vars;
   map<string, map<tuple<int, int, int>, expr_t>> AR; // AR: name -> (eqn, lag, lhs_symb_id) -> expr_t
-  map<string, map<tuple<int, int, int>, expr_t>> EC; // EC: name -> (eqn, lag, col) -> expr_t
+  map<string, map<tuple<int, int, int>, expr_t>> A0, A0star; // EC: name -> (eqn, lag, col) -> expr_t
 public:
   explicit TrendComponentModelTable(SymbolTable &symbol_table_arg);
 
@@ -84,7 +84,8 @@ public:
   void setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg);
   void setTargetVar(map<string, vector<int>> target_vars_arg);
   void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
-  void setEC(map<string, map<tuple<int, int, int>, expr_t>> EC_arg);
+  void setA0(map<string, map<tuple<int, int, int>, expr_t>> A0_arg,
+             map<string, map<tuple<int, int, int>, expr_t>> A0star_arg);
 
   //! Write output of this class
   void writeOutput(const string &basename, ostream &output) const;