From f4a7a5430c4209612bb3e8122652b77146002ff2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?=
 <stepan@adjemian.eu>
Date: Thu, 18 Nov 2021 12:24:23 +0100
Subject: [PATCH] Re-implement PAC equations.

 This commit only introduce new elements in the Dynare language (adding the
 possibility to decompose the target into stationary and non stationary
 components) and insure that all the former codes (ie without decomposition of
 the target) are still working as expected.
---
 matlab/+pac/+estimate/init.m                  |  31 +-
 matlab/+pac/+estimate/iterative_ols.m         | 174 +++++-----
 matlab/+pac/+estimate/nls.m                   |  40 +--
 matlab/+pac/+mce/parameters.m                 | 149 ++++-----
 matlab/+pac/+update/parameters.m              | 301 ++++++++----------
 matlab/+pac/initialize.m                      |  18 +-
 matlab/cherrypick.m                           |   4 +-
 matlab/pac-tools/+pac/+bgp/get.m              |   2 +-
 matlab/pac-tools/+pac/+bgp/set.m              |  36 +--
 matlab/pac-tools/geteqtag.m                   |  30 --
 matlab/pac-tools/hVectors.m                   |  47 ++-
 matlab/print_equations.m                      |  22 +-
 matlab/print_expectations.m                   |  87 ++---
 matlab/write_expectations.m                   |  22 +-
 preprocessor                                  |   2 +-
 tests/pac/trend-component-13a/example1.mod    |   6 +-
 tests/pac/trend-component-13b/example1.mod    |   6 +-
 tests/pac/trend-component-14/substitution.mod |   6 +-
 .../example1.mod                              |  12 +-
 tests/pac/trend-component-20-1/example1.mod   |   8 +-
 tests/pac/trend-component-20-2/example1.mod   |   8 +-
 tests/pac/trend-component-20-3/example1.mod   |   8 +-
 tests/pac/trend-component-20-4/example1.mod   |   8 +-
 tests/pac/trend-component-26/example1.mod     |  60 ++--
 tests/pac/trend-component-28/example4.mod     |   7 +-
 tests/pac/var-5/substitution.mod              |   6 +-
 tests/pac/var-6/substitution.mod              |   6 +-
 tests/pac/var-7/substitution.mod              |   6 +-
 28 files changed, 476 insertions(+), 636 deletions(-)
 delete mode 100644 matlab/pac-tools/geteqtag.m

diff --git a/matlab/+pac/+estimate/init.m b/matlab/+pac/+estimate/init.m
index 0dbdf70b8..cfed63b0f 100644
--- a/matlab/+pac/+estimate/init.m
+++ b/matlab/+pac/+estimate/init.m
@@ -1,7 +1,6 @@
-function [pacmodl, lhs, rhs, pnames, enames, xnames, rname, pid, eid, xid, pnames_, ipnames_, params, data, islaggedvariables, eqtag] = ...
-    init(M_, oo_, eqname, params, data, range)
+function [pacmodl, lhs, rhs, pnames, enames, xnames, rname, pid, eid, xid, pnames_, ipnames_, params, data, islaggedvariables] =  init(M_, oo_, eqname, params, data, range)
 
-% Copyright (C) 2018-2019 Dynare Team
+% Copyright © 2018-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -32,10 +31,16 @@ pacmodl = regexp(RHS, pattern, 'names');
 pacmodl = pacmodl.name;
 
 % Get the transformed equation to be estimated.
-[lhs, rhs] = get_lhs_and_rhs(eqname, M_);
+[lhs, rhs, json] = get_lhs_and_rhs(eqname, M_);
 
-% Get the equation tag (in M_.pac.(pacmodl).equations)
-eqtag = M_.pac.(pacmodl).tag_map{strcmp(M_.pac.(pacmodl).tag_map(:,1), eqname),2};
+% Get definition of aux. variable for pac expectation...
+auxrhs = json.model{strmatch(sprintf('pac_expectation_%s', pacmodl), M_.lhs, 'exact')}.rhs; 
+
+% ... and substitute in rhs.
+rhs = strrep(rhs, sprintf('pac_expectation_%s', pacmodl), auxrhs);
+
+% Get pacmodel properties
+pacmodel = M_.pac.(pacmodl);
 
 % Get the parameters and variables in the PAC equation.
 [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_);
@@ -70,15 +75,15 @@ end
 stack = dbstack;
 ipnames__ = ipnames_;                              % The user provided order.
 if isequal(stack(2).name, 'iterative_ols')
-    ipnames_ = [M_.pac.(pacmodl).equations.(eqtag).ec.params; M_.pac.(pacmodl).equations.(eqtag).ar.params'];
-    if isfield(M_.pac.(pacmodl).equations.(eqtag), 'optim_additive')
-        ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).optim_additive.params(~isnan(M_.pac.(pacmodl).equations.(eqtag).optim_additive.params))'];
+    ipnames_ = [pacmodel.ec.params; pacmodel.ar.params'];
+    if isfield(pacmodel, 'optim_additive')
+        ipnames_ = [ipnames_; pacmodel.optim_additive.params(~isnan(pacmodel.optim_additive.params))'];
     end
-    if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive')
-        ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).additive.params(~isnan(M_.pac.(pacmodl).equations.(eqtag).additive.params))'];
+    if isfield(pacmodel, 'additive')
+        ipnames_ = [ipnames_; pacmodel.additive.params(~isnan(pacmodel.additive.params))'];
     end
-    if isfield(M_.pac.(pacmodl).equations.(eqtag), 'share_of_optimizing_agents_index')
-        ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index];
+    if isfield(pacmodel, 'share_of_optimizing_agents_index')
+        ipnames_ = [ipnames_; pacmodel.share_of_optimizing_agents_index];
     end
     for i=1:length(ipnames_)
         if ~ismember(ipnames_(i), ipnames__)
diff --git a/matlab/+pac/+estimate/iterative_ols.m b/matlab/+pac/+estimate/iterative_ols.m
index 8c710c543..3b48b6645 100644
--- a/matlab/+pac/+estimate/iterative_ols.m
+++ b/matlab/+pac/+estimate/iterative_ols.m
@@ -22,7 +22,7 @@ function iterative_ols(eqname, params, data, range)
 %     equation must have NaN values in the object.
 % [4] It is assumed that the residual is additive.
 
-% Copyright (C) 2018-2021 Dynare Team
+% Copyright © 2018-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -41,31 +41,31 @@ function iterative_ols(eqname, params, data, range)
 
 global M_ oo_ options_
 
-[pacmodl, ~, rhs, ~, ~, ~, rname, ~, ~, ~, ~, ipnames_, params, data, ~, eqtag] = ...
+[pacmodl, ~, rhs, ~, ~, ~, rname, ~, ~, ~, ~, ipnames_, params, data, ~] = ...
     pac.estimate.init(M_, oo_, eqname, params, data, range);
 
 % Set initial condition.
 params0 = cell2mat(struct2cell(params));
 
 % Set flag for models with non optimizing agents.
-is_non_optimizing_agents = isfield(M_.pac.(pacmodl).equations.(eqtag), 'non_optimizing_behaviour');
+is_non_optimizing_agents = isfield(M_.pac.(pacmodl), 'non_optimizing_behaviour');
 
 % Set flag for models with exogenous variables (outside of non optimizing agents part)
-if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive')
-    is_exogenous_variables = length(M_.pac.(pacmodl).equations.(eqtag).additive.vars)>1;
+if isfield(M_.pac.(pacmodl), 'additive')
+    is_exogenous_variables = length(M_.pac.(pacmodl).additive.vars)>1;
 else
     is_exogenous_variables = false;
 end
 
 % Set flag for models with exogenous variables (in the optimizing agents part)
-if isfield(M_.pac.(pacmodl).equations.(eqtag), 'optim_additive')
-    is_optim_exogenous_variables = length(M_.pac.(pacmodl).equations.(eqtag).optim_additive.vars)>0;
+if isfield(M_.pac.(pacmodl), 'optim_additive')
+    is_optim_exogenous_variables = length(M_.pac.(pacmodl).optim_additive.vars)>0;
 else
     is_optim_exogenous_variables = false;
 end
 
 if is_non_optimizing_agents
-    non_optimizing_behaviour = M_.pac.(pacmodl).equations.(eqtag).non_optimizing_behaviour;
+    non_optimizing_behaviour = M_.pac.(pacmodl).non_optimizing_behaviour;
     non_optimizing_behaviour_params = NaN(length(non_optimizing_behaviour.params), 1);
     noparams = isnan(non_optimizing_behaviour.params);
     if ~all(noparams)
@@ -82,14 +82,14 @@ if is_non_optimizing_agents
     end
     non_optimizing_behaviour_params = non_optimizing_behaviour_params.*transpose(non_optimizing_behaviour.scaling_factor);
     % Set flag for the estimation of the share of non optimizing agents.
-    estimate_non_optimizing_agents_share = ismember(M_.param_names(M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index), fieldnames(params));
+    estimate_non_optimizing_agents_share = ismember(M_.param_names(M_.pac.(pacmodl).share_of_optimizing_agents_index), fieldnames(params));
     if ~estimate_non_optimizing_agents_share
-        share_of_optimizing_agents = M_.params(M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index);
+        share_of_optimizing_agents = M_.params(M_.pac.(pacmodl).share_of_optimizing_agents_index);
         if share_of_optimizing_agents>1 || share_of_optimizing_agents<0
             error('The share of optimizing agents shoud be in (0,1).')
         end
     end
-    share_of_optimizing_agents_index = M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index;
+    share_of_optimizing_agents_index = M_.pac.(pacmodl).share_of_optimizing_agents_index;
 else
     share_of_optimizing_agents = 1.0;
     share_of_optimizing_agents_index = [];
@@ -97,7 +97,7 @@ else
 end
 
 if is_exogenous_variables
-    additive = M_.pac.(pacmodl).equations.(eqtag).additive;
+    additive = M_.pac.(pacmodl).additive;
     residual_id = find(strcmp(rname, M_.exo_names));
     residual_jd = find(additive.vars==residual_id & ~additive.isendo);
     additive.params(residual_jd) = [];
@@ -111,55 +111,30 @@ else
 end
 
 if is_optim_exogenous_variables
-    optim_additive = M_.pac.(pacmodl).equations.(eqtag).optim_additive;
+    optim_additive = M_.pac.(pacmodl).optim_additive;
     optim_additive.estimation  = ismember(optim_additive.params, ipnames_);
 else
     optim_additive = struct('params', [], 'vars', [], 'isendo', [], 'lags', [], 'scaling_factor', [], 'estimation', []);
 end
 
-
 % Build PAC expectation matrix expression.
-dataForPACExpectation0 = dseries();
-listofvariables0 = {};
-if ~isempty(M_.pac.(pacmodl).equations.(eqtag).h0_param_indices)
-    for i=1:length(M_.pac.(pacmodl).equations.(eqtag).h0_param_indices)
-        match = regexp(rhs, sprintf('(?<var>((\\w*)|\\w*\\(-1\\)))\\*%s', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}), 'names');
-        if isempty(match)
-            match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}), 'names');
-        end
-        if isempty(strfind(match.var, '(-1)'))
-            listofvariables0{i} = match.var;
-            dataForPACExpectation0 = [dataForPACExpectation0, data{listofvariables0{i}}];
-        else
-            listofvariables0{i} = match.var(1:end-4);
-            dataForPACExpectation0 = [dataForPACExpectation0, data{match.var(1:end-4)}.lag(1)];
-        end
+dataForPACExpectation = dseries();
+listofvariables = {};
+for i=1:length(M_.pac.(pacmodl).h_param_indices)
+    match = regexp(rhs, sprintf('(?<var>((\\w*)|\\w*\\(-1\\)))\\*%s', M_.param_names{M_.pac.(pacmodl).h_param_indices(i)}), 'names');
+    if isempty(match)
+        match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).h_param_indices(i)}), 'names');
+    end
+    if isempty(strfind(match.var, '(-1)'))
+        listofvariables{i} = match.var;
+        dataForPACExpectation = [dataForPACExpectation, data{listofvariables{i}}];
+    else
+        listofvariables{i} = match.var(1:end-4);
+        dataForPACExpectation = [dataForPACExpectation, data{match.var(1:end-4)}.lag(1)];
     end
-    dataPAC0 = dataForPACExpectation0{listofvariables0{:}}(range).data;
-else
-    dataPAC0 = [];
 end
+dataPAC = dataForPACExpectation{listofvariables{:}}(range).data;
 
-dataForPACExpectation1 = dseries();
-listofvariables1 = {};
-if ~isempty(M_.pac.(pacmodl).equations.(eqtag).h1_param_indices)
-    for i=1:length(M_.pac.(pacmodl).equations.(eqtag).h1_param_indices)
-        match = regexp(rhs, sprintf('(?<var>((\\w*)|(\\w*\\(-1\\))))\\*%s', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}), 'names');
-        if isempty(match)
-            match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}), 'names');
-        end
-        if isempty(strfind(match.var, '(-1)'))
-            listofvariables1{i} = match.var;
-            dataForPACExpectation1 = [dataForPACExpectation1, data{listofvariables1{i}}];
-        else
-            listofvariables1{i} = match.var(1:end-4);
-            dataForPACExpectation1 = [dataForPACExpectation1, data{match.var(1:end-4)}.lag(1)];
-        end
-    end
-    dataPAC1 = dataForPACExpectation1{listofvariables1{:}}(range).data;
-else
-    dataPAC1 = [];
-end
 
 % Build data for non optimizing behaviour
 if is_non_optimizing_agents
@@ -256,10 +231,10 @@ end
 
 % Reorder ec.vars locally if necessary. Second variable must be the
 % endogenous variable, while the first must be the associated trend.
-if M_.pac.(pacmodl).equations.(eqtag).ec.isendo(2)
-    ecvars = M_.pac.(pacmodl).equations.(eqtag).ec.vars;
+if M_.pac.(pacmodl).ec.isendo(2)
+    ecvars = M_.pac.(pacmodl).ec.vars;
 else
-    ecvars = flip(M_.pac.(pacmodl).equations.(eqtag).ec.vars);
+    ecvars = flip(M_.pac.(pacmodl).ec.vars);
 end
 
 %% Build matrix for EC and AR terms.
@@ -268,12 +243,12 @@ DataForOLS = dseries();
 % Error correction term is trend minus the level of the endogenous variable.
 DataForOLS{'ec-term'} = data{M_.endo_names{ecvars(1)}}.lag(1)-data{M_.endo_names{ecvars(2)}}.lag(1);
 listofvariables3 = {'ec-term'};
-xparm = { M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params(1))};
-for i = 1:length(M_.pac.(pacmodl).equations.(eqtag).ar.params)
-    if islagof(M_.pac.(pacmodl).equations.(eqtag).ar.vars(i), M_.pac.(pacmodl).equations.(eqtag).lhs_var)
-        DataForOLS = [DataForOLS, data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).ar.vars(i)}}];
-        listofvariables3{i+1} = M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).ar.vars(i)};
-        xparm{i+1} = M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ar.params(i));
+xparm = { M_.param_names(M_.pac.(pacmodl).ec.params(1))};
+for i = 1:length(M_.pac.(pacmodl).ar.params)
+    if islagof(M_.pac.(pacmodl).ar.vars(i), M_.pac.(pacmodl).lhs_var)
+        DataForOLS = [DataForOLS, data{M_.endo_names{M_.pac.(pacmodl).ar.vars(i)}}];
+        listofvariables3{i+1} = M_.endo_names{M_.pac.(pacmodl).ar.vars(i)};
+        xparm{i+1} = M_.param_names(M_.pac.(pacmodl).ar.params(i));
     end
 end
 
@@ -329,7 +304,7 @@ end
 M_.params(ipnames_) = params0;
 
 % Update the reduced form PAC expectation parameters and compute the expectations.
-[PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
+[PacExpectations, M_] = UpdatePacExpectationsData(dataPAC, data, range, pacmodl, M_, oo_);
 
 noconvergence = true;
 counter = 0;
@@ -337,7 +312,7 @@ counter = 0;
 while noconvergence
     counter = counter+1;
     % Set vector for left handside variable
-    YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
+    YDATA = data{M_.endo_names{M_.pac.(pacmodl).lhs_var}}(range).data;
     if is_non_optimizing_agents
         YDATA = YDATA-share_of_optimizing_agents*PacExpectations;
         YDATA = YDATA-(1-share_of_optimizing_agents)*(dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params);
@@ -368,16 +343,16 @@ while noconvergence
     if estimate_non_optimizing_agents_share
         % First update the parameters and compute the PAC expectation reduced form parameters.
         M_.params(ipnames_(params_id_4)) = params0_;
-        [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
+        [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC, data, range, pacmodl, M_, oo_);
         % Set vector for left handside variable.
-        YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
+        YDATA = data{M_.endo_names{M_.pac.(pacmodl).lhs_var}}(range).data;
         YDATA = YDATA-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params;
         if is_exogenous_variables
             if is_any_calibrated_parameter_x
                 YDATA = YDATA-dataForExogenousVariables_(range).data;
             end
             if is_any_estimated_parameter_x
-               YDATA = YDATA-XDATA(:, length(params_id_1):end)*params0_(length(params_id_1):end);
+                YDATA = YDATA-XDATA(:, length(params_id_1):end)*params0_(length(params_id_1):end);
             end
         end
         % Set vector for regressor.
@@ -397,48 +372,41 @@ while noconvergence
         M_.params(ipnames_(params_id_0)) = share_of_optimizing_agents;
     else
         M_.params(ipnames_) = params0_;
-        [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
+        [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC, data, range, pacmodl, M_, oo_);
     end
 end
 
 % Save results
-oo_.pac.(pacmodl).equations.(eqtag).ssr = ssr;
-oo_.pac.(pacmodl).equations.(eqtag).residual = r;
-oo_.pac.(pacmodl).equations.(eqtag).estimator = params0_;
-oo_.pac.(pacmodl).equations.(eqtag).parnames = fieldnames(params);
-oo_.pac.(pacmodl).equations.(eqtag).covariance = NaN(length(params0_));
-oo_.pac.(pacmodl).equations.(eqtag).student = NaN(size(params0_));
+oo_.pac.(pacmodl).ssr = ssr;
+oo_.pac.(pacmodl).residual = r;
+oo_.pac.(pacmodl).estimator = params0_;
+oo_.pac.(pacmodl).parnames = fieldnames(params);
+oo_.pac.(pacmodl).covariance = NaN(length(params0_));
+oo_.pac.(pacmodl).student = NaN(size(params0_));
 
 
-function [PacExpectations, Model] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, Model, Output)
-    % Update PAC reduced parameters.
-    Model = pac.update.parameters(pacmodl, Model, Output, false);
-    % Compute PAC expectation.
-    if isempty(dataPAC0)
-        PacExpectations = 0;
-    else
-        PacExpectations = dataPAC0*Model.params(Model.pac.(pacmodl).equations.(eqtag).h0_param_indices);
-    end
-    if ~isempty(dataPAC1)
-        PacExpectations = PacExpectations+dataPAC1*Model.params(Model.pac.(pacmodl).equations.(eqtag).h1_param_indices);
-    end
-    % Add correction for growth neutrality if required.
-    correction = 0;
-    if isfield(Model.pac.(pacmodl), 'growth_linear_comb')
-        for iter = 1:numel(Model.pac.(pacmodl).growth_linear_comb)
-            GrowthVariable = Model.pac.(pacmodl).growth_linear_comb(iter).constant;
-            if Model.pac.(pacmodl).growth_linear_comb(iter).param_id > 0
-                GrowthVariable = GrowthVariable*Model.params(Model.pac.(pacmodl).growth_linear_comb(iter).param_id);
-            end
-            if Model.pac.(pacmodl).growth_linear_comb(iter).exo_id > 0
-                GrowthVariable = GrowthVariable*data{Model.exo_names{Model.pac.(pacmodl).growth_linear_comb(iter).exo_id}}.lag(abs(Model.pac.(pacmodl).growth_linear_comb(iter).lag));
-                GrowthVariable = GrowthVariable(range).data;
-            elseif Model.pac.(pacmodl).growth_linear_comb(iter).endo_id > 0
-                GrowthVariable = GrowthVariable*data{Model.endo_names{Model.pac.(pacmodl).growth_linear_comb(iter).endo_id}}.lag(abs(Model.pac.(pacmodl).growth_linear_comb(iter).lag));
-                GrowthVariable = GrowthVariable(range).data;
-            end
-            correction = correction + GrowthVariable;
+function [PacExpectations, Model] = UpdatePacExpectationsData(dataPAC, data, range, pacmodl, Model, Output)
+% Update PAC reduced parameters.
+Model = pac.update.parameters(pacmodl, Model, Output, false);
+% Compute PAC expectation.
+PacExpectations = dataPAC*Model.params(Model.pac.(pacmodl).h_param_indices);
+% Add correction for growth neutrality if required.
+correction = 0;
+if isfield(Model.pac.(pacmodl), 'growth_linear_comb')
+    for iter = 1:numel(Model.pac.(pacmodl).growth_linear_comb)
+        GrowthVariable = Model.pac.(pacmodl).growth_linear_comb(iter).constant;
+        if Model.pac.(pacmodl).growth_linear_comb(iter).param_id > 0
+            GrowthVariable = GrowthVariable*Model.params(Model.pac.(pacmodl).growth_linear_comb(iter).param_id);
         end
-        correction = correction*Model.params(Model.pac.(pacmodl).growth_neutrality_param_index);
+        if Model.pac.(pacmodl).growth_linear_comb(iter).exo_id > 0
+            GrowthVariable = GrowthVariable*data{Model.exo_names{Model.pac.(pacmodl).growth_linear_comb(iter).exo_id}}.lag(abs(Model.pac.(pacmodl).growth_linear_comb(iter).lag));
+            GrowthVariable = GrowthVariable(range).data;
+        elseif Model.pac.(pacmodl).growth_linear_comb(iter).endo_id > 0
+            GrowthVariable = GrowthVariable*data{Model.endo_names{Model.pac.(pacmodl).growth_linear_comb(iter).endo_id}}.lag(abs(Model.pac.(pacmodl).growth_linear_comb(iter).lag));
+            GrowthVariable = GrowthVariable(range).data;
+        end
+        correction = correction + GrowthVariable;
     end
-    PacExpectations = PacExpectations+correction;
+    correction = correction*Model.params(Model.pac.(pacmodl).growth_neutrality_param_index);
+end
+PacExpectations = PacExpectations+correction;
diff --git a/matlab/+pac/+estimate/nls.m b/matlab/+pac/+estimate/nls.m
index 8dbb66651..b9c48477f 100644
--- a/matlab/+pac/+estimate/nls.m
+++ b/matlab/+pac/+estimate/nls.m
@@ -42,7 +42,7 @@ function nls(eqname, params, data, range, optimizer, varargin)
 % is available only if the matylab optimization toolbox is installed), the
 % remaining inputs are the options (key/value) passed to the optimizers.
 
-% Copyright (C) 2018-2021 Dynare Team
+% Copyright © 2018-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -72,11 +72,11 @@ if nargin>4 && (isequal(optimizer, 'GaussNewton') || isequal(optimizer, 'lsqnonl
     end
 end
 
-[pacmodl, lhs, rhs, pnames, enames, xnames, ~, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables, eqtag] = ...
+[pacmodl, lhs, rhs, pnames, enames, xnames, ~, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables] = ...
     pac.estimate.init(M_, oo_, eqname, params, data, range);
 
 % Check that the error correction term is correct.
-if M_.pac.(pacmodl).equations.(eqtag).ec.istarget(2)
+if M_.pac.(pacmodl).ec.istarget(2)
     error(['\nThe error correction term in PAC equation (%s) is not correct.\nThe ' ...
            'error correction term should be the difference between a trend\n' ...
            'and the level of the endogenous variable.'], pacmodl);
@@ -197,8 +197,8 @@ else
         % Nothing to do here.
       case 'lsqnonlin'
         bounds = ones(length(params0),1)*[-Inf,Inf];
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),1)  = 0.0;
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),2)  = 1.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),1)  = 0.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),2)  = 1.0;
       case 'fmincon'
         if isoctave && ~user_has_octave_forge_package('optim', '1.6')
             error('Optimization algorithm ''fmincon'' requires the optim package, version 1.6 or higher')
@@ -207,8 +207,8 @@ else
         end
         minalgo = 1;
         bounds = ones(length(params0),1)*[-Inf,Inf];
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),1)  = 0.0;
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),2)  = 1.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),1)  = 0.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),2)  = 1.0;
       case 'fminunc'
         if isoctave && ~user_has_octave_forge_package('optim')
             error('Optimization algorithm ''fminunc'' requires the optim package')
@@ -230,8 +230,8 @@ else
       case 'annealing'
         minalgo = 2;
         bounds = ones(length(params0),1)*[-Inf,Inf];
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),1)  = 0.0;
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),2)  = 1.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),1)  = 0.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),2)  = 1.0;
         parameter_names = fieldnames(params);
       case 'particleswarm'
         if isoctave
@@ -241,8 +241,8 @@ else
         end
         minalgo = 12;
         bounds = ones(length(params0),1)*[-Inf,Inf];
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),1)  = 0.0;
-        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).equations.(eqtag).ec.params)),1)  = 1.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),1)  = 0.0;
+        bounds(strcmp(fieldnames(params), M_.param_names(M_.pac.(pacmodl).ec.params)),1)  = 1.0;
         parameter_names = fieldnames(params);
       otherwise
         msg = sprintf('%s is not an implemented optimization routine.\n', optimizer);
@@ -330,15 +330,15 @@ C = C/T;
 
 % Save results
 lhs = eval(strrep(lhs, 'data', 'data(range(1)-1:range(end)).data'));
-oo_.pac.(pacmodl).equations.(eqtag).lhs = dseries(lhs, range(1), sprintf('%s_lhs', eqtag));
-oo_.pac.(pacmodl).equations.(eqtag).fit = dseries(lhs-r, range(1), sprintf('%s_fit', eqtag));
-oo_.pac.(pacmodl).equations.(eqtag).residual = dseries(r, range(1), sprintf('%s_residual', eqtag));
-oo_.pac.(pacmodl).equations.(eqtag).ssr = SSR;
-oo_.pac.(pacmodl).equations.(eqtag).R2 = 1-var(r)/var(lhs);
-oo_.pac.(pacmodl).equations.(eqtag).parnames = fieldnames(params);
-oo_.pac.(pacmodl).equations.(eqtag).estimator = params1;
-oo_.pac.(pacmodl).equations.(eqtag).covariance = C;
-oo_.pac.(pacmodl).equations.(eqtag).student = params1./(sqrt(diag(C)));
+oo_.pac.(pacmodl).lhs = dseries(lhs, range(1), 'lhs');
+oo_.pac.(pacmodl).fit = dseries(lhs-r, range(1), 'fit');
+oo_.pac.(pacmodl).residual = dseries(r, range(1), 'residual');
+oo_.pac.(pacmodl).ssr = SSR;
+oo_.pac.(pacmodl).R2 = 1-var(r)/var(lhs);
+oo_.pac.(pacmodl).parnames = fieldnames(params);
+oo_.pac.(pacmodl).estimator = params1;
+oo_.pac.(pacmodl).covariance = C;
+oo_.pac.(pacmodl).student = params1./(sqrt(diag(C)));
 
 % Also save estimated parameters in M_
 M_.params(ipnames_) = params1;
diff --git a/matlab/+pac/+mce/parameters.m b/matlab/+pac/+mce/parameters.m
index 1e968496c..6098a9fdf 100644
--- a/matlab/+pac/+mce/parameters.m
+++ b/matlab/+pac/+mce/parameters.m
@@ -49,97 +49,80 @@ if ~pacmodel.model_consistent_expectations
 end
 
 % Show the equations where this PAC model is used.
-number_of_pac_eq = size(pacmodel.tag_map, 1);
-if number_of_pac_eq==1
-    fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.tag_map{1,1});
-else
-    fprintf('PAC model %s is used in %u equation(s):\n', pacname, number_of_pac_eq);
-    skipline()
-    for i=1:number_of_pac_eq
-        fprintf('    -  %s\n', pacmodel.tag_map{i,1});
-    end
-end
+fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.eq_name);
 skipline()
 
-equations = pacmodel.equations;
-
-for e=1:number_of_pac_eq
-    % Get PAC equation tag
-    eqtag = pacmodel.tag_map{e,2};
-    % Get PAC equation
-    pac_equation = equations.(eqtag);
-    % Get Error correction and autoregressive parameters in PAC equation
-    params = NaN(2+pac_equation.max_lag, 1);
-    params(1) = M_.params(pac_equation.ec.params);
-    params(1+(1:pac_equation.max_lag)) = M_.params(pac_equation.ar.params);
-    params(end) = M_.params(pacmodel.discount_index);
-    [G, alpha, beta] = buildGmatrixWithAlphaAndBeta(params);
-    M_.params(pac_equation.mce.alpha) = flip(alpha);
-    if isfield(pacmodel, 'growth_neutrality_param_index')
-        if isfield(equations.(eqtag), 'non_optimizing_behaviour')
-            gamma = M_.params(equations.(eqtag).share_of_optimizing_agents_index);
-        else
-            gamma = 1.0;
-        end
-        A = [alpha; 1];
-        A_1 = polyval(A, 1.0);
-        A_b = polyval(A, beta);
-        m = length(alpha);
-        d = A_1*A_b*(iota(m, m)'*inv((eye(m)-G)*(eye(m)-G))*iota(m, m));
-        cc = 1/gamma-(sum(params(2:end-1))+d);
-        ll = 0;
-        if isfield(equations.(eqtag), 'optim_additive')
-            % Exogenous variables are present in the λ part (optimizing agents).
-            tmp0 = 0;
-            for i=1:length(equations.(eqtag).optim_additive.params)
-                if isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
-                    tmp0 = tmp0 + equations.(eqtag).optim_additive.scaling_factor(i);
-                elseif ~isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
-                    tmp0 = tmp0 + M_.params(equations.(eqtag).optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i);
-                elseif ~islogical(equations.(eqtag).optim_additive.bgp{i})
-                    error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
-                end
-            end
-            cc = cc - tmp0;
-        end
-        if gamma<1
-            if isfield(equations.(eqtag), 'non_optimizing_behaviour') && isfield(equations.(eqtag).non_optimizing_behaviour, 'params')
-                % Exogenous variables are present in the 1-λ part (rule of thumb agents).
-                tmp0 = 0;
-                tmp1 = 0;
-                for i=1:length(equations.(eqtag).non_optimizing_behaviour.params)
-                    if isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
-                        tmp0 = tmp0 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
-                    elseif ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
-                        tmp0 = tmp0 + M_.params(equations.(eqtag).non_optimizing_behaviour.params(i))*equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
-                    elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
-                        tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
-                    elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
-                        tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.params(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
-                    end
-                end
-                cc = cc - (1-gamma)*tmp0/gamma;
-                ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
+% Get Error correction and autoregressive parameters in PAC equation
+params = NaN(2+pacmodel.max_lag, 1);
+params(1) = M_.params(pacmodel.ec.params);
+params(1+(1:pacmodel.max_lag)) = M_.params(pacmodel.ar.params);
+params(end) = M_.params(pacmodel.discount_index);
+[G, alpha, beta] = buildGmatrixWithAlphaAndBeta(params);
+M_.params(pacmodel.mce.alpha) = flip(alpha);
+if isfield(pacmodel, 'growth_neutrality_param_index')
+    if isfield(pacmodel, 'non_optimizing_behaviour')
+        gamma = M_.params(pacmodel.share_of_optimizing_agents_index);
+    else
+        gamma = 1.0;
+    end
+    A = [alpha; 1];
+    A_1 = polyval(A, 1.0);
+    A_b = polyval(A, beta);
+    m = length(alpha);
+    d = A_1*A_b*(iota(m, m)'*inv((eye(m)-G)*(eye(m)-G))*iota(m, m));
+    cc = 1/gamma-(sum(params(2:end-1))+d);
+    ll = 0;
+    if isfield(pacmodel, 'optim_additive')
+        % Exogenous variables are present in the λ part (optimizing agents).
+        tmp0 = 0;
+        for i=1:length(pacmodel.optim_additive.params)
+            if isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
+                tmp0 = tmp0 + pacmodel.optim_additive.scaling_factor(i);
+            elseif ~isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
+                tmp0 = tmp0 + M_.params(pacmodel.optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i);
+            elseif ~islogical(pacmodel.optim_additive.bgp{i})
+                error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
             end
         end
-        if isfield(equations.(eqtag), 'additive')
-            % Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
+        cc = cc - tmp0;
+    end
+    if gamma<1
+        if isfield(pacmodel, 'non_optimizing_behaviour') && isfield(pacmodel.non_optimizing_behaviour, 'params')
+            % Exogenous variables are present in the 1-λ part (rule of thumb agents).
             tmp0 = 0;
             tmp1 = 0;
-            for i=1:length(equations.(eqtag).additive.params)
-                if isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
-                    tmp0 = tmp0 + equations.(eqtag).additive.scaling_factor(i);
-                elseif ~isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
-                    tmp0 = tmp0 + M_.params(equations.(eqtag).additive.params(i))*equations.(eqtag).additive.scaling_factor(i);
-                elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && isnan(equations.(eqtag).additive.params(i))
-                    tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp{i};
-                elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && ~isnan(equations.(eqtag).additive.params(i))
-                    tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.params(i)*equations.(eqtag).additive.bgp{i};
+            for i=1:length(pacmodel.non_optimizing_behaviour.params)
+                if isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
+                    tmp0 = tmp0 + pacmodel.non_optimizing_behaviour.scaling_factor(i);
+                elseif ~isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
+                    tmp0 = tmp0 + M_.params(pacmodel.non_optimizing_behaviour.params(i))*pacmodel.non_optimizing_behaviour.scaling_factor(i);
+                elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && isnan(pacmodel.non_optimizing_behaviour.params(i))
+                    tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.bgp{i};
+                elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && ~isnan(pacmodel.non_optimizing_behaviour.params(i))
+                    tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.params(i)*pacmodel.non_optimizing_behaviour.bgp{i};
                 end
             end
-            cc = cc - tmp0/gamma;
-            ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
+            cc = cc - (1-gamma)*tmp0/gamma;
+            ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
         end
-        M_.params(pacmodel.growth_neutrality_param_index) = cc; % Multiplies the variable or expression provided though the growth option in command pac_model.
     end
+    if isfield(pacmodel, 'additive')
+        % Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
+        tmp0 = 0;
+        tmp1 = 0;
+        for i=1:length(pacmodel.additive.params)
+            if isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
+                tmp0 = tmp0 + pacmodel.additive.scaling_factor(i);
+            elseif ~isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
+                tmp0 = tmp0 + M_.params(pacmodel.additive.params(i))*pacmodel.additive.scaling_factor(i);
+            elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && isnan(pacmodel.additive.params(i))
+                tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.bgp{i};
+            elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && ~isnan(pacmodel.additive.params(i))
+                tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.params(i)*pacmodel.additive.bgp{i};
+            end
+        end
+        cc = cc - tmp0/gamma;
+        ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
+    end
+    M_.params(pacmodel.growth_neutrality_param_index) = cc; % Multiplies the variable or expression provided though the growth option in command pac_model.
 end
\ No newline at end of file
diff --git a/matlab/+pac/+update/parameters.m b/matlab/+pac/+update/parameters.m
index 8ddb17697..30b454aa0 100644
--- a/matlab/+pac/+update/parameters.m
+++ b/matlab/+pac/+update/parameters.m
@@ -66,181 +66,154 @@ if ~isfield(varcalib, 'CompanionMatrix') || any(isnan(varcalib.CompanionMatrix(:
 end
 
 % Show the equations where this PAC model is used.
-number_of_pac_eq = size(pacmodel.tag_map, 1);
-if verbose
-    fprintf('PAC model %s is used in %u equation(s):\n', pacname, number_of_pac_eq);
-    skipline()
-    for i=1:number_of_pac_eq
-        fprintf('    - %s\n', pacmodel.tag_map{i,1});
+fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.eq_name);
+skipline()
+
+% Build the vector of PAC parameters (ECM parameter + autoregressive parameters).
+pacvalues = DynareModel.params([pacmodel.ec.params; pacmodel.ar.params(1:pacmodel.max_lag)']);
+% Get the indices for the stationary/nonstationary variables in the VAR system.
+id = find(strcmp(DynareModel.endo_names{pacmodel.ec.vars(pacmodel.ec.istarget)}, varmodel.list_of_variables_in_companion_var));
+if isempty(id)
+    % Find the auxiliary variables if any
+    ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
+    if isempty(ad)
+        error('Cannot find the trend variable in the Companion VAR/VECM model.')
+    else
+        for i=1:length(ad)
+            auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)}));
+            if isequal(auxinfo.endo_index, pacmodel.ec.vars(pacmodel.ec.istarget))
+                id = ad(i);
+                break
+            end
+            if isequal(auxinfo.type, 8) && isequal(auxinfo.orig_index, pacmodel.ec.vars(pacmodel.ec.istarget))
+                id = ad(i);
+                break
+            end
+        end
+    end
+    if isempty(id)
+        error('Cannot find the trend variable in the Companion VAR/VECM model.')
     end
-    skipline()
 end
 
-equations = pacmodel.equations;
+% Infer the kind of PAC exoectation
+if isequal(pacmodel.auxiliary_model_type, 'var')
+    if varmodel.nonstationary(id)
+        kind = 'dd';
+        if varmodel.isconstant
+            id = id+1;
+        end
+    else
+        kind = 'll'
+        if varmodel.isconstant
+            id = id+1;
+        end
+    end
+else
+    % Trend component model is assumed.
+    kind = 'dd';
+end
 
-for e=1:number_of_pac_eq
-    eqtag = pacmodel.tag_map{e,2};
-    % Build the vector of PAC parameters (ECM parameter + autoregressive parameters).
-    pacvalues = DynareModel.params([equations.(eqtag).ec.params; equations.(eqtag).ar.params(1:equations.(eqtag).max_lag)']);
-    % Get the indices for the stationary/nonstationary variables in the VAR system.
-    id = find(strcmp(DynareModel.endo_names{equations.(eqtag).ec.vars(equations.(eqtag).ec.istarget)}, varmodel.list_of_variables_in_companion_var));
-    if isempty(id)
-        % Find the auxiliary variables if any
-        ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
-        if isempty(ad)
-            error('Cannot find the trend variable in the Companion VAR/VECM model.')
-        else
-            for i=1:length(ad)
-                auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)}));
-                if isequal(auxinfo.endo_index, equations.(eqtag).ec.vars(equations.(eqtag).ec.istarget))
-                    id = ad(i);
-                    break
-                end
-                if isequal(auxinfo.type, 8) && isequal(auxinfo.orig_index, equations.(eqtag).ec.vars(equations.(eqtag).ec.istarget))
-                    id = ad(i);
-                    break
-                end
-            end
-        end
-        if isempty(id)
-            error('Cannot find the trend variable in the Companion VAR/VECM model.')
-        end
-    end
-    if isequal(pacmodel.auxiliary_model_type, 'var')
-        if varmodel.nonstationary(id)
-            idns = id;
-            ids = [];
-            if varmodel.isconstant
-                idns = idns+1;
-            end
-        else
-            idns = [];
-            ids = id;
-            if varmodel.isconstant
-                ids = ids+1;
-            end
-        end
+% Override kind with the information provided by the user or update M_.pac
+
+
+% Get the value of the discount factor.
+beta = DynareModel.params(pacmodel.discount_index);
+
+% Is growth argument passed to pac_expectation?
+if isfield(pacmodel, 'growth_str')
+    growth_flag = true;
+else
+    growth_flag = false;
+end
+
+% Do we have rule of thumb agents? γ is the share of optimizing agents.
+if isfield(pacmodel, 'non_optimizing_behaviour')
+    gamma = DynareModel.params(pacmodel.share_of_optimizing_agents_index);
+else
+    gamma = 1.0;
+end
+
+% Get h vector (plus the parameter for the growth neutrality correction).
+if growth_flag
+    [h, growthneutrality] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind, id);
+else
+    h = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind, id);
+end
+
+% Update M_.params with h
+if isequal(pacmodel.auxiliary_model_type, 'var')
+    if DynareModel.var.(pacmodel.auxiliary_model_name).isconstant
+        DynareModel.params(pacmodel.h_param_indices) = h;
     else
-        % Trend component model is assumed.
-        ids = [];
-        idns = id;
+        DynareModel.params(pacmodel.h_param_indices(1)) = .0;
+        DynareModel.params(pacmodel.h_param_indices(2:end)) = h;
     end
-    % Get the value of the discount factor.
-    beta = DynareModel.params(pacmodel.discount_index);
-    % Is growth argument passed to pac_expectation?
-    if isfield(pacmodel, 'growth_str')
-        growth_flag = true;
-    else
-        growth_flag = false;
-    end
-    % Do we have rule of thumb agents? γ is the share of optimizing agents.
-    if isfield(equations.(eqtag), 'non_optimizing_behaviour')
-        gamma = DynareModel.params(equations.(eqtag).share_of_optimizing_agents_index);
-    else
-        gamma = 1.0;
-    end
-    % Get h0 and h1 vectors (plus the parameter for the growth neutrality correction).
-    if growth_flag
-        [h0, h1, growthneutrality] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, ids, idns, pacmodel.auxiliary_model_type);
-    else
-        [h0, h1] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, ids, idns, pacmodel.auxiliary_model_type);
-    end
-    % Update the parameters related to the stationary components.
-    if ~isempty(h0)
-        if isequal(pacmodel.auxiliary_model_type, 'var')
-            if DynareModel.var.(pacmodel.auxiliary_model_name).isconstant
-                DynareModel.params(equations.(eqtag).h0_param_indices) = h0;
-            else
-                DynareModel.params(equations.(eqtag).h0_param_indices(1)) = .0;
-                DynareModel.params(equations.(eqtag).h0_param_indices(2:end)) = h0;
-            end
-        else
-            DynareModel.params(equations.(eqtag).h0_param_indices) = h0;
-        end
-        DynareModel.params(pacmodel.equations.(eqtag).h0_param_indices) = h0;
-    else
-        if ~isempty(equations.(eqtag).h0_param_indices)
-            DynareModel.params(equations.(eqtag).h0_param_indices) = .0;
-        end
-    end
-    % Update the parameters related to the nonstationary components.
-    if ~isempty(h1)
-        if isequal(pacmodel.auxiliary_model_type, 'var')
-            if DynareModel.var.(pacmodel.auxiliary_model_name).isconstant
-                DynareModel.params(equations.(eqtag).h1_param_indices) = h1;
-            else
-                DynareModel.params(equations.(eqtag).h1_param_indices(1)) = .0;
-                DynareModel.params(equations.(eqtag).h1_param_indices(2:end)) = h1;
-            end
-        else
-            DynareModel.params(equations.(eqtag).h1_param_indices) = h1;
-        end
-    else
-        if ~isempty(equations.(eqtag).h1_param_indices)
-            DynareModel.params(equations.(eqtag).h1_param_indices) = .0;
-        end
-    end
-    % Update the parameter related to the growth neutrality correction.
-    if growth_flag
-        % Growth neutrality as returned by hVector is valid iff
-        % there is no exogenous variables in the model and in the
-        % absence of non optimizing agents.
-        gg = -(growthneutrality-1); % Finite sum of autoregressive parameters + infinite sum of the coefficients in the PAC expectation term.
-        cc = 1.0-gg*gamma;          % First adjustment of the growth neutrality correction (should also be divided by gamma, done below at the end of this section).
-                                    % We may have to further change the correction if we have nonzero mean exogenous variables.
-        ll = 0.0;
-        if isfield(equations.(eqtag), 'optim_additive')
-            % Exogenous variables are present in the λ part (optimizing agents).
-            tmp0 = 0;
-            for i=1:length(equations.(eqtag).optim_additive.params)
-                if isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
-                    tmp0 = tmp0 + equations.(eqtag).optim_additive.scaling_factor(i);
-                elseif ~isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
-                    tmp0 = tmp0 + DynareModel.params(equations.(eqtag).optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i);
-                elseif ~islogical(equations.(eqtag).optim_additive.bgp{i})
-                    error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
-                end
-            end
-            cc = cc - tmp0*gamma;
-        end
-        if gamma<1
-            if isfield(equations.(eqtag), 'non_optimizing_behaviour') && isfield(equations.(eqtag).non_optimizing_behaviour, 'params')
-                % Exogenous variables are present in the 1-λ part (rule of thumb agents).
-                tmp0 = 0;
-                tmp1 = 0;
-                for i=1:length(equations.(eqtag).non_optimizing_behaviour.params)
-                    if isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
-                        tmp0 = tmp0 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
-                    elseif ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
-                        tmp0 = tmp0 + DynareModel.params(equations.(eqtag).non_optimizing_behaviour.params(i))*equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
-                    elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
-                        tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
-                    elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
-                        tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.params(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
-                    end
-                end
-                cc = cc - (1.0-gamma)*tmp0;
-                ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
+else
+    DynareModel.params(pacmodel.h_param_indices) = h;
+end
+
+
+% Update the parameter related to the growth neutrality correction.
+if growth_flag
+    % Growth neutrality as returned by hVector is valid iff
+    % there is no exogenous variables in the model and in the
+    % absence of non optimizing agents.
+    gg = -(growthneutrality-1); % Finite sum of autoregressive parameters + infinite sum of the coefficients in the PAC expectation term.
+    cc = 1.0-gg*gamma;          % First adjustment of the growth neutrality correction (should also be divided by gamma, done below at the end of this section).
+                                % We may have to further change the correction if we have nonzero mean exogenous variables.
+    ll = 0.0;
+    if isfield(pacmodel, 'optim_additive')
+        % Exogenous variables are present in the λ part (optimizing agents).
+        tmp0 = 0;
+        for i=1:length(pacmodel.optim_additive.params)
+            if isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
+                tmp0 = tmp0 + pacmodel.optim_additive.scaling_factor(i);
+            elseif ~isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
+                tmp0 = tmp0 + DynareModel.params(pacmodel.optim_additive.params(i))*pacmodel.optim_additive.scaling_factor(i);
+            elseif ~islogical(pacmodel.optim_additive.bgp{i})
+                error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
             end
         end
-        if isfield(equations.(eqtag), 'additive')
-            % Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
+        cc = cc - tmp0*gamma;
+    end
+    if gamma<1
+        if isfield(pacmodel, 'non_optimizing_behaviour') && isfield(pacmodel.non_optimizing_behaviour, 'params')
+            % Exogenous variables are present in the 1-λ part (rule of thumb agents).
             tmp0 = 0;
             tmp1 = 0;
-            for i=1:length(equations.(eqtag).additive.params)
-                if isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
-                    tmp0 = tmp0 + equations.(eqtag).additive.scaling_factor(i);
-                elseif ~isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
-                    tmp0 = tmp0 + DynareModel.params(equations.(eqtag).additive.params(i))*equations.(eqtag).additive.scaling_factor(i);
-                elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && isnan(equations.(eqtag).additive.params(i))
-                    tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp{i};
-                elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && ~isnan(equations.(eqtag).additive.params(i))
-                    tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.params(i)*equations.(eqtag).additive.bgp{i};
+            for i=1:length(pacmodel.non_optimizing_behaviour.params)
+                if isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
+                    tmp0 = tmp0 + pacmodel.non_optimizing_behaviour.scaling_factor(i);
+                elseif ~isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
+                    tmp0 = tmp0 + DynareModel.params(pacmodel.non_optimizing_behaviour.params(i))*pacmodel.non_optimizing_behaviour.scaling_factor(i);
+                elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && isnan(pacmodel.non_optimizing_behaviour.params(i))
+                    tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.bgp{i};
+                elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && ~isnan(pacmodel.non_optimizing_behaviour.params(i))
+                    tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.params(i)*pacmodel.non_optimizing_behaviour.bgp{i};
                 end
             end
-            cc = cc - tmp0;
-            ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
+            cc = cc - (1.0-gamma)*tmp0;
+            ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
         end
-        DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
     end
-end
\ No newline at end of file
+    if isfield(pacmodel, 'additive')
+        % Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
+        tmp0 = 0;
+        tmp1 = 0;
+        for i=1:length(pacmodel.additive.params)
+            if isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
+                tmp0 = tmp0 + pacmodel.additive.scaling_factor(i);
+            elseif ~isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
+                tmp0 = tmp0 + DynareModel.params(pacmodel.additive.params(i))*pacmodel.additive.scaling_factor(i);
+            elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && isnan(pacmodel.additive.params(i))
+                tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.bgp{i};
+            elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && ~isnan(pacmodel.additive.params(i))
+                tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.params(i)*pacmodel.additive.bgp{i};
+            end
+        end
+        cc = cc - tmp0;
+        ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
+    end
+    DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
+end
diff --git a/matlab/+pac/initialize.m b/matlab/+pac/initialize.m
index 1a185c64b..c68086325 100644
--- a/matlab/+pac/initialize.m
+++ b/matlab/+pac/initialize.m
@@ -35,14 +35,12 @@ else
 end
 
 % Append default bgp fields if needed.
-for i=1:rows(M_.pac.(pacmodel).tag_map)
-    if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'additive')
-        M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.params));
-    end
-    if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'optim_additive')
-        M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.params));
-    end
-    if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'non_optimizing_behaviour')
-        M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.params));
-    end
+if isfield(M_.pac.(pacmodel), 'additive')
+    M_.pac.(pacmodel).additive.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).additive.params));
+end
+if isfield(M_.pac.(pacmodel), 'optim_additive')
+    M_.pac.(pacmodel).optim_additive.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).optim_additive.params));
+end
+if isfield(M_.pac.(pacmodel), 'non_optimizing_behaviour')
+    M_.pac.(pacmodel).non_optimizing_behaviour.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).non_optimizing_behaviour.params));
 end
\ No newline at end of file
diff --git a/matlab/cherrypick.m b/matlab/cherrypick.m
index 834592e5e..0f4fabed2 100644
--- a/matlab/cherrypick.m
+++ b/matlab/cherrypick.m
@@ -129,12 +129,12 @@ try
         isvar = regexp(RHS, 'var_expectation\(model_name = (?<name>\w+)\)', 'names');
         ispac = regexp(RHS, 'pac_expectation\(model_name = (?<name>\w+)\)', 'names');
         if ~isempty(isvar)
-            rhs = write_expectations(eqtags{i}, isvar.name, 'var');
+            rhs = write_expectations(isvar.name, 'var');
             lhs = sprintf('%s_VE', eqtags{i});
             RHS = strrep(RHS, sprintf('var_expectation(model_name = %s)', isvar.name), lhs);
         end
         if ~isempty(ispac)
-            [rhs, growthneutralitycorrection] = write_expectations(eqtags{i}, ispac.name, 'pac');
+            [rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac');
             if ~isempty(rhs)
                 lhs = sprintf('%s_PE', eqtags{i});
                 if isempty(growthneutralitycorrection)
diff --git a/matlab/pac-tools/+pac/+bgp/get.m b/matlab/pac-tools/+pac/+bgp/get.m
index 04037357c..3ab37a012 100644
--- a/matlab/pac-tools/+pac/+bgp/get.m
+++ b/matlab/pac-tools/+pac/+bgp/get.m
@@ -19,4 +19,4 @@ function dummy = get(pacmodel, paceq, kind, id)
 
 global M_
 
-dummy = M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{strcmp(paceq, M_.pac.(pacmodel).tag_map(:,1)),2}).(kind).bgp{id};
+dummy = M_.pac.(pacmodel).(kind).bgp{id};
diff --git a/matlab/pac-tools/+pac/+bgp/set.m b/matlab/pac-tools/+pac/+bgp/set.m
index 8a2283f73..fceb64876 100644
--- a/matlab/pac-tools/+pac/+bgp/set.m
+++ b/matlab/pac-tools/+pac/+bgp/set.m
@@ -1,4 +1,4 @@
-function dummy = set(pacmodel, paceq, variable, nonzeromean)
+function dummy = set(pacmodl, paceq, variable, nonzeromean)
 
 % Provide information about long run levels of exogenous variables in PAC equation.
 %
@@ -30,10 +30,10 @@ function dummy = set(pacmodel, paceq, variable, nonzeromean)
 
 global M_
 
-eqtag = M_.pac.(pacmodel).tag_map{strcmp(paceq, M_.pac.(pacmodel).tag_map(:,1)),2};
-
 dummy = [];
 
+pacmodel = M_.pac.(pacmodl);
+
 ide = find(strcmp(variable, M_.endo_names));
 xflag = false;
 if isempty(ide)
@@ -43,32 +43,32 @@ if isempty(ide)
 end
 
 if ~isempty(ide)
-    [M_, done] = writebgpfield('additive', pacmodel, eqtag, ide, xflag, nonzeromean, M_);
-    if done, return, end
-    [M_, done] = writebgpfield('optim_additive', pacmodel, eqtag, ide, xflag, nonzeromean, M_);
-    if done, return, end
-    [M_, done] = writebgpfield('non_optimizing_behaviour', pacmodel, eqtag, ide, xflag, nonzeromean, M_);
-    if done, return, end
+    [pacmodel, done] = writebgpfield('additive', pacmodel, ide, xflag, nonzeromean);
+    if done, M_.pac.(pacmodl) = pacmodel; return, end
+    [pacmodel, done] = writebgpfield('optim_additive', pacmodel, ide, xflag, nonzeromean);
+    if done, M_.pac.(pacmodl) = pacmodel; return, end
+    [pacmodel, done] = writebgpfield('non_optimizing_behaviour', pacmodel, ide, xflag, nonzeromean);
+    if done, M_.pac.(pacmodl) = pacmodel; return, end
     warning('%s is not an exogenous variable in equation %s.', variable, paceq)
 else
     error('Endogenous/Exogenous variable %s is unknown.', variable)
 end
 
 
-function [M_, done] = writebgpfield(type, pacmodel, eqtag, ide, xflag, nonzeromean, M_)
+function [pacmodel, done] = writebgpfield(type, pacmodel, ide, xflag, nonzeromean, M_)
 done = false;
-if isfield(M_.pac.(pacmodel).equations.(eqtag), type)
-    if ~isfield(M_.pac.(pacmodel).equations.(eqtag).additive, 'bgp')
-        M_.pac.(pacmodel).equations.(eqtag).(type).bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(eqtag).(type).params));
+if isfield(pacmodel, type)
+    if ~isfield(pacmodel.additive, 'bgp')
+        pacmodel.(type).bgp = repmat({false}, 1, length(pacmodel.(type).params));
     end
-    [isvar, ie] = ismember(ide, M_.pac.(pacmodel).equations.(eqtag).(type).vars);
+    [isvar, ie] = ismember(ide, pacmodel.(type).vars);
     if isvar
         if xflag
-            assert(~M_.pac.(pacmodel).equations.(eqtag).(type).isendo(ie), 'Variable type issue.')
+            assert(~pacmodel.(type).isendo(ie), 'Variable type issue.')
         else
-            assert(M_.pac.(pacmodel).equations.(eqtag).(type).isendo(ie), 'Variable type issue.')
+            assert(pacmodel.(type).isendo(ie), 'Variable type issue.')
         end
-        M_.pac.(pacmodel).equations.(eqtag).(type).bgp{ie} = nonzeromean;
+        pacmodel.(type).bgp{ie} = nonzeromean;
         done = true;
     end
-end
+end
\ No newline at end of file
diff --git a/matlab/pac-tools/geteqtag.m b/matlab/pac-tools/geteqtag.m
deleted file mode 100644
index 26a8eb589..000000000
--- a/matlab/pac-tools/geteqtag.m
+++ /dev/null
@@ -1,30 +0,0 @@
-function eqtag = geteqtag(eqname, pacname, Model)
-
-% Returns the equation tag, as defined by the preprocessor,
-% associated to a pac equation.
-%
-% INPUTS
-% - eqname     [char]    1×n array, name of the PAC equation.
-% - pacname    [char]    1×m array, name of the PAC model.
-%
-% OUTPUTS
-% - eqtag      [char]    equation tag associated to eqname.
-
-% Copyright © 2019 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-eqtag = Model.pac.(pacname).tag_map{strcmp(eqname, Model.pac.(pacname).tag_map(:,1)),2};
\ No newline at end of file
diff --git a/matlab/pac-tools/hVectors.m b/matlab/pac-tools/hVectors.m
index 7d44a47a2..8f9885dab 100644
--- a/matlab/pac-tools/hVectors.m
+++ b/matlab/pac-tools/hVectors.m
@@ -1,14 +1,15 @@
-function [h0, h1, longruncorrectionparameter] = hVectors(params, H, ids, idns, auxmodel)
+function [h, lrcp] = hVectors(params, H, auxmodel, kind, id)
 
 % INPUTS
 % - params          [double]     (m+1)*1 vector, PAC parameters (lag polynomial coefficients and discount factor).
 % - H               [double]     (np*np) matrix, companion matrix of the VAR(p) model.
-% - ids             [integer]    n*1 vector, selection of the stationary components of the VAR.
-% - idns            [integer]    n*1 vector, selection of the non stationary components of the VAR.
+% - auxmodel        [char]       kind of auxiliary model, possible values are: 'var' and 'trend_component'.
+% - kind            [char]       kind of expectation in PAC equation (See FRB/US doc), possible values are 'll', 'dl' and 'dd'.
+% - id              [integer]    scalar, index pointing to the target in the auxiliary model.
 %
 % OUTPUTS
-% - h0              [double]     1*n vector.
-% - h1              [double]     1*n vector.
+% - h               [double]     1*n vector of weights (used to compute the linear combination of the variables in the companion VAR representation of the auxiliary model).
+% - lrcp            [double]     scalar, parameter for the growth neutrality correction.
 %
 % REMARKS
 % The parameters are ordered as follows in the params vector:
@@ -17,7 +18,7 @@ function [h0, h1, longruncorrectionparameter] = hVectors(params, H, ids, idns, a
 %    params(2:end-1) ⟶ Autoregressive parameters.
 %    params(end)     ⟶ Discount factor.
 
-% Copyright (C) 2018 Dynare Team
+% Copyright © 2018-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -48,28 +49,24 @@ A_b = polyval(A, beta);
 m = length(alpha);
 n = length(H);
 
-tmp = eye(n*m)-kron(G, transpose(H));
+tmp = eye(n*m)-kron(G, transpose(H)); % inv(W2)
 
-if isempty(ids)
-    h0 = [];
-else
-    h0 = A_1*A_b*((kron(iota(m, m), H))'*(tmp\kron(iota(m, m), iota(n, ids))));
+switch kind
+  case 'll'
+    h = A_1*A_b*((kron(iota(m, m), H))'*(tmp\kron(iota(m, m), iota(n, id))));
+  case 'dd'
+    h = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), H')*(tmp\kron(iota(m, m), iota(n, id))));
+  case 'dl'
+    h = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), (H'-eye(length(H))))*(tmp\kron(iota(m, m), iota(n, id))));
+  otherwise
+    error('Unknown kind value in PAC model.')
 end
 
-if nargout>1
-    if isempty(idns)
-        h1 = [];
+if nargin>1
+    if isequal(kind, 'll')
+        lrcp = NaN;
     else
-        switch auxmodel
-          case {'var', 'trend_component'}
-            h1 = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), H')*(tmp\kron(iota(m, m), iota(n, idns))));
-          case 'Need to check in which case we should trigger this one...'
-            h1 = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), (H'-eye(length(H))))*(tmp\kron(iota(m, m), iota(n, idns))));
-        end
+        d = A_1*A_b*(iota(m, m)'*inv((eye(m)-G)*(eye(m)-G))*iota(m, m));
+        lrcp = (1-sum(params(2:end-1))-d);
     end
-end
-
-if nargout>2
-    d = A_1*A_b*(iota(m, m)'*inv((eye(m)-G)*(eye(m)-G))*iota(m, m));
-    longruncorrectionparameter = (1-sum(params(2:end-1))-d);
 end
\ No newline at end of file
diff --git a/matlab/print_equations.m b/matlab/print_equations.m
index 7d8069da2..7fb50b561 100644
--- a/matlab/print_equations.m
+++ b/matlab/print_equations.m
@@ -9,7 +9,7 @@ function str = print_equations(variable_name, withexpansion)
 % OUTPUTS
 % None
 
-% Copyright (C) 2019 Dynare Team
+% Copyright © 2019-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -61,14 +61,26 @@ model = jsonfile.model;
 
 % Print the equations the variable appears in.
 for it = 1:length(M_.mapping.(variable_name).eqidx)
-    rhs = model{M_.mapping.(variable_name).eqidx(it)}.rhs;
+    if M_.mapping.(variable_name).eqidx(it)>length(model)
+        % Equations appended by the preprocessor for auxiliary variables are not displayed, except if it is an aux variable
+        % for the PAC expectation.
+        % TODO Probably should do the same for VAR expectation.
+        if isfield(M_, 'lhs') && isequal(M_.lhs{M_.mapping.(variable_name).eqidx(it)}(1:16), 'pac_expectation_')
+            id = M_.mapping.(M_.lhs{M_.mapping.(variable_name).eqidx(it)}).eqidx(1);
+        else
+            continue
+        end
+    else
+        id = M_.mapping.(variable_name).eqidx(it);
+    end
+    rhs = model{id}.rhs;
     if withexpansion
         if isfield(M_, 'pac') && contains(rhs, 'pac_expectation')
             % Get the index of the equation's PAC model.
             models = fieldnames(M_.pac);
             idx = find(~cellfun('isempty',cellfun(@(s)find(contains(rhs,s)),models,'uni',0)));
             % Get the expanded PAC_EXPECTATION term.
-            [pac_expression, growthneutralitycorrection] = write_expectations(M_.pac.(models{idx}).tag_map(:,1), models{idx}, 'pac', true);
+            [pac_expression, growthneutralitycorrection] = write_expectations(models{idx}, 'pac', true);
             expression = [sprintf('\n\t + %s', growthneutralitycorrection) TransformExpandedExpr(pac_expression)];
             rhs = strrep(rhs, ['+pac_expectation(model_name = ' models{idx} ')'], expression);
         elseif isfield(M_, 'var_expectation') && contains(rhs, 'var_expectation')
@@ -76,7 +88,7 @@ for it = 1:length(M_.mapping.(variable_name).eqidx)
             models = fieldnames(M_.var_expectation);
             idx = find(~cellfun('isempty',cellfun(@(s)find(contains(rhs,s)),models,'uni',0)));
             % Get the expanded VAR_EXPECTATION term.
-            expression = write_expectations('fake', models{idx}, 'var', true);
+            expression = write_expectations(models{idx}, 'var', true);
             expression = TransformExpandedExpr(expression);
             rhs = strrep(rhs, ['+var_expectation(model_name = ' models{idx} ')'], expression);
         elseif ~isfield(M_, 'pac') && ~isfield(M_, 'var_expectation')
@@ -87,7 +99,7 @@ for it = 1:length(M_.mapping.(variable_name).eqidx)
     if nargout
         str = sprintf('%s = %s;\n', model{M_.mapping.(variable_name).eqidx(it)}.lhs, rhs);
     end
-    fprintf('%s = %s;\n', model{M_.mapping.(variable_name).eqidx(it)}.lhs, rhs);
+    fprintf('%s = %s;\n', model{id}.lhs, rhs);
 end
 
 function [transformed_expression] = TransformExpandedExpr(expression)
diff --git a/matlab/print_expectations.m b/matlab/print_expectations.m
index 3ce7c923c..783baa05f 100644
--- a/matlab/print_expectations.m
+++ b/matlab/print_expectations.m
@@ -84,11 +84,6 @@ switch expectationmodelkind
   otherwise
 end
 
-if isequal(expectationmodelkind, 'pac')
-    % Get the equation tag (in M_.pac.(pacmodl).equations)
-    eqtag = M_.pac.(expectationmodelname).tag_map{strcmp(M_.pac.(expectationmodelname).tag_map(:,1), eqname),2};
-end
-
 % Get the expectation model description
 expectationmodel = M_.(expectationmodelfield).(expectationmodelname);
 
@@ -114,11 +109,7 @@ if ~exist(sprintf('%s/model/%s', M_.fname, [expectationmodelkind '-expectations'
     mkdir(sprintf('%s/model/%s', M_.fname, [expectationmodelkind '-expectations']))
 end
 
-if isequal(expectationmodelkind, 'pac')
-    filename = sprintf('%s/model/%s/%s-%s-parameters.inc', M_.fname, [expectationmodelkind '-expectations'], eqtag, expectationmodelname);
-else
-    filename = sprintf('%s/model/%s/%s-parameters.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
-end
+filename = sprintf('%s/model/%s/%s-parameters.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
 fid = fopen(filename, 'w');
 fprintf(fid, '// This file has been generated by dynare (%s).\n\n', datestr(now));
 
@@ -135,28 +126,14 @@ switch expectationmodelkind
         end
     end
   case 'pac'
-    if ~isempty(expectationmodel.equations.(eqtag).h0_param_indices)
-        parameter_declaration = 'parameters';
-        for i=1:length(expectationmodel.equations.(eqtag).h0_param_indices)
-            parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.equations.(eqtag).h0_param_indices(i)});
-        end
-        fprintf(fid, '%s;\n\n', parameter_declaration);
-        if withcalibration
-            for i=1:length(expectationmodel.equations.(eqtag).h0_param_indices)
-                fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.equations.(eqtag).h0_param_indices(i)}, M_.params(expectationmodel.equations.(eqtag).h0_param_indices(i)));
-            end
-        end
+    parameter_declaration = 'parameters';
+    for i=1:length(expectationmodel.h_param_indices)
+        parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.h_param_indices(i)});
     end
-    if ~isempty(expectationmodel.equations.(eqtag).h1_param_indices)
-        parameter_declaration = 'parameters';
-        for i=1:length(expectationmodel.equations.(eqtag).h1_param_indices)
-            parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.equations.(eqtag).h1_param_indices(i)});
-        end
-        fprintf(fid, '%s;\n\n', parameter_declaration);
-        if withcalibration
-            for i=1:length(expectationmodel.equations.(eqtag).h1_param_indices)
-                fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.equations.(eqtag).h1_param_indices(i)}, M_.params(expectationmodel.equations.(eqtag).h1_param_indices(i)));
-            end
+    fprintf(fid, '%s;\n\n', parameter_declaration);
+    if withcalibration
+        for i=1:length(expectationmodel.h_param_indices)
+            fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.h_param_indices(i)}, M_.params(expectationmodel.h_param_indices(i)));
         end
     end
     if isfield(expectationmodel, 'growth_neutrality_param_index')
@@ -181,19 +158,15 @@ fprintf('Parameters declarations and calibrations are saved in %s.\n', filename)
 % Second print the expanded VAR_EXPECTATION/PAC_EXPECTATION term.
 %
 
-if isequal(expectationmodelkind, 'pac')
-    filename = sprintf('%s/model/%s/%s-%s-expression.inc', M_.fname, [expectationmodelkind '-expectations'], eqtag, expectationmodelname);
-else
-    filename = sprintf('%s/model/%s/%s-expression.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
-end
+filename = sprintf('%s/model/%s/%s-expression.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
 fid = fopen(filename, 'w');
 fprintf(fid, '// This file has been generated by dynare (%s).\n', datestr(now));
 
 switch expectationmodelkind
   case 'var'
-    expression = write_expectations(eqname, expectationmodelname, expectationmodelkind, true);
+    expression = write_expectations(expectationmodelname, expectationmodelkind, true);
   case 'pac'
-    [expression, growthneutralitycorrection] = write_expectations(eqname, expectationmodelname, expectationmodelkind, true);
+    [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, true);
 end
 
 fprintf(fid, '%s', expression);
@@ -206,7 +179,7 @@ fprintf('Expectation unrolled expression is saved in %s.\n', filename);
 %
 
 if isequal(expectationmodelkind, 'pac') && growth_correction
-    filename = sprintf('%s/model/%s/%s-%s-growth-neutrality-correction.inc', M_.fname, [expectationmodelkind '-expectations'], eqtag, expectationmodelname);
+    filename = sprintf('%s/model/%s/%s-growth-neutrality-correction.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
     fid = fopen(filename, 'w');
     fprintf(fid, '// This file has been generated by dynare (%s).\n', datestr(now));
     fprintf(fid, '%s', growthneutralitycorrection);
@@ -219,22 +192,11 @@ end
 %
 kind = [expectationmodelkind '_expectations'];
 mkdir(sprintf('+%s/+%s/+%s', M_.fname, kind, expectationmodelname));
-if isequal(expectationmodelkind, 'pac')
-    filename = sprintf('+%s/+%s/+%s/%s_evaluate.m', M_.fname, kind, expectationmodelname, eqtag);
-else
-    filename = sprintf('+%s/+%s/+%s/evaluate.m', M_.fname, kind, expectationmodelname);
-end
+filename = sprintf('+%s/+%s/+%s/evaluate.m', M_.fname, kind, expectationmodelname);
 fid = fopen(filename, 'w');
-if isequal(expectationmodelkind, 'pac')
-    fprintf(fid, 'function ds = %s_evaluate(dbase)\n\n', eqtag);
-else
-    fprintf(fid, 'function ds = evaluate(dbase)\n\n');
-end
-if isequal(expectationmodelkind, 'pac')
-    fprintf(fid, '%% Evaluates %s term (%s in %s).\n', kind, expectationmodelname, eqname);
-else
-    fprintf(fid, '%% Evaluates %s term (%s).\n', kind, expectationmodelname);
-end
+
+fprintf(fid, 'function ds = evaluate(dbase)\n\n');
+fprintf(fid, '%% Evaluates %s term (%s).\n', kind, expectationmodelname);
 fprintf(fid, '%%\n');
 fprintf(fid, '%% INPUTS\n');
 fprintf(fid, '%% - dbase     [dseries]  databse containing all the variables appearing in the auxiliary model for the expectation.\n');
@@ -266,8 +228,7 @@ end
 
 if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
     id = id+1;
-    expression = sprintf('%1.16f', M_.params(expectationmodel.equations.(eqtag).h0_param_indices(id))+ ...
-                         M_.params(expectationmodel.equations.(eqtag).h1_param_indices(id)));
+    expression = sprintf('%1.16f', M_.params(expectationmodel.h_param_indices(id)));
 end
 
 for i=1:maxlag
@@ -295,17 +256,7 @@ for i=1:maxlag
           case 'var'
             parameter = M_.params(expectationmodel.param_indices(id));
           case 'pac'
-            parameter = 0;
-            if ~isempty(expectationmodel.equations.(eqtag).h0_param_indices)
-                parameter = M_.params(expectationmodel.equations.(eqtag).h0_param_indices(id));
-            end
-            if ~isempty(expectationmodel.equations.(eqtag).h1_param_indices)
-                if ~parameter
-                    parameter = M_.params(expectationmodel.equations.(eqtag).h1_param_indices(id));
-                else
-                    parameter = parameter+M_.params(expectationmodel.equations.(eqtag).h1_param_indices(id));
-                end
-            end
+            parameter = M_.params(expectationmodel.h_param_indices(id));
           otherwise
         end
         switch expectationmodelkind
@@ -387,4 +338,4 @@ fprintf('Expectation dseries expression is saved in %s.\n', filename);
 
 skipline();
 
-rehash
+rehash
\ No newline at end of file
diff --git a/matlab/write_expectations.m b/matlab/write_expectations.m
index 385098931..5312e42a6 100644
--- a/matlab/write_expectations.m
+++ b/matlab/write_expectations.m
@@ -1,9 +1,8 @@
-function [expression, growthneutralitycorrection] = write_expectations(eqname, expectationmodelname, expectationmodelkind, iscrlf)
+function [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, iscrlf)
 
 % Prints the exansion of the VAR_EXPECTATION or PAC_EXPECTATION term in files.
 %
 % INPUTS
-% - eqname                      [string]    Name of the equation.
 % - epxpectationmodelname       [string]    Name of the expectation model.
 % - expectationmodelkind        [string]    Kind of the expectation model ('var' or 'pac').
 % - iscrlf                      [string]    Adds carriage return after each additive term if true.
@@ -36,8 +35,6 @@ if ismember(expectationmodelkind, {'var', 'pac'})
         expectationmodelfield = 'var_expectation';
     else
         expectationmodelfield = 'pac';
-        % Get the equation tag (in M_.pac.(pacmodl).equations)
-        eqtag = M_.pac.(expectationmodelname).tag_map{strcmp(M_.pac.(expectationmodelname).tag_map(:,1), eqname),2};
     end
 else
     error('Value of third input argument must be ''var'' or ''pac''.')
@@ -49,7 +46,7 @@ if nargout>1 && isequal(expectationmodelkind, 'var')
     error('Cannot return more than one argument if the expectation model is a VAR.')
 end
 
-if nargin<4
+if nargin<3
     iscrlf = false;
 end
 
@@ -85,8 +82,7 @@ end
 
 if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
     id = id+1;
-    expression = sprintf('%s+%s', M_.param_names{expectationmodel.equations.(eqtag).h0_param_indices(id)}, ...
-                         M_.param_names{expectationmodel.equations.(eqtag).h1_param_indices(id)});
+    expression = sprintf('%s', M_.param_names{expectationmodel.h_param_indices(id)});
 end
 
 for i=1:maxlag
@@ -114,17 +110,7 @@ for i=1:maxlag
           case 'var'
             parameter = M_.param_names{expectationmodel.param_indices(id)};
           case 'pac'
-            parameter = '';
-            if ~isempty(expectationmodel.equations.(eqtag).h0_param_indices)
-                parameter = M_.param_names{expectationmodel.equations.(eqtag).h0_param_indices(id)};
-            end
-            if ~isempty(expectationmodel.equations.(eqtag).h1_param_indices)
-                if isempty(parameter)
-                    parameter = M_.param_names{expectationmodel.equations.(eqtag).h1_param_indices(id)};
-                else
-                    parameter = sprintf('(%s+%s)', parameter, M_.param_names{expectationmodel.equations.(eqtag).h1_param_indices(id)});
-                end
-            end
+            parameter = M_.param_names{expectationmodel.h_param_indices(id)};
           otherwise
         end
         switch expectationmodelkind
diff --git a/preprocessor b/preprocessor
index 1cc512962..a210a8fd5 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 1cc512962c292b7cafecfb0550ecd7373e80fd9b
+Subproject commit a210a8fd59a052abf2c67d264f0889ba2213043d
diff --git a/tests/pac/trend-component-13a/example1.mod b/tests/pac/trend-component-13a/example1.mod
index da102a915..b980a5220 100644
--- a/tests/pac/trend-component-13a/example1.mod
+++ b/tests/pac/trend-component-13a/example1.mod
@@ -40,10 +40,10 @@ model;
 y = rho_1*y(-1) + rho_2*y(-2) + ey;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -69,6 +69,6 @@ end;
 // Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
 pac.initialize('pacman');
 
-if ~isequal(M_.pac.pacman.equations.(M_.pac.pacman.tag_map{strcmp(M_.pac.pacman.tag_map(:,1), 'zpac'),2}).ec.istarget, [true, false])
+if ~isequal(M_.pac.pacman.ec.istarget, [true, false])
    error('ec.istarget vector is wrong.')
 end
diff --git a/tests/pac/trend-component-13b/example1.mod b/tests/pac/trend-component-13b/example1.mod
index da102a915..b980a5220 100644
--- a/tests/pac/trend-component-13b/example1.mod
+++ b/tests/pac/trend-component-13b/example1.mod
@@ -40,10 +40,10 @@ model;
 y = rho_1*y(-1) + rho_2*y(-2) + ey;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -69,6 +69,6 @@ end;
 // Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
 pac.initialize('pacman');
 
-if ~isequal(M_.pac.pacman.equations.(M_.pac.pacman.tag_map{strcmp(M_.pac.pacman.tag_map(:,1), 'zpac'),2}).ec.istarget, [true, false])
+if ~isequal(M_.pac.pacman.ec.istarget, [true, false])
    error('ec.istarget vector is wrong.')
 end
diff --git a/tests/pac/trend-component-14/substitution.mod b/tests/pac/trend-component-14/substitution.mod
index 6761227e2..466807ad1 100644
--- a/tests/pac/trend-component-14/substitution.mod
+++ b/tests/pac/trend-component-14/substitution.mod
@@ -30,7 +30,7 @@ e_c_m =  .5;
 c_z_1 =  .2;
 c_z_2 = -.1;
 
-@#include "example1/model/pac-expectations/eq0-pacman-parameters.inc"
+@#include "example1/model/pac-expectations/pacman-parameters.inc"
 
 model;
 
@@ -51,7 +51,7 @@ x2bar = x2bar(-1) + ex2bar;
 
 [name='zpac']
 diff(z) = e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) +
-@#include "example1/model/pac-expectations/eq0-pacman-expression.inc"
+@#include "example1/model/pac-expectations/pacman-expression.inc"
 + ez;
 
 end;
@@ -95,4 +95,4 @@ if ~isequal(length(g1(:)), length(example1.g1(:))) || max(abs(example1.g1(:)-g1(
   warning('Jacobian matrices do not match!')
 end
 
-delete('example1.mat');
\ No newline at end of file
+delete('example1.mat');
diff --git a/tests/pac/trend-component-19-growth-lin-comb/example1.mod b/tests/pac/trend-component-19-growth-lin-comb/example1.mod
index daed6065c..5d98a310d 100644
--- a/tests/pac/trend-component-19-growth-lin-comb/example1.mod
+++ b/tests/pac/trend-component-19-growth-lin-comb/example1.mod
@@ -43,10 +43,10 @@ gg = g;
 y = rho_1*y(-1) + rho_2*y(-2) + ey;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -72,7 +72,7 @@ initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M
 TrueData = simul_backward_model(initialconditions, 5000);
 
 // NLS estimation
-// Define a structure describing the parameters to be estimated (with initial conditions). 
+// Define a structure describing the parameters to be estimated (with initial conditions).
 clear eparams
 eparams.e_c_m = .5;
 eparams.c_z_1 = .2;
@@ -85,7 +85,7 @@ pac.estimate.nls('zpac', eparams, edata, 2005Q1:2000Q1+200, 'lsqnonlin');
 e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
 c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
 c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
-resid_nls = oo_.pac.pacman.equations.eq0.residual;
+resid_nls = oo_.pac.pacman.residual;
 
 fprintf('Estimate of e_c_m: %f \n', e_c_m_nls)
 fprintf('Estimate of c_z_1: %f \n', c_z_1_nls)
@@ -94,7 +94,7 @@ fprintf('Estimate of c_z_2: %f \n', c_z_2_nls)
 skipline(2)
 
 // Iterative OLS estimation using estimates from NLS
-// Define a structure describing the parameters to be estimated (with initial conditions). 
+// Define a structure describing the parameters to be estimated (with initial conditions).
 clear eparams
 eparams.e_c_m = e_c_m_nls;
 eparams.c_z_1 = c_z_1_nls;
@@ -120,7 +120,7 @@ fprintf('\n')
 e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
 c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
 c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
-resid_iterative_ols = oo_.pac.pacman.equations.eq0.residual;
+resid_iterative_ols = oo_.pac.pacman.residual;
 
 fprintf('Estimate of e_c_m: %f \n', e_c_m_iterative_ols)
 fprintf('Estimate of c_z_1: %f \n', c_z_1_iterative_ols)
diff --git a/tests/pac/trend-component-20-1/example1.mod b/tests/pac/trend-component-20-1/example1.mod
index a0a6daaf2..eae63c7a8 100644
--- a/tests/pac/trend-component-20-1/example1.mod
+++ b/tests/pac/trend-component-20-1/example1.mod
@@ -51,10 +51,10 @@ y = rho_1*y(-1) + rho_2*y(-2) + ey;
 x = rho_3*x(-1) + rho_4*x(-2) + ex;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -67,6 +67,6 @@ diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2))
 
 end;
 
-if ~isequal(M_.endo_names(M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars), {'y'; 'x'})
+if ~isequal(M_.endo_names(M_.pac.pacman.non_optimizing_behaviour.vars), {'y'; 'x'})
    error('PAC non_optimizing_behaviour.vars field is wrong.')
-end
\ No newline at end of file
+end
diff --git a/tests/pac/trend-component-20-2/example1.mod b/tests/pac/trend-component-20-2/example1.mod
index 7406c04e0..64e16a4c2 100644
--- a/tests/pac/trend-component-20-2/example1.mod
+++ b/tests/pac/trend-component-20-2/example1.mod
@@ -51,10 +51,10 @@ y = rho_1*y(-1) + rho_2*y(-2) + ey;
 x = rho_3*x(-1) + rho_4*x(-2) + ex;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -67,6 +67,6 @@ diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2))
 
 end;
 
-if ~isequal(M_.endo_names(M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars), {'y'; 'x'})
+if ~isequal(M_.endo_names(M_.pac.pacman.non_optimizing_behaviour.vars), {'y'; 'x'})
    error('PAC non_optimizing_behaviour.vars field is wrong.')
-end
\ No newline at end of file
+end
diff --git a/tests/pac/trend-component-20-3/example1.mod b/tests/pac/trend-component-20-3/example1.mod
index 051470fdc..f7e6eb065 100644
--- a/tests/pac/trend-component-20-3/example1.mod
+++ b/tests/pac/trend-component-20-3/example1.mod
@@ -45,13 +45,13 @@ pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
 model;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:y']
 y = rho_1*y(-1) + rho_2*y(-2) + ey;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -67,6 +67,6 @@ diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2))
 
 end;
 
-if ~isequal(M_.endo_names(M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars), {'y'; 'x'})
+if ~isequal(M_.endo_names(M_.pac.pacman.non_optimizing_behaviour.vars), {'y'; 'x'})
    error('PAC non_optimizing_behaviour.vars field is wrong.')
-end
\ No newline at end of file
+end
diff --git a/tests/pac/trend-component-20-4/example1.mod b/tests/pac/trend-component-20-4/example1.mod
index 284d88b2c..aa85608ac 100644
--- a/tests/pac/trend-component-20-4/example1.mod
+++ b/tests/pac/trend-component-20-4/example1.mod
@@ -51,10 +51,10 @@ y = rho_1*y(-1) + rho_2*y(-2) + ey;
 x = rho_3*x(-1) + rho_4*x(-2) + ex;
 
 [name='eq:x1']
-diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;     
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
 
 [name='eq:x2']
-diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;     
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
 
 [name='eq:x1bar']
 x1bar = x1bar(-1) + ex1bar;
@@ -67,6 +67,6 @@ diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2))
 
 end;
 
-if ~isequal([M_.endo_names(M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars(M_.pac.pacman.equations.eq0.non_optimizing_behaviour.isendo)); M_.exo_names(M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars(~M_.pac.pacman.equations.eq0.non_optimizing_behaviour.isendo))], {'y'; 'x'; 'tt'})
+if ~isequal([M_.endo_names(M_.pac.pacman.non_optimizing_behaviour.vars(M_.pac.pacman.non_optimizing_behaviour.isendo)); M_.exo_names(M_.pac.pacman.non_optimizing_behaviour.vars(~M_.pac.pacman.non_optimizing_behaviour.isendo))], {'y'; 'x'; 'tt'})
    error('PAC non_optimizing_behaviour.vars and/or non_optimizing_behaviour.isendo fields are wrong.')
-end
\ No newline at end of file
+end
diff --git a/tests/pac/trend-component-26/example1.mod b/tests/pac/trend-component-26/example1.mod
index 2ad3694d0..ae8f3ad18 100644
--- a/tests/pac/trend-component-26/example1.mod
+++ b/tests/pac/trend-component-26/example1.mod
@@ -109,109 +109,109 @@ pac.bgp.set('pacman', 'zpac', 'dx2', .0);
 pac.update.expectation('pacman');
 
 id = find(strcmp('s', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+id = find(id==M_.pac.pacman.optim_additive.vars);
 if ~pac.bgp.get('pacman', 'zpac', 'optim_additive', id)
     error('bgp field in optim_additive for variable s is wrong.')
 end
 
 id = find(strcmp('s', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+id = find(id==M_.pac.pacman.additive.vars);
 if ~isempty(id)
-    error('Variable s should not be under M_.pac.pacman.equations.eq0.additive.')
+    error('Variable s should not be under M_.pac.pacman.additive.')
 end
 
 id = find(strcmp('s', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+id = find(id==M_.pac.pacman.non_optimizing_behaviour.vars);
 if ~isempty(id)
-    error('Variable s should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+    error('Variable s should not be under M_.pac.pacman.non_optimizing_behaviour.')
 end
 
 id = find(strcmp('dv', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+id = find(id==M_.pac.pacman.optim_additive.vars);
 if pac.bgp.get('pacman', 'zpac', 'optim_additive', id)
    error('bgp field in optim_additive for variable dv is wrong.')
 end
 
 id = find(strcmp('dv', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+id = find(id==M_.pac.pacman.additive.vars);
 if ~isempty(id)
-    error('Variable dv should not be under M_.pac.pacman.equations.eq0.additive.')
+    error('Variable dv should not be under M_.pac.pacman.additive.')
 end
 
 id = find(strcmp('dv', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+id = find(id==M_.pac.pacman.non_optimizing_behaviour.vars);
 if ~isempty(id)
-    error('Variable dv should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+    error('Variable dv should not be under M_.pac.pacman.non_optimizing_behaviour.')
 end
 
 id = find(strcmp('x', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+id = find(id==M_.pac.pacman.non_optimizing_behaviour.vars);
 if ~(abs(pac.bgp.get('pacman', 'zpac', 'non_optimizing_behaviour', id))<1e-12)
     error('bgp field in non_optimizing_behaviour for variable x is wrong.')
 end
 
 id = find(strcmp('x', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+id = find(id==M_.pac.pacman.optim_additive.vars);
 if ~isempty(id)
-    error('Variable x should not be under M_.pac.pacman.equations.eq0.optim_additive.')
+    error('Variable x should not be under M_.pac.pacman.optim_additive.')
 end
 
 id = find(strcmp('x', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+id = find(id==M_.pac.pacman.additive.vars);
 if ~isempty(id)
-    error('Variable x should not be under M_.pac.pacman.equations.eq0.additive.')
+    error('Variable x should not be under M_.pac.pacman.additive.')
 end
 
 id = find(strcmp('y', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+id = find(id==M_.pac.pacman.non_optimizing_behaviour.vars);
 if ~islogical(pac.bgp.get('pacman', 'zpac', 'non_optimizing_behaviour', id)) || pac.bgp.get('pacman', 'zpac', 'non_optimizing_behaviour', id)
    error('bgp field in non_optimizing_behaviour for variable y is wrong.')
 end
 
 id = find(strcmp('y', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+id = find(id==M_.pac.pacman.optim_additive.vars);
 if ~isempty(id)
-    error('Variable y should not be under M_.pac.pacman.equations.eq0.optim_additive.')
+    error('Variable y should not be under M_.pac.pacman.optim_additive.')
 end
 
 id = find(strcmp('y', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+id = find(id==M_.pac.pacman.additive.vars);
 if ~isempty(id)
-    error('Variable y should not be under M_.pac.pacman.equations.eq0.additive.')
+    error('Variable y should not be under M_.pac.pacman.additive.')
 end
 
 id = find(strcmp('dx2', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+id = find(id==M_.pac.pacman.additive.vars);
 if ~(abs(pac.bgp.get('pacman', 'zpac', 'additive', id))<1e-12)
     error('bgp field in additive is wrong.')
 end
 
 id = find(strcmp('dx2', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+id = find(id==M_.pac.pacman.non_optimizing_behaviour.vars);
 if ~isempty(id)
-    error('Variable y should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+    error('Variable y should not be under M_.pac.pacman.non_optimizing_behaviour.')
 end
 
 id = find(strcmp('dx2', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+id = find(id==M_.pac.pacman.optim_additive.vars);
 if ~isempty(id)
-    error('Variable dx2 should not be under M_.pac.pacman.equations.eq0.optim_additive.')
+    error('Variable dx2 should not be under M_.pac.pacman.optim_additive.')
 end
 
 id = find(strcmp('u', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+id = find(id==M_.pac.pacman.additive.vars);
 if ~islogical(pac.bgp.get('pacman', 'zpac', 'additive', id)) || pac.bgp.get('pacman', 'zpac', 'additive', id)
    error('bgp field in additive for variable u is wrong.')
 end
 
 id = find(strcmp('u', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+id = find(id==M_.pac.pacman.non_optimizing_behaviour.vars);
 if ~isempty(id)
-    error('Variable u should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+    error('Variable u should not be under M_.pac.pacman.non_optimizing_behaviour.')
 end
 
 id = find(strcmp('u', M_.endo_names));
-id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+id = find(id==M_.pac.pacman.optim_additive.vars);
 if ~isempty(id)
-    error('Variable u should not be under M_.pac.pacman.equations.eq0.optim_additive.')
+    error('Variable u should not be under M_.pac.pacman.optim_additive.')
 end
diff --git a/tests/pac/trend-component-28/example4.mod b/tests/pac/trend-component-28/example4.mod
index 5fced3966..b73c4720c 100644
--- a/tests/pac/trend-component-28/example4.mod
+++ b/tests/pac/trend-component-28/example4.mod
@@ -20,8 +20,7 @@ parameters
        a_x1_0 a_x1_1 a_x1_2 a_x1_x2_1 a_x1_x2_2
 	   a_x2_0 a_x2_1 a_x2_2 a_x2_x1_1 a_x2_x1_2
 	   e_c_m c_z_1 c_z_2 c_z_dx2 c_z_u c_z_dv c_z_s cx cy beta
-       lambda
-       px3;
+       lambda;
 
 rho_1 =  .9;
 rho_2 = -.2;
@@ -48,15 +47,13 @@ c_z_2 = -.1;
 c_z_dx2 = .3;
 c_z_u = .3;
 c_z_dv = .4;
-c_z_s  = -.2; 
+c_z_s  = -.2;
 cx = 1.0;
 cy = 1.0;
 
 
 lambda = 0.5; // Share of optimizing agents.
 
-px3 = -.1;
-
 trend_component_model(model_name=toto, eqtags=['eq:x1', 'eq:x2', 'eq:x1bar', 'eq:x2bar'], targets=['eq:x1bar', 'eq:x2bar']);
 
 pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
diff --git a/tests/pac/var-5/substitution.mod b/tests/pac/var-5/substitution.mod
index 8c2924b8f..29d837f70 100644
--- a/tests/pac/var-5/substitution.mod
+++ b/tests/pac/var-5/substitution.mod
@@ -21,7 +21,7 @@ e_c_m =  .1;
 c_z_1 =  .7;
 c_z_2 = -.3;
 
-@#include "example1/model/pac-expectations/eq0-pacman-parameters.inc"
+@#include "example1/model/pac-expectations/pacman-parameters.inc"
 
 model;
 
@@ -33,9 +33,9 @@ diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + g*(1-b_x_2)  + ex ;
 
 [name='eq:pac']
 diff(z) = e_c_m*(x(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) +
-@#include "example1/model/pac-expectations/eq0-pacman-growth-neutrality-correction.inc"
+@#include "example1/model/pac-expectations/pacman-growth-neutrality-correction.inc"
 +
-@#include "example1/model/pac-expectations/eq0-pacman-expression.inc"
+@#include "example1/model/pac-expectations/pacman-expression.inc"
 + ez;
 
 end;
diff --git a/tests/pac/var-6/substitution.mod b/tests/pac/var-6/substitution.mod
index 13978e802..f5f1c665a 100644
--- a/tests/pac/var-6/substitution.mod
+++ b/tests/pac/var-6/substitution.mod
@@ -21,7 +21,7 @@ e_c_m =  .1;
 c_z_1 =  .7;
 c_z_2 = -.3;
 
-@#include "example1/model/pac-expectations/eq0-pacman-parameters.inc"
+@#include "example1/model/pac-expectations/pacman-parameters.inc"
 
 model;
 
@@ -36,9 +36,9 @@ model;
 
   [name='eq:pac']
   diff(z) = e_c_m*(x(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) +
-  @#include "example1/model/pac-expectations/eq0-pacman-growth-neutrality-correction.inc"
+  @#include "example1/model/pac-expectations/pacman-growth-neutrality-correction.inc"
   +
-  @#include "example1/model/pac-expectations/eq0-pacman-expression.inc"
+  @#include "example1/model/pac-expectations/pacman-expression.inc"
   + ez;
 
 end;
diff --git a/tests/pac/var-7/substitution.mod b/tests/pac/var-7/substitution.mod
index 4597649c0..1e906fe72 100644
--- a/tests/pac/var-7/substitution.mod
+++ b/tests/pac/var-7/substitution.mod
@@ -24,7 +24,7 @@ c_z_2 = -.3;
 
 g = .1;
 
-@#include "example1/model/pac-expectations/eq0-pacman-parameters.inc"
+@#include "example1/model/pac-expectations/pacman-parameters.inc"
 
 model;
 
@@ -39,9 +39,9 @@ model;
 
   [name='eq:pac']
   diff(z) = e_c_m*(x(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) +
-  @#include "example1/model/pac-expectations/eq0-pacman-growth-neutrality-correction.inc"
+  @#include "example1/model/pac-expectations/pacman-growth-neutrality-correction.inc"
   +
-  @#include "example1/model/pac-expectations/eq0-pacman-expression.inc"
+  @#include "example1/model/pac-expectations/pacman-expression.inc"
   + ez;
 
 end;