From 25231f6634d6376b6a7ed93e4ee939949e5d36b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?=
 <stepan@adjemian.eu>
Date: Mon, 20 Dec 2021 20:38:58 +0100
Subject: [PATCH] Account for PAC components in evaluate routine.

The print_expectations routine was previously only considering the
aggregate expectation (for the target). Now it updates the
database (dseries) with each component of the PAC model. The growth
neutrality correction is included in the aggregate expectation but not
in the expectations of the components.
---
 matlab/print_expectations.m | 42 ++++++++++++++++++++++++++++++++-----
 matlab/write_expectations.m |  4 ++--
 2 files changed, 39 insertions(+), 7 deletions(-)

diff --git a/matlab/print_expectations.m b/matlab/print_expectations.m
index 816b0be734..4a077f2fb4 100644
--- a/matlab/print_expectations.m
+++ b/matlab/print_expectations.m
@@ -217,7 +217,10 @@ end
 % Third print a routine for evaluating VAR_EXPECTATION/PAC_EXPECTATION term (returns a dseries object).
 %
 kind = [expectationmodelkind '_expectations'];
-mkdir(sprintf('+%s/+%s/+%s', M_.fname, kind, expectationmodelname));
+ndir = sprintf('+%s/+%s/+%s', M_.fname, kind, expectationmodelname);
+if ~exist(ndir, 'dir')
+    mkdir(sprintf('+%s/+%s/+%s', M_.fname, kind, expectationmodelname));
+end
 filename = sprintf('+%s/+%s/+%s/evaluate.m', M_.fname, kind, expectationmodelname);
 fid = fopen(filename, 'w');
 
@@ -237,6 +240,12 @@ fprintf(fid, 'ds = dseries();\n\n');
 
 id = 0;
 
+if isfield(expectationmodel, 'h_param_indices')
+    decompose = false;
+else
+    decompose = true;
+end
+
 clear('expression');
 
 % Get coefficient values in the target (if any)
@@ -266,12 +275,20 @@ if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_mo
     if isfield(expectationmodel, 'h_param_indices')
         constant = M_.params(expectationmodel.h_param_indices(id));
     else
+        if decompose
+            expressions = cell(length(expectationmodel.components), 1);
+            for j=1:length(expectationmodel.components)
+                expressions{j} = sprintf('%1.16f', M_.params(expectationmodel.components(j).h_param_indices(id)));
+            end
+        end
         constant = 0;
         for j=1:length(expectationmodel.components)
             constant = constant + targetcoefficients(j)*M_.params(expectationmodel.components(j).h_param_indices(id));
         end
     end
-    expression = sprintf('%1.16f', constant);
+    if isfield(expectationmodel, 'h_param_indices')
+        expression = sprintf('%1.16f', constant);
+    end
 end
 
 for i=1:maxlag
@@ -338,6 +355,16 @@ for i=1:maxlag
                 expression = sprintf('-%1.16f*%s', -parameter, variable);
             end
         end
+        if decompose
+            for k=1:length(expectationmodel.components)
+                parameter = M_.params(expectationmodel.components(k).h_param_indices(id));
+                if parameter>=0
+                    expressions{k} = sprintf('%s+%1.16f*%s', expressions{k}, parameter, variable);;
+                else
+                    expressions{k} = sprintf('%s-%1.16f*%s', expressions{k}, -parameter, variable);
+                end
+            end
+        end
     end
 end
 
@@ -382,7 +409,7 @@ if isequal(expectationmodelkind, 'pac') && growth_correction
     else
         first = true;
         for i=1:length(expectationmodel.components)
-            if ~isequal(expectationmodel.components(i).kind, 'll') && isfield(expectationmodel.components(i), 'growth_neutrality_param_index') && isfield(expectationmodel.components(i), 'growth_linear_comb') && ~isempty(expectationmodel.components(i).growth_linear_comb) 
+            if ~isequal(expectationmodel.components(i).kind, 'll') && isfield(expectationmodel.components(i), 'growth_neutrality_param_index') && isfield(expectationmodel.components(i), 'growth_linear_comb') && ~isempty(expectationmodel.components(i).growth_linear_comb)
                 pgrowth = targetcoefficients(i)*M_.params(expectationmodel.components(i).growth_neutrality_param_index);
                 for iter = 1:numel(expectationmodel.components(i).growth_linear_comb)
                     vgrowth='';
@@ -435,11 +462,16 @@ if isequal(expectationmodelkind, 'pac') && growth_correction
     expression = sprintf('%s+%s', expression, growthcorrection);
 end % growth_correction
 
-fprintf(fid, 'ds.%s = %s;', expectationmodelname, expression);
+fprintf(fid, 'ds.%s = %s;\n', expectationmodelname, expression);
+if exist('expressions', 'var')
+    for i=1:length(expressions)
+        fprintf(fid, 'ds.%s = %s;\n', M_.lhs{expectationmodel.components(i).aux_id}, expressions{i});
+    end
+end
 fclose(fid);
 
 fprintf('Expectation dseries expression is saved in %s.\n', filename);
 
 skipline();
 
-rehash
\ No newline at end of file
+rehash
diff --git a/matlab/write_expectations.m b/matlab/write_expectations.m
index 8186527b36..06faabe44b 100644
--- a/matlab/write_expectations.m
+++ b/matlab/write_expectations.m
@@ -220,7 +220,7 @@ for i=1:maxlag
 end
 
 if aggregate
-    growthneutralitycorrection = ''
+    growthneutralitycorrection = '';
 else
     growthneutralitycorrection = {};
 end
@@ -295,4 +295,4 @@ end
 
 if nargout==1 && ~isempty(growthneutralitycorrection)
     expression = sprintf('%s + %s', expression, growthneutralitycorrection);
-end
\ No newline at end of file
+end
-- 
GitLab