diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 1bf017a7451e2ae4080454738138bf3ffc1a9ac1..755f43455b3e72845e914288568dbe1918a6dcf5 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3849,9 +3849,29 @@ DynamicModel::fillTrendComponentModelTable() const
   trend_component_model_table.setNonstationary(nonstationaryr);
 
   // Fill AR Matrix
-  map<string, map<tuple<int, int, int>, expr_t>> ARr;
+  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
+{
+  for (const auto & it : trend_component_model_table.getEqNums())
+    {
+      vector<int> nontrend_lhs;
+      vector<int> lhsv = trend_component_model_table.getLhs(it.first);
+      for (int trend_it : trend_component_model_table.getTrendEqNums(it.first))
+        nontrend_lhs.push_back(lhsv.at(distance(it.second.begin(), find(it.second.begin(), it.second.end(), trend_it))));
+
+      int i = 0;
+      map<tuple<int, int, int>, expr_t> EC;
+      for (auto eqn : it.second)
+        equations[eqn]->get_arg2()->fillErrorCorrectionRow(i++, nontrend_lhs, EC);
+      ECr[it.first] = EC;
+    }
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index b7bca9050e42bf498b91ce693e0c8570a5e38f52..936ba86539f439c4be6291f3899ca9befee7c063 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -304,9 +304,12 @@ public:
   //! Set the equations that have non-zero second derivatives
   void setNonZeroHessianEquations(map<int, string> &eqs);
 
-  //! Fill Autoregressive Matrix for var_model
+  //! Fill Autoregressive Matrix for var_model/trend_component_model
   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;
+
   //! Fill the Trend Component Model Table
   void fillTrendComponentModelTable() const;
   void fillTrendComponentModelTableFromOrigModel(StaticModel &static_model) const;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index ade8d75140dbb5ec793f4810891b0b879178e81d..f2af63bfa94a22cf3c4434f778c96775ec29171e 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -700,6 +700,11 @@ 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
+{
+}
+
 VariableNode::VariableNode(DataTree &datatree_arg, int idx_arg, int symb_id_arg, int lag_arg) :
   ExprNode(datatree_arg, idx_arg),
   symb_id(symb_id_arg),
@@ -1938,6 +1943,11 @@ 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
+{
+}
+
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, int idx_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, int param1_symb_id_arg, int param2_symb_id_arg, string adl_param_name_arg, vector<int> adl_lags_arg) :
   ExprNode(datatree_arg, idx_arg),
   arg(arg_arg),
@@ -3577,6 +3587,12 @@ UnaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<in
   arg->fillAutoregressiveRow(eqn, lhs, AR);
 }
 
+void
+UnaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &EC) const
+{
+  arg->fillErrorCorrectionRow(eqn, lhs, EC);
+}
+
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, int idx_arg, const expr_t arg1_arg,
                            BinaryOpcode op_code_arg, const expr_t arg2_arg, int powerDerivOrder_arg) :
   ExprNode(datatree_arg, idx_arg),
@@ -5481,6 +5497,70 @@ BinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<i
   arg2->fillAutoregressiveRow(eqn, lhs, AR);
 }
 
+void
+BinaryOpNode::fillErrorCorrectionRowHelper(expr_t arg1, expr_t arg2,
+                                          int eqn,
+                                          const vector<int> &lhs,
+                                          map<tuple<int, int, int>, expr_t> &EC) const
+{
+  if (op_code != BinaryOpcode::times)
+    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)
+    return;
+
+  bool encountered_trend_var = false;
+  int lhs_symb_id = -1;
+  int max_lag = 0;
+  for (const auto & it : endogs)
+    if (find(lhs.begin(), lhs.end(), it.first) != 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
+        {
+          lhs_symb_id = it.first;
+          if (it.second < max_lag)
+            max_lag = it.second;
+          encountered_trend_var = true;
+        }
+
+  if (!encountered_trend_var)
+    return;
+
+  if (EC.find(make_tuple(eqn, -max_lag, lhs_symb_id)) != 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;
+}
+
+void
+BinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &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);
+}
+
 void
 BinaryOpNode::getPacLHS(pair<int, int> &lhs)
 {
@@ -6440,6 +6520,14 @@ TrinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<
   arg3->fillAutoregressiveRow(eqn, lhs, AR);
 }
 
+void
+TrinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &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);
+}
+
 AbstractExternalFunctionNode::AbstractExternalFunctionNode(DataTree &datatree_arg,
                                                            int idx_arg,
                                                            int symb_id_arg,
@@ -7049,7 +7137,14 @@ AbstractExternalFunctionNode::substituteStaticAuxiliaryVariable() const
 void
 AbstractExternalFunctionNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const
 {
-  cerr << "External functions not supported in VARs" << endl;
+  cerr << "External functions not supported in VARs or Trend Component Models" << endl;
+  exit(EXIT_FAILURE);
+}
+
+void
+AbstractExternalFunctionNode::fillErrorCorrectionRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &EC) const
+{
+  cerr << "External functions not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
 }
 
@@ -8529,7 +8624,14 @@ VarExpectationNode::substituteStaticAuxiliaryVariable() const
 void
 VarExpectationNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const
 {
-  cerr << "Var Expectation not supported in VARs" << endl;
+  cerr << "Var Expectation not supported in VARs or Trend Component Models" << endl;
+  exit(EXIT_FAILURE);
+}
+
+void
+VarExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &EC) const
+{
+  cerr << "Var Expectation not supported in Trend Component Models" << endl;
   exit(EXIT_FAILURE);
 }
 
@@ -9020,6 +9122,13 @@ PacExpectationNode::fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<t
   exit(EXIT_FAILURE);
 }
 
+void
+PacExpectationNode::fillErrorCorrectionRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &EC) const
+{
+  cerr << "Pac Expectation not supported in Trend Component Models" << endl;
+  exit(EXIT_FAILURE);
+}
+
 void
 PacExpectationNode::writeJsonOutput(ostream &output,
                                     const temporary_terms_t &temporary_terms,
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 4f1f5f1991701b28504a0d49840f320a65fa853e..b0a33132f6bdbf55759877548eee00b69332e96f 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -571,6 +571,9 @@ class ExprNode
       //! Fills the AR matrix structure
       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;
+
       //! Returns true if PacExpectationNode encountered
       virtual bool containsPacExpectation(const string &pac_model_name = "") const = 0;
 
@@ -658,6 +661,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;
   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;
@@ -758,6 +762,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;
   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;
@@ -882,6 +887,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;
   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;
@@ -1036,6 +1042,9 @@ public:
   void fillAutoregressiveRowHelper(expr_t arg1, expr_t arg2,
                                    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;
   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;
@@ -1144,6 +1153,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;
   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;
@@ -1266,6 +1276,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;
   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;
@@ -1475,6 +1486,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;
   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;
@@ -1572,6 +1584,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;
   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/SubModel.cc b/src/SubModel.cc
index 917eaf84b84da1becfb7f7f661b801f469a3bde5..a4da2b3ea02710131706c47843845ab5118b7ae6 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -125,6 +125,19 @@ TrendComponentModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> A
   AR = move(AR_arg);
 }
 
+void
+TrendComponentModelTable::setEC(map<string, map<tuple<int, int, int>, expr_t>> EC_arg)
+{
+  for (const auto & itmap1 : EC_arg)
+    for (const auto &itmap2 : itmap1.second)
+      if (get<1>(itmap2.first) != 1)
+        {
+          cerr << "Error in EC component: lag must be equal to 1" << endl;
+          exit(EXIT_FAILURE);
+        }
+  EC = move(EC_arg);
+}
+
 map<string, vector<string>>
 TrendComponentModelTable::getEqTags() const
 {
@@ -188,6 +201,19 @@ TrendComponentModelTable::getTrendEqNums() const
   return trend_eqnums;
 }
 
+vector<int>
+TrendComponentModelTable::getTrendEqNums(const string &name_arg) const
+{
+  checkModelName(name_arg);
+  return trend_eqnums.find(name_arg)->second;
+}
+
+map<string, vector<int>>
+TrendComponentModelTable::getNonTrendEqNums() const
+{
+  return nontrend_eqnums;
+}
+
 vector<int>
 TrendComponentModelTable::getNonTrendEqNums(const string &name_arg) const
 {
@@ -235,16 +261,16 @@ TrendComponentModelTable::getOrigDiffVar(const string &name_arg) const
 void
 TrendComponentModelTable::writeOutput(const string &basename, ostream &output) const
 {
-  string filename = "+" + basename + "/trend_component_ar.m";
-  ofstream ar_output;
-  ar_output.open(filename, ios::out | ios::binary);
-  if (!ar_output.is_open())
+  string filename = "+" + basename + "/trend_component_ar_ec.m";
+  ofstream ar_ec_output;
+  ar_ec_output.open(filename, ios::out | ios::binary);
+  if (!ar_ec_output.is_open())
     {
       cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  ar_output << "function ar = trend_component_ar(model_name, params)" << endl
-            << "%function ar = trend_component_ar(model_name, params)" << endl
+  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
             << "% File automatically generated by the Dynare preprocessor" << endl << endl;
 
   for (const auto &name : names)
@@ -309,29 +335,45 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
         output << (it >= 0 ? symbol_table.getTypeSpecificID(it) + 1 : -1) << " ";
       output << "];" << endl;
 
-      vector<int> nontrend_lhs;
+      vector<int> nontrend_lhs, trend_lhs;
       vector<int> lhsv = getLhs(name);
       vector<int> eqnumsv = getEqNums(name);
       for (int nontrend_it : getNonTrendEqNums(name))
         nontrend_lhs.push_back(lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), nontrend_it))));
 
-      ar_output << "if strcmp(model_name, '" << name << "')" << endl
+      for (int trend_it : getTrendEqNums(name))
+        trend_lhs.push_back(lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), trend_it))));
+
+      ar_ec_output << "if strcmp(model_name, '" << name << "')" << endl
+                << "    % AR" << endl
                 << "    ar = zeros(" << nontrend_lhs.size() << ", " << nontrend_lhs.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(nontrend_lhs.begin(), find(nontrend_lhs.begin(), nontrend_lhs.end(), lhs_symb_id));
-          ar_output << "    ar(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
-          it.second->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
-          ar_output << ";" << endl;
+          ar_ec_output << "    ar(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
+          it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_ec_output << ";" << endl;
         }
-      ar_output << "    return" << endl
+      ar_ec_output << endl
+                   << "    % EC" << endl
+                   << "    ec = zeros(" << nontrend_lhs.size() << ", " << nontrend_lhs.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.begin(), find(trend_lhs.begin(), trend_lhs.end(), lhs_symb_id));
+          ar_ec_output << "    ec(" << eqn + 1 << ", " << colidx + 1 << ", 1) = ";
+          it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_ec_output << ";" << endl;
+        }
+      ar_ec_output << "    return" << endl
                 << "end" << endl << endl;
     }
-  ar_output << "error([model_name ' is not a valid trend_component_model name'])" << endl
+  ar_ec_output << "error([model_name ' is not a valid trend_component_model name'])" << endl
             << "end" << endl;
-  ar_output.close();
+  ar_ec_output.close();
 }
 
 void
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 649b389698159322fb260118cea53da345a2bd16..c3a0bc57b360cea52f8d8c6da2bf51e6b90d0a4b 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -45,7 +45,7 @@ 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; // AR: name -> (eqn, lag, lhs_symb_id) -> param_expr_t
+  map<string, map<tuple<int, int, int>, expr_t>> AR, EC; // AR/EC: name -> (eqn, lag, lhs_symb_id) -> expr_t
 public:
   TrendComponentModelTable(SymbolTable &symbol_table_arg);
 
@@ -61,6 +61,7 @@ public:
   map<string, vector<string>> getTrendEqTags() const;
   map<string, vector<int>> getEqNums() const;
   map<string, vector<int>> getTrendEqNums() const;
+  vector<int> getTrendEqNums(const string &name_arg) const;
   vector<int> getEqNums(const string &name_arg) const;
   vector<int> getMaxLags(const string &name_arg) const;
   int getMaxLag(const string &name_arg) const;
@@ -68,6 +69,7 @@ public:
   vector<expr_t> getLhsExprT(const string &name_arg) const;
   vector<bool> getDiff(const string &name_arg) const;
   vector<int> getOrigDiffVar(const string &name_arg) const;
+  map<string, vector<int>> getNonTrendEqNums() const;
   vector<int> getNonTrendEqNums(const string &name_arg) const;
   vector<bool> getNonstationary(const string &name_arg) const;
 
@@ -82,6 +84,7 @@ public:
   void setNonstationary(map<string, vector<bool>> nonstationary_arg);
   void setTrendVar(map<string, vector<int>> trend_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 setNonTrendEqNums(map<string, vector<int>> trend_eqnums_arg);
 
   //! Write output of this class