diff --git a/matlab/+var_expectation/+update/parameters.m b/matlab/+var_expectation/+update/parameters.m
index b3582e294..8052109ec 100644
--- a/matlab/+var_expectation/+update/parameters.m
+++ b/matlab/+var_expectation/+update/parameters.m
@@ -84,8 +84,17 @@ if discountfactor>1
     error('The discount cannot be greater than one.')
 end
 
-% Set variable_id in VAR model
-variable_id_in_var = strcmp(auxmodel.list_of_variables_in_companion_var, varexpectationmodel.variable);
+% Set variables_id in VAR model
+m = length(varexpectationmodel.expr.vars);
+variables_id_in_var = NaN(m,1);
+for i = 1:m
+    j = find(strcmp(auxmodel.list_of_variables_in_companion_var, DynareModel.endo_names{varexpectationmodel.expr.vars(i)}));
+    if isempty(j)
+        error('Cannot find variable %s in the companion VAR', DynareModel.endo_names{varexpectationmodel.expr.vars(i)})
+    else
+        variables_id_in_var(i) = find(strcmp(auxmodel.list_of_variables_in_companion_var, DynareModel.endo_names{varexpectationmodel.expr.vars(i)}));
+    end
+end
 
 % Get the horizon parameter.
 horizon = varexpectationmodel.horizon;
@@ -127,8 +136,10 @@ CompanionMatrix = auxcalib.CompanionMatrix;
 n = length(CompanionMatrix);
 
 % Set the selection vector
-alpha = zeros(1, length(CompanionMatrix));
-alpha(variable_id_in_var) = 1;
+alpha = zeros(1, n);
+alpha(variables_id_in_var) = varexpectationmodel.expr.constants;
+params_id_in_var = ~isnan(varexpectationmodel.expr.params);
+alpha(variables_id_in_var(params_id_in_var)) = varexpectationmodel.expr.params(params_id_in_var);
 
 if length(horizon)==1
     % Compute the reduced form parameters of the (discounted) forecast in period t+horizon(1)
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a621aca03..9794e1ec7 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -368,6 +368,7 @@ MODFILES = \
     var-expectations/7/substitution.mod \
     var-expectations/8/example.mod \
     var-expectations/8/substitution.mod \
+    var-expectations/9/example.mod \
 	trend-component-and-var-models/vm1.mod \
 	trend-component-and-var-models/vm2.mod \
 	trend-component-and-var-models/vm3.mod \
diff --git a/tests/var-expectations/1/example.mod b/tests/var-expectations/1/example.mod
index 6b78d9acf..0a607e56a 100644
--- a/tests/var-expectations/1/example.mod
+++ b/tests/var-expectations/1/example.mod
@@ -36,7 +36,7 @@ var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 ** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
 */
 
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = x, auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
 
 
 model;
diff --git a/tests/var-expectations/2/example.mod b/tests/var-expectations/2/example.mod
index 236a55085..8becf8d6c 100644
--- a/tests/var-expectations/2/example.mod
+++ b/tests/var-expectations/2/example.mod
@@ -36,7 +36,7 @@ var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 ** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
 */
 
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 2, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = x, auxiliary_model_name = toto, horizon = 2, discount = beta)  ;
 
 
 model;
diff --git a/tests/var-expectations/3/example.mod b/tests/var-expectations/3/example.mod
index 1b65b6df6..6298ab64b 100644
--- a/tests/var-expectations/3/example.mod
+++ b/tests/var-expectations/3/example.mod
@@ -36,7 +36,7 @@ var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 ** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
 */
 
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 0:Inf, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = x, auxiliary_model_name = toto, horizon = 0:Inf, discount = beta)  ;
 
 
 model;
diff --git a/tests/var-expectations/4/example.mod b/tests/var-expectations/4/example.mod
index b05049f9f..2ec4fe30a 100644
--- a/tests/var-expectations/4/example.mod
+++ b/tests/var-expectations/4/example.mod
@@ -36,7 +36,7 @@ var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 ** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
 */
 
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 0:500, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = x, auxiliary_model_name = toto, horizon = 0:500, discount = beta)  ;
 
 
 model;
diff --git a/tests/var-expectations/5/example.mod b/tests/var-expectations/5/example.mod
index 010eceb92..8f088186b 100644
--- a/tests/var-expectations/5/example.mod
+++ b/tests/var-expectations/5/example.mod
@@ -36,7 +36,7 @@ var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 ** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
 */
 
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 15:50, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = x, auxiliary_model_name = toto, horizon = 15:50, discount = beta)  ;
 
 model;
 [ name = 'X' ]
diff --git a/tests/var-expectations/6/example.mod b/tests/var-expectations/6/example.mod
index 3d64bc559..65ba9eff6 100644
--- a/tests/var-expectations/6/example.mod
+++ b/tests/var-expectations/6/example.mod
@@ -47,7 +47,7 @@ trend_component_model(model_name=toto, eqtags=['eq:x1', 'eq:x2', 'eq:x1bar', 'eq
 ** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
 */
 
-var_expectation_model(model_name = varexp, variable = x1, auxiliary_model_name = toto, horizon = 15:50, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = x1, auxiliary_model_name = toto, horizon = 15:50, discount = beta)  ;
 
 
 model;
@@ -81,14 +81,29 @@ end
 
 var_expectation.print('varexp');
 
+shocks;
+  var ex1 = .01;
+  var ex2 = .01;
+  var ex1bar = .02;
+  var ex2bar = .02;
+end;
+
+
 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 = example.set_auxiliary_variables(y, x, M_.params);
-  y = [y(find(M_.lead_lag_incidence(1,:))); y];
-  [residual, g1] = example.dynamic(y, x', M_.params, oo_.steady_state, 1);
-  ynames = M_.endo_names;
-  save('example.mat', 'residual', 'g1');
+  initialconditions =zeros(3,5);
+  initialconditions(3,1) = .10; % foo(-1)
+  initialconditions(3,2) = .20; % x1(-1)
+  initialconditions(2,2) = .22; % x1(-2)
+  initialconditions(1,2) = .24; % x1(-3)
+  initialconditions(3,3) = .30; % x2(-1)
+  initialconditions(2,3) = .32; % x2(-2)
+  initialconditions(1,3) = .34; % x2(-3)
+  initialconditions(3,4) = .25; % x1bar(-1)
+  initialconditions(3,5) = .25; % x2bar(-1)
+  initialconditions = ...
+  dseries(initialconditions, dates('2000Q1'), {'foo', 'x1','x2', 'x1bar', 'x2bar'});
+    set_dynare_seed('default');
+  ts = simul_backward_model(initialconditions, 100);
+  foo = ts.foo.data;
+  save('example.mat', 'foo');
 end;
\ No newline at end of file
diff --git a/tests/var-expectations/6/substitution.mod b/tests/var-expectations/6/substitution.mod
index daf189169..6deae093f 100644
--- a/tests/var-expectations/6/substitution.mod
+++ b/tests/var-expectations/6/substitution.mod
@@ -45,23 +45,32 @@ foo = .5*foo(-1) +
 ;
 end;
 
-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];
-  example = load('example.mat');
-  [residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
+shocks;
+  var ex1 = .01;
+  var ex2 = .01;
+  var ex1bar = .02;
+  var ex2bar = .02;
 end;
 
-if max(abs(example.residual-residual))>1e-8
-  error('Residuals do not match!')
-end
 
-if max(max(abs(example.g1-g1)))>1e-8
-  error('Jacobian matrices do not match!')
-end
-
-delete('example.mat');
\ No newline at end of file
+verbatim;
+  initialconditions =zeros(3,5);
+  initialconditions(3,1) = .10; % foo(-1)
+  initialconditions(3,2) = .20; % x1(-1)
+  initialconditions(2,2) = .22; % x1(-2)
+  initialconditions(1,2) = .24; % x1(-3)
+  initialconditions(3,3) = .30; % x2(-1)
+  initialconditions(2,3) = .32; % x2(-2)
+  initialconditions(1,3) = .34; % x2(-3)
+  initialconditions(3,4) = .25; % x1bar(-1)
+  initialconditions(3,5) = .25; % x2bar(-1)
+  initialconditions = ...
+  dseries(initialconditions, dates('2000Q1'), {'foo', 'x1','x2', 'x1bar', 'x2bar'});
+    set_dynare_seed('default');
+  ts = simul_backward_model(initialconditions, 100);
+  ex = load('example.mat');
+  delete('example.mat')
+  if max(abs(ex.foo-ts.foo.data))>1e-12
+     error('Simulations do not match!')
+  end
+end;
\ No newline at end of file
diff --git a/tests/var-expectations/7/example.mod b/tests/var-expectations/7/example.mod
index ea01f0da8..4aa7b4374 100644
--- a/tests/var-expectations/7/example.mod
+++ b/tests/var-expectations/7/example.mod
@@ -17,7 +17,7 @@ beta = 1/(1+.02);
 var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 
 // Define a VAR_EXPECTATION_MODEL
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = diff(x), auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
 
 
 model;
@@ -40,13 +40,26 @@ var_expectation.update('varexp');
 // Print expanded VAR_EXPECTATION expression in a file (to be included in substitution.mod).
 var_expectation.print('varexp');
 
-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 = example.set_auxiliary_variables(y, x, M_.params);
-  y = [y(find(M_.lead_lag_incidence(1,:))); y];
-  [residual, g1] = example.dynamic(y, x', M_.params, oo_.steady_state, 1);
-  save('example.mat', 'residual', 'g1');
+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, 100);
+  foo = ts.foo.data;
+  save('example.mat', 'foo');
 end;
diff --git a/tests/var-expectations/7/substitution.mod b/tests/var-expectations/7/substitution.mod
index 43002de67..2cfb853f9 100644
--- a/tests/var-expectations/7/substitution.mod
+++ b/tests/var-expectations/7/substitution.mod
@@ -31,23 +31,31 @@ foo = .5*foo(-1) +
 ;
 end;
 
-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];
-  [residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
-  example = load('example.mat');
+shocks;
+  var e_x = .01;
+  var e_y = .01;
+  var e_z = .01;
 end;
 
-if max(abs(example.residual-residual))>1e-8
-  error('Residuals do not match!')
+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, 100);
+  ex = load('example.mat');
+end;
+
+delete('example.mat')
+
+if max(abs(ex.foo-ts.foo.data))>1e-12
+   error('Simulations do not match!')
 end
 
-if max(max(abs(example.g1-g1)))>1e-8
-  error('Jacobian matrices do not match!')
-end
-
-delete('example.mat')
\ No newline at end of file
diff --git a/tests/var-expectations/8/example.mod b/tests/var-expectations/8/example.mod
index 44b81430f..cc14a3e9e 100644
--- a/tests/var-expectations/8/example.mod
+++ b/tests/var-expectations/8/example.mod
@@ -17,7 +17,7 @@ beta = 1/(1+.02);
 var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
 
 // Define a VAR_EXPECTATION_MODEL
-var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+var_expectation_model(model_name = varexp, expression = diff(log(x)), auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
 
 model;
 [ name = 'X' ]
@@ -30,22 +30,34 @@ log(y) = d*log(y(-2)) + e*diff(z(-1)) + e_y;
 foo = .5*foo(-1) + 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;
-  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 = example.set_auxiliary_variables(y, x, M_.params);
-  y = [y(find(M_.lead_lag_incidence(1,:))); y];
-  [residual, g1] = example.dynamic(y, x', M_.params, oo_.steady_state, 1);
-  save('example.mat', 'residual', 'g1');
+  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, 100);
+  foo = ts.foo.data;
+  save('example.mat', 'foo');
 end;
\ No newline at end of file
diff --git a/tests/var-expectations/8/substitution.mod b/tests/var-expectations/8/substitution.mod
index 561d48034..96b9bbb61 100644
--- a/tests/var-expectations/8/substitution.mod
+++ b/tests/var-expectations/8/substitution.mod
@@ -28,23 +28,30 @@ foo = .5*foo(-1) +
 ;
 end;
 
-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];
-  [residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
-  example = load('example.mat');
+shocks;
+  var e_x = .01;
+  var e_y = .01;
+  var e_z = .01;
 end;
 
-if max(abs(example.residual-residual))>1e-8
-  error('Residuals do not match!')
-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, 100);
+  ex = load('example.mat');
+end;
 
-if max(max(abs(example.g1-g1)))>1e-8
-  error('Jacobian matrices do not match!')
-end
+delete('example.mat')
 
-delete('example.mat')
\ No newline at end of file
+if max(abs(ex.foo-ts.foo.data))>1e-12
+   error('Simulations do not match!')
+end
\ No newline at end of file
diff --git a/tests/var-expectations/9/example.mod b/tests/var-expectations/9/example.mod
new file mode 100644
index 000000000..b39bf4e90
--- /dev/null
+++ b/tests/var-expectations/9/example.mod
@@ -0,0 +1,82 @@
+// --+ options: stochastic,transform_unary_ops,json=compute +--
+
+var foo x1 x2 x1bar x2bar;
+
+varexo ex1 ex2 ex1bar ex2bar;
+
+parameters a_x1_0 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
+       beta
+       alpha ;
+
+a_x1_0 =  -.9;
+a_x1_0_ =  -.8;
+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 = 1/(1+.02);
+alpha = .5;
+
+// Define a TREND_COMPONENT model from a subset of equations in the model block.
+trend_component_model(model_name=toto, eqtags=['eq:x1', 'eq:x2', 'eq:x1bar', 'eq:x2bar'], targets=['eq:x1bar', 'eq:x2bar']);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:             the name of the VAR_EXPECTATION_MODEL (mandatory).
+** auxiliary_model_name:   the name of the VAR model used for the expectations (mandatory).
+** variable:               the name of the variable to be forecasted (mandatory).
+** horizon:                the horizon forecast (mandatory).
+** discount:               the discount factor, which can be a value or a declared parameter (default is 1.0, no discounting).
+**
+**
+** The `horizon` parameter can be an integer in which case the (discounted) `horizon` step ahead forecast
+** is computed using the VAR model `var_model_name`. Alternatively, `horizon` can be a range. In this case
+** VAR_EXPECTATION_MODEL returns a discounted sum of expected values. If `horizon` is set equal to the range
+** 0:Inf, then VAR_EXPECTATION_MODEL computes:
+**
+**                                       ∑ βʰ Eₜ[yₜ₊ₕ]
+**
+** where the sum is over h=0,…,∞ and the conditional expectations are computed with VAR model `var_model_name`.
+*/
+
+var_expectation_model(model_name = varexp, expression = x1/100-alpha*x2, auxiliary_model_name = toto, horizon = 15:50, discount = beta)  ;
+
+
+model;
+
+[name='eq:x1', data_type='nonstationary']
+diff(x1) = a_x1_0*(x1(-1)-x1bar(-1))+a_x1_0_*(x2(-1)-x2bar(-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', data_type='nonstationary']
+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', data_type='nonstationary']
+x1bar = x1bar(-1) + ex1bar;
+
+[name='eq:x2bar', data_type='nonstationary']
+x2bar = x2bar(-1) + ex2bar;
+
+foo = .5*foo(-1) + 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');
+
+weights = M_.params(M_.var_expectation.varexp.param_indices);
+
+if ~all(weights(1:6)) || ~all(weights(9:10)) || weights(7) || weights(8) || weights(11) || weights(12)
+   error('Wrong reduced form parameter for VAR_EXPECTATION_MODEL')
+end
\ No newline at end of file