diff --git a/matlab/+var_expectation/update_parameters.m b/matlab/+var_expectation/update_parameters.m
index e947c493cd8c164dbfcf3e898f3a272fa6501e3c..c57d5eed382620a4996daa23422d421213850b8d 100644
--- a/matlab/+var_expectation/update_parameters.m
+++ b/matlab/+var_expectation/update_parameters.m
@@ -84,6 +84,18 @@ if discountfactor>1
     error('The discount cannot be greater than one.')
 end
 
+% Set timeshift
+if isfield(varexpectationmodel, 'time_shift')
+    timeshift = varexpectationmodel.time_shift;
+else
+    timeshift = 0;
+end
+
+% Throw an error if timeshift is positive
+if timeshift>0
+    error('Option time_shift of the var_expectation command cannot be positive.')
+end
+
 % Set variables_id in VAR model
 m = length(varexpectationmodel.expr.vars);
 variables_id_in_var = NaN(m,1);
@@ -148,11 +160,17 @@ if length(horizon)==1
     elseif horizon>1 
         parameters = alpha*mpower(discountfactor*CompanionMatrix, varexpectationmodel.horizon);
     end
+    if timeshift<0
+        parameters = parameters*mpower(CompanionMatrix, -timeshift);
+    end
 else
     % Compute the reduced form parameters of the discounted sum of forecasts between t+horizon(1) and
     % t+horizon(2). Not that horzizon(2) need not be finite.
     if horizon(1)==0 && isinf(horizon(2))
         parameters = alpha/(eye(n)-discountfactor*CompanionMatrix);
+        if timeshift<0
+            parameters = parameters*mpower(CompanionMatrix, -timeshift);
+        end
     elseif horizon(1)>0 && isinf(horizon(2))
         % Define the discounted companion matrix
         DiscountedCompanionMatrix = discountfactor*CompanionMatrix;
@@ -162,8 +180,14 @@ else
             tmp1 = tmp1 + mpower(DiscountedCompanionMatrix, h); 
         end
         tmp1 = alpha*tmp1;
+        if timeshift<0
+            tmp1 = tmp1*mpower(CompanionMatrix, -timeshift);
+        end
         % Second compute the parameters implied by the discounted sum from h=0 to h=Inf 
         tmp2 = alpha/(eye(n)-DiscountedCompanionMatrix);
+        if timeshift<0
+            tmp2 = tmp2*mpower(CompanionMatrix, -timeshift);
+        end
         % Finally
         parameters = tmp2-tmp1;
     elseif isfinite(horizon(2))
@@ -174,6 +198,9 @@ else
             tmp = tmp + mpower(DiscountedCompanionMatrix, h);
         end
         parameters = alpha*tmp;
+        if timeshift<0
+            parameters = parameters*mpower(CompanionMatrix, -timeshift);
+        end
     end
 end
 
diff --git a/matlab/print_equations.m b/matlab/print_equations.m
index 0bf9b60f8117d0101ca77dc1d80dcbca80c2cef3..7d8069da2b7b0e56948b2a56932dc0deca71cff5 100644
--- a/matlab/print_equations.m
+++ b/matlab/print_equations.m
@@ -1,4 +1,4 @@
-function print_equations(variable_name, withexpansion)
+function str = print_equations(variable_name, withexpansion)
 
 % Prints equations where the variable appears in.
 %
@@ -84,6 +84,9 @@ for it = 1:length(M_.mapping.(variable_name).eqidx)
             withexpansion = false;
         end
     end
+    if nargout
+        str = sprintf('%s = %s;\n', model{M_.mapping.(variable_name).eqidx(it)}.lhs, rhs);
+    end
     fprintf('%s = %s;\n', model{M_.mapping.(variable_name).eqidx(it)}.lhs, rhs);
 end
 
diff --git a/matlab/write_expectations.m b/matlab/write_expectations.m
index 0b04e3237b7fb9014b559313fe5a0f68ed782759..831bea27db4508640d562e6072ccbb48ef644bb0 100644
--- a/matlab/write_expectations.m
+++ b/matlab/write_expectations.m
@@ -74,6 +74,10 @@ end
 
 id = 0;
 
+if isequal(expectationmodelkind, 'var')
+    timeindices = (0:(maxlag-1))+abs(expectationmodel.time_shift);
+end
+
 for i=1:maxlag
     for j=1:length(auxmodel.list_of_variables_in_companion_var)
         id = id+1;
@@ -114,8 +118,8 @@ for i=1:maxlag
         end
         switch expectationmodelkind
           case 'var'
-            if i>1
-                variable = sprintf('%s(-%d)', variable, i-1);
+            if timeindices(i)
+                variable = sprintf('%s(-%d)', variable, timeindices(i));
             end
           case 'pac'
             variable = sprintf('%s(-%d)', variable, i);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index dfaee4eeab93f59fff54c3f6d1ee6880b79e8a8f..3ca8a269467432560b895b5476d4861896856146 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -440,12 +440,19 @@ MODFILES = \
 
 ECB_MODFILES = \
 	var-expectations/1/example1.mod \
+    var-expectations/1-with-time-shift-a/example1.mod \
 	var-expectations/2/example1.mod \
+	var-expectations/2-with-time-shift/example1.mod \
 	var-expectations/3/example1.mod \
+	var-expectations/3-with-time-shift/example1.mod \
 	var-expectations/4/example1.mod \
+	var-expectations/4-with-time-shift/example1.mod \
 	var-expectations/5/example1.mod \
+	var-expectations/5-with-time-shift/example1.mod \
 	var-expectations/6/example1.mod \
 	var-expectations/6/substitution.mod \
+	var-expectations/6-with-time-shift/example1.mod \
+	var-expectations/6-with-time-shift/substitution.mod \
 	var-expectations/7/example1.mod \
 	var-expectations/7/substitution.mod \
 	var-expectations/8/example1.mod \
@@ -588,7 +595,8 @@ XFAIL_MODFILES = ramst_xfail.mod \
 	estimation/tune_mh_jscale/fs2000_1_xfail.mod \
 	estimation/tune_mh_jscale/fs2000_2_xfail.mod \
 	estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL.mod \
-	estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL2.mod
+	estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL2.mod \
+	var-expectations/1-with-time-shift-b/example1.mod
 
 MFILES = histval_initval_file/ramst_initval_file_data.m
 
@@ -821,6 +829,8 @@ kalman/likelihood_from_dynare/fs2000ns_corr_ME_missing.o.trs: kalman/likelihood_
 
 var-expectations/6/substitution.m.trs: var-expectations/6/example1.m.trs
 var-expectations/6/substitution.o.trs: var-expectations/6/example1.o.trs
+var-expectations/6-with-time-shift/substitution.m.trs: var-expectations/6-with-time-shift/example1.m.trs
+var-expectations/6-with-time-shift/substitution.o.trs: var-expectations/6-with-time-shift/example1.o.trs
 var-expectations/7/substitution.m.trs: var-expectations/7/example1.m.trs
 var-expectations/7/substitution.o.trs: var-expectations/7/example1.o.trs
 var-expectations/8/substitution.m.trs: var-expectations/8/example1.m.trs
@@ -839,6 +849,9 @@ lmmcp/sw_newton.o.trs: lmmcp/sw_lmmcp.o.trs
 var-expectations/4/example1.m.trs: var-expectations/3/example1.m.trs
 var-expectations/4/example1.o.trs: var-expectations/3/example1.o.trs
 
+var-expectations/4-with-time-shift/example1.m.trs: var-expectations/3-with-time-shift/example1.m.trs
+var-expectations/4-with-time-shift/example1.o.trs: var-expectations/3-with-time-shift/example1.o.trs
+
 trend-component-and-var-models/tcm7.m.trs: trend-component-and-var-models/tcm6.m.trs
 trend-component-and-var-models/tcm7.o.trs: trend-component-and-var-models/tcm6.o.trs
 trend-component-and-var-models/tcm8.m.trs: trend-component-and-var-models/tcm6.m.trs
diff --git a/tests/var-expectations/1-with-time-shift-a/example1.mod b/tests/var-expectations/1-with-time-shift-a/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..c13808424be9fb696fc44a40cd7b97dcdbca94a9
--- /dev/null
+++ b/tests/var-expectations/1-with-time-shift-a/example1.mod
@@ -0,0 +1,64 @@
+// --+ 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(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:       the name of the VAR_EXPECTATION_MODEL (mandatory).
+** var_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).
+** time_shift:       shifts the information set to the past, must be a non positive scalar. By default, expectations
+**                   about `variable` in period `t+horizon` are formed in period t (time_shift=0)
+**
+**
+** 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 = x, auxiliary_model_name = toto, horizon = 1, time_shift=-1, discount = beta)  ;
+
+
+model;
+[ name = 'X' ]
+x = a*x(-1) + b*x(-2) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+y = d*y(-2) + e*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');
+
+if ~isfield(M_.var_expectation.varexp, 'time_shift') || ~isequal(M_.var_expectation.varexp.time_shift, -1)
+    error('Preprocessor does not honour time_shift option as expected.')
+end
diff --git a/tests/var-expectations/1-with-time-shift-b/example1.mod b/tests/var-expectations/1-with-time-shift-b/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..198d520f64808e9e52c5a59b46bb2a3ad25b4e81
--- /dev/null
+++ b/tests/var-expectations/1-with-time-shift-b/example1.mod
@@ -0,0 +1,66 @@
+// --+ 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(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:       the name of the VAR_EXPECTATION_MODEL (mandatory).
+** var_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).
+** time_shift:       shifts the information set to the past, must be a non positive scalar. By default, expectations 
+**                   about `variable` in period `t+horizon` are formed in period t (time_shift=0) 
+**
+**
+** 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 = x, auxiliary_model_name = toto, horizon = 1, time_shift=1, discount = beta)  ;
+
+
+model;
+[ name = 'X' ]
+x = a*x(-1) + b*x(-2) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+y = d*y(-2) + e*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
+try
+    var_expectation.update('varexp');
+    if isfield(M_.var_expectation.varexp, 'time_shift')
+        error('Preprocessor should not allow positive values for time_shift option.')
+    end
+catch
+end
diff --git a/tests/var-expectations/2-with-time-shift/example1.mod b/tests/var-expectations/2-with-time-shift/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..c2c91d02c108295b81e01d4ed48cf96821cb22bd
--- /dev/null
+++ b/tests/var-expectations/2-with-time-shift/example1.mod
@@ -0,0 +1,86 @@
+// --+ 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(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:       the name of the VAR_EXPECTATION_MODEL (mandatory).
+** var_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).
+** time_shift:       shifts the information set to the past, must be a non positive scalar. By default, expectations
+**                   about `variable` in period `t+horizon` are formed in period t (time_shift=0)
+**
+**
+** 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 = x, auxiliary_model_name = toto, time_shift=-2, horizon = 2, discount = beta)  ;
+
+
+model;
+[ name = 'X' ]
+x = a*x(-1) + b*x(-2) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+y = d*y(-2) + e*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 equations where the variable appears in
+fprintf('x is in: \n')
+print_equations('x')
+fprintf('\n')
+
+fprintf('y is in: \n')
+str = print_equations('y', true);
+fprintf('\n')
+
+
+if ~isfield(M_.var_expectation.varexp, 'time_shift') || ~isequal(M_.var_expectation.varexp.time_shift, -2)
+    error('Preprocessor does not honour time_shift option as expected.')
+end
+
+str = strrep(str, 'foo = .5*foo(-1)', '');
+str = strrep(str, '+ var_expectation_model_varexp_x_0*x(-2)', '');
+str = strrep(str, '+ var_expectation_model_varexp_y_0*y(-2)', '');
+str = strrep(str, '+ var_expectation_model_varexp_z_0*z(-2)', '');
+str = strrep(str, '+ var_expectation_model_varexp_x_1*x(-3)', '');
+str = strrep(str, '+ var_expectation_model_varexp_y_1*y(-3)', '');
+str = strrep(str, '+ var_expectation_model_varexp_z_1*z(-3)', '');
+str = strrep(str, ';', '');
+if ~isempty(strtrim(str))
+    error('Printed equation is wrong.')
+end
diff --git a/tests/var-expectations/3-with-time-shift/example1.mod b/tests/var-expectations/3-with-time-shift/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..952599ec1bdfc07ffd4b61aa9a3a409bd76f96ca
--- /dev/null
+++ b/tests/var-expectations/3-with-time-shift/example1.mod
@@ -0,0 +1,72 @@
+// --+ 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(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:       the name of the VAR_EXPECTATION_MODEL (mandatory).
+** var_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).
+** time_shift:       shifts the information set to the past, must be a non positive scalar. By default, expectations
+**                   about `variable` in period `t+horizon` are formed in period t (time_shift=0)
+**
+** 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 = x, auxiliary_model_name = toto, time_shift=-1, horizon = 0:Inf, discount = beta)  ;
+
+
+model;
+[ name = 'X' ]
+x = a*x(-1) + b*x(-2) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+y = d*y(-2) + e*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');
+
+/*
+** REMARK The VAR model is such that x depends on past values of x
+** (xₜ₋₁ and xₜ₋₂) and on zₜ₋₂. Consequently the reduced
+** form parameters associated to yₜ₋₁, yₜ₋₂ have to be zero.
+*/
+
+weights = M_.params(M_.var_expectation.varexp.param_indices);
+
+if weights(2) || ~weights(3) || weights(5) || ~weights(1) || ~weights(4) || ~weights(6)
+   error('Wrong reduced form parameter for VAR_EXPECTATION_MODEL')
+end
+
+save('weights.mat','weights');
diff --git a/tests/var-expectations/3/example1.mod b/tests/var-expectations/3/example1.mod
index 6298ab64bb138190b8d36e346c4f6f07aeda39df..e85c31f23a01778ca6a657d5065a00e932b25b45 100644
--- a/tests/var-expectations/3/example1.mod
+++ b/tests/var-expectations/3/example1.mod
@@ -68,4 +68,4 @@ if weights(2) || ~weights(3) || weights(5) || ~weights(1) || ~weights(4) || ~wei
    error('Wrong reduced form parameter for VAR_EXPECTATION_MODEL')
 end
 
-save('weights.mat','weights');
\ No newline at end of file
+save('weights.mat','weights');
diff --git a/tests/var-expectations/4-with-time-shift/example1.mod b/tests/var-expectations/4-with-time-shift/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..b78c8a22f512193f74fbce9cfff9a1b3579d80a9
--- /dev/null
+++ b/tests/var-expectations/4-with-time-shift/example1.mod
@@ -0,0 +1,78 @@
+// --+ 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(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:       the name of the VAR_EXPECTATION_MODEL (mandatory).
+** var_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).
+** time_shift:       shifts the information set to the past, must be a non positive scalar. By default, expectations
+**                   about `variable` in period `t+horizon` are formed in period t (time_shift=0)
+**
+** 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 = x, auxiliary_model_name = toto, time_shift=-1, horizon = 0:500, discount = beta)  ;
+
+
+model;
+[ name = 'X' ]
+x = a*x(-1) + b*x(-2) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+y = d*y(-2) + e*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');
+
+/*
+** REMARK The VAR model is such that x depends on past values of x
+** (xₜ₋₁ and xₜ₋₂) and on zₜ₋₂. Consequently the reduced
+** form parameters associated to yₜ₋₁, yₜ₋₂ have to be zero.
+*/
+
+weights = M_.params(M_.var_expectation.varexp.param_indices);
+
+if weights(2) || ~weights(3) || weights(5) || ~weights(1) || ~weights(4) || ~weights(6)
+   error('Wrong reduced form parameter for VAR_EXPECTATION_MODEL')
+end
+
+WEIGHTS = load('../3-with-time-shift/weights.mat');
+
+if max(abs(weights-WEIGHTS.weights))>1e-12
+   error('Inconsistent results in var-expectations/3-with-time-shift and var-expectations/4-with-time-shift.')
+end
+
+delete('../3-with-time-shift/weights.mat')
diff --git a/tests/var-expectations/5-with-time-shift/example1.mod b/tests/var-expectations/5-with-time-shift/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..c527e9da6ec302f89697240caf18d295138067d4
--- /dev/null
+++ b/tests/var-expectations/5-with-time-shift/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(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+/* Define a VAR_EXPECTATION_MODEL
+** ------------------------------
+**
+** model_name:       the name of the VAR_EXPECTATION_MODEL (mandatory).
+** var_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).
+** time_shift:       shifts the information set to the past, must be a non positive scalar. By default, expectations
+**                   about `variable` in period `t+horizon` are formed in period t (time_shift=0)
+**
+** 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 = x, auxiliary_model_name = toto, time_shift=-2, horizon = 15:50, discount = beta)  ;
+
+model;
+[ name = 'X' ]
+x = a*x(-1) + b*x(-2) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+y = d*y(-2) + e*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');
+
+/*
+** REMARK The VAR model is such that x depends on past values of x
+** (xₜ₋₁ and xₜ₋₂) and on zₜ₋₂. Consequently the reduced
+** form parameters associated to yₜ₋₁, yₜ₋₂ have to be zero.
+*/
+
+weights = M_.params(M_.var_expectation.varexp.param_indices);
+
+if weights(2) || ~weights(3) || weights(5) || ~weights(1) || ~weights(4) || ~weights(6)
+   error('Wrong reduced form parameter for VAR_EXPECTATION_MODEL')
+end
diff --git a/tests/var-expectations/6-with-time-shift/example1.mod b/tests/var-expectations/6-with-time-shift/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..359c3efc6b7953c86a48516c44b0442ec31b1b3f
--- /dev/null
+++ b/tests/var-expectations/6-with-time-shift/example1.mod
@@ -0,0 +1,109 @@
+// --+ 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 ;
+
+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);
+
+// 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, auxiliary_model_name = toto, time_shift=-1, 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
+
+var_expectation.print('varexp');
+
+shocks;
+  var ex1 = .01;
+  var ex2 = .01;
+  var ex1bar = .02;
+  var ex2bar = .02;
+end;
+
+
+verbatim;
+  initialconditions =zeros(4,5);
+  initialconditions(4,1) = .10; % foo(-1)
+  initialconditions(4,2) = .20; % x1(-1)
+  initialconditions(3,2) = .22; % x1(-2)
+  initialconditions(2,2) = .24; % x1(-3)
+  initialconditions(4,3) = .30; % x2(-1)
+  initialconditions(3,3) = .32; % x2(-2)
+  initialconditions(2,3) = .34; % x2(-3)
+  initialconditions(4,4) = .25; % x1bar(-1)
+  initialconditions(4,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('example1.mat', 'foo');
+end;
diff --git a/tests/var-expectations/6-with-time-shift/substitution.mod b/tests/var-expectations/6-with-time-shift/substitution.mod
new file mode 100644
index 0000000000000000000000000000000000000000..277819822ba98c9d28ded365b16c6d5581254b85
--- /dev/null
+++ b/tests/var-expectations/6-with-time-shift/substitution.mod
@@ -0,0 +1,76 @@
+// --+ 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 ;
+
+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;
+
+@#include "example1/model/var-expectations/varexp-parameters.inc"
+
+beta = 1/(1+.02);
+
+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) +
+@#include "example1/model/var-expectations/varexp-expression.inc"
+;
+end;
+
+shocks;
+  var ex1 = .01;
+  var ex2 = .01;
+  var ex1bar = .02;
+  var ex2bar = .02;
+end;
+
+
+verbatim;
+  initialconditions =zeros(4,5);
+  initialconditions(4,1) = .10; % foo(-1)
+  initialconditions(4,2) = .20; % x1(-1)
+  initialconditions(3,2) = .22; % x1(-2)
+  initialconditions(2,2) = .24; % x1(-3)
+  initialconditions(4,3) = .30; % x2(-1)
+  initialconditions(3,3) = .32; % x2(-2)
+  initialconditions(2,3) = .34; % x2(-3)
+  initialconditions(4,4) = .25; % x1bar(-1)
+  initialconditions(4,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('example1.mat');
+  delete('example1.mat')
+  if max(abs(ex.foo-ts.foo.data))>1e-12
+     error('Simulations do not match!')
+  end
+end;