diff --git a/matlab/+pac/+mce/parameters.m b/matlab/+pac/+mce/parameters.m
index 9c51cbfc792dace250007eb73e6388eb776940f6..69c5b2189f972740520a062ae5183744dcfbdc33 100644
--- a/matlab/+pac/+mce/parameters.m
+++ b/matlab/+pac/+mce/parameters.m
@@ -76,11 +76,72 @@ for e=1:number_of_pac_eq
     [G, alpha, beta] = buildGmatrixWithAlphaAndBeta(params);
     M_.params(pac_equation.mce.alpha) = 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));
-        M_.params(pacmodel.growth_neutrality_param_index) = 1-sum(params(2:end-1))-d;
+        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
+            equations.(eqtag).('non_optimizing_behaviour')
+            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.
+            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.
+            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};
+                end
+            end
+            cc
+            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
 end
\ No newline at end of file
diff --git a/matlab/+pac/+update/parameters.m b/matlab/+pac/+update/parameters.m
index 226ddfc1f9efdc082d0874ae59cbbe87bb5b0c05..4f50f06cc3ff3ed09a01c1f35581523066c83b9e 100644
--- a/matlab/+pac/+update/parameters.m
+++ b/matlab/+pac/+update/parameters.m
@@ -173,7 +173,7 @@ for e=1:number_of_pac_eq
         if ~isempty(equations.(eqtag).h1_param_indices)
             DynareModel.params(equations.(eqtag).h1_param_indices) = .0;
         end
-    end 
+    end
     % Update the parameter related to the growth neutrality correction.
     if isfield(equations.(eqtag), 'non_optimizing_behaviour')
         gamma = DynareModel.params(equations.(eqtag).share_of_optimizing_agents_index);
@@ -185,46 +185,61 @@ for e=1:number_of_pac_eq
         % 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-gamma*gg;          % 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.
+        cc = 1.0/gamma-gg;          % 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).
-            tmp = 0;
+            tmp0 = 0;
             for i=1:length(equations.(eqtag).optim_additive.params)
-                if isnan(equations.(eqtag).optim_additive.params(i)) && equations.(eqtag).optim_additive.bgp(i)
-                    tmp = tmp + equations.(eqtag).optim_additive.scaling_factor(i)*equations.(eqtag).optim_additive.bgp(i);
-                elseif ~isnan(equations.(eqtag).optim_additive.params(i)) && equations.(eqtag).optim_additive.bgp(i)
-                    tmp = tmp + DynareModel.params(equations.(eqtag).optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i)*equations.(eqtag).optim_additive.bgp(i);
+                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 - gamma*tmp;
+            cc = cc - tmp0;
         end
         if gamma<1
             if isfield(equations.(eqtag), 'non_optimizing_behaviour.params')
                 % Exogenous variables are present in the 1-λ part (rule of thumb agents).
-                tmp = 0;
+                tmp0 = 0;
+                tmp1 = 0;
                 for i=1:length(equations.(eqtag).non_optimizing_behaviour.params)
-                    if isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && equations.(eqtag).non_optimizing_behaviour.bgp(i)
-                        tmp = tmp + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp(i);
-                    elseif ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && equations.(eqtag).non_optimizing_behaviour.bgp(i)
-                        tmp = tmp + DynareModel.params(equations.(eqtag).non_optimizing_behaviour.params(i))*equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp(i);
+                    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)*tmp;
+                cc = cc - (1.0-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
         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.
-            tmp = 0;
+            tmp0 = 0;
+            tmp1 = 0;
             for i=1:length(equations.(eqtag).additive.params)
-                if isnan(equations.(eqtag).additive.params(i)) && equations.(eqtag).additive.bgp(i)
-                    tmp = tmp + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp(i);
-                elseif ~isnan(equations.(eqtag).additive.params(i)) && equations.(eqtag).additive.bgp(i)
-                    tmp = tmp + DynareModel.params(equations.(eqtag).additive.params(i))*equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp(i);
+                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};
                 end
             end
-            cc = cc - tmp;
+            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
-        DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma;
+        DynareModel.params(pacmodel.growth_neutrality_param_index) = cc; % Multiplies the variable or expression provided though the growth option in command pac_model.
     end
 end
\ No newline at end of file
diff --git a/matlab/+pac/initialize.m b/matlab/+pac/initialize.m
index 5ace9c916e733536ed246bec34e70aa1cd28bfef..1a185c64bba24bdea3268c3585d388b1b5e62d3a 100644
--- a/matlab/+pac/initialize.m
+++ b/matlab/+pac/initialize.m
@@ -8,7 +8,7 @@ function initialize(pacmodel)
 % OUTPUTS
 % None
 
-% Copyright (C) 2018-2019 Dynare Team
+% Copyright © 2018-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,12 +37,12 @@ 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 = false(1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.params));
+        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 = false(1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.params));
+        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 = false(1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.params));
+        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
 end
\ No newline at end of file
diff --git a/matlab/pac-tools/+pac/+bgp/get.m b/matlab/pac-tools/+pac/+bgp/get.m
index d2342e61e4a89d60b3a1e61f6d7ffa51c62d0ab9..04037357cc4c6c74b75c29f9a3dadec731724662 100644
--- a/matlab/pac-tools/+pac/+bgp/get.m
+++ b/matlab/pac-tools/+pac/+bgp/get.m
@@ -1,6 +1,6 @@
 function dummy = get(pacmodel, paceq, kind, id)
 
-% Copyright (C) 2019 Dynare Team
+% Copyright © 2019-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -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).equations.(M_.pac.(pacmodel).tag_map{strcmp(paceq, M_.pac.(pacmodel).tag_map(:,1)),2}).(kind).bgp{id};
diff --git a/matlab/pac-tools/+pac/+bgp/set.m b/matlab/pac-tools/+pac/+bgp/set.m
index d09c5acb7388fbd2a7fea9900c41edfdc05d3ec3..8a2283f736d25267015a67b34861ec3c2ae2fa00 100644
--- a/matlab/pac-tools/+pac/+bgp/set.m
+++ b/matlab/pac-tools/+pac/+bgp/set.m
@@ -59,7 +59,7 @@ function [M_, done] = writebgpfield(type, pacmodel, eqtag, ide, xflag, nonzerome
 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 = zeros(1, length(M_.pac.(pacmodel).equations.(eqtag).(type).params));
+        M_.pac.(pacmodel).equations.(eqtag).(type).bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(eqtag).(type).params));
     end
     [isvar, ie] = ismember(ide, M_.pac.(pacmodel).equations.(eqtag).(type).vars);
     if isvar
@@ -68,7 +68,7 @@ if isfield(M_.pac.(pacmodel).equations.(eqtag), type)
         else
             assert(M_.pac.(pacmodel).equations.(eqtag).(type).isendo(ie), 'Variable type issue.')
         end
-        M_.pac.(pacmodel).equations.(eqtag).(type).bgp(ie) = nonzeromean;
+        M_.pac.(pacmodel).equations.(eqtag).(type).bgp{ie} = nonzeromean;
         done = true;
     end
 end
diff --git a/tests/pac/trend-component-26/example1.mod b/tests/pac/trend-component-26/example1.mod
index d07be333b5e4b172e38b3a0d80044c7711ef7042..2ad3694d07679f0eb8f29a11947b47931ae495a0 100644
--- a/tests/pac/trend-component-26/example1.mod
+++ b/tests/pac/trend-component-26/example1.mod
@@ -36,7 +36,7 @@ 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;
 
@@ -101,45 +101,117 @@ end;
 pac.initialize('pacman');
 
 // Exogenous variables with non zero mean:
-pac.bgp.set('pacman', 'zpac', 's', .2);
-pac.bgp.set('pacman', 'zpac', 'x', .2);
-pac.bgp.set('pacman', 'zpac', 'dx2', .2);
+pac.bgp.set('pacman', 'zpac', 's', true);
+pac.bgp.set('pacman', 'zpac', 'x', .0);
+pac.bgp.set('pacman', 'zpac', 'dx2', .0);
 
 // Update the parameters of the PAC expectation model (h0 and h1 vectors, growth neutrality correction).
 pac.update.expectation('pacman');
 
+id = find(strcmp('s', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.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);
-if ~pac.bgp.get('pacman', 'zpac', 'additive', id)
-   error('bgp field in additive is wrong.')
+if ~isempty(id)
+    error('Variable s should not be under M_.pac.pacman.equations.eq0.additive.')
+end
+
+id = find(strcmp('s', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+if ~isempty(id)
+    error('Variable s should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+end
+
+id = find(strcmp('dv', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.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);
-if pac.bgp.get('pacman', 'zpac', 'additive', id)
-   error('bgp field in additive is wrong.')
+if ~isempty(id)
+    error('Variable dv should not be under M_.pac.pacman.equations.eq0.additive.')
+end
+
+id = find(strcmp('dv', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
+if ~isempty(id)
+    error('Variable dv should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+end
+
+id = find(strcmp('x', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.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);
+if ~isempty(id)
+    error('Variable x should not be under M_.pac.pacman.equations.eq0.optim_additive.')
 end
 
 id = find(strcmp('x', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+if ~isempty(id)
+    error('Variable x should not be under M_.pac.pacman.equations.eq0.additive.')
+end
+
+id = find(strcmp('y', M_.endo_names));
 id = find(id==M_.pac.pacman.equations.eq0.non_optimizing_behaviour.vars);
-if ~pac.bgp.get('pacman', 'zpac', 'non_optimizing_behaviour', id)
-   error('bgp field in non_optimizing_behaviour is wrong.')
+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);
+if ~isempty(id)
+    error('Variable y should not be under M_.pac.pacman.equations.eq0.optim_additive.')
+end
+
+id = find(strcmp('y', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
+if ~isempty(id)
+    error('Variable y should not be under M_.pac.pacman.equations.eq0.additive.')
+end
+
+id = find(strcmp('dx2', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.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);
-if pac.bgp.get('pacman', 'zpac', 'non_optimizing_behaviour', id)
-   error('bgp field in non_optimizing_behaviour is wrong.')
+if ~isempty(id)
+    error('Variable y should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
 end
 
 id = find(strcmp('dx2', M_.endo_names));
 id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
-if ~pac.bgp.get('pacman', 'zpac', 'optim_additive', id)
-   error('bgp field in optim_additive is wrong.')
+if ~isempty(id)
+    error('Variable dx2 should not be under M_.pac.pacman.equations.eq0.optim_additive.')
 end
 
 id = find(strcmp('u', M_.endo_names));
 id = find(id==M_.pac.pacman.equations.eq0.additive.vars);
-if pac.bgp.get('pacman', 'zpac', 'optim_additive', id)
-   error('bgp field in optim_additive is wrong.')
+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);
+if ~isempty(id)
+    error('Variable u should not be under M_.pac.pacman.equations.eq0.non_optimizing_behaviour.')
+end
+
+id = find(strcmp('u', M_.endo_names));
+id = find(id==M_.pac.pacman.equations.eq0.optim_additive.vars);
+if ~isempty(id)
+    error('Variable u should not be under M_.pac.pacman.equations.eq0.optim_additive.')
 end