diff --git a/src/CodeInterpreter.hh b/src/CodeInterpreter.hh
index 899fb0142adbb4d9261812905e2c544ae80ccc1b..efdb0367b4deb99645c00d13f7384b4009384b7e 100644
--- a/src/CodeInterpreter.hh
+++ b/src/CodeInterpreter.hh
@@ -241,6 +241,14 @@ enum class PriorDistributions
    weibull = 8
   };
 
+enum class PacTargetKind
+  {
+    unspecified, // Must be the first one, because it’s the default initializer
+    ll,
+    dl,
+    dd
+  };
+
 struct Block_contain_type
 {
   int Equation, Variable, Own_Derivative;
diff --git a/src/DataTree.cc b/src/DataTree.cc
index d460ad08731891439d611dec81c03327800499ce..654ff9d5950112b5265fec7d4d470957aadc9b13 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -93,6 +93,7 @@ DataTree::operator=(const DataTree &d)
   external_function_node_map.clear();
   var_expectation_node_map.clear();
   pac_expectation_node_map.clear();
+  pac_target_nonstationary_node_map.clear();
   first_deriv_external_function_node_map.clear();
   second_deriv_external_function_node_map.clear();
 
@@ -704,6 +705,20 @@ DataTree::AddPacExpectation(const string &model_name)
   return p;
 }
 
+expr_t
+DataTree::AddPacTargetNonstationary(const string &model_name)
+{
+  if (auto it = pac_target_nonstationary_node_map.find(model_name);
+      it != pac_target_nonstationary_node_map.end())
+    return it->second;
+
+  auto sp = make_unique<PacTargetNonstationaryNode>(*this, node_list.size(), model_name);
+  auto p = sp.get();
+  node_list.push_back(move(sp));
+  pac_target_nonstationary_node_map[model_name] = p;
+  return p;
+}
+
 BinaryOpNode *
 DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
 {
diff --git a/src/DataTree.hh b/src/DataTree.hh
index e7baa641bbccfb765afad51a8e7643d2558ed82e..565a0189fd31cce1d6bea632562c9e7a68e1acb3 100644
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -81,6 +81,10 @@ private:
   using pac_expectation_node_map_t = map<string, PacExpectationNode *>;
   pac_expectation_node_map_t pac_expectation_node_map;
 
+  // model_name -> PacTargetNonstationaryNode
+  using pac_target_nonstationary_node_map_t = map<string, PacTargetNonstationaryNode *>;
+  pac_target_nonstationary_node_map_t pac_target_nonstationary_node_map;
+
   // (arguments, deriv_idx, symb_id) -> FirstDerivExternalFunctionNode
   using first_deriv_external_function_node_map_t = map<tuple<vector<expr_t>, int, int>, FirstDerivExternalFunctionNode *>;
   first_deriv_external_function_node_map_t first_deriv_external_function_node_map;
@@ -101,13 +105,6 @@ protected:
   //! Internal implementation of ParamUsedWithLeadLag()
   bool ParamUsedWithLeadLagInternal() const;
 
-  /*! Takes a MATLAB/Octave package name (possibly with several levels nested using dots),
-    and returns the name of the corresponding filesystem directory (which
-    is created by the function if it does not exist).
-    In practice the package nesting is used for the planner_objective (stored
-    inside +objective subdir). */
-  static string packageDir(const string &package);
-
   /* Writes the contents of “new_contents” to the file “filename”. However, if
      the file already exists and would not be modified by this operation, then do
      nothing. */
@@ -260,6 +257,8 @@ public:
   expr_t AddVarExpectation(const string &model_name);
   //! Adds pac_expectation command to model tree
   expr_t AddPacExpectation(const string &model_name);
+  //! Adds a pac_target_nonstationary node to model tree
+  expr_t AddPacTargetNonstationary(const string &model_name);
   //! Adds a model local variable with its value
   void AddLocalVariable(int symb_id, expr_t value) noexcept(false);
   //! Adds an external function node
@@ -344,6 +343,13 @@ public:
   {
     no_commutativity = true;
   }
+
+  /*! Takes a MATLAB/Octave package name (possibly with several levels nested using dots),
+    and returns the name of the corresponding filesystem directory (which
+    is created by the function if it does not exist).
+    In practice the package nesting is used for the planner_objective (stored
+    inside +objective subdir). */
+  static string packageDir(const string &package);
 };
 
 inline expr_t
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 39a6ce520ac9c57c3b25e242f8b1791364b1853a..6a4885da4638a2680562d7dce0f798682690b012 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4080,6 +4080,80 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
   pac_expectation_substitution[name] = AddVariable(expect_var_id);
 }
 
+void
+DynamicModel::computePacBackwardExpectationSubstitutionWithComponents(const string &name,
+                                                                      const vector<int> &lhs,
+                                                                      int max_lag,
+                                                                      const string &aux_model_type,
+                                                                      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);
+      }
+  };
+
+  expr_t substexpr = Zero;
+  int component_idx = 1;
+
+  for (auto &[component, growth, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth, growth_info] : pac_target_components)
+    {
+      string name_component = name + "_component" + to_string(component_idx);
+
+      // Create the linear combination of the variables from the auxiliary model
+      expr_t auxdef = Zero;
+      if (aux_model_type == "var")
+        {
+          /* If the auxiliary model is a VAR, add a parameter corresponding
+             to the constant. */
+          int new_param_symb_id = create_aux_param("h_" + name_component + "_constant");
+          h_indices.push_back(new_param_symb_id);
+          auxdef = AddPlus(auxdef, AddVariable(new_param_symb_id));
+        }
+      for (int i = 1; i < max_lag + 1; i++)
+        for (auto lhsit : lhs)
+          {
+            int new_param_symb_id = create_aux_param("h_" + name_component + "_var_" + symbol_table.getName(lhsit) + "_lag_" + to_string(i));
+            h_indices.push_back(new_param_symb_id);
+            auxdef = AddPlus(auxdef, AddTimes(AddVariable(new_param_symb_id),
+                                              AddVariable(lhsit, -i)));
+          }
+
+      // If needed, add the growth neutrality correction for this component
+      if (growth)
+        {
+          growth_neutrality_param = create_aux_param(name_component + "_pac_growth_neutrality_correction");
+          auxdef = AddPlus(auxdef, AddTimes(growth, AddVariable(growth_neutrality_param)));
+        }
+      else
+        growth_neutrality_param = -1;
+
+      // Create the auxiliary variable for this component
+      int aux_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, auxdef);
+      expr_t auxvar = AddVariable(aux_id);
+
+      // Add the equation defining the auxiliary variable for this component
+      expr_t neweq = AddEqual(auxvar, auxdef);
+      addEquation(neweq, -1);
+      addAuxEquation(neweq);
+
+      // Update the expression to be substituted for the pac_expectation operator
+      substexpr = AddPlus(substexpr, AddTimes(coeff, auxvar));
+
+      component_idx++;
+    }
+
+  pac_expectation_substitution[name] = substexpr;
+}
+
 void
 DynamicModel::substitutePacExpectation(const map<string, expr_t> &pac_expectation_substitution,
                                        const map<string, string> &pac_eq_name)
@@ -4093,6 +4167,13 @@ DynamicModel::substitutePacExpectation(const map<string, expr_t> &pac_expectatio
     }
 }
 
+void
+DynamicModel::substitutePacTargetNonstationary(const string &pac_model_name, expr_t substexpr)
+{
+  for (auto &eq : equations)
+    eq = dynamic_cast<BinaryOpNode *>(eq->substitutePacTargetNonstationary(pac_model_name, substexpr));
+}
+
 void
 DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
@@ -6482,3 +6563,15 @@ DynamicModel::simplifyEquations()
       findConstantEquationsWithoutMcpTag(subst_table);
     }
 }
+
+void
+DynamicModel::checkNoRemainingPacTargetNonstationary() const
+{
+  for (size_t eq = 0; eq < equations.size(); eq++)
+    if (equations[eq]->containsPacTargetNonstationary())
+      {
+        cerr << "ERROR: in equation " << equation_tags.getTagValueByEqnAndKey(eq, "name")
+             << ", the pac_target_nonstationary operator does not match a corresponding 'pac_target_info' block" << endl;
+        exit(EXIT_FAILURE);
+      }
+}
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 05584e4270bb6d81ffb5480277606df1f5f72a43..6907bcda9104bfa01226e12ea02448807ca55080 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -582,10 +582,28 @@ public:
                                                  map<string, vector<int>> &pac_aux_param_symb_ids,
                                                  map<string, expr_t> &pac_expectation_substitution);
 
+  /* Same as above, but for PAC models which have an associated
+     pac_target_info.
+     Contrary to the above routine, this one will create the growth correction
+     parameters as needed.
+     Those parameter IDs, as well as the IDs for the h parameters, are stored
+     in target_components.
+     The routine also creates the auxiliary variables for the components, and
+     adds the corresponding equations. */
+  void computePacBackwardExpectationSubstitutionWithComponents(const string &name,
+                                                               const vector<int> &lhs,
+                                                               int max_lag,
+                                                               const string &aux_model_type,
+                                                               vector<PacModelTable::target_component_t> &pac_target_components,
+                                                               map<string, expr_t> &pac_expectation_substitution);
+
   //! Substitutes pac_expectation operator with expectation based on auxiliary model
   void substitutePacExpectation(const map<string, expr_t> &pac_expectation_substitution,
                                 const map<string, string> &pac_eq_name);
 
+  //! Substitutes the pac_target_nonstationary operator of a given pac_model
+  void substitutePacTargetNonstationary(const string &pac_model_name, expr_t substexpr);
+
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
                                  const ExprNode::subst_table_t &diff_subst_table) const;
@@ -605,6 +623,10 @@ public:
     out otherwise */
   void checkNoRemainingPacExpectation() const;
 
+  /*! Checks that all pac_target_nonstationary operators have been substituted, error
+    out otherwise */
+  void checkNoRemainingPacTargetNonstationary() const;
+
   auto
   getStaticOnlyEquationsInfo() const
   {
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 1d2dde48801292428cfcc53848fd0ed6e0a90cc9..87c0f9ba428b9867215ef031cc66389e3b54683f 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -174,6 +174,8 @@ class ParsingDriver;
 %token RANDOM_FUNCTION_CONVERGENCE_CRITERION RANDOM_PARAMETER_CONVERGENCE_CRITERION NO_INIT_ESTIMATION_CHECK_FIRST_OBS
 %token HETEROSKEDASTIC_FILTER TIME_SHIFT STRUCTURAL TERMINAL_STEADY_STATE_AS_GUESS_VALUE
 %token SURPRISE OCCBIN_CONSTRAINTS
+%token PAC_TARGET_INFO COMPONENT TARGET AUXNAME AUXNAME_TARGET_NONSTATIONARY PAC_TARGET_NONSTATIONARY
+%token <string> KIND LL DL DD
 /* Method of Moments */
 %token METHOD_OF_MOMENTS MOM_METHOD
 %token BARTLETT_KERNEL_LAG WEIGHTING_MATRIX WEIGHTING_MATRIX_SCALING_FACTOR ANALYTIC_STANDARD_ERRORS ANALYTIC_JACOBIAN PENALIZED_ESTIMATOR VERBOSE 
@@ -214,6 +216,7 @@ class ParsingDriver;
 %type <vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>>> occbin_constraints_regimes_list
 %type <map<string, expr_t>> occbin_constraints_regime_options_list
 %type <pair<string, expr_t>> occbin_constraints_regime_option
+%type <PacTargetKind> pac_target_kind
 %%
 
 %start statement_list;
@@ -348,6 +351,7 @@ statement : parameters
           | model_replace
           | model_options
           | var_remove
+          | pac_target_info
           ;
 
 dsample : DSAMPLE INT_NUMBER ';'
@@ -945,6 +949,49 @@ occbin_constraints_regime_option : BIND hand_side ';'
                                    { $$ = { "error_relax", $2 }; }
                                  ;
 
+pac_target_info : PAC_TARGET_INFO '(' symbol ')' ';'
+                  { driver.begin_pac_target_info($3); }
+                  pac_target_info_statement_list
+                  END ';'
+                  { driver.end_pac_target_info(); }
+                  ;
+
+pac_target_info_statement_list : pac_target_info_statement
+                               | pac_target_info_statement_list pac_target_info_statement
+                               ;
+
+pac_target_info_statement : TARGET hand_side ';'
+                            { driver.set_pac_target_info_target($2); }
+                          | AUXNAME_TARGET_NONSTATIONARY symbol ';'
+                            { driver.set_pac_target_info_auxname_target_nonstationary($2); }
+                          | pac_target_info_component
+                          ;
+
+pac_target_info_component : COMPONENT hand_side ';'
+                            pac_target_info_component_list
+                            { driver.add_pac_target_info_component($2); }
+                          ;
+
+pac_target_info_component_list : pac_target_info_component_elem
+                               | pac_target_info_component_list pac_target_info_component_elem
+                               ;
+
+pac_target_info_component_elem : GROWTH hand_side ';'
+                                 { driver.set_pac_target_info_component_growth($2); }
+                               | AUXNAME symbol ';'
+                                 { driver.set_pac_target_info_component_auxname($2); }
+                               | KIND pac_target_kind ';'
+                                 { driver.set_pac_target_info_component_kind($2); }
+                               ;
+
+pac_target_kind : LL
+                  { $$ = PacTargetKind::ll; }
+                | DL
+                  { $$ = PacTargetKind::dl; }
+                | DD
+                  { $$ = PacTargetKind::dd; }
+                ;
+
 /* The tokens below must be accepted in both DYNARE_STATEMENT and DYNARE_BLOCK
    states in the lexer, because of model block and model_options statement */
 model_option : BLOCK { driver.block(); }
@@ -1036,6 +1083,8 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_var_expectation($3); }
           | PAC_EXPECTATION '(' symbol ')'
             { $$ = driver.add_pac_expectation($3); }
+          | PAC_TARGET_NONSTATIONARY '(' symbol ')'
+            { $$ = driver.add_pac_target_nonstationary($3); }
           | MINUS hand_side %prec UNARY
             { $$ = driver.add_uminus($2); }
           | PLUS hand_side %prec UNARY
@@ -4301,6 +4350,10 @@ symbol : NAME
        | RELAX
        | ERROR_BIND
        | ERROR_RELAX
+       | KIND
+       | LL
+       | DL
+       | DD
        ;
 
 %%
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 1bc54b01d05f013aaba48f655759d29f581b4f79..f9be86f995da44c407fbd87d02ed90c0d78caa95 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -236,6 +236,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>matched_moments {BEGIN DYNARE_BLOCK; return token::MATCHED_MOMENTS;}
 <INITIAL>occbin_constraints {BEGIN DYNARE_BLOCK; return token::OCCBIN_CONSTRAINTS;}
 <INITIAL>model_replace {BEGIN DYNARE_BLOCK; return token::MODEL_REPLACE;}
+<INITIAL>pac_target_info {BEGIN DYNARE_BLOCK; return token::PAC_TARGET_INFO;}
 
  /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
@@ -684,7 +685,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>epilogue {return token::EPILOGUE;}
 <DYNARE_STATEMENT>growth_factor {return token::GROWTH_FACTOR;}
 <DYNARE_STATEMENT>log_growth_factor {return token::LOG_GROWTH_FACTOR;}
-<DYNARE_STATEMENT>growth {return token::GROWTH;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>growth {return token::GROWTH;}
 <DYNARE_STATEMENT>cova_compute {return token::COVA_COMPUTE;}
 <DYNARE_STATEMENT>discretionary_tol {return token::DISCRETIONARY_TOL;}
 <DYNARE_STATEMENT>analytic_derivation {return token::ANALYTIC_DERIVATION;}
@@ -800,6 +801,27 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK># {return Dynare::parser::token_type (yytext[0]);}
 
 <DYNARE_BLOCK>restriction {return token::RESTRICTION;}
+<DYNARE_BLOCK>component {return token::COMPONENT;}
+<DYNARE_BLOCK>target {return token::TARGET;}
+<DYNARE_BLOCK>auxname {return token::AUXNAME;}
+<DYNARE_BLOCK>auxname_target_nonstationary {return token::AUXNAME_TARGET_NONSTATIONARY;}
+<DYNARE_BLOCK>kind {
+  yylval->build<string>(yytext);
+  return token::KIND;
+}
+<DYNARE_BLOCK>ll {
+  yylval->build<string>(yytext);
+  return token::LL;
+}
+<DYNARE_BLOCK>dl {
+  yylval->build<string>(yytext);
+  return token::DL;
+}
+<DYNARE_BLOCK>dd {
+  yylval->build<string>(yytext);
+  return token::DD;
+}
+
 
  /* Inside Dynare statement */
 <DYNARE_STATEMENT>solve_algo {return token::SOLVE_ALGO;}
@@ -945,6 +967,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
 <DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
+<DYNARE_BLOCK>pac_target_nonstationary {return token::PAC_TARGET_NONSTATIONARY;}
 <DYNARE_STATEMENT>discount {return token::DISCOUNT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 6ce01c80ceaa9481cf82399521135f3c0f7125e7..cc7b73fb0888004be8c85aeda5ff941893407168 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -667,6 +667,12 @@ NumConstNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return const_cast<NumConstNode *>(this);
 }
 
+expr_t
+NumConstNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  return const_cast<NumConstNode *>(this);
+}
+
 expr_t
 NumConstNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -694,6 +700,12 @@ NumConstNode::containsPacExpectation(const string &pac_model_name) const
   return false;
 }
 
+bool
+NumConstNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  return false;
+}
+
 expr_t
 NumConstNode::replaceTrendVar() const
 {
@@ -1588,6 +1600,15 @@ VariableNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return const_cast<VariableNode *>(this);
 }
 
+expr_t
+VariableNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->substitutePacTargetNonstationary(name, subexpr);
+
+  return const_cast<VariableNode *>(this);
+}
+
 expr_t
 VariableNode::decreaseLeadsLags(int n) const
 {
@@ -1820,6 +1841,15 @@ VariableNode::containsPacExpectation(const string &pac_model_name) const
   return false;
 }
 
+bool
+VariableNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->containsPacTargetNonstationary(pac_model_name);
+
+  return false;
+}
+
 expr_t
 VariableNode::replaceTrendVar() const
 {
@@ -3538,6 +3568,13 @@ UnaryOpNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
+expr_t
+UnaryOpNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  expr_t argsubst = arg->substitutePacTargetNonstationary(name, subexpr);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 expr_t
 UnaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -3665,6 +3702,12 @@ UnaryOpNode::containsPacExpectation(const string &pac_model_name) const
   return arg->containsPacExpectation(pac_model_name);
 }
 
+bool
+UnaryOpNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  return arg->containsPacTargetNonstationary(pac_model_name);
+}
+
 expr_t
 UnaryOpNode::replaceTrendVar() const
 {
@@ -5086,6 +5129,14 @@ BinaryOpNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+expr_t
+BinaryOpNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  expr_t arg1subst = arg1->substitutePacTargetNonstationary(name, subexpr);
+  expr_t arg2subst = arg2->substitutePacTargetNonstationary(name, subexpr);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
 expr_t
 BinaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -5120,6 +5171,12 @@ BinaryOpNode::containsPacExpectation(const string &pac_model_name) const
   return arg1->containsPacExpectation(pac_model_name) || arg2->containsPacExpectation(pac_model_name);
 }
 
+bool
+BinaryOpNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  return arg1->containsPacTargetNonstationary(pac_model_name) || arg2->containsPacTargetNonstationary(pac_model_name);
+}
+
 expr_t
 BinaryOpNode::replaceTrendVar() const
 {
@@ -6341,6 +6398,15 @@ TrinaryOpNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+expr_t
+TrinaryOpNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  expr_t arg1subst = arg1->substitutePacTargetNonstationary(name, subexpr);
+  expr_t arg2subst = arg2->substitutePacTargetNonstationary(name, subexpr);
+  expr_t arg3subst = arg3->substitutePacTargetNonstationary(name, subexpr);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 expr_t
 TrinaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -6368,6 +6434,14 @@ TrinaryOpNode::containsPacExpectation(const string &pac_model_name) const
   return (arg1->containsPacExpectation(pac_model_name) || arg2->containsPacExpectation(pac_model_name) || arg3->containsPacExpectation(pac_model_name));
 }
 
+bool
+TrinaryOpNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  return arg1->containsPacTargetNonstationary(pac_model_name)
+    || arg2->containsPacTargetNonstationary(pac_model_name)
+    || arg3->containsPacTargetNonstationary(pac_model_name);
+}
+
 expr_t
 TrinaryOpNode::replaceTrendVar() const
 {
@@ -6741,6 +6815,15 @@ AbstractExternalFunctionNode::substitutePacExpectation(const string &name, expr_
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+expr_t
+AbstractExternalFunctionNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  vector<expr_t> arguments_subst;
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substitutePacTargetNonstationary(name, subexpr));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 expr_t
 AbstractExternalFunctionNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -6835,6 +6918,15 @@ AbstractExternalFunctionNode::containsPacExpectation(const string &pac_model_nam
   return false;
 }
 
+bool
+AbstractExternalFunctionNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  for (auto argument : arguments)
+    if (argument->containsPacTargetNonstationary(pac_model_name))
+      return true;
+  return false;
+}
+
 expr_t
 AbstractExternalFunctionNode::replaceTrendVar() const
 {
@@ -8344,12 +8436,24 @@ VarExpectationNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return const_cast<VarExpectationNode *>(this);
 }
 
+expr_t
+VarExpectationNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
 bool
 VarExpectationNode::containsPacExpectation(const string &pac_model_name) const
 {
   return false;
 }
 
+bool
+VarExpectationNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  return false;
+}
+
 void
 VarExpectationNode::writeJsonAST(ostream &output) const
 {
@@ -8417,6 +8521,12 @@ PacExpectationNode::containsPacExpectation(const string &pac_model_name) const
     return pac_model_name == model_name;
 }
 
+bool
+PacExpectationNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  return false;
+}
+
 void
 PacExpectationNode::writeJsonAST(ostream &output) const
 {
@@ -8443,6 +8553,99 @@ PacExpectationNode::substitutePacExpectation(const string &name, expr_t subexpr)
   return subexpr;
 }
 
+expr_t
+PacExpectationNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+PacTargetNonstationaryNode::PacTargetNonstationaryNode(DataTree &datatree_arg,
+                                                       int idx_arg,
+                                                       string model_name_arg) :
+  SubModelNode{datatree_arg, idx_arg, move(model_name_arg)}
+{
+}
+
+expr_t
+PacTargetNonstationaryNode::clone(DataTree &datatree) const
+{
+  return datatree.AddPacTargetNonstationary(model_name);
+}
+
+void
+PacTargetNonstationaryNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                                        const temporary_terms_t &temporary_terms,
+                                        const temporary_terms_idxs_t &temporary_terms_idxs,
+                                        const deriv_node_temp_terms_t &tef_terms) const
+{
+  assert(output_type != ExprNodeOutputType::matlabOutsideModel);
+  if (isLatexOutput(output_type))
+    {
+      output << "PAC_TARGET_NONSTATIONARY" << LEFT_PAR(output_type) << model_name << RIGHT_PAR(output_type);
+      return;
+    }
+}
+
+int
+PacTargetNonstationaryNode::maxLagWithDiffsExpanded() const
+{
+  // This node will be replaced by the target lagged by one
+  return 1;
+}
+
+expr_t
+PacTargetNonstationaryNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  return const_cast<PacTargetNonstationaryNode *>(this);
+}
+
+bool
+PacTargetNonstationaryNode::containsPacExpectation(const string &pac_model_name) const
+{
+  return false;
+}
+
+bool
+PacTargetNonstationaryNode::containsPacTargetNonstationary(const string &pac_model_name) const
+{
+  if (pac_model_name.empty())
+    return true;
+  else
+    return pac_model_name == model_name;
+}
+
+void
+PacTargetNonstationaryNode::writeJsonAST(ostream &output) const
+{
+  output << R"({"node_type" : "PacTargetNonstationaryNode", )"
+         << R"("name" : ")" << model_name << R"("})";
+}
+
+void
+PacTargetNonstationaryNode::writeJsonOutput(ostream &output,
+                                            const temporary_terms_t &temporary_terms,
+                                            const deriv_node_temp_terms_t &tef_terms,
+                                            bool isdynamic) const
+{
+  output << "pac_target_nonstationary("
+         << "model_name = " << model_name
+         << ")";
+}
+
+expr_t
+PacTargetNonstationaryNode::substitutePacExpectation(const string &name, expr_t subexpr)
+{
+  return const_cast<PacTargetNonstationaryNode *>(this);
+}
+
+expr_t
+PacTargetNonstationaryNode::substitutePacTargetNonstationary(const string &name, expr_t subexpr)
+{
+  if (model_name != name)
+    return const_cast<PacTargetNonstationaryNode *>(this);
+  return subexpr;
+}
+
 void
 ExprNode::decomposeAdditiveTerms(vector<pair<expr_t, int>> &terms, int current_sign) const
 {
@@ -8649,6 +8852,8 @@ ExprNode::matchParamTimesTargetMinusVariable(int symb_id) const
       if (datatree.symbol_table.isAuxiliaryVariable(target->symb_id))
         {
           auto avi = datatree.symbol_table.getAuxVarInfo(target->symb_id);
+          if (avi.get_type() == AuxVarType::pacTargetNonstationary && target->lag == -1)
+            return true;
           return (avi.get_type() == AuxVarType::unaryOp
                   && avi.get_unary_op() == "log"
                   && avi.get_orig_symb_id() != -1
@@ -8664,3 +8869,59 @@ ExprNode::matchParamTimesTargetMinusVariable(int symb_id) const
   else
     throw MatchFailureException{"Neither factor is of the form (target-variable) where target is endo or exo (possibly logged), and has one lag"};
 }
+
+pair<int, expr_t>
+ExprNode::matchEndogenousTimesConstant() const
+{
+  throw MatchFailureException{"This expression is not of the form endogenous*constant"};
+}
+
+pair<int, expr_t>
+VariableNode::matchEndogenousTimesConstant() const
+{
+  if (get_type() == SymbolType::endogenous)
+    return { symb_id, datatree.One };
+  else
+    throw MatchFailureException{"This expression is not of the form endogenous*constant"};
+}
+
+pair<int, expr_t>
+BinaryOpNode::matchEndogenousTimesConstant() const
+{
+  if (op_code == BinaryOpcode::times)
+    {
+      if (auto varg1 = dynamic_cast<VariableNode *>(arg1);
+          varg1 && varg1->get_type() == SymbolType::endogenous && arg2->isConstant())
+        return { varg1->symb_id, arg2 };
+      if (auto varg2 = dynamic_cast<VariableNode *>(arg2);
+          varg2 && varg2->get_type() == SymbolType::endogenous && arg1->isConstant())
+        return { varg2->symb_id, arg1 };
+    }
+  throw MatchFailureException{"This expression is not of the form endogenous*constant"};
+}
+
+pair<vector<pair<int, expr_t>>, expr_t>
+ExprNode::matchLinearCombinationOfEndogenousWithConstant() const
+{
+  vector<pair<expr_t, int>> all_terms;
+  decomposeAdditiveTerms(all_terms);
+
+  vector<pair<int, expr_t>> endo_terms;
+  expr_t intercept = datatree.Zero;
+  for (auto [term, sign] : all_terms)
+    if (term->isConstant())
+      {
+        if (sign == -1)
+          intercept = datatree.AddMinus(intercept, term);
+        else
+          intercept = datatree.AddPlus(intercept, term);
+      }
+    else
+      {
+        auto [endo_id, constant] = term->matchEndogenousTimesConstant();
+        if (sign == -1)
+          constant = datatree.AddUMinus(constant);
+        endo_terms.emplace_back(endo_id, constant);
+      }
+  return { endo_terms, intercept };
+}
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 17bf0405dd333105956e19cd3a0ed9a475d21d9d..c4b20990cb2c1088c073dd0273c9ba6908815ade 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -612,6 +612,9 @@ public:
   //! Substitute pac_expectation operator
   virtual expr_t substitutePacExpectation(const string &name, expr_t subexpr) = 0;
 
+  //! Substitute pac_target_nonstationary operator
+  virtual expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) = 0;
+
   virtual int findTargetVariable(int lhs_symb_id) const = 0;
 
   //! Add ExprNodes to the provided datatree
@@ -626,7 +629,7 @@ public:
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const = 0;
 
-  //! Matches a linear combination of variables, where scalars can be constant*parameter
+  //! Matches a linear combination of variables (endo or exo), where scalars can be constant*parameter
   /*! Returns a list of (variable_id, lag, param_id, constant)
     corresponding to the terms in the expression. When there is no
     parameter in a term, param_id == -1.
@@ -638,6 +641,15 @@ public:
   */
   vector<tuple<int, int, int, double>> matchLinearCombinationOfVariables(bool variable_obligatory_in_each_term = true) const;
 
+  /* Matches a linear combination of endogenous, where scalars can be any
+     constant expression (i.e. containing no endogenous, no exogenous and no
+     exogenous deterministic). The linear combination can contain constant
+     terms (intercept).
+     Returns a pair composed of:
+     – the terms of the form endogenous*scalar, as a list of (endo_id, constant expr);
+     – the sum of all constant (intercept) terms */
+  pair<vector<pair<int, expr_t>>, expr_t> matchLinearCombinationOfEndogenousWithConstant() const;
+
   pair<int, vector<tuple<int, int, int, double>>> matchParamTimesLinearCombinationOfVariables() const;
 
   /* Matches an expression of the form parameter*(var1-endo2).
@@ -663,6 +675,9 @@ public:
   //! Returns true if PacExpectationNode encountered
   virtual bool containsPacExpectation(const string &pac_model_name = "") const = 0;
 
+  //! Returns true if PacTargetNonstationaryNode encountered
+  virtual bool containsPacTargetNonstationary(const string &pac_model_name = "") const = 0;
+
   //! Decompose an expression into its additive terms
   /*! Returns a list of terms, with their sign (either 1 or -1, depending
     on whether the terms appears with a plus or a minus).
@@ -691,7 +706,14 @@ public:
   */
   tuple<int, int, int, double> matchVariableTimesConstantTimesParam(bool variable_obligatory = true) const;
 
-  //! Exception thrown by matchVariableTimesConstantTimesParam when matching fails
+  /* Matches an expression of the form endogenous*constant where constant is an
+     expression containing no endogenous, no exogenous and no exogenous deterministic.
+     Returns (endo_id, constant expr).
+     Note that it will also match a simple endogenous (in which case the
+     constant will of course be equal to one). */
+  virtual pair<int, expr_t> matchEndogenousTimesConstant() const;
+
+  //! Exception thrown when matching fails
   class MatchFailureException
   {
   public:
@@ -780,6 +802,7 @@ public:
   expr_t substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -792,6 +815,7 @@ public:
   bool isInStaticForm() const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   bool isParamTimesEndogExpr() const override;
   expr_t substituteStaticAuxiliaryVariable() const override;
 };
@@ -850,6 +874,7 @@ public:
   expr_t substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -862,10 +887,12 @@ public:
   bool isInStaticForm() const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   bool isParamTimesEndogExpr() const override;
   //! Substitute auxiliary variables by their expression in static model
   expr_t substituteStaticAuxiliaryVariable() const override;
   void matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const override;
+  pair<int, expr_t> matchEndogenousTimesConstant() const override;
 };
 
 //! Unary operator node
@@ -951,6 +978,7 @@ public:
   expr_t substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -963,6 +991,7 @@ public:
   bool isInStaticForm() const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   bool isParamTimesEndogExpr() const override;
   //! Substitute auxiliary variables by their expression in static model
   expr_t substituteStaticAuxiliaryVariable() const override;
@@ -1056,6 +1085,7 @@ public:
   expr_t substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -1077,6 +1107,7 @@ public:
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   /*
     ec_params_and_vars:
     - 1st element = feedback force parameter
@@ -1108,6 +1139,7 @@ public:
   void decomposeAdditiveTerms(vector<pair<expr_t, int>> &terms, int current_sign) const override;
   void decomposeMultiplicativeFactors(vector<pair<expr_t, int>> &factors, int current_exponent = 1) const override;
   void matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const override;
+  pair<int, expr_t> matchEndogenousTimesConstant() const override;
 };
 
 //! Trinary operator node
@@ -1187,6 +1219,7 @@ public:
   expr_t substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -1199,6 +1232,7 @@ public:
   bool isInStaticForm() const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   bool isParamTimesEndogExpr() const override;
   //! Substitute auxiliary variables by their expression in static model
   expr_t substituteStaticAuxiliaryVariable() const override;
@@ -1294,6 +1328,7 @@ public:
   expr_t substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1308,6 +1343,7 @@ public:
   bool isInStaticForm() const override;
   expr_t replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   bool isParamTimesEndogExpr() const override;
   //! Substitute auxiliary variables by their expression in static model
   expr_t substituteStaticAuxiliaryVariable() const override;
@@ -1498,7 +1534,9 @@ public:
   int maxLagWithDiffsExpanded() const override;
   expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
 };
@@ -1512,7 +1550,25 @@ public:
   int maxLagWithDiffsExpanded() const override;
   expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
+  void writeJsonAST(ostream &output) const override;
+  void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
+};
+
+class PacTargetNonstationaryNode : public SubModelNode
+{
+public:
+  PacTargetNonstationaryNode(DataTree &datatree_arg, int idx_arg, string model_name);
+  void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
+  expr_t clone(DataTree &datatree) const override;
+  int maxLagWithDiffsExpanded() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
+  expr_t substitutePacExpectation(const string &name, expr_t subexpr) override;
+  expr_t substitutePacTargetNonstationary(const string &name, expr_t subexpr) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
+  bool containsPacTargetNonstationary(const string &pac_model_name = "") const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
 };
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 902f05634034cba4e1e67061249d4427d08ac32c..916d471769d95813424b40bac430d0dd77d4baf3 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -442,7 +442,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
   original_model.fillVarModelTableFromOrigModel();
 
   // PAC model
-  pac_model_table.transformPass(diff_subst_table, dynamic_model, var_model_table,
+  pac_model_table.transformPass(unary_ops_nodes, unary_ops_subst_table,
+                                diff_nodes, diff_subst_table,
+                                dynamic_model, var_model_table,
                                 trend_component_model_table);
 
   // Create auxiliary vars for Expectation operator
@@ -1119,6 +1121,8 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
 
       // Create epilogue file
       epilogue.writeEpilogueFile(basename);
+
+      pac_model_table.writeTargetCoefficientsFile(basename);
     }
 
   cout << "done" << endl;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index bbbde0c2573b289dde98422d19c63d954bb0c7c3..3337eb7221b9575f917af3fc5e90963cb3c32aca 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -2047,3 +2047,12 @@ ModelTree::updateReverseVariableEquationOrderings()
       eq_idx_orig2block[eq_idx_block2orig[i]] = i;
     }
 }
+
+expr_t
+ModelTree::getRHSFromLHS(expr_t lhs) const
+{
+  for (auto eq : equations)
+    if (eq->arg1 == lhs)
+      return eq->arg2;
+  throw ExprNode::MatchFailureException{"Cannot find an equation with the requested LHS"};
+}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 373273e6b11708fc339c5934f5d4078e2b8658de..5421792cd1425ae8dae781d99a4ce07223d32bfc 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -444,6 +444,10 @@ public:
   //! Helper for writing the Jacobian elements in MATLAB and C
   /*! Writes either (i+1,j+1) or [i+j*no_eq] */
   void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
+  /* Given an expression, searches for the first equation that has exactly this
+     expression on the LHS, and returns the RHS of that equation.
+     If no such equation can be found, throws an ExprNode::MatchFailureExpression */
+  expr_t getRHSFromLHS(expr_t lhs) const;
 
   //! Returns all the equation tags associated to an equation
   inline map<string, string>
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 3094863436761eaa8f8a6ab84792947d7de96114..215148a4266c762c9e65b60ac339e0916ba5efeb 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2615,6 +2615,15 @@ ParsingDriver::add_pac_expectation(const string &model_name)
   return data_tree->AddPacExpectation(model_name);
 }
 
+expr_t
+ParsingDriver::add_pac_target_nonstationary(const string &model_name)
+{
+  if (data_tree == occbin_constraints_tree.get())
+    error("The 'pac_target_nonstationary' operator is forbidden in 'occbin_constraints'.");
+
+  return data_tree->AddPacTargetNonstationary(model_name);
+}
+
 void
 ParsingDriver::begin_pac_growth()
 {
@@ -3502,6 +3511,57 @@ ParsingDriver::end_occbin_constraints(const vector<tuple<string, BinaryOpNode *,
   reset_data_tree();
 }
 
+void
+ParsingDriver::begin_pac_target_info(string name)
+{
+  pac_target_info_name = move(name);
+  set_current_data_tree(&mod_file->dynamic_model);
+}
+
+void
+ParsingDriver::end_pac_target_info()
+{
+  reset_data_tree();
+}
+
+void
+ParsingDriver::set_pac_target_info_target(expr_t target)
+{
+  mod_file->pac_model_table.setTargetExpr(pac_target_info_name, target);
+}
+
+void
+ParsingDriver::set_pac_target_info_auxname_target_nonstationary(string auxname)
+{
+  mod_file->pac_model_table.setTargetAuxnameNonstationary(pac_target_info_name, move(auxname));
+}
+
+void
+ParsingDriver::add_pac_target_info_component(expr_t component_expr)
+{
+  get<0>(pac_target_info_component) = component_expr;
+  mod_file->pac_model_table.addTargetComponent(pac_target_info_name, pac_target_info_component);
+  pac_target_info_component = {};
+}
+
+void
+ParsingDriver::set_pac_target_info_component_growth(expr_t growth)
+{
+  get<1>(pac_target_info_component) = growth;
+}
+
+void
+ParsingDriver::set_pac_target_info_component_auxname(string auxname)
+{
+  get<2>(pac_target_info_component) = move(auxname);
+}
+
+void
+ParsingDriver::set_pac_target_info_component_kind(PacTargetKind kind)
+{
+  get<3>(pac_target_info_component) = kind;
+}
+
 vector<string>
 ParsingDriver::strsplit(const string &str, char delim)
 {
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index b6cd328e9ae5379fce6e40805212f56617db3286..8c9d08fb68ba045683f622dab511ebf4889efa75 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -249,6 +249,9 @@ private:
   SymbolList graph_formats;
   //! Temporary storage for equation tags
   map<string, string> eq_tags;
+  // Temporary storages for pac_target_info
+  string pac_target_info_name;
+  PacModelTable::target_component_t pac_target_info_component;
 
   //! The mod file representation constructed by this ParsingDriver
   unique_ptr<ModFile> mod_file;
@@ -741,6 +744,8 @@ public:
   expr_t add_var_expectation(const string &model_name);
   //! Writes token "PAC_EXPECTATION(model_name, discount, growth)" to model tree
   expr_t add_pac_expectation(const string &model_name);
+  //! Adds a pac_target_nonstationary(model_name, discount, growth) node to model tree
+  expr_t add_pac_target_nonstationary(const string &model_name);
   //! Creates pac_model statement
   void begin_pac_growth();
   void begin_pac_model();
@@ -897,6 +902,14 @@ public:
   void begin_model_replace(const vector<pair<string, string>> &listed_eqs_by_tags);
   // Add a var_remove statement
   void var_remove();
+  void begin_pac_target_info(string name);
+  void end_pac_target_info();
+  void set_pac_target_info_target(expr_t target);
+  void set_pac_target_info_auxname_target_nonstationary(string auxname);
+  void add_pac_target_info_component(expr_t component_expr);
+  void set_pac_target_info_component_growth(expr_t growth);
+  void set_pac_target_info_component_auxname(string auxname);
+  void set_pac_target_info_component_kind(PacTargetKind kind);
   // Equivalent of MATLAB’s strsplit. Returns an empty vector given an empty string.
   static vector<string> strsplit(const string &str, char delim);
   // Returns true iff the string is a legal symbol identifier (see NAME token in lexer)
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 28b09016c2dfd2b17b073fcc4e696510512a0d1a..73d2ce57cb5ac18a34efff77405867ef7f2cd725 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -18,6 +18,7 @@
  */
 
 #include <algorithm>
+#include <cassert>
 
 #include "SubModel.hh"
 #include "DynamicModel.hh"
@@ -791,7 +792,60 @@ PacModelTable::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation
 {
   for (auto &[name, gv] : growth)
     if (gv)
-      gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
+      {
+        if (target_info.find(name) != target_info.end())
+          {
+            cerr << "ERROR: for PAC model '" << name << "', it is not possible to declare a 'growth' option in the 'pac_model' command when there is also a 'pac_target_info' block" << endl;
+            exit(EXIT_FAILURE);
+          }
+        gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
+      }
+
+  for (const auto &[name, ti] : target_info)
+    for (auto &[expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv, gv_info] : get<2>(ti))
+      if (gv)
+        gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
+
+  for (const auto &[name, ti] : target_info)
+    {
+      auto &[target, auxname_target_nonstationary, components] = ti;
+      if (!target)
+        {
+          cerr << "ERROR: the block 'pac_target_info(" << name << ")' is missing the 'target' statement" << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (auxname_target_nonstationary.empty())
+        {
+          cerr << "ERROR: the block 'pac_target_info(" << name << ")' is missing the 'auxname_target_nonstationary' statement" << endl;
+          exit(EXIT_FAILURE);
+        }
+      int nonstationary_nb = 0;
+      for (auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : components)
+        {
+          if (auxname.empty())
+            {
+              cerr << "ERROR: the block 'pac_target_info(" << name << ")' is missing the 'auxname' statement in some 'component'" << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (kind == PacTargetKind::unspecified)
+            {
+              cerr << "ERROR: the block 'pac_target_info(" << name << ")' is missing the 'kind' statement in some 'component'" << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (kind == PacTargetKind::ll && growth_component)
+            {
+              cerr << "ERROR: in the block 'pac_target_info(" << name << ")', a component of 'kind ll' (i.e. stationary) has a 'growth' option. This is not permitted." << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (kind == PacTargetKind::dd || kind == PacTargetKind::dl)
+            nonstationary_nb++;
+        }
+      if (!nonstationary_nb)
+        {
+          cerr << "ERROR: the block 'pac_target_info(" << name << ")' must contain at least one nonstationary component (i.e. of 'kind' equal to either 'dd' or 'dl')." << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
 }
 
 void
@@ -800,6 +854,11 @@ PacModelTable::findDiffNodesInGrowth(lag_equivalence_table_t &diff_nodes) const
   for (auto &[name, gv] : growth)
     if (gv)
       gv->findDiffNodes(diff_nodes);
+
+  for (const auto &[name, ti] : target_info)
+    for (auto &[expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv, gv_info] : get<2>(ti))
+      if (gv)
+        gv->findDiffNodes(diff_nodes);
 }
 
 void
@@ -808,10 +867,18 @@ PacModelTable::substituteDiffNodesInGrowth(const lag_equivalence_table_t &diff_n
   for (auto &[name, gv] : growth)
     if (gv)
       gv = gv->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+
+  for (auto &[name, ti] : target_info)
+    for (auto &[expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv, gv_info] : get<2>(ti))
+      if (gv)
+        gv = gv->substituteDiff(diff_nodes, diff_subst_table, neweqs);
 }
 
 void
-PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
+PacModelTable::transformPass(const lag_equivalence_table_t &unary_ops_nodes,
+                             ExprNode::subst_table_t &unary_ops_subst_table,
+                             const lag_equivalence_table_t &diff_nodes,
+                             ExprNode::subst_table_t &diff_subst_table,
                              DynamicModel &dynamic_model, const VarModelTable &var_model_table,
                              const TrendComponentModelTable &trend_component_model_table)
 {
@@ -834,6 +901,123 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
             exit(EXIT_FAILURE);
           }
 
+      // Perform transformations for the pac_target_info block (if any)
+      if (target_info.find(name) != target_info.end())
+        {
+          // Substitute unary ops and diffs in the target…
+          expr_t &target = get<0>(target_info[name]);
+          vector<BinaryOpNode *> neweqs;
+          target = target->substituteUnaryOpNodes(unary_ops_nodes, unary_ops_subst_table, neweqs);
+          if (neweqs.size() > 0)
+            {
+              cerr << "ERROR: the 'target' expression of 'pac_target_info(" << name << ")' contains a variable with a unary operator that is not present in the model" << endl;
+              exit(EXIT_FAILURE);
+            }
+          target = target->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+          if (neweqs.size() > 0)
+            {
+              cerr << "ERROR: the 'target' expression of 'pac_target_info(" << name << ")' contains a diff'd variable that is not present in the model" << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // …and in component expressions
+          auto &components = get<2>(target_info[name]);
+          for (auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : components)
+            {
+              component = component->substituteUnaryOpNodes(unary_ops_nodes, unary_ops_subst_table, neweqs);
+              if (neweqs.size() > 0)
+                {
+                  cerr << "ERROR: a 'component' expression of 'pac_target_info(" << name << ")' contains a variable with a unary operator that is not present in the model" << endl;
+                  exit(EXIT_FAILURE);
+                }
+              component = component->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+              if (neweqs.size() > 0)
+                {
+                  cerr << "ERROR: a 'component' expression of 'pac_target_info(" << name << ")' contains a diff'd variable that is not present in the model" << endl;
+                  exit(EXIT_FAILURE);
+                }
+            }
+
+          /* Fill the growth_info structure.
+             Cannot be done in an earlier pass since growth terms can be
+             transformed by DynamicModel::substituteDiff(). */
+          for (auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : components)
+            {
+              if (growth_component)
+                try
+                  {
+                    growth_component_info = growth_component->matchLinearCombinationOfVariables(false);
+                  }
+                catch (ExprNode::MatchFailureException &e)
+                  {
+                    cerr << "ERROR: PAC growth must be a linear combination of variables" << endl;
+                    exit(EXIT_FAILURE);
+                  }
+            }
+
+          // Identify the model equation defining the target
+          expr_t target_expr;
+          try
+            {
+              target_expr = dynamic_model.getRHSFromLHS(target);
+            }
+          catch (ExprNode::MatchFailureException)
+            {
+              cerr << "ERROR: there is no equation whose LHS is equal to the 'target' of 'pac_target_info(" << name << ")'" << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // Parse that model equation
+          vector<pair<int, expr_t>> terms;
+          expr_t constant;
+          try
+            {
+              tie(terms, constant) = target_expr->matchLinearCombinationOfEndogenousWithConstant();
+            }
+          catch (ExprNode::MatchFailureException)
+            {
+              cerr << "ERROR: the model equation defining the 'target' of 'pac_target_info(" << name << ")' is not of the right form (should be a linear combination of endogenous variables)" << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // Associate the coefficients of the linear combination with the right components
+          for (auto [var, coeff] : terms)
+            if (auto it = find_if(components.begin(), components.end(),
+                                  [&](const auto &v) { return get<0>(v) == dynamic_model.AddVariable(var); });
+                it != components.end())
+              get<4>(*it) = coeff;
+            else
+              {
+                cerr << "ERROR: the model equation defining the 'target' of 'pac_target_info(" << name << ")' contains a variable (" << symbol_table.getName(var) << ") that is not declared as a 'component'" << endl;
+                exit(EXIT_FAILURE);
+              }
+
+          // Verify that all declared components appear in that equation
+          for (const auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : components)
+            if (!coeff)
+              {
+                cerr << "ERROR: a 'component' of 'pac_target_info(" << name << ")' does not appear in the model equation defining the 'target'" << endl;
+                exit(EXIT_FAILURE);
+              }
+
+          /* Add the variable and equation defining the stationary part of the
+             target. Note that it includes the constant. */
+          expr_t yns = constant;
+          for (const auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : components)
+            if (kind != PacTargetKind::ll)
+              yns = dynamic_model.AddPlus(yns, dynamic_model.AddTimes(coeff, component));
+          int target_nonstationary_id = symbol_table.addPacTargetNonstationaryAuxiliaryVar(get<1>(target_info[name]), yns);
+          expr_t neweq = dynamic_model.AddEqual(dynamic_model.AddVariable(target_nonstationary_id), yns);
+          dynamic_model.addEquation(neweq, -1);
+          dynamic_model.addAuxEquation(neweq);
+
+          /* Perform the substitution of the pac_target_nonstationary operator.
+             This needs to be done here, otherwise
+             DynamicModel::analyzePacEquationStructure() will not be able to
+             identify the error-correction part */
+          dynamic_model.substitutePacTargetNonstationary(name, dynamic_model.AddVariable(target_nonstationary_id, -1));
+        }
+
       // Collect some information about PAC models
       int max_lag;
       if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name[name]))
@@ -878,31 +1062,106 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
       if (growth[name])
         growth_correction_term = dynamic_model.AddTimes(growth[name], dynamic_model.AddVariable(growth_neutrality_params[name]));
       if (aux_model_name[name].empty())
-        dynamic_model.computePacModelConsistentExpectationSubstitution(name,
-                                                                       symbol_table.getID(discount[name]),
-                                                                       pacEquationMaxLag(name),
-                                                                       growth_correction_term,
-                                                                       diff_subst_table,
-                                                                       aux_var_symb_ids,
-                                                                       aux_param_symb_ids,
-                                                                       pac_expectation_substitution);
+        {
+          if (target_info.find(name) != target_info.end())
+            {
+              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);
+            }
+          else
+            dynamic_model.computePacModelConsistentExpectationSubstitution(name,
+                                                                           symbol_table.getID(discount[name]),
+                                                                           pacEquationMaxLag(name),
+                                                                           growth_correction_term,
+                                                                           diff_subst_table,
+                                                                           aux_var_symb_ids,
+                                                                           aux_param_symb_ids,
+                                                                           pac_expectation_substitution);
+        }
       else
-        dynamic_model.computePacBackwardExpectationSubstitution(name, lhs[name], max_lag,
-                                                                aux_model_type[name],
-                                                                growth_correction_term,
-                                                                aux_var_symb_ids,
-                                                                aux_param_symb_ids,
-                                                                pac_expectation_substitution);
+        {
+          if (target_info.find(name) != target_info.end())
+            {
+              assert(growth_correction_term == dynamic_model.Zero);
+              dynamic_model.computePacBackwardExpectationSubstitutionWithComponents(name, lhs[name],
+                                                                                    max_lag,
+                                                                                    aux_model_type[name],
+                                                                                    get<2>(target_info[name]),
+                                                                                    pac_expectation_substitution);
+            }
+          else
+            dynamic_model.computePacBackwardExpectationSubstitution(name, lhs[name], max_lag,
+                                                                    aux_model_type[name],
+                                                                    growth_correction_term,
+                                                                    aux_var_symb_ids,
+                                                                    aux_param_symb_ids,
+                                                                    pac_expectation_substitution);
+        }
     }
 
   // Actually perform the substitution of pac_expectation
   dynamic_model.substitutePacExpectation(pac_expectation_substitution, eq_name);
   dynamic_model.checkNoRemainingPacExpectation();
+
+  // Check that there is no remaining pac_target_nonstationary operator
+  dynamic_model.checkNoRemainingPacTargetNonstationary();
 }
 
 void
 PacModelTable::writeOutput(const string &basename, ostream &output) const
 {
+  // Helper to print the “growth_info” structure (linear decomposition of growth)
+  auto growth_info_helper = [&](const string &fieldname, const growth_info_t &gi)
+  {
+    int i = 1;
+    for (auto [growth_symb_id, growth_lag, param_id, constant] : gi)
+      {
+        string structname = fieldname + "(" + to_string(i++) + ").";
+        if (growth_symb_id >= 0)
+          {
+            string var_field = "endo_id";
+            if (symbol_table.getType(growth_symb_id) == SymbolType::exogenous)
+              {
+                var_field = "exo_id";
+                output << structname << "endo_id = 0;" << endl;
+              }
+            else
+              output << structname << "exo_id = 0;" << endl;
+            try
+              {
+                // case when this is not the highest lag of the growth variable
+                int aux_symb_id = symbol_table.searchAuxiliaryVars(growth_symb_id, growth_lag);
+                output << structname << var_field << " = " << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
+                       << structname << "lag = 0;" << endl;
+              }
+            catch (...)
+              {
+                try
+                  {
+                    // case when this is the highest lag of the growth variable
+                    int tmp_growth_lag = growth_lag + 1;
+                    int aux_symb_id = symbol_table.searchAuxiliaryVars(growth_symb_id, tmp_growth_lag);
+                    output << structname << var_field << " = " << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
+                           << structname << "lag = -1;" << endl;
+                  }
+                catch (...)
+                  {
+                    // case when there is no aux var for the variable
+                    output << structname << var_field << " = "<< symbol_table.getTypeSpecificID(growth_symb_id) + 1 << ";" << endl
+                           << structname << "lag = " << growth_lag << ";" << endl;
+                  }
+              }
+          }
+        else
+          output << structname << "endo_id = 0;" << endl
+                 << structname << "exo_id = 0;" << endl
+                 << structname << "lag = 0;" << endl;
+        output << structname << "param_id = "
+               << (param_id == -1 ? 0 : symbol_table.getTypeSpecificID(param_id) + 1) << ";" << endl
+               << structname << "constant = " << constant << ";" << endl;
+      }
+  };
+
   for (const auto &name : names)
     {
       output << "M_.pac." << name << ".auxiliary_model_name = '" << aux_model_name.at(name) << "';" << endl
@@ -913,53 +1172,7 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const
           output << "M_.pac." << name << ".growth_str = '";
           original_growth.at(name)->writeJsonOutput(output, {}, {}, true);
           output << "';" << endl;
-          int i = 0;
-          for (auto [growth_symb_id, growth_lag, param_id, constant] : growth_info.at(name))
-            {
-              string structname = "M_.pac." + name + ".growth_linear_comb(" + to_string(++i) + ").";
-              if (growth_symb_id >= 0)
-                {
-                  string var_field = "endo_id";
-                  if (symbol_table.getType(growth_symb_id) == SymbolType::exogenous)
-                    {
-                      var_field = "exo_id";
-                      output << structname << "endo_id = 0;" << endl;
-                    }
-                  else
-                    output << structname << "exo_id = 0;" << endl;
-                  try
-                    {
-                      // case when this is not the highest lag of the growth variable
-                      int aux_symb_id = symbol_table.searchAuxiliaryVars(growth_symb_id, growth_lag);
-                      output << structname << var_field << " = " << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
-                             << structname << "lag = 0;" << endl;
-                    }
-                  catch (...)
-                    {
-                      try
-                        {
-                          // case when this is the highest lag of the growth variable
-                          int tmp_growth_lag = growth_lag + 1;
-                          int aux_symb_id = symbol_table.searchAuxiliaryVars(growth_symb_id, tmp_growth_lag);
-                          output << structname << var_field << " = " << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
-                                 << structname << "lag = -1;" << endl;
-                        }
-                      catch (...)
-                        {
-                          // case when there is no aux var for the variable
-                          output << structname << var_field << " = "<< symbol_table.getTypeSpecificID(growth_symb_id) + 1 << ";" << endl
-                                 << structname << "lag = " << growth_lag << ";" << endl;
-                        }
-                    }
-                }
-              else
-                output << structname << "endo_id = 0;" << endl
-                       << structname << "exo_id = 0;" << endl
-                       << structname << "lag = 0;" << endl;
-              output << structname << "param_id = "
-                     << (param_id == -1 ? 0 : symbol_table.getTypeSpecificID(param_id) + 1) << ";" << endl
-                     << structname << "constant = " << constant << ";" << endl;
-            }
+          growth_info_helper("M_.pac." + name + ".growth_linear_comb", growth_info.at(name));
         }
     }
 
@@ -1160,6 +1373,34 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const
           output << "];" << endl;
         }
     }
+
+  for (auto &[name, val] : target_info)
+    {
+      int component_idx = 1;
+      for (auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : get<2>(val))
+        {
+          string fieldname = "M_.pac." + name + ".components(" + to_string(component_idx) + ")";
+          output << fieldname << ".aux_id = " << symbol_table.getTypeSpecificID(auxname) + 1 << ";" << endl
+                 << fieldname << ".endo_var = " << symbol_table.getTypeSpecificID(dynamic_cast<VariableNode *>(component)->symb_id) + 1 << ";" << endl
+                 << fieldname << ".kind = '" << kindToString(kind) << "';" << endl
+                 << fieldname << ".h_param_indices = [";
+          for (int id : h_indices)
+            output << symbol_table.getTypeSpecificID(id) + 1 << " ";
+          output << "];" << endl
+                 << fieldname << ".coeff_str = '";
+          coeff->writeJsonOutput(output, {}, {}, true);
+          output << "';" << endl;
+          if (growth_component)
+            {
+              output << fieldname << ".growth_neutrality_param_index = " << symbol_table.getTypeSpecificID(growth_neutrality_param) + 1 << ";" << endl
+                     << fieldname << ".growth_str = '";
+              original_growth_component->writeJsonOutput(output, {}, {}, true);
+              output << "';" << endl;
+              growth_info_helper(fieldname + ".growth_linear_comb", growth_component_info);
+            }
+          component_idx++;
+        }
+    }
 }
 
 void
@@ -1167,6 +1408,8 @@ PacModelTable::writeJsonOutput(ostream &output) const
 {
   for (const auto &name : names)
     {
+      /* The calling method has already added a comma, so don’t output one for
+         the first statement */
       if (name != *names.begin())
         output << ", ";
       output << R"({"statementName": "pac_model",)"
@@ -1181,6 +1424,31 @@ PacModelTable::writeJsonOutput(ostream &output) const
         }
       output << "}" << endl;
     }
+
+  for (auto &[name, val] : target_info)
+    {
+      output << R"(, {"statementName": "pac_target_info", "model_name": ")" << name
+             << R"(", "target": ")";
+      get<0>(val)->writeJsonOutput(output, {}, {}, true);
+      output << R"(", "auxname_target_nonstationary": ")" << get<1>(val)
+             << R"(", "components": [)";
+      for (auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : get<2>(val))
+        {
+          if (component != get<0>(get<2>(val).front()))
+            output << ", ";
+          output << R"({"component": ")";
+          component->writeJsonOutput(output, {}, {}, true);
+          output << R"(", "auxname": ")" << auxname
+                 << R"(", "kind": ")" << kindToString(kind);
+          if (growth_component)
+            {
+              output << R"(", "growth_str": ")";
+              original_growth_component->writeJsonOutput(output, {}, {}, true);
+            }
+          output << R"("})";
+        }
+      output << "]}" << endl;
+    }
 }
 
 int
@@ -1188,3 +1456,75 @@ PacModelTable::pacEquationMaxLag(const string &name_arg) const
 {
   return get<2>(equation_info.at(name_arg)).size();
 }
+
+string
+PacModelTable::kindToString(PacTargetKind kind)
+{
+  switch (kind)
+    {
+    case PacTargetKind::unspecified:
+      cerr << "Internal error: kind should not be unspecified" << endl;
+      exit(EXIT_FAILURE);
+    case PacTargetKind::ll:
+      return "ll";
+    case PacTargetKind::dl:
+      return "dl";
+    case PacTargetKind::dd:
+      return "dd";
+    }
+  // Silent GCC warning
+  assert(false);
+}
+
+void
+PacModelTable::setTargetExpr(const string &name_arg, expr_t target)
+{
+  get<0>(target_info[name_arg]) = target;
+}
+
+void
+PacModelTable::setTargetAuxnameNonstationary(const string &name_arg, string auxname)
+{
+  get<1>(target_info[name_arg]) = move(auxname);
+}
+
+void
+PacModelTable::addTargetComponent(const string &name_arg, target_component_t component)
+{
+  get<7>(component) = get<1>(component); // original_growth = growth
+  get<2>(target_info[name_arg]).emplace_back(move(component));
+}
+
+void
+PacModelTable::writeTargetCoefficientsFile(const string &basename) const
+{
+  if (target_info.empty())
+    return;
+
+  string filename = DataTree::packageDir(basename) + "/pac_target_coefficients.m";
+  ofstream output;
+  output.open(filename, ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  output << "function coeffs = pac_target_coefficients(model_name, params)" << endl;
+  for (auto &[model_name, val] : target_info)
+    {
+      output << "  if strcmp(model_name, '" << model_name << "')" << endl
+             << "    coeffs = NaN(" << get<2>(val).size() << ",1);" << endl;
+      int i = 1;
+      for (auto &[component, growth_component, auxname, kind, coeff, growth_neutrality_param, h_indices, original_growth_component, growth_component_info] : get<2>(val))
+        {
+          output << "    coeffs(" << i++ << ") = ";
+          coeff->writeOutput(output, ExprNodeOutputType::matlabDynamicModel);
+          output << ";" << endl;
+        }
+      output << "    return" << endl
+             << "  end" << endl;
+    }
+  output << "  error([ 'Unknown PAC model: ' model_name ])" << endl
+         << "end" << endl;
+  output.close();
+}
diff --git a/src/SubModel.hh b/src/SubModel.hh
index d45f29c707e392f041899c933f4b50ca9b74dc1d..74d924c74045473bfa7dfd3af8921af26ecc79c5 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -197,12 +197,15 @@ private:
   set<string> names;
   map<string, string> aux_model_name;
   map<string, string> discount;
-  // The growth expressions belong to the main dynamic_model from the ModFile instance
+  /* The growth expressions belong to the main dynamic_model from the ModFile
+     instance. The growth expression is necessarily nullptr for a model with a
+     pac_target_info block. */
   map<string, expr_t> growth, original_growth;
   /* Information about the structure of growth expressions (which must be a
      linear combination of variables).
      Each tuple represents a term: (endo_id, lag, param_id, constant) */
-  map<string, vector<tuple<int, int, int, double>>> growth_info;
+  using growth_info_t = vector<tuple<int, int, int, double>>;
+  map<string, growth_info_t> growth_info;
 
   /* Stores the name of the PAC equation associated to the model.
      pac_model_name → eq_name */
@@ -214,6 +217,8 @@ private:
      pac_expectation value
      - in the MCE case, this auxiliary represents Z₁ (i.e. without the growth
      correction term)
+     Note that this structure is not used in the presence of the
+     pac_target_info block.
      pac_model_name → symb_id */
   map<string, int> aux_var_symb_ids;
   /* Stores symb_ids for auxiliary parameters created for the expression
@@ -221,10 +226,13 @@ private:
      neutrality correction):
      - in the backward case, contains the “h” parameters
      - in the MCE case, contains the “α” parameters
+     Note that this structure is not used in the presence of the
+     pac_target_info block.
      pac_model_name → symb_ids */
   map<string, vector<int>> aux_param_symb_ids;
   /* Stores indices for growth neutrality parameters
-     pac_model_name → growth_neutrality_param_index */
+     pac_model_name → growth_neutrality_param_index.
+     This map is not used for PAC models with a pac_target_info block. */
   map<string, int> growth_neutrality_params;
 
   // Stores LHS vars (only for backward PAC models)
@@ -243,8 +251,21 @@ public:
 private:
   equation_info_t equation_info;
 
+public:
+  /* (component variable/expr, growth, auxname, kind, coeff. in the linear
+     combination, growth_param ID possibly equal to -1, vector of h parameters,
+     original_growth, growth_info) */
+  using target_component_t = tuple<expr_t, expr_t, string, PacTargetKind, expr_t, int, vector<int>, expr_t, growth_info_t>;
+
+private:
+  // pac_model_name → (target variable/expr, auxname_target_nonstationary, target components)
+  map<string, tuple<expr_t, string, vector<target_component_t>>> target_info;
+
   int pacEquationMaxLag(const string &name_arg) const;
 
+  // Return a text representation of a kind (but fails on “unspecified” kind value)
+  static string kindToString(PacTargetKind kind);
+
 public:
   explicit PacModelTable(SymbolTable &symbol_table_arg);
   void addPacModel(string name_arg, string aux_model_name_arg, string discount_arg, expr_t growth_arg);
@@ -255,11 +276,20 @@ public:
   // Called by DynamicModel::substituteDiff()
   void substituteDiffNodesInGrowth(const lag_equivalence_table_t &diff_nodes, ExprNode::subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs);
   // Must be called after substituteDiffNodesInGrowth()
-  void transformPass(ExprNode::subst_table_t &diff_subst_table,
+  void transformPass(const lag_equivalence_table_t &unary_ops_nodes,
+                     ExprNode::subst_table_t &unary_ops_subst_table,
+                     const lag_equivalence_table_t &diff_nodes,
+                     ExprNode::subst_table_t &diff_subst_table,
                      DynamicModel &dynamic_model, const VarModelTable &var_model_table,
                      const TrendComponentModelTable &trend_component_model_table);
   void writeOutput(const string &basename, ostream &output) const;
   void writeJsonOutput(ostream &output) const;
+  void setTargetExpr(const string &name_arg, expr_t target);
+  void setTargetAuxnameNonstationary(const string &name_arg, string auxname);
+  /* Only the first four elements of the tuple are expected to be set by the
+     caller. The other ones will be filled by this class. */
+  void addTargetComponent(const string &name_arg, target_component_t component);
+  void writeTargetCoefficientsFile(const string &basename) const;
 };
 
 
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 23c028ea585f9e979c63a9fc7a34b43457e8cbb8..32ea5670013005722d631780c1b2c6ecb280a84d 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -377,6 +377,7 @@ SymbolTable::writeOutput(ostream &output) const noexcept(false)
             break;
           case AuxVarType::expectation:
           case AuxVarType::pacExpectation:
+          case AuxVarType::pacTargetNonstationary:
             break;
           case AuxVarType::diff:
           case AuxVarType::diffLag:
@@ -687,6 +688,24 @@ SymbolTable::addPacExpectationAuxiliaryVar(const string &name, expr_t expr_arg)
   return symb_id;
 }
 
+int
+SymbolTable::addPacTargetNonstationaryAuxiliaryVar(const string &name, expr_t expr_arg)
+{
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(name, SymbolType::endogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: the variable/parameter '" << name << "' conflicts with a variable that will be generated for a 'pac_target_nonstationary' expression. Please rename it." << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_vars.emplace_back(symb_id, AuxVarType::pacTargetNonstationary, 0, 0, 0, 0, expr_arg, "");
+  return symb_id;
+}
+
 int
 SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const noexcept(false)
 {
@@ -1004,6 +1023,7 @@ SymbolTable::writeJsonOutput(ostream &output) const
 	      break;
 	    case AuxVarType::expectation:
 	    case AuxVarType::pacExpectation:
+            case AuxVarType::pacTargetNonstationary:
 	      break;
 	    case AuxVarType::diff:
 	    case AuxVarType::diffLag:
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index b3fd7be48db353386486e8e1bc5487669fd62c0e..9d3579ac4bb2a30c8c0c5c9736fb53f54d3fc1d6 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -53,7 +53,8 @@ enum class AuxVarType
    diffLag = 9, //!< Variable for timing between Diff operators (lag)
    unaryOp = 10, //!< Variable for allowing the undiff operator to work when diff was taken of unary op, eg diff(log(x))
    diffLead = 11, //!< Variable for timing between Diff operators (lead)
-   pacExpectation = 12 //!< Variable created for the substitution of the pac_expectation operator
+   pacExpectation = 12, //!< Variable created for the substitution of the pac_expectation operator
+   pacTargetNonstationary = 13 //!< Variable created for the substitution of the pac_target_nonstationary operator
   };
 
 //! Information on some auxiliary variables
@@ -337,6 +338,8 @@ public:
   int addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, string unary_op, int orig_symb_id = -1, int orig_lag = 0) noexcept(false);
   //! An auxiliary variable for a pac_expectation operator
   int addPacExpectationAuxiliaryVar(const string &name, expr_t expr_arg);
+  //! An auxiliary variable for a pac_target_nonstationary operator
+  int addPacTargetNonstationaryAuxiliaryVar(const string &name, expr_t expr_arg);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const