diff --git a/matlab/+pac/+estimate/init.m b/matlab/+pac/+estimate/init.m
index ceadf809e..b6caf2561 100644
--- a/matlab/+pac/+estimate/init.m
+++ b/matlab/+pac/+estimate/init.m
@@ -1,4 +1,4 @@
-function [pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, pnames_, ipnames_, params, data, islaggedvariables, eqtag] = ...
+function [pacmodl, lhs, rhs, pnames, enames, xnames, rname, pid, eid, xid, pnames_, ipnames_, params, data, islaggedvariables, eqtag] = ...
     init(M_, oo_, eqname, params, data, range)
 
 % Copyright (C) 2018-2019 Dynare Team
@@ -70,7 +70,10 @@ end
 stack = dbstack;
 ipnames__ = ipnames_;                              % The user provided order.
 if isequal(stack(2).name, 'iterative_ols')
-    ipnames_  = [M_.pac.(pacmodl).equations.(eqtag).ec.params; M_.pac.(pacmodl).equations.(eqtag).ar.params']; % The correct order.
+    ipnames_ = [M_.pac.(pacmodl).equations.(eqtag).ec.params; M_.pac.(pacmodl).equations.(eqtag).ar.params'];
+    if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive')
+        ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).additive.params(~isnan(M_.pac.(pacmodl).equations.(eqtag).additive.params))'];
+    end
     if isfield(M_.pac.(pacmodl).equations.(eqtag), 'share_of_optimizing_agents_index')
         ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index];
     end
diff --git a/matlab/+pac/+estimate/iterative_ols.m b/matlab/+pac/+estimate/iterative_ols.m
index fdf95fa82..af5a2425c 100644
--- a/matlab/+pac/+estimate/iterative_ols.m
+++ b/matlab/+pac/+estimate/iterative_ols.m
@@ -43,7 +43,7 @@ global M_ oo_
 
 debug = false; % If true, prints the value of SSR for the initial guess.
 
-[pacmodl, ~, rhs, ~, ~, ~, ~, ~, ~, ~, ipnames_, params, data, ~, eqtag] = ...
+[pacmodl, ~, rhs, ~, ~, ~, rname, ~, ~, ~, ~, ipnames_, params, data, ~, eqtag] = ...
     pac.estimate.init(M_, oo_, eqname, params, data, range);
 
 % Set initial condition.
@@ -52,12 +52,19 @@ params0 = cell2mat(struct2cell(params));
 % Set flag for models with non optimizing agents.
 is_non_optimizing_agents = isfield(M_.pac.(pacmodl).equations.(eqtag), 'non_optimizing_behaviour');
 
+% Set flag for models with exogenous variables (outside of non optimizing agents part)
+if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive')
+    is_exogenous_variables = length(M_.pac.(pacmodl).equations.(eqtag).additive.vars)>1;
+else
+    is_exogenous_variables = false;
+end
+
 if is_non_optimizing_agents
     non_optimizing_behaviour = M_.pac.(pacmodl).equations.(eqtag).non_optimizing_behaviour;
     non_optimizing_behaviour_params = NaN(length(non_optimizing_behaviour.params), 1);
     noparams = isnan(non_optimizing_behaviour.params);
     if ~all(noparams)
-        % Find estimated non optimizing behaviour parameters (if any.
+        % Find estimated non optimizing behaviour parameters (if any).
         non_optimizing_behaviour_estimated_params = ismember(M_.param_names(non_optimizing_behaviour.params), fieldnames(params));
         if any(non_optimizing_behaviour_estimated_params)
             error('The estimation of non optimizing behaviour parameters is not yet allowed.')
@@ -77,11 +84,32 @@ if is_non_optimizing_agents
             error('The share of optimizing agents shoud be in (0,1).')
         end
     end
+    share_of_optimizing_agents_index = M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index;
 else
     share_of_optimizing_agents = 1.0;
+    share_of_optimizing_agents_index = [];
     estimate_non_optimizing_agents_share = false;
 end
 
+if is_exogenous_variables
+    additive = M_.pac.(pacmodl).equations.(eqtag).additive;
+    residual_id = find(strcmp(rname, M_.exo_names));
+    residual_jd = find(additive.vars==residual_id & ~additive.isendo);
+    additive.params(residual_jd) = [];
+    additive.vars(residual_jd) = [];
+    additive.isendo(residual_jd) = [];
+    additive.lags(residual_jd) = [];
+    additive.scaling_factor(residual_jd) = [];
+    additive.estimation  = ismember(additive.params, ipnames_);
+else
+    additive.params = [];
+    additive.vars = [];
+    additive.isendo = [];
+    additive.lags = [];
+    additive.scaling_factor = [];
+    additive.estimation = [];
+end
+
 % Build PAC expectation matrix expression.
 dataForPACExpectation0 = dseries();
 listofvariables0 = {};
@@ -144,6 +172,42 @@ else
     dataForNonOptimizingBehaviour = dseries();
 end
 
+% Build data for exogenous variables (out of non optimizing behaviour term).
+if is_exogenous_variables
+    listofvariables2 = {}; j = 0;
+    dataForExogenousVariables = dseries();  % Estimated parameters
+    dataForExogenousVariables_ = 0;         % Calibrated parameters
+    is_any_calibrated_parameter_x = false;
+    is_any_estimated_parameter_x = false;
+    for i=1:length(additive.vars)
+        if additive.isendo(i)
+            variable = M_.endo_names{additive.vars(i)};
+        else
+            variable = M_.exo_names{additive.vars(i)};
+        end
+        if additive.estimation(i)
+            j = j+1;
+            is_any_estimated_parameter_x = true;
+            listofvariables2{j} = variable;
+            dataForExogenousVariables = [dataForExogenousVariables, additive.scaling_factor(i)*data{variable}.lag(additive.lags(i))];
+        else
+            is_any_calibrated_parameter_x = true;
+            tmp = data{variable}.lag(additive.lags(i)).data;
+            if ~isnan(additive.params(i))
+                tmp = M_.params(additive.params(i))*tmp;
+            end
+            tmp = additive.scaling_factor(i)*tmp;
+            dataForExogenousVariables_ = dataForExogenousVariables_+tmp;
+        end
+    end
+    if is_any_calibrated_parameter_x
+        dataForExogenousVariables_ = dseries(dataForExogenousVariables_, data.dates(1), 'exogenous_variables_associated_with_calibrated_parameters');
+    end
+else
+    dataForExogenousVariables = dseries();
+    dataForExogenousVariables_ = dseries();
+end
+
 % Reorder ec.vars locally if necessary. Second variable must be the
 % endogenous variable, while the first must be the associated trend.
 if M_.pac.(pacmodl).equations.(eqtag).ec.isendo(2)
@@ -152,7 +216,7 @@ else
     ecvars = flip(M_.pac.(pacmodl).equations.(eqtag).ec.vars);
 end
 
-% Build matrix for EC and AR terms.
+%% Build matrix for EC and AR terms.
 DataForOLS = dseries();
 
 % Error correction term is trend minus the level of the endogenous variable.
@@ -169,16 +233,39 @@ end
 
 XDATA = DataForOLS{listofvariables3{:}}(range).data;
 
-if estimate_non_optimizing_agents_share
-    [~, ecm_params_id] = setdiff(ipnames_, M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index);
-    [~, share_param_id] = setdiff(1:length(ipnames_), ecm_params_id);
-    params0_ = params0(ecm_params_id);
-    share_of_optimizing_agents = params0(share_param_id);
-    if share_of_optimizing_agents>1 || share_of_optimizing_agents<0
-        error('Initial value for the share of optimizing agents shoud be in (0,1).')
-    end
+% Get index in params0 for share of optimizing agents parameter (if
+% not estimated, params_id_0 is empty).
+if is_non_optimizing_agents
+    params_id_0 = find(ipnames_==share_of_optimizing_agents_index);
 else
-    params0_ = params0;
+    params_id_0 = [];
+end
+
+% Get indices in params0 for EC and AR parameters
+[~, params_id_1] = setdiff(ipnames_, [share_of_optimizing_agents_index, additive.params]);
+
+% Get indices in params0 for other parameters (optimizing agents
+% share plus parameters related to exogenous variables).
+[~, params_id_2] = setdiff(1:length(ipnames_), params_id_1);
+
+% Get indices in params0 for the parameters associated to the
+% exogenous variables.
+params_id_3 = setdiff(params_id_2, params_id_0);
+
+% Get values for EC and AR parameters
+params0_ = params0(params_id_1);
+
+% Get values for parameters associated to the exogenous variables.
+params0__ = params0(params_id_3);
+
+% Get value of the share of optimizing agents.
+if estimate_non_optimizing_agents_share
+    share_of_optimizing_agents = params0(params_id_0);
+end
+
+% Check that the share is in (0,1)
+if share_of_optimizing_agents>1 || share_of_optimizing_agents<0
+    error('Initial value for the share of optimizing agents shoud be in (0,1).')
 end
 
 % Update the vector of parameters.
@@ -187,19 +274,6 @@ M_.params(ipnames_) = params0;
 % Update the reduced form PAC expectation parameters and compute the expectations.
 [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
 
-if debug
-    YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
-    if is_non_optimizing_agents
-        YDATA = YDATA-share_of_optimizing_agents*PacExpectations;
-        YDATA = YDATA-(1-share_of_optimizing_agents)*(dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params);
-    else
-        YDATA = YDATA-PacExpectations;
-    end
-    r = YDATA-XDATA*(params0_*share_of_optimizing_agents);
-    ssr = r'*r;
-    fprintf('\nInitial value of the objective (SSR) is %s.\n\n', num2str(ssr));
-end
-
 noconvergence = true;
 counter = 0;
 
@@ -213,6 +287,14 @@ while noconvergence
     else
         YDATA = YDATA-PacExpectations;
     end
+    if is_exogenous_variables
+        if is_any_calibrated_parameter_x
+            YDATA = YDATA-dataForExogenousVariables_(range).data;
+        end
+        if is_any_estimated_parameter_x
+            YDATA = YDATA-dataForExogenousVariables{listofvariables2{:}}(range).data*params0__;
+        end
+    end
     % Run OLS to estimate PAC parameters (autoregressive parameters and error correction parameter).
     params1_ = (XDATA\YDATA)/share_of_optimizing_agents;
     % Compute residuals and sum of squareed residuals.
@@ -226,18 +308,53 @@ while noconvergence
     % Update the value of the share of non optimizing agents (if estimated)
     if estimate_non_optimizing_agents_share
         % First update the parameters and compute the PAC expectation reduced form parameters.
-        M_.params(ipnames_(ecm_params_id)) = params0_;
+        M_.params(ipnames_(params_id_1)) = params0_;
         [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
         % Set vector for left handside variable.
         YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
         YDATA = YDATA-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params;
+        if is_exogenous_variables && is_any_calibrated_parameter_x
+            YDATA = YDATA-dataForExogenousVariables_(range).data;
+        end
         % Set vector for regressor.
         ZDATA = XDATA*params0_+PacExpectations-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params;
+        if is_exogenous_variables && is_any_estimated_parameter_x
+            ZDATA = [ZDATA, dataForExogenousVariables{listofvariables2{:}}(range).data];
+        end
         % Update the (estimated) share of optimizing agents by running OLS
-        share_of_optimizing_agents = (ZDATA\YDATA);
+        beta = (ZDATA\YDATA);
+        share_of_optimizing_agents = beta(1);
+        if is_exogenous_variables && is_any_estimated_parameter_x
+            params0__ = beta(2:end);
+        end
         % Force the share of optimizing agents to be in [0,1].
         share_of_optimizing_agents = max(min(share_of_optimizing_agents, 1.0), 0.0);
-        M_.params(ipnames_(share_param_id)) = share_of_optimizing_agents;
+        M_.params(ipnames_(params_id_0)) = share_of_optimizing_agents;
+        if is_exogenous_variables && is_any_estimated_parameter_x
+            M_.params(ipnames_(params_id_3)) = params0__;
+        end
+    elseif is_exogenous_variables && is_any_estimated_parameter_x
+        % First update the parameters and compute the PAC expectation reduced form parameters.
+        M_.params(ipnames_(params_id_1)) = params0_;
+        [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
+        % Set vector for left handside variable.
+        YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
+        if is_any_calibrated_parameter_x
+            YDATA = YDATA-dataForExogenousVariables_(range).data;
+        end
+        if is_non_optimizing_agents
+            YDATA = YDATA-share_of_optimizing_agents*(XDATA*params0_+PacExpectations) - ...
+                    (1-share_of_optimizing_agents)*(dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params);
+        else
+            YDATA = YDATA-XDATA*params0_-PacExpectations;
+        end
+        % Set vector for regressor
+        ZDATA = dataForExogenousVariables{listofvariables2{:}}(range).data;
+        % Update the (estimated) parameters associated to the
+        % exogenous variables.
+        beta = (ZDATA\YDATA);
+        params0__ = beta;
+        M_.params(ipnames_(params_id_3)) = params0__;
     else
         M_.params(ipnames_) = params0_;
         [PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
diff --git a/matlab/+pac/+estimate/nls.m b/matlab/+pac/+estimate/nls.m
index 66721a709..e71b341a0 100644
--- a/matlab/+pac/+estimate/nls.m
+++ b/matlab/+pac/+estimate/nls.m
@@ -76,7 +76,7 @@ if nargin>4 && (isequal(optimizer, 'GaussNewton') || isequal(optimizer, 'lsqnonl
     objective = 'r_';
 end
 
-[pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables, eqtag] = ...
+[pacmodl, lhs, rhs, pnames, enames, xnames, ~, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables, eqtag] = ...
     pac.estimate.init(M_, oo_, eqname, params, data, range);
 
 % Check that the error correction term is correct.
diff --git a/preprocessor b/preprocessor
index 2312ce13d..5f013756f 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 2312ce13dc58d3645e0f37f1a757fbe50b6932fd
+Subproject commit 5f013756f23fe89ec946fcff2cfdad7627693e12
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 82f1b764d..48cabfd98 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -426,6 +426,10 @@ MODFILES = \
     pac/trend-component-20-2/example.mod \
     pac/trend-component-20-3/example.mod \
     pac/trend-component-20-4/example.mod \
+    pac/trend-component-21/example.mod \
+    pac/trend-component-22/example.mod \
+    pac/trend-component-23/example.mod \
+    pac/trend-component-24/example.mod \
     ecb/backward-models/irf/solow_1.mod \
     ecb/backward-models/irf/solow_2.mod \
 	dynare-command-options/ramst.mod
diff --git a/tests/pac/trend-component-21/clean b/tests/pac/trend-component-21/clean
new file mode 100755
index 000000000..cf3492a28
--- /dev/null
+++ b/tests/pac/trend-component-21/clean
@@ -0,0 +1,8 @@
+#!/bin/sh
+
+rm -rf example
+rm  -rf +example
+rm -f example.log
+rm -f *.mat
+rm -f *.m
+rm -f *.dat
diff --git a/tests/pac/trend-component-21/example.mod b/tests/pac/trend-component-21/example.mod
new file mode 100644
index 000000000..d127a0db8
--- /dev/null
+++ b/tests/pac/trend-component-21/example.mod
@@ -0,0 +1,160 @@
+// --+ options: json=compute, stochastic +--
+
+var x1 x2 x1bar x2bar z y x;
+
+varexo ex1 ex2 ex1bar ex2bar ez ey ex;
+
+parameters
+       rho_1 rho_2 rho_3 rho_4
+       a_x1_0 a_x1_1 a_x1_2 a_x1_x2_1 a_x1_x2_2
+	   a_x2_0 a_x2_1 a_x2_2 a_x2_x1_1 a_x2_x1_2
+	   e_c_m c_z_1 c_z_2 c_z_dx2 beta
+       lambda;
+
+rho_1 =  .9;
+rho_2 = -.2;
+rho_3 =  .4;
+rho_4 = -.3;
+
+
+a_x1_0 =  -.9;
+a_x1_1 =  .4;
+a_x1_2 =  .3;
+a_x1_x2_1 = .1;
+a_x1_x2_2 = .2;
+
+a_x2_0 =  -.9;
+a_x2_1 =   .2;
+a_x2_2 =  -.1;
+a_x2_x1_1 = -.1;
+a_x2_x1_2 = .2;
+
+beta  =  .2;
+e_c_m =  .5;
+c_z_1 =  .2;
+c_z_2 = -.1;
+c_z_dx2 = .3;
+
+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);
+
+model;
+
+[name='eq:y']
+y = rho_1*y(-1) + rho_2*y(-2) + ey;
+
+[name='eq:x']
+x = rho_3*x(-1) + rho_4*x(-2) + ex;
+
+[name='eq:x1']
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
+
+[name='eq:x2']
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
+
+[name='eq:x1bar']
+x1bar = x1bar(-1) + ex1bar;
+
+[name='eq:x2bar']
+x2bar = x2bar(-1) + ex2bar;
+
+[name='zpac']
+diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + ez;
+
+end;
+
+shocks;
+    var ex1 = 1.0;
+    var ex2 = 1.0;
+    var ex1bar = 1.0;
+    var ex2bar = 1.0;
+    var ez = 1.0;
+    var ey = 0.1;
+    var ex = 0.1;
+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);
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.c_z_dx2 =  .5;
+eparams.lambda  =  .7;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
+
+e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
+c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
+c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
+c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
+lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
+disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols))
+disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols))
+
+skipline(2)
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.c_z_dx2 =  .5;
+eparams.lambda  =  .7;
+
+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'));
+c_z_dx2_nls = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
+lambda_nls = M_.params(strmatch('lambda', 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))
+disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls))
+disp(sprintf('Estimate of lambda: %f', lambda_nls))
+
+if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
+end
+
+if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
+end
+
+if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
+end
+
+if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)')
+end
+
+if abs(lambda_nls-lambda_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (lambda)')
+end
\ No newline at end of file
diff --git a/tests/pac/trend-component-22/clean b/tests/pac/trend-component-22/clean
new file mode 100755
index 000000000..cf3492a28
--- /dev/null
+++ b/tests/pac/trend-component-22/clean
@@ -0,0 +1,8 @@
+#!/bin/sh
+
+rm -rf example
+rm  -rf +example
+rm -f example.log
+rm -f *.mat
+rm -f *.m
+rm -f *.dat
diff --git a/tests/pac/trend-component-22/example.mod b/tests/pac/trend-component-22/example.mod
new file mode 100644
index 000000000..1f8d36b47
--- /dev/null
+++ b/tests/pac/trend-component-22/example.mod
@@ -0,0 +1,150 @@
+// --+ options: json=compute, stochastic +--
+
+var x1 x2 x1bar x2bar z y x;
+
+varexo ex1 ex2 ex1bar ex2bar ez ey ex;
+
+parameters
+       rho_1 rho_2 rho_3 rho_4
+       a_x1_0 a_x1_1 a_x1_2 a_x1_x2_1 a_x1_x2_2
+	   a_x2_0 a_x2_1 a_x2_2 a_x2_x1_1 a_x2_x1_2
+	   e_c_m c_z_1 c_z_2 c_z_dx2 beta
+       lambda;
+
+rho_1 =  .9;
+rho_2 = -.2;
+rho_3 =  .4;
+rho_4 = -.3;
+
+
+a_x1_0 =  -.9;
+a_x1_1 =  .4;
+a_x1_2 =  .3;
+a_x1_x2_1 = .1;
+a_x1_x2_2 = .2;
+
+a_x2_0 =  -.9;
+a_x2_1 =   .2;
+a_x2_2 =  -.1;
+a_x2_x1_1 = -.1;
+a_x2_x1_2 = .2;
+
+beta  =  .2;
+e_c_m =  .5;
+c_z_1 =  .2;
+c_z_2 = -.1;
+c_z_dx2 = .3;
+
+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);
+
+model;
+
+[name='eq:y']
+y = rho_1*y(-1) + rho_2*y(-2) + ey;
+
+[name='eq:x']
+x = rho_3*x(-1) + rho_4*x(-2) + ex;
+
+[name='eq:x1']
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
+
+[name='eq:x2']
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
+
+[name='eq:x1bar']
+x1bar = x1bar(-1) + ex1bar;
+
+[name='eq:x2bar']
+x2bar = x2bar(-1) + ex2bar;
+
+[name='zpac']
+diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + ez;
+
+end;
+
+shocks;
+    var ex1 = 1.0;
+    var ex2 = 1.0;
+    var ex1bar = 1.0;
+    var ex2bar = 1.0;
+    var ez = 1.0;
+    var ey = 0.1;
+    var ex = 0.1;
+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);
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.c_z_dx2 =  .5;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
+
+e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
+c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
+c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
+c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
+disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols))
+
+skipline(2)
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.c_z_dx2 =  .5;
+
+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'));
+c_z_dx2_nls = M_.params(strmatch('c_z_dx2', 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))
+disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls))
+
+if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
+end
+
+if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
+end
+
+if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
+end
+
+if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)')
+end
diff --git a/tests/pac/trend-component-23/clean b/tests/pac/trend-component-23/clean
new file mode 100755
index 000000000..cf3492a28
--- /dev/null
+++ b/tests/pac/trend-component-23/clean
@@ -0,0 +1,8 @@
+#!/bin/sh
+
+rm -rf example
+rm  -rf +example
+rm -f example.log
+rm -f *.mat
+rm -f *.m
+rm -f *.dat
diff --git a/tests/pac/trend-component-23/example.mod b/tests/pac/trend-component-23/example.mod
new file mode 100644
index 000000000..c299ab6ce
--- /dev/null
+++ b/tests/pac/trend-component-23/example.mod
@@ -0,0 +1,165 @@
+// --+ options: json=compute, stochastic +--
+
+var x1 x2 x1bar x2bar z y x u;
+
+varexo ex1 ex2 ex1bar ex2bar ez ey ex eu;
+
+parameters
+       rho_1 rho_2 rho_3 rho_4
+       a_x1_0 a_x1_1 a_x1_2 a_x1_x2_1 a_x1_x2_2
+	   a_x2_0 a_x2_1 a_x2_2 a_x2_x1_1 a_x2_x1_2
+	   e_c_m c_z_1 c_z_2 c_z_dx2 c_z_u beta
+       lambda;
+
+rho_1 =  .9;
+rho_2 = -.2;
+rho_3 =  .4;
+rho_4 = -.3;
+
+
+a_x1_0 =  -.9;
+a_x1_1 =  .4;
+a_x1_2 =  .3;
+a_x1_x2_1 = .1;
+a_x1_x2_2 = .2;
+
+a_x2_0 =  -.9;
+a_x2_1 =   .2;
+a_x2_2 =  -.1;
+a_x2_x1_1 = -.1;
+a_x2_x1_2 = .2;
+
+beta  =  .2;
+e_c_m =  .5;
+c_z_1 =  .2;
+c_z_2 = -.1;
+c_z_dx2 = .3;
+c_z_u = .3;
+
+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);
+
+model;
+
+[name='eq:u']
+u = .5*u(-1) - .2*u(-2) + eu;
+
+[name='eq:y']
+y = rho_1*y(-1) + rho_2*y(-2) + ey;
+
+[name='eq:x']
+x = rho_3*x(-1) + rho_4*x(-2) + ex;
+
+[name='eq:x1']
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
+
+[name='eq:x2']
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
+
+[name='eq:x1bar']
+x1bar = x1bar(-1) + ex1bar;
+
+[name='eq:x2bar']
+x2bar = x2bar(-1) + ex2bar;
+
+[name='zpac']
+diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + c_z_u*u + ez;
+
+end;
+
+shocks;
+    var ex1 = 1.0;
+    var ex2 = 1.0;
+    var ex1bar = 1.0;
+    var ex2bar = 1.0;
+    var ez = 1.0;
+    var ey = 0.1;
+    var ex = 0.1;
+    var eu = 0.05;
+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);
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.c_z_dx2 =  .5;
+eparams.lambda  =  .7;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
+
+e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
+c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
+c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
+c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
+lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
+disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols))
+disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols))
+
+skipline(2)
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.c_z_dx2 =  .5;
+eparams.lambda  =  .7;
+
+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'));
+c_z_dx2_nls = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
+lambda_nls = M_.params(strmatch('lambda', 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))
+disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls))
+disp(sprintf('Estimate of lambda: %f', lambda_nls))
+
+if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
+end
+
+if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
+end
+
+if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
+end
+
+if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)')
+end
+
+if abs(lambda_nls-lambda_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (lambda)')
+end
\ No newline at end of file
diff --git a/tests/pac/trend-component-24/clean b/tests/pac/trend-component-24/clean
new file mode 100755
index 000000000..cf3492a28
--- /dev/null
+++ b/tests/pac/trend-component-24/clean
@@ -0,0 +1,8 @@
+#!/bin/sh
+
+rm -rf example
+rm  -rf +example
+rm -f example.log
+rm -f *.mat
+rm -f *.m
+rm -f *.dat
diff --git a/tests/pac/trend-component-24/example.mod b/tests/pac/trend-component-24/example.mod
new file mode 100644
index 000000000..e3ff87b5f
--- /dev/null
+++ b/tests/pac/trend-component-24/example.mod
@@ -0,0 +1,155 @@
+// --+ options: json=compute, stochastic +--
+
+var x1 x2 x1bar x2bar z y x u;
+
+varexo ex1 ex2 ex1bar ex2bar ez ey ex eu;
+
+parameters
+       rho_1 rho_2 rho_3 rho_4
+       a_x1_0 a_x1_1 a_x1_2 a_x1_x2_1 a_x1_x2_2
+	   a_x2_0 a_x2_1 a_x2_2 a_x2_x1_1 a_x2_x1_2
+	   e_c_m c_z_1 c_z_2 c_z_dx2 c_z_u beta
+       lambda;
+
+rho_1 =  .9;
+rho_2 = -.2;
+rho_3 =  .4;
+rho_4 = -.3;
+
+
+a_x1_0 =  -.9;
+a_x1_1 =  .4;
+a_x1_2 =  .3;
+a_x1_x2_1 = .1;
+a_x1_x2_2 = .2;
+
+a_x2_0 =  -.9;
+a_x2_1 =   .2;
+a_x2_2 =  -.1;
+a_x2_x1_1 = -.1;
+a_x2_x1_2 = .2;
+
+beta  =  .2;
+e_c_m =  .5;
+c_z_1 =  .2;
+c_z_2 = -.1;
+c_z_dx2 = .3;
+c_z_u = .3;
+
+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);
+
+model;
+
+[name='eq:u']
+u = .5*u(-1) - .2*u(-2) + eu;
+
+[name='eq:y']
+y = rho_1*y(-1) + rho_2*y(-2) + ey;
+
+[name='eq:x']
+x = rho_3*x(-1) + rho_4*x(-2) + ex;
+
+[name='eq:x1']
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1;
+
+[name='eq:x2']
+diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2;
+
+[name='eq:x1bar']
+x1bar = x1bar(-1) + ex1bar;
+
+[name='eq:x2bar']
+x2bar = x2bar(-1) + ex2bar;
+
+[name='zpac']
+diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + c_z_u*u + ez;
+
+end;
+
+shocks;
+    var ex1 = 1.0;
+    var ex2 = 1.0;
+    var ex1bar = 1.0;
+    var ex2bar = 1.0;
+    var ez = 1.0;
+    var ey = 0.1;
+    var ex = 0.1;
+    var eu = 0.05;
+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);
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.lambda  =  .7;
+
+// Define the dataset used for estimation
+edata = TrueData;
+edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
+pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
+
+e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
+c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
+c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
+lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact'));
+
+disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
+disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
+disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
+disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols))
+
+skipline(2)
+
+// Define a structure describing the parameters to be estimated (with initial conditions).
+clear eparams
+eparams.e_c_m   =  .9;
+eparams.c_z_1   =  .5;
+eparams.c_z_2   =  .2;
+eparams.lambda  =  .7;
+
+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'));
+lambda_nls = M_.params(strmatch('lambda', 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))
+disp(sprintf('Estimate of lambda: %f', lambda_nls))
+
+if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
+end
+
+if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
+end
+
+if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
+end
+
+if abs(lambda_nls-lambda_iterative_ols)>.01
+   error('Iterative OLS and NLS do not provide consistent estimates (lambda)')
+end
\ No newline at end of file