diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 9d79ef1d7fe9efd3f25c7a57b26b5bfb71844daa..8a3eebcd397cbbc038278e5e4bd4fc2e75daf6c9 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3268,12 +3268,13 @@ DynamicModel::fillVarModelTable() const
 void
 DynamicModel::fillVarModelTableFromOrigModel() const
 {
-  map<string, vector<int>> lags, orig_diff_var;
+  map<string, vector<int>> lags;
+  map<string, vector<optional<int>>> orig_diff_var;
   map<string, vector<bool>> diff;
   for (const auto &[model_name, eqns] : var_model_table.getEqNums())
     {
       set<expr_t> lhs;
-      vector<int> orig_diff_var_vec;
+      vector<optional<int>> orig_diff_var_vec;
       vector<bool> diff_vec;
       for (auto eqn : eqns)
         {
@@ -3323,7 +3324,7 @@ DynamicModel::fillVarModelTableFromOrigModel() const
               orig_diff_var_vec.push_back(diff_set.begin()->first);
             }
           else
-            orig_diff_var_vec.push_back(-1);
+            orig_diff_var_vec.push_back(nullopt);
 
         }
 
@@ -3616,12 +3617,13 @@ DynamicModel::computeErrorComponentMatrices(const ExprNode::subst_table_t &diff_
 void
 DynamicModel::fillTrendComponentModelTableFromOrigModel() const
 {
-  map<string, vector<int>> lags, orig_diff_var;
+  map<string, vector<int>> lags;
+  map<string, vector<optional<int>>> orig_diff_var;
   map<string, vector<bool>> diff;
   for (const auto &[model_name, eqns] : trend_component_model_table.getEqNums())
     {
       set<expr_t> lhs;
-      vector<int> orig_diff_var_vec;
+      vector<optional<int>> orig_diff_var_vec;
       vector<bool> diff_vec;
       for (auto eqn : eqns)
         {
@@ -3666,7 +3668,7 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel() const
               orig_diff_var_vec.push_back(diff_set.begin()->first);
             }
           else
-            orig_diff_var_vec.push_back(-1);
+            orig_diff_var_vec.push_back(nullopt);
         }
 
       if (eqns.size() != lhs.size())
@@ -3710,7 +3712,6 @@ DynamicModel::getUndiffLHSForPac(const string &aux_model_name,
   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);
   vector<bool> diff = trend_component_model_table.getDiff(aux_model_name);
-  vector<int> orig_diff_var = trend_component_model_table.getOrigDiffVar(aux_model_name);
   vector<int> eqnumber = trend_component_model_table.getEqNums(aux_model_name);
   vector<int> nontrend_eqnums = trend_component_model_table.getNonTargetEqNums(aux_model_name);
 
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 0fd6f5ec4024292b5486a0356762ca4fe033775a..5d46d3d4ae3cce4169d1463f4c5be99d9014f771 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -101,7 +101,7 @@ TrendComponentModelTable::setDiff(map<string, vector<bool>> diff_arg)
 }
 
 void
-TrendComponentModelTable::setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg)
+TrendComponentModelTable::setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg)
 {
   orig_diff_var = move(orig_diff_var_arg);
 }
@@ -240,13 +240,6 @@ TrendComponentModelTable::getDiff(const string &name_arg) const
   return diff.find(name_arg)->second;
 }
 
-const vector<int> &
-TrendComponentModelTable::getOrigDiffVar(const string &name_arg) const
-{
-  checkModelName(name_arg);
-  return orig_diff_var.find(name_arg)->second;
-}
-
 void
 TrendComponentModelTable::writeOutput(const string &basename, ostream &output) const
 {
@@ -297,8 +290,8 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
         output << boolalpha << it << " ";
       output << "];" << endl
              << "M_.trend_component." << name << ".orig_diff_var = [";
-      for (auto it : orig_diff_var.at(name))
-        output << (it >= 0 ? symbol_table.getTypeSpecificID(it) + 1 : -1) << " ";
+      for (const auto &it : orig_diff_var.at(name))
+        output << (it ? symbol_table.getTypeSpecificID(*it) + 1 : -1) << " ";
       output << "];" << endl
              << "M_.trend_component." << name << ".nonstationary = [";
       for (size_t i = 0; i < diff.at(name).size(); i++)
@@ -483,8 +476,8 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
       output << "];" << endl
              << "M_.var." << name << ".nonstationary = M_.var." << name << ".diff;" << endl
              << "M_.var." << name << ".orig_diff_var = [";
-      for (auto it : orig_diff_var.at(name))
-        output << (it >= 0 ? symbol_table.getTypeSpecificID(it) + 1 : -1) << " ";
+      for (const auto &it : orig_diff_var.at(name))
+        output << (it ? symbol_table.getTypeSpecificID(*it) + 1 : -1) << " ";
       output << "];" << endl;
       int i = 1;
       for (const auto &it : rhs.at(name))
@@ -681,7 +674,7 @@ VarModelTable::setDiff(map<string, vector<bool>> diff_arg)
 }
 
 void
-VarModelTable::setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg)
+VarModelTable::setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg)
 {
   orig_diff_var = move(orig_diff_var_arg);
 }
diff --git a/src/SubModel.hh b/src/SubModel.hh
index c2d2fa6a91b00b0410eee0e24553d0f4116d6046..96cda36113d4e4ac7a5efce6d582021c489b0c40 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -46,7 +46,8 @@ private:
   SymbolTable &symbol_table;
   set<string> names;
   map<string, vector<string>> eqtags, target_eqtags;
-  map<string, vector<int>> eqnums, target_eqnums, nontarget_eqnums, max_lags, lhs, target_lhs, nontarget_lhs, orig_diff_var;
+  map<string, vector<int>> eqnums, target_eqnums, nontarget_eqnums, max_lags, lhs, target_lhs, nontarget_lhs;
+  map<string, vector<optional<int>>> orig_diff_var;
   map<string, vector<set<pair<int, int>>>> rhs;
   map<string, vector<bool>> diff;
   map<string, vector<expr_t>> lhs_expr_t;
@@ -77,7 +78,6 @@ public:
   const vector<int> &getLhs(const string &name_arg) const;
   const vector<expr_t> &getLhsExprT(const string &name_arg) const;
   const vector<bool> &getDiff(const string &name_arg) const;
-  const vector<int> &getOrigDiffVar(const string &name_arg) const;
   const map<string, vector<int>> &getNonTargetEqNums() const;
   const vector<int> &getNonTargetEqNums(const string &name_arg) const;
   const vector<int> &getNonTargetLhs(const string &name_arg) const;
@@ -89,7 +89,7 @@ public:
   void setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg);
   void setMaxLags(map<string, vector<int>> max_lags_arg);
   void setDiff(map<string, vector<bool>> diff_arg);
-  void setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg);
+  void setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg);
   void setTargetVar(map<string, vector<optional<int>>> target_vars_arg);
   void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
   void setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg,
@@ -125,7 +125,8 @@ private:
   set<string> names;
   map<string, bool> structural; // Whether VARs are structural or reduced-form
   map<string, vector<string>> eqtags;
-  map<string, vector<int>> eqnums, max_lags, lhs, lhs_orig_symb_ids, orig_diff_var;
+  map<string, vector<int>> eqnums, max_lags, lhs, lhs_orig_symb_ids;
+  map<string, vector<optional<int>>> orig_diff_var;
   map<string, vector<set<pair<int, int>>>> rhs; // name -> for each equation: set of pairs (var, lag)
   map<string, vector<bool>> diff;
   map<string, vector<expr_t>> lhs_expr_t;
@@ -164,7 +165,7 @@ public:
   void setLhsExprT(map<string, vector<expr_t>> lhs_expr_t_arg);
   void setDiff(map<string, vector<bool>> diff_arg);
   void setMaxLags(map<string, vector<int>> max_lags_arg);
-  void setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg);
+  void setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg);
   void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
   void setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg);
   void setConstants(map<string, map<int, expr_t>> constants_arg);