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))