src/Bytecode.hh

@@ -20,6 +20,7 @@
 #ifndef BYTECODE_HH
 #define BYTECODE_HH
 
+#include <concepts>
 #include <filesystem>
 #include <fstream>
 #include <ios>
@@ -139,6 +140,9 @@ protected:
   ~Instruction() = default;
 };
 
+template<typename T>
+concept IsInstruction = derived_from<T, Instruction>;
+
 struct FLDZ final : public Instruction
 {
   FLDZ() : Instruction {Tag::FLDZ}
@@ -426,7 +430,7 @@ struct FSTPV final : public Instruction
 
 class FCALL final : public Instruction
 {
-  template<typename B>
+  template<IsInstruction B>
   friend Writer& operator<<(Writer& code_file, const B& instr);
 
 private:
@@ -602,7 +606,7 @@ public:
 
 class FBEGINBLOCK final : public Instruction
 {
-  template<typename B>
+  template<IsInstruction B>
   friend Writer& operator<<(Writer& code_file, const B& instr);
 
 private:
@@ -760,7 +764,7 @@ public:
 // Superclass of std::ofstream for writing a sequence of bytecode instructions
 class Writer : private ofstream
 {
-  template<typename B>
+  template<IsInstruction B>
   friend Writer& operator<<(Writer& code_file, const B& instr);
 
 private:
@@ -779,7 +783,7 @@ public:
      It is the responsibility of the caller to ensure that the new instruction
      occupies exactly as many bytes as the former one. */
   void
-  overwriteInstruction(int instruction_number, const auto& new_instruction)
+  overwriteInstruction(int instruction_number, const IsInstruction auto& new_instruction)
   {
     seekp(instructions_positions.at(instruction_number));
     *this << new_instruction;
@@ -790,7 +794,7 @@ public:
 
 // Overloads of operator<< for writing bytecode instructions
 
-template<typename B>
+template<IsInstruction B>
 Writer&
 operator<<(Writer& code_file, const B& instr)
 {

src/DynamicModel.cc

@@ -2084,6 +2084,203 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(
   pac_expectation_substitution[name] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction_term);
 }
 
+void
+DynamicModel::computePacModelConsistentExpectationSubstitutionWithComponents(
+    const string& name, int discount_symb_id, int pac_eq_max_lag,
+    ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+    vector<PacModelTable::target_component_t>& pac_target_components,
+    map<string, expr_t>& pac_expectation_substitution)
+{
+  auto create_aux_param = [&](const string& param_name) {
+    try
+      {
+        return symbol_table.addSymbol(param_name, SymbolType::parameter);
+      }
+    catch (SymbolTable::AlreadyDeclaredException)
+      {
+        cerr << "ERROR: the variable/parameter '" << param_name
+             << "' conflicts with some auxiliary parameter that will be generated for the '" << name
+             << "' PAC model. Please rename that parameter." << endl;
+        exit(EXIT_FAILURE);
+      }
+  };
+  int neqs = 0;
+
+  // At the end of this loop:
+  //         ₘ                          ₘ
+  // A_1 = 1+∑ αᵢ = A(1) and A_beta = 1+∑ αᵢβⁱ
+  //        ᵢ₌₁                        ᵢ₌₁
+  expr_t A_1 = One;
+  expr_t A_beta = One;
+  expr_t beta = AddVariable(discount_symb_id);
+  for (int i = 1; i <= pac_eq_max_lag + 1; i++)
+    {
+      string param_name = "mce_alpha_" + name + "_" + to_string(i);
+      try
+        {
+          int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+          pac_aux_param_symb_ids[name].push_back(alpha_i_symb_id);
+          A_1 = AddPlus(A_1, AddVariable(alpha_i_symb_id));
+          A_beta = AddPlus(A_beta, AddTimes(AddVariable(alpha_i_symb_id),
+                                            AddPower(beta, AddPossiblyNegativeConstant(i))));
+        }
+      catch (SymbolTable::AlreadyDeclaredException& e)
+        {
+          cerr << "The variable/parameter '" << param_name
+               << "' conflicts with a parameter that will be generated for the '" << name
+               << "' PAC model. Please rename it." << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
+
+  auto create_target_lag = [&](int variable, int lag) {
+    if (symbol_table.isAuxiliaryVariable(variable))
+      return AddVariable(symbol_table.getOrigSymbIdForAuxVar(variable), lag);
+    else
+      return AddVariable(variable, lag);
+  };
+
+  expr_t substexpr = Zero;
+  for (int component_idx {1};
+       auto& [component, growth, auxname, kind, coeff, growth_neutrality_param, h_indices,
+              original_growth, growth_info] : pac_target_components)
+    {
+      // Create the auxiliary variable for this component
+      int aux_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, component);
+
+      // Get the component variable id
+      int component_id = dynamic_cast<VariableNode*>(component)->symb_id;
+
+      //      ₘ
+      // fp = ∑ αᵢβⁱZₜ₊ᵢ
+      //     ᵢ₌₁
+      expr_t fp = Zero;
+      for (int i = 1; i <= pac_eq_max_lag + 1; i++)
+        {
+          int alpha_i_symb_id = -1;
+          string param_name = "mce_alpha_" + name + "_" + to_string(i);
+          try
+            {
+              alpha_i_symb_id = symbol_table.getID(param_name);
+            }
+          catch (SymbolTable::UnknownSymbolNameException& e)
+            {
+              alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+            }
+          fp = AddPlus(fp, AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
+                                             AddPower(beta, AddPossiblyNegativeConstant(i))),
+                                    AddVariable(aux_id, i)));
+        }
+      if (kind != PacTargetKind::ll) // non-stationary component y¹ₜ
+        {
+          // Add diff nodes and eqs for the non-stationary component
+          const VariableNode* aux_target_ns_var_diff_node;
+          expr_t diff_node_to_search = AddDiff(create_target_lag(component_id, 0));
+          if (auto sit = diff_subst_table.find(diff_node_to_search); sit != diff_subst_table.end())
+            aux_target_ns_var_diff_node = sit->second;
+          else
+            {
+              int symb_id
+                  = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
+              aux_target_ns_var_diff_node = AddVariable(symb_id);
+              auto neweq = AddEqual(const_cast<VariableNode*>(aux_target_ns_var_diff_node),
+                                    AddMinus(create_target_lag(component_id, 0),
+                                             create_target_lag(component_id, -1)));
+              addEquation(neweq, nullopt);
+              addAuxEquation(neweq);
+              neqs++;
+            }
+
+          map<int, VariableNode*> target_ns_aux_var_to_add;
+          const VariableNode* last_aux_var = aux_target_ns_var_diff_node;
+          for (int i = 1; i <= pac_eq_max_lag; i++, neqs++)
+            {
+              expr_t this_diff_node = AddDiff(create_target_lag(component_id, i));
+              int symb_id = symbol_table.addDiffLeadAuxiliaryVar(
+                  this_diff_node->idx, this_diff_node, last_aux_var->symb_id, 1);
+              VariableNode* current_aux_var = AddVariable(symb_id);
+              auto neweq = AddEqual(current_aux_var, AddVariable(last_aux_var->symb_id, 1));
+              addEquation(neweq, nullopt);
+              addAuxEquation(neweq);
+              last_aux_var = current_aux_var;
+              target_ns_aux_var_to_add[i] = current_aux_var;
+            }
+
+          // At the end of this loop,
+          //      ₘ₋₁ ⎛ ₘ₋₁     ⎞
+          //  fs = ∑  ⎢  ∑  αᵢβⁱ⎥ Δy¹ₜ₊ᵢ
+          //      ₖ₌₁ ⎝ ᵢ₌ₖ     ⎠
+          expr_t fs = Zero;
+          for (int k = 1; k <= pac_eq_max_lag; k++)
+            {
+              expr_t ssum = Zero;
+              for (int j = k + 1; j <= pac_eq_max_lag + 1; j++)
+                {
+                  int alpha_j_symb_id = -1;
+                  string param_name = "mce_alpha_" + name + "_" + to_string(j);
+                  try
+                    {
+                      alpha_j_symb_id = symbol_table.getID(param_name);
+                    }
+                  catch (SymbolTable::UnknownSymbolNameException& e)
+                    {
+                      alpha_j_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+                    }
+                  ssum = AddPlus(ssum, AddTimes(AddVariable(alpha_j_symb_id),
+                                                AddPower(beta, AddPossiblyNegativeConstant(j))));
+                }
+              fs = AddPlus(fs, AddTimes(ssum, target_ns_aux_var_to_add[k]));
+            }
+
+          // Assembling the equation Zₜ¹ = A_1 ( Δyₜ¹ - fs ) - fp_1
+          auto neweq_1 = AddEqual(
+              AddVariable(aux_id),
+              AddMinus(
+                  AddTimes(A_1,
+                           AddMinus(const_cast<VariableNode*>(aux_target_ns_var_diff_node), fs)),
+                  fp));
+          addEquation(neweq_1, nullopt);
+          neqs++;
+          /* This equation is not added to the list of auxiliary equations, because it
+             is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
+        }
+      else // Stationary component yₜ⁰
+        {
+          // Assembling the equation Zₜ⁰ = A_beta A_1 yₜ⁰ - fp
+          auto neweq_0 = AddEqual(AddVariable(aux_id),
+                                  AddMinus(AddTimes(A_beta, AddTimes(A_1, component)), fp));
+          addEquation(neweq_0, nullopt);
+          neqs++;
+          /* This equation is not added to the list of auxiliary equations, because it
+             is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
+        }
+
+      // If needed, add the growth neutrality correction for this component
+      expr_t growth_correction_term = Zero;
+      string name_component = name + "_component" + to_string(component_idx);
+      if (growth)
+        {
+          growth_neutrality_param
+              = create_aux_param(name_component + "_pac_growth_neutrality_correction");
+          growth_correction_term = AddTimes(growth, AddVariable(growth_neutrality_param));
+        }
+      else
+        growth_neutrality_param = -1;
+
+      substexpr = AddPlus(substexpr,
+                          AddTimes(coeff, AddPlus(AddVariable(aux_id), growth_correction_term)));
+      component_idx++;
+    }
+
+  cout << "PAC Model Consistent Expectation: added " << neqs
+       << " auxiliary variables and equations for model " << name << "." << endl;
+
+  /* The growth correction term is not added to the definition of Z₁
+     because the latter is recursive. Rather put it at the level of the
+     substition of pac_expectation operator. */
+  pac_expectation_substitution[name] = substexpr;
+}
+
 void
 DynamicModel::computePacBackwardExpectationSubstitution(
     const string& name, const vector<int>& lhs, int max_lag, const string& aux_model_type,

src/DynamicModel.hh

@@ -638,6 +638,18 @@ public:
       map<string, int>& pac_aux_var_symb_ids, map<string, vector<int>>& pac_aux_param_symb_ids,
       map<string, expr_t>& pac_expectation_substitution);
 
+  /* For a PAC MCE model with an associated pac_target_info, fill pac_expectation_substitution with
+     the expression that will be substituted for the pac_expectation operator. In the process, add
+     the variables and the equations defining Z₁ and Z₀ for each component. The new auxiliary
+     parameters are added to pac_mce_alpha_symb_ids. The routine also creates the auxiliary
+     variables for the components, and adds the corresponding equations.
+  */
+  void computePacModelConsistentExpectationSubstitutionWithComponents(
+      const string& name, int discount_symb_id, int pac_eq_max_lag,
+      ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+      vector<PacModelTable::target_component_t>& pac_target_components,
+      map<string, expr_t>& pac_expectation_substitution);
+
   /* For a PAC backward model, fill pac_expectation_substitution with the
      expression that will be substituted for the pac_expectation operator.
      The symbol IDs of the new parameters are also added to pac_aux_param_symb_ids.

src/DynareBison.yy

@@ -215,7 +215,7 @@ str_tolower(string s)
 %token ENDVAL_STEADY STEADY_SOLVE_ALGO STEADY_MAXIT STEADY_TOLF STEADY_TOLX STEADY_MARKOWITZ
 %token HOMOTOPY_MAX_COMPLETION_SHARE HOMOTOPY_MIN_STEP_SIZE HOMOTOPY_INITIAL_STEP_SIZE HOMOTOPY_STEP_SIZE_INCREASE_SUCCESS_COUNT
 %token HOMOTOPY_LINEARIZATION_FALLBACK HOMOTOPY_MARGINAL_LINEARIZATION_FALLBACK HOMOTOPY_EXCLUDE_VAREXO FROM_INITVAL_TO_ENDVAL
-%token STATIC_MFS RELATIVE_TO_INITVAL MATCHED_IRFS MATCHED_IRFS_WEIGHTS WEIGHTS
+%token STATIC_MFS RELATIVE_TO_INITVAL MATCHED_IRFS MATCHED_IRFS_WEIGHTS WEIGHTS PERPENDICULAR
 
 %token <vector<string>> SYMBOL_VEC
 
@@ -1011,6 +1011,14 @@ equation : hand_side EQUAL hand_side ';'
            { $$ = driver.add_model_equal($4, $6, $2); }
          | '[' tag_pair_list ']' hand_side ';'
            { $$ = driver.add_model_equal_with_zero_rhs($4, $2); }
+         | hand_side EQUAL hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal($1, $3, {}, $5); }
+         | hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal_with_zero_rhs($1, {}, $3); }
+         | '[' tag_pair_list ']' hand_side EQUAL hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal($4, $6, $2, $8); }
+         | '[' tag_pair_list ']' hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal_with_zero_rhs($4, $2, $6); }
          ;
 
 tag_pair_list : tag_pair_list COMMA tag_pair

src/DynareFlex.ll

@@ -1016,6 +1016,8 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT,DYNARE_BLOCK>nan {return token::NAN_CONSTANT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>inf {return token::INF_CONSTANT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>constants {return token::CONSTANTS;}
+<DYNARE_BLOCK>⟂ {return token::PERPENDICULAR;}
+<DYNARE_BLOCK>_\|_ {return token::PERPENDICULAR;}
 
  /* options for GSA module by Marco Ratto */
 <DYNARE_STATEMENT>identification {return token::IDENTIFICATION;}

src/ExprNode.hh

 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -20,6 +20,7 @@
 #ifndef EXPR_NODE_HH
 #define EXPR_NODE_HH
 
+#include <concepts>
 #include <functional>
 #include <map>
 #include <optional>
@@ -1173,6 +1174,7 @@ protected:
   // Returns the node obtained by applying a transformation recursively on the argument (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1339,6 +1341,7 @@ private:
   // Returns the node obtained by applying a transformation recursively on the arguments (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1532,6 +1535,7 @@ private:
   // Returns the node obtained by applying a transformation recursively on the arguments (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1660,6 +1664,7 @@ private:
   // Returns the node obtained by applying a transformation recursively on the arguments (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {

src/ParsingDriver.cc

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -2626,7 +2626,8 @@ ParsingDriver::extended_path()
 }
 
 expr_t
-ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_tags)
+ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_tags,
+                               expr_t complementarity_condition)
 {
   expr_t id = model_tree->AddEqual(arg1, arg2);
 
@@ -2634,6 +2635,28 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
     if (key == "endogenous")
       declare_or_change_type(SymbolType::endogenous, value);
 
+  if (eq_tags.contains("mcp"))
+    {
+      if (complementarity_condition)
+        error("Can't have both an 'mcp' tag and a complementarity condition after the "
+              "perpendicular symbol");
+      else
+        warning("Specifying complementarity conditions with the 'mcp' tag is obsolete. Please "
+                "consider switching to the new syntax using the perpendicular symbol.");
+    }
+
+  if (complementarity_condition)
+    {
+      if (auto bcomp = dynamic_cast<BinaryOpNode*>(complementarity_condition);
+          !(bcomp
+            && (bcomp->op_code == BinaryOpcode::less || bcomp->op_code == BinaryOpcode::lessEqual
+                || bcomp->op_code == BinaryOpcode::greater
+                || bcomp->op_code == BinaryOpcode::greaterEqual)))
+        error("The complementarity constraint must be an inequality.");
+
+      eq_tags.emplace("mcp", complementarity_condition->toString());
+    }
+
   if (eq_tags.contains("static"))
     {
       // If the equation is tagged [static]
@@ -2723,9 +2746,10 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
 }
 
 expr_t
-ParsingDriver::add_model_equal_with_zero_rhs(expr_t arg, map<string, string> eq_tags)
+ParsingDriver::add_model_equal_with_zero_rhs(expr_t arg, map<string, string> eq_tags,
+                                             expr_t complementarity_condition)
 {
-  return add_model_equal(arg, model_tree->Zero, move(eq_tags));
+  return add_model_equal(arg, model_tree->Zero, move(eq_tags), complementarity_condition);
 }
 
 void

src/ParsingDriver.hh

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -748,9 +748,11 @@ public:
   //! Extended path
   void extended_path();
   //! Writes token "arg1=arg2" to model tree
-  expr_t add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_tags);
+  expr_t add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_tags,
+                         expr_t complementarity_condition = nullptr);
   //! Writes token "arg=0" to model tree
-  expr_t add_model_equal_with_zero_rhs(expr_t arg, map<string, string> eq_tags);
+  expr_t add_model_equal_with_zero_rhs(expr_t arg, map<string, string> eq_tags,
+                                       expr_t complementarity_condition = nullptr);
   //! Writes token "arg1+arg2" to model tree
   expr_t add_plus(expr_t arg1, expr_t arg2);
   //! Writes token "arg1-arg2" to model tree

src/Statement.hh

 /*
- * Copyright © 2006-2023 Dynare Team
+ * Copyright © 2006-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -254,6 +254,12 @@ public:
   {
   };
 
+  // pair<string, string> corresponds to a pair of numerical values
+  // vector<vector<string>> corresponds to a vector of vectors of numerical values
+  using OptionValue
+      = variant<NumVal, pair<string, string>, StringVal, DateVal, SymbolListVal, vector<int>,
+                VecStrVal, VecCellStrVal, VecValueVal, vector<vector<string>>>;
+
   [[nodiscard]] bool empty() const;
   void clear();
   // Whether there is an option with that name that has been given a value
@@ -281,7 +287,8 @@ public:
   /* Retrieves the value of the option with that name.
      Throws UnknownOptionException if there is no option with that name.
      Throws bad_variant_access if the option has a value of a different type. */
-  template<class T>
+  template<typename T>
+    requires requires(T p) { std::get<T>(OptionValue {}); }
   T
   get(const string& name) const
   {
@@ -295,7 +302,8 @@ public:
   /* Retrieves the value of the option with that name.
      Returns nullopt if there is no option with that name.
      Throws bad_variant_access if the option has a value of a different type. */
-  template<class T>
+  template<typename T>
+    requires requires(T p) { std::get<T>(OptionValue {}); }
   optional<T>
   get_if(const string& name) const
   {
@@ -309,7 +317,7 @@ public:
   /* Applies a variant visitor to the value of the option with that name.
      Throws UnknownOptionException if there is no option with that name. */
   template<class Visitor>
-  decltype(auto)
+    requires invocable<Visitor, OptionValue> decltype(auto)
   visit(const string& name, Visitor&& vis) const
   {
     auto it = options.find(name);
@@ -324,11 +332,7 @@ public:
   void writeJsonOutput(ostream& output) const;
 
 private:
-  // pair<string, string> corresponds to a pair of numerical values
-  // vector<vector<string>> corresponds to a vector of vectors of numerical values
-  map<string, variant<NumVal, pair<string, string>, StringVal, DateVal, SymbolListVal, vector<int>,
-                      VecStrVal, VecCellStrVal, VecValueVal, vector<vector<string>>>>
-      options;
+  map<string, OptionValue> options;
   void writeOutputCommon(ostream& output, const string& option_group) const;
   // Helper constant for visitors
   template<class>

src/SubModel.cc

@@ -1384,11 +1384,11 @@ PacModelTable::transformPass(const lag_equivalence_table_t& unary_ops_nodes,
         {
           if (target_info.contains(name))
             {
-              cerr << "ERROR: the block 'pac_target_info(" << name
-                   << ")' is not supported in the context of a PAC model with model-consistent "
-                      "expectations (MCE)."
-                   << endl;
-              exit(EXIT_FAILURE);
+              assert(growth_correction_term == dynamic_model.Zero);
+              dynamic_model.computePacModelConsistentExpectationSubstitutionWithComponents(
+                  name, symbol_table.getID(discount[name]), pacEquationMaxLag(name),
+                  diff_subst_table, aux_param_symb_ids, get<2>(target_info[name]),
+                  pac_expectation_substitution);
             }
           else
             dynamic_model.computePacModelConsistentExpectationSubstitution(
@@ -1507,7 +1507,7 @@ PacModelTable::writeOutput(ostream& output) const
 
   // Write the auxiliary variable IDs created for the pac_expectation operator
   for (auto& [name, id] : aux_var_symb_ids)
-    output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.z1" : "aux_id")
+    output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.z" : "aux_id")
            << " = " << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
 
   // Write PAC equation name info
@@ -1709,6 +1709,21 @@ PacModelTable::writeOutput(ostream& output) const
           }
         component_idx++;
       }
+
+  for (auto& [name, val] : target_info)
+    if (aux_model_name.at(name).empty())
+      {
+        string pac_model_name = "M_.pac." + name + ".mce.";
+        output << pac_model_name << "z = NaN(" << get<2>(val).size() << ",1);" << endl;
+        for (int component_idx {1};
+             auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                    h_indices, original_growth_component, growth_component_info] : get<2>(val))
+          {
+            output << pac_model_name << "z(" << component_idx
+                   << ") = " << symbol_table.getTypeSpecificID(auxname) + 1 << ";" << endl;
+            component_idx++;
+          }
+      }
 }
 
 void

src/SymbolTable.cc

@@ -731,7 +731,7 @@ SymbolTable::getOrigSymbIdForAuxVar(int aux_var_symb_id_arg) const noexcept(fals
     if ((aux_var.type == AuxVarType::endoLag || aux_var.type == AuxVarType::exoLag
          || aux_var.type == AuxVarType::diff || aux_var.type == AuxVarType::diffLag
          || aux_var.type == AuxVarType::diffLead || aux_var.type == AuxVarType::diffForward
-         || aux_var.type == AuxVarType::unaryOp)
+         || aux_var.type == AuxVarType::unaryOp || aux_var.type == AuxVarType::pacExpectation)
         && aux_var.symb_id == aux_var_symb_id_arg)
       {
         if (optional<int> r = aux_var.orig_symb_id; r)