diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 27a4f3a8495edfe07e5534c62f267db1f7397f52..5d89412aaccfa50853312ad2b7588baccb152da6 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3755,6 +3755,7 @@ DynamicModel::walkPacParameters()
 
       if (equation->containsPacExpectation())
         {
+          int optim_share_index = -1;
           set<int> optim_share;
           expr_t optim_part = nullptr;
           expr_t non_optim_part = nullptr;
@@ -3766,10 +3767,12 @@ DynamicModel::walkPacParameters()
             equation->get_arg2()->getPacOptimizingPart(ec_params_and_vars, ar_params_and_vars);
           else
             {
+              optim_share_index = *(optim_share.begin());
               optim_part->getPacOptimizingPart(ec_params_and_vars, ar_params_and_vars);
               non_optim_part->getPacNonOptimizingPart(non_optim_params_vars_and_scaling_factor);
             }
           equation->addParamInfoToPac(lhs,
+                                      optim_share_index,
                                       ec_params_and_vars, ar_params_and_vars,
                                       non_optim_params_vars_and_scaling_factor);
         }
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index a3f17923ed53e4d6b7acfb7e78bef95f1e92fa68..39dce9ac8750011c8623de95a40495f78d6cb02d 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -675,7 +675,7 @@ NumConstNode::getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>
 }
 
 void
-NumConstNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+NumConstNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
 }
 
@@ -1859,7 +1859,7 @@ VariableNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-VariableNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+VariableNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
 }
 
@@ -3481,9 +3481,9 @@ UnaryOpNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-UnaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+UnaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
-  arg->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+  arg->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
 }
 
 void
@@ -5312,10 +5312,10 @@ BinaryOpNode::getPacLHS(pair<int, int> &lhs)
 }
 
 void
-BinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+BinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
-  arg1->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
-  arg2->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+  arg1->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+  arg2->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
 }
 
 void
@@ -6203,11 +6203,11 @@ TrinaryOpNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-TrinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+TrinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
-  arg1->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
-  arg2->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
-  arg3->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+  arg1->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+  arg2->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+  arg3->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
 }
 
 void
@@ -6725,10 +6725,10 @@ AbstractExternalFunctionNode::getPacOptimizingShareAndExprNodes(set<int> &optim_
 }
 
 void
-AbstractExternalFunctionNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+AbstractExternalFunctionNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
   for (auto argument : arguments)
-    argument->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
+    argument->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, params_vars_and_scaling_factor_arg);
 }
 
 void
@@ -8290,7 +8290,7 @@ VarExpectationNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-VarExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+VarExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
 }
 
@@ -8378,6 +8378,10 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     output << "M_.pac." << model_name << ".growth_neutrality_param_index = "
            << datatree.symbol_table.getTypeSpecificID(growth_param_index) + 1 << ";" << endl;
 
+  if (optim_share_index >= 0)
+    output << "M_.pac." << model_name << ".share_of_optimizing_agents_value = "
+           << datatree.symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
+
   output << "M_.pac." << model_name << ".ec.params = ";
   auto it = ec_params_and_vars.begin();
   output << datatree.symbol_table.getTypeSpecificID(it->first) + 1
@@ -8812,7 +8816,7 @@ PacExpectationNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
+PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg)
 {
   if (lhs_arg.first == -1)
     {
@@ -8827,6 +8831,7 @@ PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pai
     }
 
   lhs_pac_var = lhs_arg;
+  optim_share_index = optim_share_arg;
   ar_params_and_vars = ar_params_and_vars_arg;
   ec_params_and_vars = ec_params_and_vars_arg;
   params_vars_and_scaling_factor = params_vars_and_scaling_factor_arg;
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 875bd16273aadb855c11ec662bb239c0651f817c..1c72f15e6781dca9d12e1beb2d77a7cc48c0bc68 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -544,7 +544,7 @@ class ExprNode
                                                      expr_t &optim_part,
                                                      expr_t &non_optim_part) const = 0;
       //! Adds PAC equation param info to pac_expectation
-      virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) = 0;
+      virtual void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) = 0;
 
       //! Fills var_model info for pac_expectation node
       virtual void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) = 0;
@@ -632,7 +632,7 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
@@ -730,7 +730,7 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
@@ -852,7 +852,7 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
@@ -1001,7 +1001,7 @@ public:
   //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
   expr_t getNonZeroPartofEquation() const;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacOptimizingPart(set<pair<int, pair<int, int>>> &ec_params_and_vars,
@@ -1107,7 +1107,7 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
@@ -1227,7 +1227,7 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override = 0;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
@@ -1434,7 +1434,7 @@ public:
   expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>
@@ -1463,6 +1463,7 @@ private:
   int max_lag, pac_max_lag;
   vector<int> h0_indices, h1_indices;
   int growth_param_index, equation_number;
+  int optim_share_index;
   set<pair<int, pair<int, int>>> ec_params_and_vars;
   set<pair<int, pair<int, int>>> ar_params_and_vars;
   set<pair<int, pair<pair<int, int>, double>>> params_vars_and_scaling_factor;
@@ -1528,7 +1529,7 @@ public:
   expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
+  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, set<pair<int, pair<pair<int, int>, double>>> &params_vars_and_scaling_factor_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
   bool containsPacExpectation(const string &pac_model_name = "") const override;
   void getPacNonOptimizingPart(set<pair<int, pair<pair<int, int>, double>>>