From 6417520c09c6888b0820551c3288145c4c77ba62 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?=
 <stepan@adjemian.eu>
Date: Sat, 23 Sep 2023 17:43:19 +0200
Subject: [PATCH] Allow cherrypicking and aggregation with PAC/MCE.

---
 matlab/write_expectations.m         |  10 ++-
 matlab/write_pac_mce_expectations.m | 121 ++++++++++++++++++++++++++++
 2 files changed, 130 insertions(+), 1 deletion(-)
 create mode 100644 matlab/write_pac_mce_expectations.m

diff --git a/matlab/write_expectations.m b/matlab/write_expectations.m
index 06faabe44..e0eb51d4e 100644
--- a/matlab/write_expectations.m
+++ b/matlab/write_expectations.m
@@ -11,7 +11,7 @@ function [expression, growthneutralitycorrection] = write_expectations(expectati
 % - expression                  [string]    Unrolled expectation expression.
 % - growthneutralitycorrection  [string]
 
-% Copyright © 2019-2021 Dynare Team
+% Copyright © 2019-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -42,6 +42,14 @@ end
 
 expectationmodel = M_.(expectationmodelfield).(expectationmodelname);
 
+if isfield(expectationmodel, 'model_consistent_expectations') && expectationmodel.model_consistent_expectations
+    expression = '';
+    if nargout>1
+        growthneutralitycorrection = 0;
+    end
+    return
+end
+
 if nargout>1 && isequal(expectationmodelkind, 'var')
     error('Cannot return more than one argument if the expectation model is a VAR.')
 end
diff --git a/matlab/write_pac_mce_expectations.m b/matlab/write_pac_mce_expectations.m
new file mode 100644
index 000000000..05fe8e025
--- /dev/null
+++ b/matlab/write_pac_mce_expectations.m
@@ -0,0 +1,121 @@
+function [expression, growthneutralitycorrection] = write_pac_mce_expectations(eqname, expectationmodelname, iscrlf)
+
+% Prints the expansion of the PAC_EXPECTATION term in files.
+%
+% INPUTS
+% - eqname                         [char]       Name of the equation.
+% - epxpectationmodelname          [char]       Name of the expectation model.
+% - iscrlf                         [logical]    Adds carriage return after each additive term if true (default is false).
+%
+% OUTPUTS
+% - expression                     [char]       Unrolled expectation expression.
+% - growthneutralitycorrection     [char]
+
+
+% Copyright © 2019-2023 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/>.
+
+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;
+betaid = expectationmodel.discount_index;
+
+target = M_.endo_names{targetid};
+transformations = {};
+
+if isauxiliary(targetid)
+    ida = get_aux_variable_id(target);
+    op = 0;
+    while ida
+        op = op+1;
+        if isequal(M_.aux_vars(ida).type, 8)
+            transformations(op) = {'diff'};
+            target = M_.endo_names{M_.aux_vars(ida).orig_index};
+            ida = get_aux_variable_id(target);
+        elseif isequal(M_.aux_vars(ida).type, 10)
+            transformations(op) = {M_.aux_vars(ida).unary_op};
+            target = M_.endo_names{M_.aux_vars(ida).orig_index};
+            ida = get_aux_variable_id(target);
+        else
+            error('This case is not implemented.')
+        end
+    end
+end
+
+expression = '';
+A1 = '1';
+
+for i=1:length(alphaid)
+    expression = sprintf('%s-%s*(%s^%i)*%s(%i)', expression, ...
+                         M_.param_names{alphaid(i)}, ...
+                         M_.param_names{betaid}, ...
+                         i, auxname, i);
+    A1 = sprintf('%s+%s', A1, M_.param_names{alphaid(i)});
+end
+
+expression = sprintf('%s+(%s)', expression, A1);
+
+if isempty(transformations)
+    expression = sprintf('%s*(diff(%s)', expression, target);
+else
+    variable = target;
+    for k=length(transformations):-1:1
+        variable = sprintf('%s(%s)', transformations{k}, variable);
+    end
+    expression = sprintf('%s*(diff(%s)', expression, variable);
+end
+
+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);
+    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);
+    end
+end
+
+expression = sprintf('%s)', expression);
+
+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);
+    else
+        growthneutralitycorrection = sprintf('%s*(%s)', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str);
+    end
+else
+    growthneutralitycorrection = '';
+end