diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index ebb3d71c2b7bf60b06d23edda89df4b1538eb9e9..b42db040cf78b34fc9310edfb3c4955d2595309f 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3844,25 +3844,32 @@ DynamicModel::fillTrendComponentModelTable() const
     }
   trend_component_model_table.setRhs(rhsr);
   trend_component_model_table.setVals(eqnums, trend_eqnums, lhsr, lhs_expr_tr, nonstationaryr);
-
-  // Fill AR Matrix
-  map<string, map<tuple<int, int, int>, expr_t>> ARr, ECr;
-  fillAutoregressiveMatrix(ARr, true);
-  trend_component_model_table.setAR(ARr);
-  fillErrorComponentMatrix(ECr);
-  trend_component_model_table.setEC(ECr);
 }
 
 void
-DynamicModel::fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &ECr) const
+DynamicModel::fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &ECr,
+                                       ExprNode::subst_table_t &diff_subst_table) const
 {
-  for (const auto & it : trend_component_model_table.getNonTrendEqNums())
+  for (const auto & it : trend_component_model_table.getEqNums())
     {
       int i = 0;
       map<tuple<int, int, int>, expr_t> EC;
       vector<int> trend_lhs = trend_component_model_table.getTrendLhs(it.first);
+      vector<int> nontrend_eqnums = trend_component_model_table.getNonTrendEqNums(it.first);
+      vector<int> undiff_nontrend_lhs = getUndiffLHSForPac(it.first, diff_subst_table);
+      vector<int> parsed_undiff_nontrend_lhs;
+
+      for (auto eqn : it.second)
+        {
+          if (find(nontrend_eqnums.begin(), nontrend_eqnums.end(), eqn) != nontrend_eqnums.end())
+            parsed_undiff_nontrend_lhs.push_back(undiff_nontrend_lhs.at(i));
+          i++;
+        }
+
+      i = 0;
       for (auto eqn : it.second)
-        equations[eqn]->get_arg2()->fillErrorCorrectionRow(i++, trend_lhs, EC);
+        if (find(nontrend_eqnums.begin(), nontrend_eqnums.end(), eqn) != nontrend_eqnums.end())
+          equations[eqn]->get_arg2()->fillErrorCorrectionRow(i++, parsed_undiff_nontrend_lhs, trend_lhs, EC);
       ECr[it.first] = EC;
     }
 }
@@ -3934,6 +3941,16 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel(StaticModel &static_mode
   trend_component_model_table.setOrigDiffVar(orig_diff_var);
 }
 
+void
+DynamicModel::fillTrendComponentmodelTableAREC(ExprNode::subst_table_t &diff_subst_table) const
+{
+  map<string, map<tuple<int, int, int>, expr_t>> ARr, ECr;
+  fillAutoregressiveMatrix(ARr, true);
+  trend_component_model_table.setAR(ARr);
+  fillErrorComponentMatrix(ECr, diff_subst_table);
+  trend_component_model_table.setEC(ECr);
+}
+
 void
 DynamicModel::addEquationsForVar()
 {
@@ -3994,7 +4011,7 @@ DynamicModel::addEquationsForVar()
 
 vector<int>
 DynamicModel::getUndiffLHSForPac(const string &aux_model_name,
-                                 ExprNode::subst_table_t &diff_subst_table)
+                                 ExprNode::subst_table_t &diff_subst_table) const
 {
   vector<expr_t> lhs_expr_t = trend_component_model_table.getLhsExprT(aux_model_name);
   vector<int> lhs = trend_component_model_table.getLhs(aux_model_name);
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 936ba86539f439c4be6291f3899ca9befee7c063..5f5c67fe6c12bcfc4afb82370f4a8460d9fb379e 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -308,11 +308,12 @@ 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) const;
+  void fillErrorComponentMatrix(map<string, map<tuple<int, int, int>, expr_t>> &ECr, ExprNode::subst_table_t &diff_subst_table) const;
 
   //! Fill the Trend Component Model Table
   void fillTrendComponentModelTable() const;
   void fillTrendComponentModelTableFromOrigModel(StaticModel &static_model) const;
+  void fillTrendComponentmodelTableAREC(ExprNode::subst_table_t &diff_subst_table) const;
 
   //! Fill the Var Model Table
   void fillVarModelTable() const;
@@ -447,7 +448,7 @@ public:
 
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
-                                 ExprNode::subst_table_t &diff_subst_table);
+                                 ExprNode::subst_table_t &diff_subst_table) const;
 
   //! Transforms the model by replacing trend variables with a 1
   void removeTrendVariableFromEquations();
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 6f4878bf54e0828b78c6d3cf4bede91004c0fcff..63ec686d371b43b3cc497a59bd8ef499216fd1d5 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -701,7 +701,7 @@ NumConstNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<i
 }
 
 void
-NumConstNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
 }
 
@@ -1944,7 +1944,7 @@ VariableNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<i
 }
 
 void
-VariableNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
 }
 
@@ -3588,9 +3588,9 @@ UnaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<in
 }
 
 void
-UnaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
-  arg->fillErrorCorrectionRow(eqn, lhs, EC);
+  arg->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
 }
 
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, int idx_arg, const expr_t arg1_arg,
@@ -5500,75 +5500,87 @@ BinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<i
 
 void
 BinaryOpNode::fillErrorCorrectionRowHelper(expr_t arg1, expr_t arg2,
-                                          int eqn,
-                                          const vector<int> &lhs,
-                                          map<tuple<int, int, int>, expr_t> &EC) const
+                                           int eqn,
+                                           const vector<int> &nontrend_lhs,
+                                           const vector<int> &trend_lhs,
+                                           map<tuple<int, int, int>, expr_t> &EC) const
 {
   if (op_code != BinaryOpcode::times)
     return;
 
+  set<pair<int, int>> endogs, tmp;
+  arg1->collectDynamicVariables(SymbolType::endogenous, tmp);
+  arg1->collectDynamicVariables(SymbolType::exogenous, tmp);
+  if (tmp.size() != 0)
+    return;
+
   BinaryOpNode *multiplicandr = dynamic_cast<BinaryOpNode *>(arg2);
   if (multiplicandr == nullptr
       || multiplicandr->get_op_code() != BinaryOpcode::minus)
     return;
 
-  set<pair<int, int>> endogs, tmp;
   arg2->collectDynamicVariables(SymbolType::endogenous, endogs);
   if (endogs.size() != 2)
     return;
 
-  arg1->collectDynamicVariables(SymbolType::endogenous, tmp);
-  arg1->collectDynamicVariables(SymbolType::exogenous, tmp);
-  if (tmp.size() != 0)
+  arg2->collectDynamicVariables(SymbolType::exogenous, endogs);
+  arg2->collectDynamicVariables(SymbolType::parameter, endogs);
+  if (endogs.size() != 2)
     return;
 
-  bool encountered_trend_var = false;
-  int lhs_symb_id = -1;
+  int endog1, lag1, endog2, lag2;
+  tie(endog1, lag1) = *endogs.begin();
+  tie(endog2, lag2) = *next(endogs.begin(), 1);
+  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 max_lag = 0;
-  for (const auto & it : endogs)
+  int colidx = -1;
+  if (find(nontrend_lhs.begin(), nontrend_lhs.end(), endog1) != nontrend_lhs.end())
     {
-      bool isauxvar = datatree.symbol_table.isAuxiliaryVariable(it.first);
-      int tmp_symb_id = isauxvar ?
-        datatree.symbol_table.getOrigSymbIdForDiffAuxVar(it.first) : it.first;
-
-      if (find(lhs.begin(), lhs.end(), tmp_symb_id) != lhs.end())
-        if (encountered_trend_var)
-          {
-            cerr << "BinaryOpNode::fillErrorCorrectionRowHelper: Error filling EC matrix: "
-                 << "two trend variables encountered in EC term" << endl;
-            exit(EXIT_FAILURE);
-          }
-        else
-          {
-            encountered_trend_var = true;
-            lhs_symb_id = tmp_symb_id;
-            int tmp_lag = it.second;
-            if (isauxvar)
-              tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(it.first);
-            if (tmp_lag < max_lag)
-              max_lag = tmp_lag;
-          }
+      colidx = (int) distance(nontrend_lhs.begin(), find(nontrend_lhs.begin(), nontrend_lhs.end(), endog1));
+      int tmp_lag = lag2;
+      if (isauxvar2)
+        tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(orig_endog2);
+      if (tmp_lag < max_lag)
+        max_lag = tmp_lag;
     }
-
-  if (!encountered_trend_var)
+  else if (find(nontrend_lhs.begin(), nontrend_lhs.end(), endog2) != nontrend_lhs.end())
+    {
+      colidx = (int) distance(nontrend_lhs.begin(), find(nontrend_lhs.begin(), nontrend_lhs.end(), endog2));
+      int tmp_lag = lag1;
+      if (isauxvar1)
+        tmp_lag = -1 * datatree.symbol_table.getOrigLeadLagForDiffAuxVar(orig_endog1);
+      if (tmp_lag < max_lag)
+        max_lag = tmp_lag;
+    }
+  else
     return;
 
-  if (EC.find(make_tuple(eqn, -max_lag, lhs_symb_id)) != EC.end())
+  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);
     }
-  EC[make_tuple(eqn, -max_lag, lhs_symb_id)] = arg1;
+  EC[make_tuple(eqn, -max_lag, colidx)] = arg1;
 }
 
 void
-BinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
-  fillErrorCorrectionRowHelper(arg1, arg2, eqn, lhs, EC);
-  fillErrorCorrectionRowHelper(arg2, arg1, eqn, lhs, EC);
-  arg1->fillErrorCorrectionRow(eqn, lhs, EC);
-  arg2->fillErrorCorrectionRow(eqn, lhs, EC);
+  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);
 }
 
 void
@@ -6531,11 +6543,11 @@ TrinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<
 }
 
 void
-TrinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
-  arg1->fillErrorCorrectionRow(eqn, lhs, EC);
-  arg2->fillErrorCorrectionRow(eqn, lhs, EC);
-  arg3->fillErrorCorrectionRow(eqn, lhs, EC);
+  arg1->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
+  arg2->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
+  arg3->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, EC);
 }
 
 AbstractExternalFunctionNode::AbstractExternalFunctionNode(DataTree &datatree_arg,
@@ -7152,7 +7164,7 @@ AbstractExternalFunctionNode::fillAutoregressiveRow(int eqn, const vector<int> &
 }
 
 void
-AbstractExternalFunctionNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
   cerr << "External functions not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
@@ -8639,7 +8651,7 @@ VarExpectationNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<t
 }
 
 void
-VarExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) const
 {
   cerr << "Var Expectation not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
@@ -9133,7 +9145,7 @@ PacExpectationNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<t
 }
 
 void
-PacExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &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> &EC) 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 b0a33132f6bdbf55759877548eee00b69332e96f..679221d0e2a230b29b7da72f60892e4ad37181aa 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -572,7 +572,8 @@ class ExprNode
       virtual void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const = 0;
 
       //! Fills the EC matrix structure
-      virtual void fillErrorCorrectionRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &EC) const = 0;
+      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;
 
       //! Returns true if PacExpectationNode encountered
       virtual bool containsPacExpectation(const string &pac_model_name = "") const = 0;
@@ -661,7 +662,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
@@ -762,7 +763,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
@@ -887,7 +888,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
@@ -1043,8 +1044,9 @@ public:
                                    int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const;
   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> &lhs, map<tuple<int, int, int>, expr_t> &AR) const;
-  void fillErrorCorrectionRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
+                                    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;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacOptimizingPart(int lhs_orig_symb_id, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars,
                             set<pair<int, pair<int, int>>> &params_and_vars) const override;
@@ -1153,7 +1155,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
@@ -1276,7 +1278,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
@@ -1486,7 +1488,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
@@ -1584,7 +1586,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, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(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 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> &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> &EC) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
                                &params_vars_and_scaling_factor) const override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index c94cb80f3010757788accc50bb83ac567653c942..2aacae8d784ed61609e8599b71846c19a8ca6c81 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -396,6 +396,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   // Fill Trend Component Model Table
   dynamic_model.fillTrendComponentModelTable();
   original_model.fillTrendComponentModelTableFromOrigModel(diff_static_model);
+  dynamic_model.fillTrendComponentmodelTableAREC(diff_subst_table);
   dynamic_model.fillVarModelTable();
   original_model.fillVarModelTableFromOrigModel(diff_static_model);
 
diff --git a/src/SubModel.cc b/src/SubModel.cc
index ed6f43ff5a20f678c796037337ca7a140efeb4bc..9a1def8b5d15a5ca63cac4d76720806463d22118 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -356,9 +356,8 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
                    << "    ec = zeros(" << nontrend_lhs_vec.size() << ", " << nontrend_lhs_vec.size() << ", 1);" << endl;
       for (const auto & it : EC.at(name))
         {
-          int eqn, lag, lhs_symb_id;
-          tie (eqn, lag, lhs_symb_id) = it.first;
-          int colidx = (int) distance(trend_lhs_vec.begin(), find(trend_lhs_vec.begin(), trend_lhs_vec.end(), lhs_symb_id));
+          int eqn, lag, colidx;
+          tie (eqn, lag, colidx) = it.first;
           ar_ec_output << "    ec(" << eqn + 1 << ", " << colidx + 1 << ", 1) = ";
           it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 5d3c0a5b713bca84f224a8a499f8241f82be25eb..65016a7f1f68b9bda4eacbba620565994017bee6 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -45,7 +45,8 @@ private:
   map<string, vector<bool>> diff, nonstationary;
   map<string, vector<expr_t>> lhs_expr_t;
   map<string, vector<int>> trend_vars;
-  map<string, map<tuple<int, int, int>, expr_t>> AR, EC; // AR/EC: name -> (eqn, lag, lhs_symb_id) -> expr_t
+  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
 public:
   TrendComponentModelTable(SymbolTable &symbol_table_arg);