diff --git a/matlab/get_companion_matrix.m b/matlab/get_companion_matrix.m
index 04dfcdf93..6fc137b8d 100644
--- a/matlab/get_companion_matrix.m
+++ b/matlab/get_companion_matrix.m
@@ -42,9 +42,9 @@ if nargin < 2
 end
 
 if strcmp(auxiliary_model_type, 'var')
-    AR = feval([M_.fname '.var_ar'], auxiliary_model_name, M_.params);
+    [AR, ~] = feval(sprintf('%s.varmatrices', M_.fname), auxiliary_model_name, M_.params, M_.var.(auxiliary_model_name).structural);
 elseif strcmp(auxiliary_model_type, 'trend_component')
-    [AR, A0, A0star] = feval([M_.fname '.trend_component_ar_a0'], auxiliary_model_name, M_.params);
+    [AR, A0, A0star] = feval(sprintf('%s.trend_component_ar_a0', M_.fname), auxiliary_model_name, M_.params);
 else
     error('Unknown type of auxiliary model.')
 end
diff --git a/preprocessor b/preprocessor
index c9c36c037..6b9d94405 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit c9c36c037aa6d6c67f8fbb2236b03729a22e5878
+Subproject commit 6b9d94405c44fe2533a2d8bc377173a2a70f3a54
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 8866b9645..34922721b 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -461,6 +461,8 @@ ECB_MODFILES = \
 	var-expectations/8/substitution.mod \
 	var-expectations/9/example1.mod \
 	var-expectations/10/example1.mod \
+	var-expectations/11/example1.mod \
+	var-expectations/12/example1.mod \
 	trend-component-and-var-models/vm1.mod \
 	trend-component-and-var-models/vm2.mod \
 	trend-component-and-var-models/vm3.mod \
@@ -494,6 +496,8 @@ ECB_MODFILES = \
 	pac/var-5/substitution.mod \
 	pac/var-6/example1.mod \
 	pac/var-6/substitution.mod \
+	pac/var-7/example1.mod \
+	pac/var-7/substitution.mod \
 	pac/trend-component-1/example1.mod \
 	pac/trend-component-2/example1.mod \
 	pac/trend-component-3/example1.mod \
@@ -842,6 +846,8 @@ pac/var-5/substitution.m.trs: pac/var-5/example1.m.trs
 pac/var-5/substitution.o.trs: pac/var-5/example1.o.trs
 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/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/pac/run_all_tests.m b/tests/pac/run_all_tests.m
index 32717f04b..f98e12ddc 100644
--- a/tests/pac/run_all_tests.m
+++ b/tests/pac/run_all_tests.m
@@ -1,6 +1,6 @@
 function run_all_tests()
 
-% Copyright (C) 2018 Dynare Team
+% Copyright © 2018-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -23,6 +23,9 @@ r = [r; run_this_test('var-1')];
 r = [r; run_this_test('var-2')];
 r = [r; run_this_test('var-3')];
 r = [r; run_this_test('var-4')];
+r = [r; run_this_test('var-5')];
+r = [r; run_this_test('var-6')];
+r = [r; run_this_test('var-7')];
 
 r = [r; run_this_test('trend-component-1')];
 r = [r; run_this_test('trend-component-2')];
@@ -31,6 +34,8 @@ r = [r; run_this_test('trend-component-4')];
 r = [r; run_this_test('trend-component-5')];
 r = [r; run_this_test('trend-component-6')];
 r = [r; run_this_test('trend-component-7')];
+r = [r; run_this_test('trend-component-8-mc-iterative-ols')];
+r = [r; run_this_test('trend-component-8-mc-nls')];
 r = [r; run_this_test('trend-component-9')];
 r = [r; run_this_test('trend-component-10')];
 r = [r; run_this_test('trend-component-11')];
@@ -40,7 +45,25 @@ r = [r; run_this_test('trend-component-13b')];
 r = [r; run_this_test('trend-component-14')];
 r = [r; run_this_test('trend-component-15')];
 r = [r; run_this_test('trend-component-16')];
+r = [r; run_this_test('trend-component-17')];
+r = [r; run_this_test('trend-component-18')];
+r = [r; run_this_test('trend-component-19')];
 r = [r; run_this_test('trend-component-19-growth-lin-comb')];
+r = [r; run_this_test('trend-component-20-1')];
+r = [r; run_this_test('trend-component-20-2')];
+r = [r; run_this_test('trend-component-20-3')];
+r = [r; run_this_test('trend-component-20-4')];
+r = [r; run_this_test('trend-component-21')];
+r = [r; run_this_test('trend-component-22')];
+r = [r; run_this_test('trend-component-23')];
+r = [r; run_this_test('trend-component-24')];
+r = [r; run_this_test('trend-component-25')];
+r = [r; run_this_test('trend-component-26')];
+r = [r; run_this_test('trend-component-27')];
+r = [r; run_this_test('trend-component-28')];
+r = [r; run_this_test('trend-component-29')];
+r = [r; run_this_test('trend-component-30')];
+r = [r; run_this_test('trend-component-31')];
 
 print_results(r);
 
diff --git a/tests/pac/var-7/clean b/tests/pac/var-7/clean
new file mode 100755
index 000000000..1c4fb0aab
--- /dev/null
+++ b/tests/pac/var-7/clean
@@ -0,0 +1,10 @@
+#!/bin/sh
+
+rm -rf example1
+rm -rf +example1
+rm -f example1*.mat
+rm -f example1.log
+rm -rf substitution
+rm -rf +substitution
+rm -f substitution*.mat
+rm -f substitution.log
diff --git a/tests/pac/var-7/example1.mod b/tests/pac/var-7/example1.mod
new file mode 100644
index 000000000..ccda33575
--- /dev/null
+++ b/tests/pac/var-7/example1.mod
@@ -0,0 +1,79 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z u;
+
+varexo ex ey ez eu;
+
+parameters a_y_0 a_y_1 a_y_2 b_y_1 b_y_2 b_x_0 b_x_1 b_x_2 ; // VAR parameters
+
+parameters g beta e_c_m c_z_1 c_z_2;                         // PAC equation parameters
+
+a_y_0 =  .4;
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_0 = -.1;
+b_x_1 = -.1;
+b_x_2 = -.2;
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+g = .1;
+
+var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman, growth=diff(u(-1)));
+
+model;
+
+[name='eq:u']
+diff(u) = g + eu;
+
+[name='eq:y']
+y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + a_y_0*diff(x) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+[name='eq:x']
+diff(x) = b_x_0*y + b_x_1*y(-2) + b_x_2*diff(x(-1)) + g*(1-b_x_2)  + ex ;
+
+[name='eq:pac']
+diff(z) = e_c_m*(x(-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;
+    var eu = 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 20 periods
+set_dynare_seed('default');
+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-7/substitution.mod b/tests/pac/var-7/substitution.mod
new file mode 100644
index 000000000..4597649c0
--- /dev/null
+++ b/tests/pac/var-7/substitution.mod
@@ -0,0 +1,88 @@
+// --+ options: transform_unary_ops, json=compute, stochastic +--
+
+var y x z u;
+
+varexo ex ey ez eu;
+
+parameters a_y_0 a_y_1 a_y_2 b_y_1 b_y_2 b_x_0 b_x_1 b_x_2 g; // VAR parameters
+
+parameters e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_0 =  .4;
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_0 = -.1;
+b_x_1 = -.1;
+b_x_2 = -.2;
+
+beta  =  .9;
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+g = .1;
+
+@#include "example1/model/pac-expectations/eq0-pacman-parameters.inc"
+
+model;
+
+  [name='eq:u']
+  diff(u) = g + eu;
+
+  [name='eq:y']
+  y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + a_y_0*diff(x) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+  [name='eq:x']
+  diff(x) = b_x_0*y + b_x_1*y(-2) + b_x_2*diff(x(-1)) + g*(1-b_x_2)  + ex ;
+
+  [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/eq0-pacman-growth-neutrality-correction.inc"
+  +
+  @#include "example1/model/pac-expectations/eq0-pacman-expression.inc"
+  + ez;
+
+end;
+
+shocks;
+  var ex = 1.0;
+  var ey = 1.0;
+  var ez = 1.0;
+  var eu = 0.1;
+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
+
+  delete('example1.mat');
+
+end;
diff --git a/tests/var-expectations/11/example1.mod b/tests/var-expectations/11/example1.mod
new file mode 100644
index 000000000..d5c21ae90
--- /dev/null
+++ b/tests/var-expectations/11/example1.mod
@@ -0,0 +1,46 @@
+// --+ options: stochastic,json=compute +--
+
+var foo z x y;
+varexo e_x e_y e_z;
+parameters a b c d e f beta ;
+
+a =  .9;
+b = -.2;
+c =  .3;
+f =  .8;
+d =  .5;
+e =  .4;
+
+beta = 1/(1+.02);
+
+// Define a VAR model from a subset of equations in the model block.
+var_model(structural, model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+// Define a VAR_EXPECTATION_MODEL
+var_expectation_model(model_name = varexp, expression = diff(log(x)), auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+
+model;
+[ name = 'X' ]
+diff(log(x)) = b*diff(z) + a*diff(log(x(-1))) + (1-a)*diff(log(x(-2))) + c*diff(z(-2)) + e_x;
+[ name = 'Z' ]
+diff(z) = f*(diff(z(-1))-diff(log(x)))+c*diff(z(-2)) + e_z;
+[ name = 'Y' ]
+log(y) = diff(log(x)) + d*log(y(-2)) + e*diff(z(-1)) + e_y;
+
+foo = var_expectation(varexp);
+end;
+
+[ar, a0] = example1.varmatrices('toto', M_.params);
+
+assert(isequal(diag(a0), ones(3,1)), 'Diagonal of a0 is wrong.')
+
+assert(a0(1,3)==-b, 'Element (1,3) in A0 is wrong.')
+assert(a0(1,2)==0, 'Element (1,2) in A0 is wrong.')
+assert(a0(2,1)==-1, 'Element (2,1) in A0 is wrong.')
+assert(a0(2,3)==0, 'Element (2,3) in A0 is wrong.')
+assert(a0(3,1)==f, 'Element (3,1) in A0 is wrong.')
+assert(a0(3,2)==0, 'Element (3,1) in A0 is wrong.')
+
+assert(isequal(ar(:,:,1), [a 0 0; 0 0 e; 0 0 f]), 'First autoregressive matrix is wrong');
+assert(isequal(ar(:,:,2), [1-a 0 c; 0 d 0; 0 0 c]), 'Second autoregressive matrix is wrong');
+
diff --git a/tests/var-expectations/12/example1.mod b/tests/var-expectations/12/example1.mod
new file mode 100644
index 000000000..a46ab788b
--- /dev/null
+++ b/tests/var-expectations/12/example1.mod
@@ -0,0 +1,69 @@
+// --+ options: stochastic,json=compute +--
+
+var foo z x y;
+varexo e_x e_y e_z;
+parameters a b c d e f beta ;
+
+a =  .9;
+b = -.2;
+c =  .3;
+f =  .8;
+d =  .5;
+e =  .4;
+
+beta = 1/(1+.02);
+
+// Define a VAR model from a subset of equations in the model block.
+var_model(structural, model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+// Define a VAR_EXPECTATION_MODEL
+var_expectation_model(model_name = varexp, expression = diff(log(x)), auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+
+model;
+[ name = 'X' ]
+diff(log(x)) = b*diff(z) + a*diff(log(x(-1))) + (1-a)*diff(log(x(-2))) + c*diff(z(-2)) + e_x;
+[ name = 'Z' ]
+diff(z) = f*(diff(z(-1))-diff(log(x)))+c*diff(z(-2)) + e_z;
+[ name = 'Y' ]
+log(y) = diff(log(x)) + d*log(y(-2)) + e*diff(z(-1)) + e_y;
+
+foo = var_expectation(varexp);
+end;
+
+// Initialize the VAR expectation model, will build the companion matrix of the VAR.
+var_expectation.initialize('varexp')
+
+// Update VAR_EXPECTATION reduced form parameters
+var_expectation.update('varexp');
+
+// Print expanded VAR_EXPECTATION expression in a file (to be included in substitution.mod).
+var_expectation.print('varexp');
+
+shocks;
+  var e_x = .01;
+  var e_y = .01;
+  var e_z = .01;
+end;
+
+verbatim;
+  initialconditions =zeros(3,4);
+  initialconditions(3,1) = .1; % foo(-1)
+  initialconditions(:,2) = .2; % y(-1)
+  initialconditions(3,3) = .3; % z(-1)
+  initialconditions(2,3) = .4; % z(-2)
+  initialconditions(3,4) = .5; % x(-1)
+  initialconditions(2,4) = .6; % x(-2)
+  initialconditions(1,4) = .7; % x(-3)
+  initialconditions = ...
+  dseries(initialconditions, dates('2000Q1'), {'foo', 'y','z', 'x'});
+  set_dynare_seed('default');
+  ts = simul_backward_model(initialconditions, 15);
+  foo = ts.foo.data;
+  % Evaluate the (VAR) expectation term
+  ts{'toto'} = example1.var_expectations.varexp.evaluate(ts);
+  % Check that the evaluation is correct.
+  range = dates('2000Q4'):dates('2004Q2');
+  if max(abs(ts(range).foo.data-ts(range).toto.data))>1e-5
+     error('Expectation term evaluations do not match!')
+  end
+end;