Bug fix in the definition of the PAC expectation (MCE version).

Replaced

         ᵐ                              ᵐ⁻¹ ᵐ⁻¹
  Zₜ = - ∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1) [ 𝛥 yₜ + ∑   ∑   𝛼ⱼ₊₁𝛽ʲ⁺¹𝛥 yₜ₊ₖ ]
        ᵢ₌₁                             ₖ₌₁ ⱼ₌ₖ

by

         ᵐ                              ᵐ⁻¹ ᵐ⁻¹
  Zₜ = - ∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1) [ 𝛥 yₜ - ∑   ∑   𝛼ⱼ₊₁𝛽ʲ⁺¹𝛥 yₜ₊ₖ ]
        ᵢ₌₁                             ₖ₌₁ ⱼ₌ₖ

Also:

 - Added comments and comparison with the formula in Brayton et alii (2000),
 - Removed unused input argument,
 - Removed unused variables.

matlab/write_pac_mce_expectations.m

-function [expression, growthneutralitycorrection] = write_pac_mce_expectations(eqname, expectationmodelname, iscrlf)
+function [expression, growthneutralitycorrection] = write_pac_mce_expectations(eqname, expectationmodelname)
 
 % Prints the expansion of the PAC_EXPECTATION term in files.
 %
@@ -11,8 +11,7 @@ function [expression, growthneutralitycorrection] = write_pac_mce_expectations(e
 % - expression                     [char]       Unrolled expectation expression.
 % - growthneutralitycorrection     [char]
 
-
-% Copyright © 2019-2023 Dynare Team
+% Copyright © 2019-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,20 +30,12 @@ function [expression, growthneutralitycorrection] = write_pac_mce_expectations(e
 
 global M_
 
-if nargin<4
-    iscrlf = false;
-end
-
 expectationmodel = M_.pac.(expectationmodelname);
 
-pacequation = expectationmodel;
-
-params = M_.params([pacequation.ec.params; pacequation.ar.params(:); expectationmodel.discount_index]);
-
 auxname = sprintf('%s_Z', eqname);
 
-targetid = pacequation.ec.vars((pacequation.ec.istarget==true));
-alphaid = pacequation.mce.alpha;
+targetid = expectationmodel.ec.vars((expectationmodel.ec.istarget==true));
+alphaid = expectationmodel.mce.alpha;
 betaid = expectationmodel.discount_index;
 
 target = M_.endo_names{targetid};
@@ -69,9 +60,33 @@ if isauxiliary(targetid)
     end
 end
 
+
+% In Brayton, Davis and Tulip (2000) the formula for the PAC expectation under Model Consistent Expectations is given as:
+%
+%       ᵐ                              ᵐ⁻¹ ᵐ⁻¹
+%  Zₜ =  ∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1) [ 𝛥 yₜ − ∑   ∑   𝛼ⱼ₊₁𝛽ʲ⁺¹𝛥 yₜ₊ₖ ]
+%       ᵢ₌₁                            ₖ₌₁ ⱼ₌ₖ
+%
+% where yₜ is the target in the PAC equation, see equation 10 or A.92 in the appendix. The sign before the sum of expected
+% Z is incorrect, in the follwing code we define Z as:
+%
+%        ᵐ                              ᵐ⁻¹ ᵐ⁻¹
+%  Zₜ =  -∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1) [ 𝛥 yₜ − ∑   ∑   𝛼ⱼ₊₁𝛽ʲ⁺¹𝛥 yₜ₊ₖ ]
+%        ᵢ₌₁                            ₖ₌₁ ⱼ₌ₖ
+
 expression = '';
 A1 = '1';
 
+% First loop to build
+%
+%     ᵐ
+%    -∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ
+%    ᵢ₌₁
+%
+%              ᵐ
+% and A(1)=1 + ∑ αᵢ
+%             ᵢ₌₁
+
 for i=1:length(alphaid)
     expression = sprintf('%s-%s*(%s^%i)*%s(%i)', expression, ...
                          M_.param_names{alphaid(i)}, ...
@@ -80,11 +95,25 @@ for i=1:length(alphaid)
     A1 = sprintf('%s+%s', A1, M_.param_names{alphaid(i)});
 end
 
+% Write
+%
+%     ᵐ
+%    -∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1)
+%    ᵢ₌₁
+
+
 expression = sprintf('%s+(%s)', expression, A1);
 
+% Write
+%
+%     ᵐ
+%    -∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1) [Δyₜ
+%    ᵢ₌₁
+
 if isempty(transformations)
     expression = sprintf('%s*(diff(%s)', expression, target);
 else
+    % Typically the target may be the log of a variable or the first difference of the log of a variable.
     variable = target;
     for k=length(transformations):-1:1
         variable = sprintf('%s(%s)', transformations{k}, variable);
@@ -92,24 +121,32 @@ else
     expression = sprintf('%s*(diff(%s)', expression, variable);
 end
 
+% Second loop (for the double sum) to write
+%
+%        ᵐ                              ᵐ⁻¹ ᵐ⁻¹
+%  Zₜ =  -∑  𝛼ᵢ 𝛽ⁱ⁺¹ Zₜ₊ᵢ + A(1) [ 𝛥 yₜ − ∑   ∑   𝛼ⱼ₊₁𝛽ʲ⁺¹𝛥 yₜ₊ₖ
+%        ᵢ₌₁                            ₖ₌₁ ⱼ₌ₖ
+
 for i=1:length(alphaid)-1
     tmp = sprintf('%s*%s^%i', M_.param_names{alphaid(i+1)}, M_.param_names{betaid}, i+1);
     for j=i+1:length(alphaid)-1
         tmp = sprintf('%s+%s*%s^%i', tmp, M_.param_names{alphaid(j+1)}, M_.param_names{betaid}, j+1);
     end
     if isempty(transformations)
-        expression = sprintf('%s+(%s)*diff(%s(%i))', expression, tmp, target, i);
+        expression = sprintf('%s-(%s)*diff(%s(%i))', expression, tmp, target, i);
     else
         variable = sprintf('%s(%i)', target, i);
         for k=length(transformations):-1:1
             variable = sprintf('%s(%s)', transformations{k}, variable);
         end
-        expression = sprintf('%s+(%s)*diff(%s(%i))', expression, tmp, variable);
+        expression = sprintf('%s-(%s)*diff(%s(%i))', expression, tmp, variable);
     end
 end
 
+% Close brackets.
 expression = sprintf('%s)', expression);
 
+% Add growth neutrality correction if required.
 if isfield(expectationmodel, 'growth_neutrality_param_index')
     if numel(expectationmodel.growth_linear_comb) == 1
         growthneutralitycorrection = sprintf('%s*%s', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str);