From b297353b06c9c35223500439c0bc63321675768b 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, 25 Nov 2021 16:16:42 +0100
Subject: [PATCH] PAC decomposition between stationary and non-stationary
 components.

The taget in PAC equation can be decomposed into an arbitrary number of components (variables
in the VAR auxiliary model).

TODO Iterative OLS estimation (which is not the preferred estimation routine).
TODO Decomposition in the routine evaluating the forecasts for each component.
---
 matlab/+pac/+estimate/init.m                  |  18 +-
 matlab/+pac/+update/parameters.m              | 296 ++++++++++++------
 matlab/aggregate.m                            |   2 +-
 matlab/cherrypick.m                           |  95 +++++-
 ...t_variables_and_parameters_in_expression.m |  29 +-
 matlab/print_expectations.m                   | 186 ++++++++---
 matlab/write_expectations.m                   | 160 +++++++++-
 preprocessor                                  |   2 +-
 scripts/dynare.el                             |   2 +-
 tests/Makefile.am                             |  15 +
 tests/ecb/cherrypick/test1.mod                |  10 +-
 tests/ecb/cherrypick/test2.mod                | 110 +++++--
 tests/pac/var-10e/example1.mod                | 103 ++++++
 tests/pac/var-10e/example2.mod                | 104 ++++++
 tests/pac/var-11e/example1.mod                | 100 ++++++
 tests/pac/var-8/example1.mod                  |  90 ++++++
 tests/pac/var-8/substitution.mod              |  83 +++++
 tests/pac/var-8e/example1.mod                 |  94 ++++++
 tests/pac/var-9/example1.mod                  |  94 ++++++
 tests/pac/var-9/substitution.mod              |  84 +++++
 tests/pac/var-9e/example1.mod                 |  98 ++++++
 21 files changed, 1568 insertions(+), 207 deletions(-)
 create mode 100644 tests/pac/var-10e/example1.mod
 create mode 100644 tests/pac/var-10e/example2.mod
 create mode 100644 tests/pac/var-11e/example1.mod
 create mode 100644 tests/pac/var-8/example1.mod
 create mode 100644 tests/pac/var-8/substitution.mod
 create mode 100644 tests/pac/var-8e/example1.mod
 create mode 100644 tests/pac/var-9/example1.mod
 create mode 100644 tests/pac/var-9/substitution.mod
 create mode 100644 tests/pac/var-9e/example1.mod

diff --git a/matlab/+pac/+estimate/init.m b/matlab/+pac/+estimate/init.m
index cfed63b0f..8a63a24ef 100644
--- a/matlab/+pac/+estimate/init.m
+++ b/matlab/+pac/+estimate/init.m
@@ -30,14 +30,28 @@ pattern = '(\(model_name\s*=\s*)(?<name>\w+)\)';
 pacmodl = regexp(RHS, pattern, 'names');
 pacmodl = pacmodl.name;
 
+pacmodel = M_.pac.(pacmodl);
+
 % Get the transformed equation to be estimated.
 [lhs, rhs, json] = get_lhs_and_rhs(eqname, M_);
 
 % Get definition of aux. variable for pac expectation...
-auxrhs = json.model{strmatch(sprintf('pac_expectation_%s', pacmodl), M_.lhs, 'exact')}.rhs; 
+if isfield(pacmodel, 'aux_id')
+    auxrhs = {M_.lhs{pacmodel.aux_id}, json.model{pacmodel.aux_id}.rhs};
+elseif isfield(pacmodel, 'components')
+    auxrhs = cell(length(pacmodel.components), 2);
+    for i=1:length(pacmodel.components)
+        auxrhs{i,1} = M_.lhs{pacmodel.components(i).aux_id};
+        auxrhs{i,2} = sprintf('(%s)', json.model{pacmodel.components(i).aux_id}.rhs);
+    end
+else
+    error('Cannot find auxiliary variables for PAC expectation.')
+end
 
 % ... and substitute in rhs.
-rhs = strrep(rhs, sprintf('pac_expectation_%s', pacmodl), auxrhs);
+for i=1:rows(auxrhs)
+    rhs = strrep(rhs, auxrhs{i,1}, auxrhs{i,2});
+end
 
 % Get pacmodel properties
 pacmodel = M_.pac.(pacmodl);
diff --git a/matlab/+pac/+update/parameters.m b/matlab/+pac/+update/parameters.m
index 30b454aa0..09a5e4ffb 100644
--- a/matlab/+pac/+update/parameters.m
+++ b/matlab/+pac/+update/parameters.m
@@ -66,56 +66,118 @@ if ~isfield(varcalib, 'CompanionMatrix') || any(isnan(varcalib.CompanionMatrix(:
 end
 
 % Show the equations where this PAC model is used.
-fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.eq_name);
-skipline()
+if verbose
+    fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.eq_name);
+    skipline()
+end
+
+% Do we need to decompose the PAC expectation?
+if isfield(pacmodel, 'components')
+    numberofcomponents = length(pacmodel.components);
+else
+    numberofcomponents = 0;
+end
 
 % 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
-end
 
-% 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;
+% Get the indices for the stationary/nonstationary variables in the VAR system.
+if numberofcomponents
+    id = cell(numberofcomponents, 1);
+    for i=1:numberofcomponents
+        id(i) = {find(strcmp(DynareModel.endo_names{pacmodel.components(i).endo_var}, varmodel.list_of_variables_in_companion_var))};
+        if isempty(id{i})
+            % 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 j=1:length(ad)
+                    auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(j)}));
+                    if isequal(auxinfo.endo_index, pacmodel.components(i).endo_var)
+                        id(i) = {ad(j)};
+                        break
+                    end
+                    if isequal(auxinfo.type, 8) && isequal(auxinfo.orig_index, pacmodel.components(i).endo_var)
+                        id(i) = {ad(j)};
+                        break
+                    end
+                end
+            end
+            if isempty(id{i})
+                error('Cannot find the trend variable in the Companion VAR/VECM model.')
+            end
         end
     end
 else
-    % Trend component model is assumed.
-    kind = 'dd';
+    id = {find(strcmp(DynareModel.endo_names{pacmodel.ec.vars(pacmodel.ec.istarget)}, varmodel.list_of_variables_in_companion_var))};
+    if isempty(id{1})
+        % 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{1})
+            error('Cannot find the trend variable in the Companion VAR/VECM model.')
+        end
+    end
+end
+
+if ~numberofcomponents
+    % Infer the kind of PAC exoectation
+    if isequal(pacmodel.auxiliary_model_type, 'var')
+        if varmodel.nonstationary(id{1})
+            kind = {'dd'};
+            if varmodel.isconstant
+                id{1} = id{1}+1;
+            end
+        else
+            kind = {'ll'};
+            if varmodel.isconstant
+                id{1} = id{1}+1;
+            end
+        end
+    else
+        % Trend component model is assumed.
+        kind = {'dd'};
+    end
+else
+    if varmodel.isconstant
+        for i=1:numberofcomponents
+            id{i} = id{i}+1;
+        end
+    end
 end
 
 % Override kind with the information provided by the user or update M_.pac
-
+if ~numberofcomponents
+    if ~isempty(pacmodel.kind)
+        kind = {pacmodel.kind};
+    else
+        pacmodel.kind = kind{1};
+    end
+else
+    kind = cell(numberofcomponents,1);
+    for i=1:numberofcomponents
+        if isempty(pacmodel.components(i).kind)
+            error('kind declaration is mandatory for each component in pac_target_info.')
+        else
+            kind{i} = pacmodel.components(i).kind;
+        end
+    end
+end
 
 % Get the value of the discount factor.
 beta = DynareModel.params(pacmodel.discount_index);
@@ -125,6 +187,12 @@ if isfield(pacmodel, 'growth_str')
     growth_flag = true;
 else
     growth_flag = false;
+    for i=1:numberofcomponents
+        if isfield(pacmodel.components(i), 'growth_str')
+            growth_flag = true;
+            break
+        end
+    end
 end
 
 % Do we have rule of thumb agents? γ is the share of optimizing agents.
@@ -136,84 +204,114 @@ 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);
+    h = cell(1,length(id));
+    growthneutrality = cell(1,length(id));
+    for i=1:length(id)
+        [h{i}, growthneutrality{i}] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind{i}, id{i});
+    end
 else
-    h = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind, id);
+    h = cell(1,length(id));
+    for i=1:length(id)
+        h(i) = {hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind{i}, id{i})};
+    end
 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;
+        if isfield(pacmodel, 'h_param_indices')
+            % No decomposition
+            DynareModel.params(pacmodel.h_param_indices) = h{1};
+        else
+            for i=1:numberofcomponents
+                DynareModel.params(pacmodel.components(i).h_param_indices) = h{i};
+            end
+        end
     else
-        DynareModel.params(pacmodel.h_param_indices(1)) = .0;
-        DynareModel.params(pacmodel.h_param_indices(2:end)) = h;
-    end
+        if isfield(pacmodel, 'h_param_indices')
+            % No decomposition
+            DynareModel.params(pacmodel.h_param_indices(1)) = .0;
+            DynareModel.params(pacmodel.h_param_indices(2:end)) = h{1};
+        else
+            for i=1:numberofcomponents
+                DynareModel.params(pacmodel.components(i).h_param_indices(1)) = .0;
+                DynareModel.params(pacmodel.components(i).h_param_indices(2:end)) = h{i};
+            end
+        end
+    end % If the auxiliary model (VAR) has no constant.
 else
-    DynareModel.params(pacmodel.h_param_indices) = h;
-end
-
+    DynareModel.params(pacmodel.h_param_indices) = h{1};
+end % if auxiliary model is a VAR
 
 % 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
+    for j=1:length(id)
+        if isnan(growthneutrality{j})
+            continue
         end
-        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).
+        gg = -(growthneutrality{j}-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;
-            tmp1 = 0;
-            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};
+            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
-            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.
+            cc = cc - tmp0*gamma;
         end
-    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 + 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};
+        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(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 - (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
         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.
+        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
+        if isfield(pacmodel, 'growth_neutrality_param_index')
+            DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
+        else
+            DynareModel.params(pacmodel.components(j).growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
+        end
     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
diff --git a/matlab/aggregate.m b/matlab/aggregate.m
index 38ff663a0..9fdd67670 100644
--- a/matlab/aggregate.m
+++ b/matlab/aggregate.m
@@ -90,7 +90,7 @@ for i=1:length(varargin)
                     model{j} = strrep(model{j}, eqtagname{1}, eqtagname_);
                 end
             else
-                model{j} = eqtagname_;
+                model{j} = strcat('[', eqtagname_, ']');
             end
             % Add equation tag with block name.
             if ~isempty(rootfolder)
diff --git a/matlab/cherrypick.m b/matlab/cherrypick.m
index 0f4fabed2..c1f4d983b 100644
--- a/matlab/cherrypick.m
+++ b/matlab/cherrypick.m
@@ -82,10 +82,10 @@ try
         % Get the original equation.
         [LHS, RHS] = get_lhs_and_rhs(eqtags{i}, M_, true, json);
         % Get the parameters, endogenous and exogenous variables in the current equation.
-        [pnames, ~, xnames] = get_variables_and_parameters_in_equation(LHS, RHS, M_);
+        [pnames, enames, xnames] = get_variables_and_parameters_in_equation(LHS, RHS, M_);
         lhs_expression = LHS;
         LHS = get_variables_and_parameters_in_expression(LHS);
-        enames = LHS;
+        enames = union(enames, LHS);
         if length(LHS)>1
             error('Expressions with more than one variable on the LHS are not allowed.')
         end
@@ -111,7 +111,11 @@ try
         end
         % Remove residual from equation if required.
         if noresids
-            exogenous_variables_to_be_removed = ~ismember(xnames, M_.simulation_exo_names);
+            if isfield(M_, 'simulation_exo_names')
+                exogenous_variables_to_be_removed = ~ismember(xnames, M_.simulation_exo_names);
+            else
+                exogenous_variables_to_be_removed = false(size(xnames));
+            end
             if any(exogenous_variables_to_be_removed)
                 switch sum(exogenous_variables_to_be_removed)
                   case 1
@@ -128,31 +132,77 @@ try
         % Unroll expectation terms if any.
         isvar = regexp(RHS, 'var_expectation\(model_name = (?<name>\w+)\)', 'names');
         ispac = regexp(RHS, 'pac_expectation\(model_name = (?<name>\w+)\)', 'names');
+        istar = regexp(RHS, 'pac_target_nonstationary\(model_name = (?<name>\w+)\)', 'names');
         if ~isempty(isvar)
             rhs = write_expectations(isvar.name, 'var');
-            lhs = sprintf('%s_VE', eqtags{i});
-            RHS = strrep(RHS, sprintf('var_expectation(model_name = %s)', isvar.name), lhs);
+            auxlhs = sprintf('%s_VE', eqtags{i});
+            RHS = strrep(RHS, sprintf('var_expectation(model_name = %s)', isvar.name), auxlhs);
         end
         if ~isempty(ispac)
-            [rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac');
+            if isfield(M_.pac.(ispac.name), 'components')
+                [rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac', false, false);
+            else
+                [rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac');
+            end
             if ~isempty(rhs)
-                lhs = sprintf('%s_PE', eqtags{i});
-                if isempty(growthneutralitycorrection)
+                if iscell(rhs) % PAC expectations are decomposed.
+                    auxlhs = cell(size(rhs));
+                    for k=1:length(M_.pac.(ispac.name).components)
+                        if isequal(k, 1)
+                            lhs = M_.lhs{M_.pac.(ispac.name).components(k).aux_id};
+                            auxlhs{k} = lhs;
+                            if ~isempty(growthneutralitycorrection{k})
+                                lhs = sprintf('%s+%s', lhs, growthneutralitycorrection{k});
+                            end
+                            if  ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
+                                if isempty(growthneutralitycorrection{k})
+                                    lhs = sprintf('%s*%s', M_.pac.(ispac.name).components(k).coeff_str, lhs);
+                                else
+                                    lhs = sprintf('%s*(%s)', M_.pac.(ispac.name).components(k).coeff_str, lhs);
+                                end
+                            end
+                        else
+                            lhs_ = M_.lhs{M_.pac.(ispac.name).components(k).aux_id};
+                            auxlhs{k} = lhs_;
+                            if ~isempty(growthneutralitycorrection{k})
+                                lhs_ = sprintf('%s+%s', lhs_, growthneutralitycorrection{k});
+                            end
+                            if  ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
+                                if isempty(growthneutralitycorrection{k})
+                                    lhs_ = sprintf('%s*%s', M_.pac.(ispac.name).components(k).coeff_str, lhs_);
+                                else
+                                    lhs_ = sprintf('%s*(%s)', M_.pac.(ispac.name).components(k).coeff_str, lhs_);
+                                end
+                            end
+                            lhs = sprintf('%s+%s', lhs, lhs_);
+                        end
+                    end
                     RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), lhs);
                 else
-                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', lhs, growthneutralitycorrection));
+                    auxlhs = M_.lhs{M_.pac.(ispac.name).aux_id};
+                    if isempty(growthneutralitycorrection)
+                        RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), auxlhs);
+                    else
+                        RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', auxlhs, growthneutralitycorrection));
+                    end
                 end
             else
                 % MCE version of the PAC equation.
                 [rhs, growthneutralitycorrection] = write_pac_mce_expectations(eqtags{i}, ispac.name);
-                lhs = sprintf('%s_Z', eqtags{i});
+                auxlhs = sprintf('%s_Z', eqtags{i});
                 if isempty(growthneutralitycorrection)
-                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), lhs);
+                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), auxlhs);
                 else
-                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', lhs, growthneutralitycorrection));
+                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', auxlhs, growthneutralitycorrection));
                 end
             end
         end
+        if ~isempty(istar)
+            % Note that auxlhs and rhs are already defined in the previous block, it is not possible to be here if isempty(ispac) is true.
+            auxlhs{end+1} = M_.endo_names{M_.pac.(ispac.name).ec.vars(M_.pac.(ispac.name).ec.istarget)};
+            rhs{end+1} = M_.aux_vars(strmatch(auxlhs{end}, M_.endo_names, 'exact')==[M_.aux_vars(:).endo_index]).orig_expr;
+            RHS = strrep(RHS, sprintf('pac_target_nonstationary(model_name = %s)', ispac.name), sprintf('%s(-1)', auxlhs{end}));
+        end
         % Print equation for unrolled PAC/VAR-expectation and update
         % list of parameters and endogenous variables (if any).
         if ~isempty(rhs)
@@ -160,13 +210,28 @@ try
             % will not return the lhs variable in expectation_enames since
             % the name is created on the fly and is not a  member of M_.endo_names.
             expectation_pnames = get_variables_and_parameters_in_equation('', rhs, M_);
-            expectation_enames = get_variables_and_parameters_in_expression(lhs);
+            expectation_enames = get_variables_and_parameters_in_expression(auxlhs);
             expectation_xnames = get_variables_and_parameters_in_expression(rhs);
             pnames = union(pnames, expectation_pnames);
             xnames = union(xnames, setdiff(expectation_xnames, expectation_pnames));
             enames = union(enames, expectation_enames);
-            fprintf(fid, '[name=''%s'']\n', lhs);
-            fprintf(fid, '%s = %s;\n\n', lhs, rhs);
+            if ischar(rhs)
+                fprintf(fid, '[name=''%s'']\n', auxlhs);
+                fprintf(fid, '%s = %s;\n\n', auxlhs, rhs);
+            else
+                for k=1:length(rhs)
+                    fprintf(fid, '[name=''%s'']\n', auxlhs{k});
+                    fprintf(fid, '%s = %s;\n\n', auxlhs{k}, rhs{k});
+                end
+            end
+            if isfield(M_.pac.(ispac.name), 'components')
+                for k=1:length(M_.pac.(ispac.name).components)
+                    if ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
+                        params = get_variables_and_parameters_in_expression(M_.pac.(ispac.name).components(k).coeff_str);
+                        pnames = union(pnames, params);
+                    end
+                end
+            end
         else
             pRHS = get_variables_and_parameters_in_equation('', RHS, M_);
             xRHS = get_variables_and_parameters_in_expression(RHS);
diff --git a/matlab/get_variables_and_parameters_in_expression.m b/matlab/get_variables_and_parameters_in_expression.m
index c12d69b24..26af543c0 100644
--- a/matlab/get_variables_and_parameters_in_expression.m
+++ b/matlab/get_variables_and_parameters_in_expression.m
@@ -8,7 +8,7 @@ function objects = get_variables_and_parameters_in_expression(expr)
 % OUTPUTS
 % - objects    [cell]             cell of row char arrays, names of the variables and parameters in expr.
 
-% Copyright (C) 2020 Dynare Team
+% Copyright © 2020-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -25,16 +25,27 @@ function objects = get_variables_and_parameters_in_expression(expr)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-objects = strsplit(expr, {'+','-','*','/','^', ...
+if iscell(expr)
+    objects = splitexpr(expr{1});
+    for i=2:length(expr)
+        objects = [objects, splitexpr(expr{i})];
+    end
+else
+    objects = splitexpr(expr);
+end
+
+% Filter out the numbers, punctuation.
+objects(cellfun(@(x) all(isstrprop(x, 'digit')+isstrprop(x, 'punct')), objects)) = [];
+
+                        % Filter out empty elements.
+                        objects(cellfun(@(x) all(isempty(x)), objects)) = [];
+
+
+                        function o = splitexpr(expr)
+                        o = strsplit(expr, {'+','-','*','/','^', ...
                     'log(', 'log10(', 'ln(', 'exp(', ...
                     'sqrt(', 'abs(', 'sign(', ...
                     'sin(', 'cos(', 'tan(', 'asin(', 'acos(', 'atan(', ...
                     'min(', 'max(', ...
                     'normcdf(', 'normpdf(', 'erf(', ...
-                    'diff(', 'adl(', '(', ')', '\n', '\t', ' '});
-
-% Filter out the numbers, punctuation.
-objects(cellfun(@(x) all(isstrprop(x, 'digit')+isstrprop(x, 'punct')), objects)) = [];
-
-% Filter out empty elements.
-objects(cellfun(@(x) all(isempty(x)), objects)) = [];
\ No newline at end of file
+                    'diff(', 'adl(', '(', ')', '\n', '\t', ' '});
\ No newline at end of file
diff --git a/matlab/print_expectations.m b/matlab/print_expectations.m
index 783baa05f..816b0be73 100644
--- a/matlab/print_expectations.m
+++ b/matlab/print_expectations.m
@@ -127,13 +127,29 @@ switch expectationmodelkind
     end
   case 'pac'
     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)});
+    if isfield(expectationmodel, 'h_param_indices')
+        for i=1:length(expectationmodel.h_param_indices)
+            parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.h_param_indices(i)});
+        end
+    else
+        for j=1:length(expectationmodel.components)
+            for i=1:length(expectationmodel.components(j).h_param_indices)
+                parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.components(j).h_param_indices(i)});
+            end
+        end
     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)));
+        if isfield(expectationmodel, 'h_param_indices')
+            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
+        else
+            for j=1:length(expectationmodel.components)
+                for i=1:length(expectationmodel.components(j).h_param_indices)
+                    fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.components(j).h_param_indices(i)}, M_.params(expectationmodel.components(j).h_param_indices(i)));
+                end
+            end
         end
     end
     if isfield(expectationmodel, 'growth_neutrality_param_index')
@@ -145,10 +161,20 @@ switch expectationmodelkind
         growth_correction = true;
     else
         growth_correction = false;
+        if isfield(expectationmodel, 'components')
+            for j=1:length(expectationmodel.components)
+                if isfield(expectationmodel.components(j), 'growth_neutrality_param_index') && ~isempty(expectationmodel.components(j).growth_neutrality_param_index)
+                    fprintf(fid, '\n');
+                    fprintf(fid, 'parameters %s;\n\n', M_.param_names{expectationmodel.components(j).growth_neutrality_param_index});
+                    if withcalibration
+                        fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.components(j).growth_neutrality_param_index}, M_.params(expectationmodel.components(j).growth_neutrality_param_index));
+                    end
+                    growth_correction = true;
+                end
+            end
+        end
     end
-  otherwise
 end
-
 fclose(fid);
 
 skipline()
@@ -211,6 +237,13 @@ fprintf(fid, 'ds = dseries();\n\n');
 
 id = 0;
 
+clear('expression');
+
+% Get coefficient values in the target (if any)
+if exist(sprintf('+%s/pac_target_coefficients.m', M_.fname), 'file')
+    targetcoefficients = feval(sprintf('%s.pac_target_coefficients', M_.fname), expectationmodelname, M_.params);
+end
+
 maxlag = max(auxmodel.max_lag);
 if isequal(expectationmodel.auxiliary_model_type, 'trend_component')
     % Need to add a lag since the error correction equations are rewritten in levels.
@@ -222,13 +255,23 @@ if isequal(expectationmodelkind, 'var')
 end
 
 if isequal(expectationmodelkind, 'var') && isequal(expectationmodel.auxiliary_model_type, 'var')
+    % Constant in the VAR auxiliary model
     id = id+1;
     expression = sprintf('%1.16f', M_.params(expectationmodel.param_indices(id)));
 end
 
 if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
+    % Constant in the VAR auxiliary model
     id = id+1;
-    expression = sprintf('%1.16f', M_.params(expectationmodel.h_param_indices(id)));
+    if isfield(expectationmodel, 'h_param_indices')
+        constant = M_.params(expectationmodel.h_param_indices(id));
+    else
+        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);
 end
 
 for i=1:maxlag
@@ -256,7 +299,14 @@ for i=1:maxlag
           case 'var'
             parameter = M_.params(expectationmodel.param_indices(id));
           case 'pac'
-            parameter = M_.params(expectationmodel.h_param_indices(id));
+            if isfield(expectationmodel, 'h_param_indices')
+                parameter = M_.params(expectationmodel.h_param_indices(id));
+            else
+                parameter = 0;
+                for k=1:length(expectationmodel.components)
+                    parameter = parameter+targetcoefficients(k)*M_.params(expectationmodel.components(k).h_param_indices(id));
+                end
+            end
           otherwise
         end
         switch expectationmodelkind
@@ -275,61 +325,115 @@ for i=1:maxlag
                 variable = sprintf('%s.%s()', variable, transformations{k});
             end
         end
-        if isequal(id, 1)
-            if isequal(expectationmodelkind, 'pac') && growth_correction
-                pgrowth = M_.params(expectationmodel.growth_neutrality_param_index);
-                linearCombination = '';
-                for iter = 1:numel(expectationmodel.growth_linear_comb)
+        if exist('expression','var')
+            if parameter>=0
+                expression = sprintf('%s+%1.16f*%s', expression, parameter, variable);
+            elseif parameter<0
+                expression = sprintf('%s-%1.16f*%s', expression, -parameter, variable);
+            end
+        else
+            if parameter>=0
+                expression = sprintf('%1.16f*%s', parameter, variable);
+            elseif parameter<0
+                expression = sprintf('-%1.16f*%s', -parameter, variable);
+            end
+        end
+    end
+end
+
+if isequal(expectationmodelkind, 'pac') && growth_correction
+    if isfield(expectationmodel, 'growth_neutrality_param_index')
+        pgrowth = M_.params(expectationmodel.growth_neutrality_param_index);
+        for iter = 1:numel(expectationmodel.growth_linear_comb)
+            vgrowth='';
+            if expectationmodel.growth_linear_comb(iter).exo_id > 0
+                vgrowth = strcat('dbase.', M_.exo_names{expectationmodel.growth_linear_comb(iter).exo_id});
+            elseif expectationmodel.growth_linear_comb(iter).endo_id > 0
+                vgrowth = strcat('dbase.', M_.endo_names{expectationmodel.growth_linear_comb(iter).endo_id});
+            end
+            if expectationmodel.growth_linear_comb(iter).lag ~= 0
+                vgrowth = sprintf('%s(%d)', vgrowth, expectationmodel.growth_linear_comb(iter).lag);
+            end
+            if expectationmodel.growth_linear_comb(iter).param_id > 0
+                if ~isempty(vgrowth)
+                    vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.growth_linear_comb(iter).param_id), vgrowth);
+                else
+                    vgrowth = num2str(M_.params(expectationmodel.growth_linear_comb(iter).param_id), '%1.16f');
+                end
+            end
+            if abs(expectationmodel.growth_linear_comb(iter).constant) ~= 1
+                if ~isempty(vgrowth)
+                    vgrowth = sprintf('%1.16f*%s', expectationmodel.growth_linear_comb(iter).constant, vgrowth);
+                else
+                    vgrowth = num2str(expectationmodel.growth_linear_comb(iter).constant, '%1.16f');
+                end
+            end
+            if iter > 1
+                if expectationmodel.growth_linear_comb(iter).constant > 0
+                    linearCombination = sprintf('%s+%s', linearCombination, vgrowth);
+                else
+                    linearCombination = sprintf('%s-%s', linearCombination, vgrowth);
+                end
+            else
+                linearCombination = vgrowth;
+            end
+        end % loop over growth linear combination elements
+        growthcorrection = sprintf('%1.16f*(%s)', pgrowth, linearCombination);
+    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) 
+                pgrowth = targetcoefficients(i)*M_.params(expectationmodel.components(i).growth_neutrality_param_index);
+                for iter = 1:numel(expectationmodel.components(i).growth_linear_comb)
                     vgrowth='';
-                    if expectationmodel.growth_linear_comb(iter).exo_id > 0
-                        vgrowth = strcat('dbase.', M_.exo_names{expectationmodel.growth_linear_comb(iter).exo_id});
-                    elseif expectationmodel.growth_linear_comb(iter).endo_id > 0
-                        vgrowth = strcat('dbase.', M_.endo_names{expectationmodel.growth_linear_comb(iter).endo_id});
+                    if expectationmodel.components(i).growth_linear_comb(iter).exo_id > 0
+                        vgrowth = strcat('dbase.', M_.exo_names{expectationmodel.components(i).growth_linear_comb(iter).exo_id});
+                    elseif expectationmodel.components(i).growth_linear_comb(iter).endo_id > 0
+                        vgrowth = strcat('dbase.', M_.endo_names{expectationmodel.components(i).growth_linear_comb(iter).endo_id});
+                        % TODO Check if we should not substitute auxiliary variables with original transformed variables here.
                     end
-                    if expectationmodel.growth_linear_comb(iter).lag ~= 0
-                        vgrowth = sprintf('%s(%d)', vgrowth, expectationmodel.growth_linear_comb(iter).lag);
+                    if expectationmodel.components(i).growth_linear_comb(iter).lag ~= 0
+                        vgrowth = sprintf('%s(%d)', vgrowth, expectationmodel.components(i).growth_linear_comb(iter).lag);
                     end
-                    if expectationmodel.growth_linear_comb(iter).param_id > 0
+                    if expectationmodel.components(i).growth_linear_comb(iter).param_id > 0
                         if ~isempty(vgrowth)
-                            vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.growth_linear_comb(iter).param_id), vgrowth);
+                            vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.components(i).growth_linear_comb(iter).param_id), vgrowth);
                         else
-                            vgrowth = num2str(M_.params(expectationmodel.growth_linear_comb(iter).param_id), '%1.16f');
+                            vgrowth = num2str(M_.params(expectationmodel.components(i).growth_linear_comb(iter).param_id), '%1.16f');
                         end
                     end
-                    if abs(expectationmodel.growth_linear_comb(iter).constant) ~= 1
+                    if abs(expectationmodel.components(i).growth_linear_comb(iter).constant) ~= 1
                         if ~isempty(vgrowth)
-                            vgrowth = sprintf('%1.16f*%s', expectationmodel.growth_linear_comb(iter).constant, vgrowth);
+                            vgrowth = sprintf('%1.16f*%s', expectationmodel.components(i).growth_linear_comb(iter).constant, vgrowth);
                         else
-                            vgrowth = num2str(expectationmodel.growth_linear_comb(iter).constant, '%1.16f');
+                            vgrowth = num2str(expectationmodel.components(i).growth_linear_comb(iter).constant, '%1.16f');
                         end
                     end
                     if iter > 1
-                        if expectationmodel.growth_linear_comb(iter).constant > 0
+                        if expectationmodel.components(i).growth_linear_comb(iter).constant > 0
                             linearCombination = sprintf('%s+%s', linearCombination, vgrowth);
                         else
                             linearCombination = sprintf('%s-%s', linearCombination, vgrowth);
                         end
                     else
-                        linearCombination=vgrowth;
+                        linearCombination = vgrowth;
+                    end
+                end % loop over growth linear combination elements
+                if first
+                    growthcorrection = sprintf('%1.16f*(%s)', pgrowth, linearCombination);
+                    first = false;
+                else
+                    if pgrowth>0
+                        growthcorrection = sprintf('%s+%1.16f*(%s)', growthcorrection, pgrowth, linearCombination);
+                    elseif pgrowth<0
+                        growthcorrection = sprintf('%s-%1.16f*(%s)', growthcorrection, -pgrowth, linearCombination);
                     end
                 end
-                if parameter >= 0
-                    expression = sprintf('%1.16f*(%s)+%1.16f*%s', pgrowth, linearCombination, parameter, variable);
-                else
-                    expression = sprintf('%1.16f*(%s)-%1.16f*%s', pgrowth, linearCombination, -parameter, variable);
-                end
-            else
-                expression = sprintf('%1.16f*%s', parameter, variable);
-            end
-        else
-            if parameter>=0
-                expression = sprintf('%s+%1.16f*%s', expression, parameter, variable);
-            else
-                expression = sprintf('%s-%1.16f*%s', expression, -parameter, variable);
             end
         end
     end
-end
+    expression = sprintf('%s+%s', expression, growthcorrection);
+end % growth_correction
 
 fprintf(fid, 'ds.%s = %s;', expectationmodelname, expression);
 fclose(fid);
diff --git a/matlab/write_expectations.m b/matlab/write_expectations.m
index 5312e42a6..8186527b3 100644
--- a/matlab/write_expectations.m
+++ b/matlab/write_expectations.m
@@ -1,4 +1,4 @@
-function [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, iscrlf)
+function [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, iscrlf, aggregate)
 
 % Prints the exansion of the VAR_EXPECTATION or PAC_EXPECTATION term in files.
 %
@@ -48,6 +48,16 @@ end
 
 if nargin<3
     iscrlf = false;
+    aggregate = true;
+end
+
+if nargin<4
+    aggregate = true;
+end
+
+if isfield(expectationmodel, 'h_param_indices')
+    % Disaggregation requires components...
+    aggregate = true;
 end
 
 % Get the name of the associated VAR model and test its existence.
@@ -82,7 +92,29 @@ end
 
 if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
     id = id+1;
-    expression = sprintf('%s', M_.param_names{expectationmodel.h_param_indices(id)});
+    if isfield(expectationmodel, 'h_param_indices')
+        expression = sprintf('%s', M_.param_names{expectationmodel.h_param_indices(id)});
+    else
+        if aggregate
+            if isequal(expectationmodel.components(1).coeff_str, '1')
+                expression = sprintf('%s', M_.param_names{expectationmodel.components(1).h_param_indices(id)});
+            else
+                expression = sprintf('%s*%s', expectationmodel.components(1).coeff_str, M_.param_names{expectationmodel.components(1).h_param_indices(id)});
+            end
+            for i=2:length(expectationmodel.components)
+                if isequal(expectationmodel.components(i).coeff_str, '1')
+                    expression = sprintf('%s+%s', expression, M_.param_names{expectationmodel.components(i).h_param_indices(id)});
+                else
+                    expression = sprintf('%s+%s*%s', expression, expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).h_param_indices(id)});
+                end
+            end
+        else
+            expression = cell(length(expectationmodel.components), 1);
+            for i=1:length(expectationmodel.components)
+                expression(i) = {M_.param_names{expectationmodel.components(i).h_param_indices(id)}};
+            end
+        end
+    end
 end
 
 for i=1:maxlag
@@ -110,7 +142,31 @@ for i=1:maxlag
           case 'var'
             parameter = M_.param_names{expectationmodel.param_indices(id)};
           case 'pac'
-            parameter = M_.param_names{expectationmodel.h_param_indices(id)};
+            if isfield(expectationmodel, 'h_param_indices')
+                parameter = M_.param_names{expectationmodel.h_param_indices(id)};
+            else
+                if aggregate
+                    % TODO Check if we can have parameters entering with a minus sign in the linear combination defining the target.
+                    if isequal(expectationmodel.components(1).coeff_str, '1')
+                        parameter = M_.param_names{expectationmodel.components(1).h_param_indices(id)};
+                    else
+                        parameter = sprintf('%s*%s', expectationmodel.components(1).coeff_str, M_.param_names{expectationmodel.components(1).h_param_indices(id)});
+                    end
+                    for k=2:length(expectationmodel.components)
+                        if isequal(expectationmodel.components(k).coeff_str, '1')
+                            parameter = sprintf('%s+%s', parameter, M_.param_names{expectationmodel.components(k).h_param_indices(id)});
+                        else
+                            parameter = sprintf('%s+%s*%s', parameter, expectationmodel.components(k).coeff_str, M_.param_names{expectationmodel.components(k).h_param_indices(id)});
+                        end
+                    end
+                    parameter = sprintf('(%s)', parameter);
+                else
+                    parameter = cell(length(expectationmodel.components), 1);
+                    for k=1:length(expectationmodel.components)
+                        parameter(k) = {M_.param_names{expectationmodel.components(k).h_param_indices(id)}};
+                    end
+                end
+            end
           otherwise
         end
         switch expectationmodelkind
@@ -128,21 +184,47 @@ for i=1:maxlag
             end
         end
         if isequal(id, 1)
-            if iscrlf
-                expression = sprintf('%s*%s\n', parameter, variable);
+            if aggregate
+                if iscrlf
+                    expression = sprintf('%s*%s\n', parameter, variable);
+                else
+                    expression = sprintf('%s*%s', parameter, variable);
+                end
             else
-                expression = sprintf('%s*%s', parameter, variable);
+                for k=1:length(expectationmodel.components)
+                    if iscrlf
+                        expression(k) = {sprintf('%s*%s\n', parameter{k}, variable)};
+                    else
+                        expression(k) = {sprintf('%s*%s', parameter{k}, variable)};
+                    end
+                end
             end
         else
-            if iscrlf
-                expression = sprintf('%s + %s*%s\n', expression, parameter, variable);
+            if aggregate
+                if iscrlf
+                    expression = sprintf('%s + %s*%s\n', expression, parameter, variable);
+                else
+                    expression = sprintf('%s + %s*%s', expression, parameter, variable);
+                end
             else
-                expression = sprintf('%s + %s*%s', expression, parameter, variable);
+                for k=1:length(expectationmodel.components)
+                    if iscrlf
+                        expression(k) = {sprintf('%s + %s*%s\n', expression{k}, parameter{k}, variable)};
+                    else
+                        expression(k) = {sprintf('%s + %s*%s', expression{k}, parameter{k}, variable)};
+                    end
+                end
             end
         end
     end
 end
 
+if aggregate
+    growthneutralitycorrection = ''
+else
+    growthneutralitycorrection = {};
+end
+
 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);
@@ -150,7 +232,65 @@ if isfield(expectationmodel, 'growth_neutrality_param_index')
         growthneutralitycorrection = sprintf('%s*(%s)', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str);
     end
 else
-    growthneutralitycorrection = '';
+    if isfield(expectationmodel, 'components')
+        if aggregate
+            growthneutralitycorrection = '';
+            for i=1:length(expectationmodel.components)
+                if ~isequal(expectationmodel.components(i).kind, 'll')
+                    if isfield(expectationmodel.components(i), 'growth_neutrality_param_index')
+                        if isempty(growthneutralitycorrection)
+                            if ~isempty(expectationmodel.components(i).growth_str)
+                                if isequal(expectationmodel.components(i).coeff_str, '1')
+                                    if numel(expectationmodel.components(i).growth_linear_comb) == 1
+                                        growthneutralitycorrection = sprintf('%s*%s', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    else
+                                        growthneutralitycorrection = sprintf('%s*(%s)', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    end
+                                else
+                                    if numel(expectationmodel.components(i).growth_linear_comb) == 1
+                                        growthneutralitycorrection = sprintf('%s*%s*%s', expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    else
+                                        growthneutralitycorrection = sprintf('%s*%s*(%s)', expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    end
+                                end
+                            end
+                        else
+                            if ~isempty(expectationmodel.components(i).growth_str)
+                                if isequal(expectationmodel.components(i).coeff_str, '1')
+                                    if numel(expectationmodel.components(i).growth_linear_comb) == 1
+                                        growthneutralitycorrection = sprintf('%s+%s*%s', growthneutralitycorrection, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    else
+                                        growthneutralitycorrection = sprintf('%s+%s*(%s)', growthneutralitycorrection, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    end
+                                else
+                                    if numel(expectationmodel.components(i).growth_linear_comb) == 1
+                                        growthneutralitycorrection = sprintf('%s+%s*%s*%s', growthneutralitycorrection, expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    else
+                                        growthneutralitycorrection = sprintf('%s+%s*%s*(%s)', growthneutralitycorrection, expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
+                                    end
+                                end
+                            end
+                        end
+                    end % if growth neutrality correction for this component
+                end % if non stationary component
+            end
+        else
+            growthneutralitycorrection = repmat({''}, length(expectationmodel.components), 1);
+            for i=1:length(growthneutralitycorrection)
+                if ~isequal(expectationmodel.components(i).kind, 'll')
+                    if isfield(expectationmodel.components(i), 'growth_neutrality_param_index')
+                        if ~isempty(expectationmodel.components(i).growth_str)
+                            if numel(expectationmodel.components(i).growth_linear_comb) == 1
+                                growthneutralitycorrection(i) = {sprintf('%s*%s', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str)};
+                            else
+                                growthneutralitycorrection(i) = {sprintf('%s*(%s)', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str)};
+                            end
+                        end
+                    end % if growth neutrality correction for this component
+                end % if non stationary component
+            end
+        end % if aggregate
+    end
 end
 
 if nargout==1 && ~isempty(growthneutralitycorrection)
diff --git a/preprocessor b/preprocessor
index 7dde09169..a8fce06dc 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 7dde09169ecd4b5e32b3cd35bef88958a262ec11
+Subproject commit a8fce06dc46ca09329e6899e1fde47f2cd81809b
diff --git a/scripts/dynare.el b/scripts/dynare.el
index 39ae6dd11..c3fe64a61 100644
--- a/scripts/dynare.el
+++ b/scripts/dynare.el
@@ -102,7 +102,7 @@
       "osr_params_bounds" "observation_trends" "deterministic_trends" "optim_weights"
       "homotopy_setup" "conditional_forecast_paths" "svar_identification"
       "moment_calibration" "irf_calibration" "ramsey_constraints" "generate_irfs"
-      "matched_moments" "occbin_constraints" "model_replace" "verbatim")
+      "matched_moments" "occbin_constraints" "model_replace" "verbatim" "pac_target_info")
     "Dynare block keywords."))
 
 ;; Mathematical functions and operators used in model equations (see "hand_side" in Bison file)
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9bff02c43..37fecb676 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -518,6 +518,15 @@ ECB_MODFILES = \
 	pac/var-6/substitution.mod \
 	pac/var-7/example1.mod \
 	pac/var-7/substitution.mod \
+	pac/var-8/example1.mod \
+	pac/var-8/substitution.mod \
+	pac/var-8e/example1.mod \
+	pac/var-9/example1.mod \
+	pac/var-9/substitution.mod \
+	pac/var-9e/example1.mod \
+	pac/var-10e/example1.mod \
+	pac/var-10e/example2.mod \
+	pac/var-11e/example1.mod \
 	pac/trend-component-1/example1.mod \
 	pac/trend-component-2/example1.mod \
 	pac/trend-component-3/example1.mod \
@@ -878,6 +887,12 @@ pac/var-6/substitution.m.trs: pac/var-6/example1.m.trs
 pac/var-6/substitution.o.trs: pac/var-6/example1.o.trs
 pac/var-7/substitution.m.trs: pac/var-7/example1.m.trs
 pac/var-7/substitution.o.trs: pac/var-7/example1.o.trs
+pac/var-8/substitution.m.trs: pac/var-8/example1.m.trs
+pac/var-8/substitution.o.trs: pac/var-8/example1.o.trs
+pac/var-9/substitution.m.trs: pac/var-9/example1.m.trs
+pac/var-9/substitution.o.trs: pac/var-9/example1.o.trs
+pac/var-10e/example2.m.trs: pac/var-10e/example1.m.trs
+pac/var-10e/example2.o.trs: pac/var-10e/example1.o.trs
 pac/trend-component-14/substitution.m.trs: pac/trend-component-14/example1.m.trs
 pac/trend-component-14/substitution.o.trs: pac/trend-component-14/example1.o.trs
 
diff --git a/tests/ecb/cherrypick/test1.mod b/tests/ecb/cherrypick/test1.mod
index be64d23f2..d8c063a37 100644
--- a/tests/ecb/cherrypick/test1.mod
+++ b/tests/ecb/cherrypick/test1.mod
@@ -45,7 +45,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;
 
@@ -54,7 +54,7 @@ lambda = 0.5; // Share of optimizing agents.
 
 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);
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman, auxname=rototo);
 
 model;
 
@@ -125,7 +125,7 @@ pac.update.expectation('pacman');
 // variables and exogenous variables in .inc files under ./simulation-files folder. Note that
 // innovations ex1bar and ex2bar will not appear in the equations.
 verbatim;
-  cherrypick('test1', 'simulation-files1', {'z1', 'z2', 'z3'}, true);
-  cherrypick('test1', 'simulation-files2', {'zpac', 'eq:x1', 'eq:x', 'eq:y'}, true);
-  aggregate('toto.mod', {}, '', 'simulation-files1', 'simulation-files2');
+  cherrypick('test1', 'simulation-files-1-a', {'z1', 'z2', 'z3'}, true);
+  cherrypick('test1', 'simulation-files-1-b', {'zpac', 'eq:x1', 'eq:x', 'eq:y'}, true);
+  aggregate('toto1.mod', {}, '', 'simulation-files-1-a', 'simulation-files-1-b');
 end;
diff --git a/tests/ecb/cherrypick/test2.mod b/tests/ecb/cherrypick/test2.mod
index 57ace0d96..d78621d6c 100644
--- a/tests/ecb/cherrypick/test2.mod
+++ b/tests/ecb/cherrypick/test2.mod
@@ -1,36 +1,100 @@
-var x y z;
+// --+ options: json=compute, stochastic +--
 
-model;
+var y x z v u w z1 z2 z3;
 
-x = y(-1)+z(1);
+varexo ex ey ez eu ew;
 
-y = x-z;
+parameters a_y_1 a_y_2 b_y_1 b_y_2 a_u_1 a_u_2 b_u_1 b_u_2 b_x_1 b_x_2 d_y;
 
-diff(z) = 0;
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+a_u_1 =  .2;
+a_u_2 =  .3;
+b_u_1 =  .1;
+b_u_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  auxname pv_y_;
+  kind ll;
+
+  component x;
+  auxname pv_dx_;
+  growth diff(x(-1));
+  kind dl;
 
 end;
 
+model;
 
-if ~isfield(M_, 'exo_names')
-   error('M_.exo_names not defined.')
-end
+  [name='eq:u']
+  u = a_u_1*u(-1) + a_u_2*diff(x(-1)) + b_u_1*y(-2) + b_u_2*diff(x(-2)) + eu ;
 
-if ~iscell(M_.exo_names)
-   error('M_.exo_names not a cell array.')
-end
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
 
-if ~isempty(M_.exo_names)
-   error('M_.exo_names not an empty cell array.')
-end
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
 
-if ~isfield(M_, 'param_names')
-   error('M_.param_names not defined.')
-end
+  [name='eq:v']
+  v = x + d_y*y ;
 
-if ~iscell(M_.param_names)
-   error('M_.param_names not a cell array.')
-end
+  [name='eq:pac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
 
-if ~isempty(M_.param_names)
-   error('M_.param_names not an empty cell array.')
-end
\ No newline at end of file
+  [name='eq:w']
+  diff(w) = -.1*w + ew ;
+
+  [name='eq:z1']
+  z1 = z+y-x+u;
+
+  [name='eq:z2']
+  z2 = z-y+x-u;
+
+  [name='eq:z3']
+  z3 = u-diff(v);
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+  var ew = 1.0;
+  var eu = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Select a subset of the equations and print the equations, the list of parameters, endogenous
+// variables and exogenous variables in .inc files under ./simulation-files folder. Note that
+// innovations ex1bar and ex2bar will not appear in the equations.
+verbatim;
+  cherrypick('test2', 'simulation-files-2-a', {'eq:z1', 'eq:z2', 'eq:z3'}, true);
+  cherrypick('test2', 'simulation-files-2-b', {'eq:pac', 'eq:x', 'eq:y'}, true);
+  aggregate('toto2.mod', {}, '', 'simulation-files-2-a', 'simulation-files-2-b');
+end;
diff --git a/tests/pac/var-10e/example1.mod b/tests/pac/var-10e/example1.mod
new file mode 100644
index 000000000..b185d7345
--- /dev/null
+++ b/tests/pac/var-10e/example1.mod
@@ -0,0 +1,103 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  growth diff(y(-1));
+  auxname pv_y_;
+  kind dd;
+
+  component x;
+  growth diff(x(-1));
+  auxname pv_dx_;
+  kind dd;
+
+end;
+
+model(use_dll);
+
+  [name='eq:y']
+  diff(y) = .01 + a_y_1*diff(y(-1)) + a_y_2*diff(x(-1)) + b_y_1*diff(y(-2)) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = .05 + b_x_1*diff(y(-2)) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [name='zpac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 500);
+
+TrueData.save('example1.data')
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+400, 'fmincon');
+
+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'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
+
+save('example1.estimation.mat', 'e_c_m_nls', 'c_z_1_nls', 'c_z_2_nls')
diff --git a/tests/pac/var-10e/example2.mod b/tests/pac/var-10e/example2.mod
new file mode 100644
index 000000000..ddc56ad8c
--- /dev/null
+++ b/tests/pac/var-10e/example2.mod
@@ -0,0 +1,104 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y', 'eq:v']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, kind=dd, growth=diff(v(-1)), model_name=pacman);
+
+model(use_dll);
+
+  [name='eq:y']
+  diff(y) = .01 + a_y_1*diff(y(-1)) + a_y_2*diff(x(-1)) + b_y_1*diff(y(-2)) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = .05 + b_x_1*diff(y(-2)) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  diff(v) = diff(x) + d_y*diff(y) ;
+
+  [name='zpac']
+  diff(z) = e_c_m*(v(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 500);
+
+TrueData.save('example2.data')
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+400, 'fmincon');
+
+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'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
+
+% Check consistency with disaggregated target
+ts1 = dseries('example1.data.mat');
+ts2 = dseries('example2.data.mat');
+
+if max(abs(ts1.z.data-ts2.z.data))>1e-12
+    error('Simulations in example1 and example2 are not consistent.')
+end
+
+e1 = load('example1.estimation.mat');
+
+if abs(e_c_m_nls-e1.e_c_m_nls)>1e-6 || abs(c_z_1_nls-e1.c_z_1_nls)>1e-6 || abs(c_z_2_nls-e1.c_z_2_nls)>1e-6
+    error('Estimations in example1 and example2 are not consistent.')
+end
+
+delete example1.data.mat
+delete example2.data.mat
+
+delete example1.estimation.mat
diff --git a/tests/pac/var-11e/example1.mod b/tests/pac/var-11e/example1.mod
new file mode 100644
index 000000000..89a56773d
--- /dev/null
+++ b/tests/pac/var-11e/example1.mod
@@ -0,0 +1,100 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y cst; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+cst = .1;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  auxname pv_y_;
+  kind ll;
+
+  component x;
+  growth diff(x(-1));
+  auxname pv_dx_;
+  kind dd;
+
+end;
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = .01 + x + d_y*y ;
+
+  [name='zpac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 5000);
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
+
+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'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
diff --git a/tests/pac/var-8/example1.mod b/tests/pac/var-8/example1.mod
new file mode 100644
index 000000000..8885873b5
--- /dev/null
+++ b/tests/pac/var-8/example1.mod
@@ -0,0 +1,90 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  auxname pv_y_;
+  kind ll;
+
+  component x;
+  auxname pv_dx_;
+  kind dd;
+
+end;
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [name='eq:pac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 20);
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+verbatim;
+  set_dynare_seed('default');
+  y = zeros(M_.endo_nbr,1);
+  y(1:M_.orig_endo_nbr) = rand(M_.orig_endo_nbr, 1);
+  x = randn(M_.exo_nbr,1);
+  y = example1.set_auxiliary_variables(y, x, M_.params);
+  y = [y(find(M_.lead_lag_incidence(1,:))); y];
+  [residual, g1] = example1.dynamic(y, x', M_.params, oo_.steady_state, 1);
+  save('example1.mat', 'residual', 'g1', 'TrueData');
+end;
diff --git a/tests/pac/var-8/substitution.mod b/tests/pac/var-8/substitution.mod
new file mode 100644
index 000000000..9f909dbb3
--- /dev/null
+++ b/tests/pac/var-8/substitution.mod
@@ -0,0 +1,83 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+
+@#include "example1/model/pac-expectations/pacman-parameters.inc"
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [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/pacman-expression.inc"
+  + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 20 periods
+set_dynare_seed('default');
+TrueData = simul_backward_model(initialconditions, 20);
+
+verbatim;
+
+  set_dynare_seed('default');
+  y = zeros(M_.endo_nbr,1);
+  y(1:M_.orig_endo_nbr) = rand(M_.orig_endo_nbr, 1);
+  x = randn(M_.exo_nbr,1);
+  y = substitution.set_auxiliary_variables(y, x, M_.params);
+  y = [y(find(M_.lead_lag_incidence(1,:))); y];
+  example1 = load('example1.mat');
+  [residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
+
+  if max(abs(example1.TrueData.data(:)-TrueData.data(:)))>1e-9
+    error('Simulations do not match.')
+  end
+
+  if ~isequal(length(residual), length(example1.residual)) || max(abs(example1.residual-residual))>1e-8
+    warning('Residuals do not match!')
+  end
+
+  if ~isequal(length(g1(:)), length(example1.g1(:))) || max(abs(example1.g1(:)-g1(:)))>1e-8
+    warning('Jacobian matrices do not match!')
+  end
+
+end;
diff --git a/tests/pac/var-8e/example1.mod b/tests/pac/var-8e/example1.mod
new file mode 100644
index 000000000..9fb1bbc9b
--- /dev/null
+++ b/tests/pac/var-8e/example1.mod
@@ -0,0 +1,94 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  auxname pv_y_;
+  kind ll;
+
+  component x;
+  auxname pv_dx_;
+  kind dd;
+
+end;
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [name='zpac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 5000);
+
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
+
+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'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
diff --git a/tests/pac/var-9/example1.mod b/tests/pac/var-9/example1.mod
new file mode 100644
index 000000000..0a5fcda4a
--- /dev/null
+++ b/tests/pac/var-9/example1.mod
@@ -0,0 +1,94 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  auxname pv_y_;
+  kind ll;
+
+  component x;
+  growth diff(x(-1));
+  auxname pv_dx_;
+  kind dd;
+
+end;
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [name='eq:pac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 20);
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+verbatim;
+  set_dynare_seed('default');
+  y = zeros(M_.endo_nbr,1);
+  y(1:M_.orig_endo_nbr) = rand(M_.orig_endo_nbr, 1);
+  x = randn(M_.exo_nbr,1);
+  y = example1.set_auxiliary_variables(y, x, M_.params);
+  y = [y(find(M_.lead_lag_incidence(1,:))); y];
+  [residual, g1] = example1.dynamic(y, x', M_.params, oo_.steady_state, 1);
+  save('example1.mat', 'residual', 'g1', 'TrueData');
+end;
diff --git a/tests/pac/var-9/substitution.mod b/tests/pac/var-9/substitution.mod
new file mode 100644
index 000000000..d9e70378b
--- /dev/null
+++ b/tests/pac/var-9/substitution.mod
@@ -0,0 +1,84 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+@#include "example1/model/pac-expectations/pacman-parameters.inc"
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [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/pacman-expression.inc"
+  +
+  @#include "example1/model/pac-expectations/pacman-growth-neutrality-correction.inc"
+  + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 20 periods
+set_dynare_seed('default');
+TrueData = simul_backward_model(initialconditions, 20);
+
+verbatim;
+
+  set_dynare_seed('default');
+  y = zeros(M_.endo_nbr,1);
+  y(1:M_.orig_endo_nbr) = rand(M_.orig_endo_nbr, 1);
+  x = randn(M_.exo_nbr,1);
+  y = substitution.set_auxiliary_variables(y, x, M_.params);
+  y = [y(find(M_.lead_lag_incidence(1,:))); y];
+  example1 = load('example1.mat');
+  [residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
+
+  if max(abs(example1.TrueData.data(:)-TrueData.data(:)))>1e-9
+    error('Simulations do not match.')
+  end
+
+  if ~isequal(length(residual), length(example1.residual)) || max(abs(example1.residual-residual))>1e-8
+    warning('Residuals do not match!')
+  end
+
+  if ~isequal(length(g1(:)), length(example1.g1(:))) || max(abs(example1.g1(:)-g1(:)))>1e-8
+    warning('Jacobian matrices do not match!')
+  end
+
+end;
diff --git a/tests/pac/var-9e/example1.mod b/tests/pac/var-9e/example1.mod
new file mode 100644
index 000000000..4d0169f33
--- /dev/null
+++ b/tests/pac/var-9e/example1.mod
@@ -0,0 +1,98 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z v;
+
+varexo ex ey ez ;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
+
+parameters beta e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+d_y = .5;
+
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
+
+pac_target_info(pacman);
+  target v;
+  auxname_target_nonstationary vns;
+
+  component y;
+  auxname pv_y_;
+  kind ll;
+
+  component x;
+  growth diff(x(-1));
+  auxname pv_dx_;
+  kind dd;
+
+end;
+
+model;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+
+  [name='eq:x']
+  diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
+
+  [name='eq:v']
+  v = x + d_y*y ;
+
+  [name='zpac']
+  diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+end;
+
+// Initialize the PAC model (build the Companion VAR representation for the auxiliary model).
+pac.initialize('pacman');
+
+// Update the parameters of the PAC expectation model (h0 and h1 vectors).
+pac.update.expectation('pacman');
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 500 periods
+TrueData = simul_backward_model(initialconditions, 5000);
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman', 'eq:pac');
+
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
+
+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'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))