From d96740c2ef3167b2e4ab6f913bd6efaf782be8a0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemia=20=28Scylla=29?= <stepan@dynare.org>
Date: Sun, 14 Oct 2018 17:01:02 +0200
Subject: [PATCH] Added integration tests (var and pac expectations).

---
 tests/Makefile.am                             | 22 +++++
 tests/pac/trend-component-14/clean            |  9 ++
 tests/pac/trend-component-14/example.mod      | 93 ++++++++++++++++++
 tests/pac/trend-component-14/substitution.mod | 98 +++++++++++++++++++
 tests/pac/var-5/clean                         | 10 ++
 tests/pac/var-5/example.mod                   | 73 ++++++++++++++
 tests/pac/var-5/substitution.mod              | 77 +++++++++++++++
 tests/var-expectations/6/example.mod          | 18 +++-
 tests/var-expectations/6/substitution.mod     | 67 +++++++++++++
 tests/var-expectations/7/example.mod          | 52 ++++++++++
 tests/var-expectations/7/substitution.mod     | 53 ++++++++++
 tests/var-expectations/8/example.mod          | 51 ++++++++++
 tests/var-expectations/8/substitution.mod     | 50 ++++++++++
 13 files changed, 671 insertions(+), 2 deletions(-)
 create mode 100755 tests/pac/trend-component-14/clean
 create mode 100644 tests/pac/trend-component-14/example.mod
 create mode 100644 tests/pac/trend-component-14/substitution.mod
 create mode 100755 tests/pac/var-5/clean
 create mode 100644 tests/pac/var-5/example.mod
 create mode 100644 tests/pac/var-5/substitution.mod
 create mode 100644 tests/var-expectations/6/substitution.mod
 create mode 100644 tests/var-expectations/7/example.mod
 create mode 100644 tests/var-expectations/7/substitution.mod
 create mode 100644 tests/var-expectations/8/example.mod
 create mode 100644 tests/var-expectations/8/substitution.mod

diff --git a/tests/Makefile.am b/tests/Makefile.am
index 11794a2bc..b2658f0a3 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -358,6 +358,12 @@ MODFILES = \
 	var-expectations/3/example.mod \
 	var-expectations/4/example.mod \
 	var-expectations/5/example.mod \
+    var-expectations/6/example.mod \
+    var-expectations/6/substitution.mod \
+    var-expectations/7/example.mod \
+    var-expectations/7/substitution.mod \
+    var-expectations/8/example.mod \
+    var-expectations/8/substitution.mod \
 	trend-component-and-var-models/vm1.mod \
 	trend-component-and-var-models/vm2.mod \
 	trend-component-and-var-models/vm3.mod \
@@ -381,6 +387,8 @@ MODFILES = \
 	pac/var-2/example.mod \
 	pac/var-3/example.mod \
 	pac/var-4/example.mod \
+    pac/var-5/example.mod \
+    pac/var-5/substitution.mod \
 	pac/trend-component-1/example.mod \
 	pac/trend-component-2/example.mod \
 	pac/trend-component-3/example.mod \
@@ -395,6 +403,8 @@ MODFILES = \
 	pac/trend-component-13a/example.mod \
 	pac/trend-component-13b/example.mod \
     estimation/univariate/bayesian.mod \
+    pac/trend-component-14/example.mod \
+    pac/trend-component-14/substitution.mod \
     dynare-command-options/ramst.mod
 
 PARTICLEFILES = \
@@ -591,6 +601,18 @@ kalman/likelihood_from_dynare/fs2000ns_uncorr_ME_missing.o.trs: kalman/likelihoo
 kalman/likelihood_from_dynare/fs2000ns_corr_ME_missing.m.trs: kalman/likelihood_from_dynare/fs2000ns_uncorr_ME.m.trs
 kalman/likelihood_from_dynare/fs2000ns_corr_ME_missing.o.trs: kalman/likelihood_from_dynare/fs2000ns_uncorr_ME.o.trs
 
+var-expectations/6/substitution.m.trs: var-expectations/6/example.m.trs
+var-expectations/6/substitution.o.trs: var-expectations/6/example.o.trs
+var-expectations/7/substitution.m.trs: var-expectations/7/example.m.trs
+var-expectations/7/substitution.o.trs: var-expectations/7/example.o.trs
+var-expectations/8/substitution.m.trs: var-expectations/8/example.m.trs
+var-expectations/8/substitution.o.trs: var-expectations/8/example.o.trs
+
+pac/var-5/substitution.m.trs: pac/var-5/example.m.trs
+pac/var-5/substitution.o.trs: pac/var-5/example.o.trs
+pac/trend-component-14/substitution.m.trs: pac/trend-component-14/example.m.trs
+pac/trend-component-14/substitution.o.trs: pac/trend-component-14/example.o.trs
+
 lmmcp/sw_newton.m.trs: lmmcp/sw_lmmcp.m.trs
 lmmcp/sw_newton.o.trs: lmmcp/sw_lmmcp.o.trs
 
diff --git a/tests/pac/trend-component-14/clean b/tests/pac/trend-component-14/clean
new file mode 100755
index 000000000..4ca44cac8
--- /dev/null
+++ b/tests/pac/trend-component-14/clean
@@ -0,0 +1,9 @@
+#!/bin/sh
+
+rm -rf example
+rm  -rf +example
+rm -f example.log
+rm -f *.mat
+rm -rf substitution
+rm  -rf +substitution
+rm -f substitution.log
diff --git a/tests/pac/trend-component-14/example.mod b/tests/pac/trend-component-14/example.mod
new file mode 100644
index 000000000..09919248c
--- /dev/null
+++ b/tests/pac/trend-component-14/example.mod
@@ -0,0 +1,93 @@
+// --+ options: transform_unary_ops, json=compute, stochastic +--
+
+var x1 x2 x1bar x2bar z y;
+
+varexo ex1 ex2 ex1bar ex2bar ez ey;
+
+parameters
+       rho_1 rho_2
+       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 beta;
+
+rho_1 = .9;
+rho_2 = -.2;
+
+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;
+
+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:x1', data_type='nonstationary']
+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', 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;
+
+[name='zpac']
+diff(z) = e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) + pac_expectation(pacman) + 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;
+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
+set_dynare_seed('default');
+TrueData = simul_backward_model(initialconditions, 300);
+
+// Print expanded PAC_EXPECTATION term.
+pac.print('pacman');
+
+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', 'TrueData');
+end;
\ No newline at end of file
diff --git a/tests/pac/trend-component-14/substitution.mod b/tests/pac/trend-component-14/substitution.mod
new file mode 100644
index 000000000..db5d12cd7
--- /dev/null
+++ b/tests/pac/trend-component-14/substitution.mod
@@ -0,0 +1,98 @@
+// --+ options: transform_unary_ops, json=compute, stochastic +--
+
+var x1 x2 x1bar x2bar z y;
+
+varexo ex1 ex2 ex1bar ex2bar ez ey;
+
+parameters
+       rho_1 rho_2
+       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 beta;
+
+rho_1 = .9;
+rho_2 = -.2;
+
+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;
+
+@#include "example/model/pac-expectations/pacman-parameters.inc"
+
+model;
+
+[name='eq:y']
+y = rho_1*y(-1) + rho_2*y(-2) + ey;
+
+[name='eq:x1', data_type='nonstationary']
+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', 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;
+
+[name='zpac']
+diff(z) = e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1))  + c_z_2*diff(z(-2)) +
+@#include "example/model/pac-expectations/pacman-expression.inc"
++ 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;
+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 500 periods
+set_dynare_seed('default');
+TrueData = simul_backward_model(initialconditions, 300);
+
+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);
+end;
+
+if max(abs(example.TrueData.data(:)-TrueData.data(:)))>1e-9
+   error('Simulations do not match.')
+end
+
+if ~isequal(length(residual), length(example.residual)) || max(abs(example.residual-residual))>1e-8
+  warning('Residuals do not match!')
+end
+
+if ~isequal(length(g1(:)), length(example.g1(:))) || max(abs(example.g1(:)-g1(:)))>1e-8
+  warning('Jacobian matrices do not match!')
+end
+
+delete('example.mat');
\ No newline at end of file
diff --git a/tests/pac/var-5/clean b/tests/pac/var-5/clean
new file mode 100755
index 000000000..20de477e3
--- /dev/null
+++ b/tests/pac/var-5/clean
@@ -0,0 +1,10 @@
+#!/bin/sh
+
+rm -rf example
+rm -rf +example
+rm -f example*.mat
+rm -f example.log
+rm -rf substitution
+rm -rf +substitution
+rm -f substitution*.mat
+rm -f substitution.log
diff --git a/tests/pac/var-5/example.mod b/tests/pac/var-5/example.mod
new file mode 100644
index 000000000..80f570d07
--- /dev/null
+++ b/tests/pac/var-5/example.mod
@@ -0,0 +1,73 @@
+// --+ options: json=compute, stochastic +--
+
+var y x z;
+
+varexo ex ey ez;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 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_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+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, eqtags=['eq:x', 'eq:y']);
+
+pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman, growth=g);
+
+model;
+
+[name='eq:y']
+y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+[name='eq:x', data_type='nonstationary']
+diff(x) = 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;
+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');
+
+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', 'TrueData');
+end;
diff --git a/tests/pac/var-5/substitution.mod b/tests/pac/var-5/substitution.mod
new file mode 100644
index 000000000..41c55fbfb
--- /dev/null
+++ b/tests/pac/var-5/substitution.mod
@@ -0,0 +1,77 @@
+// --+ options: transform_unary_ops, json=compute, stochastic +--
+
+var y x z;
+
+varexo ex ey ez;
+
+parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 g; // VAR parameters
+
+parameters e_c_m c_z_1 c_z_2;               // PAC equation parameters
+
+a_y_1 =  .2;
+a_y_2 =  .3;
+b_y_1 =  .1;
+b_y_2 =  .4;
+b_x_1 = -.1;
+b_x_2 = -.2;
+
+g = .1;
+
+e_c_m =  .1;
+c_z_1 =  .7;
+c_z_2 = -.3;
+
+@#include "example/model/pac-expectations/pacman-parameters.inc"
+
+model;
+
+[name='eq:y']
+y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
+
+[name='eq:x', data_type='nonstationary']
+diff(x) = 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 "example/model/pac-expectations/pacman-expression.inc"
++ ez;
+
+end;
+
+shocks;
+    var ex = 1.0;
+    var ey = 1.0;
+    var ez = 1.0;
+end;
+
+// Set initial conditions to zero. Please use more sensible values if any...
+initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names));
+
+// Simulate the model for 20 periods
+set_dynare_seed('default');
+TrueData = simul_backward_model(initialconditions, 20);
+
+verbatim;
+  set_dynare_seed('default');
+  y = zeros(M_.endo_nbr,1);
+  y(1:M_.orig_endo_nbr) = rand(M_.orig_endo_nbr, 1);
+  x = randn(M_.exo_nbr,1);
+  y = substitution.set_auxiliary_variables(y, x, M_.params);
+  y = [y(find(M_.lead_lag_incidence(1,:))); y];
+  example = load('example.mat');
+  [residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
+end;
+
+if max(abs(example.TrueData.data(:)-TrueData.data(:)))>1e-9
+   error('Simulations do not match.')
+end
+
+if ~isequal(length(residual), length(example.residual)) || max(abs(example.residual-residual))>1e-8
+  warning('Residuals do not match!')
+end
+
+if ~isequal(length(g1(:)), length(example.g1(:))) || max(abs(example.g1(:)-g1(:)))>1e-8
+  warning('Jacobian matrices do not match!')
+end
+
+delete('example.mat');
\ No newline at end of file
diff --git a/tests/var-expectations/6/example.mod b/tests/var-expectations/6/example.mod
index e50d0f4ec..3d64bc559 100644
--- a/tests/var-expectations/6/example.mod
+++ b/tests/var-expectations/6/example.mod
@@ -1,4 +1,4 @@
-// --+ options: stochastic,json=compute +--
+// --+ options: stochastic,transform_unary_ops,json=compute +--
 
 var foo x1 x2 x1bar x2bar;
 
@@ -77,4 +77,18 @@ 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
+end
+
+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);
+  ynames = M_.endo_names;
+  save('example.mat', 'residual', 'g1');
+end;
\ No newline at end of file
diff --git a/tests/var-expectations/6/substitution.mod b/tests/var-expectations/6/substitution.mod
new file mode 100644
index 000000000..daf189169
--- /dev/null
+++ b/tests/var-expectations/6/substitution.mod
@@ -0,0 +1,67 @@
+// --+ 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 "example/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 "example/model/var-expectations/varexp-expression.inc"
+;
+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);
+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
diff --git a/tests/var-expectations/7/example.mod b/tests/var-expectations/7/example.mod
new file mode 100644
index 000000000..ea01f0da8
--- /dev/null
+++ b/tests/var-expectations/7/example.mod
@@ -0,0 +1,52 @@
+// --+ 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
+var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+
+
+model;
+[ name = 'X' ]
+diff(x) = a*diff(x(-1)) + b*diff(x(-2)) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+log(y) = d*log(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 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');
+end;
diff --git a/tests/var-expectations/7/substitution.mod b/tests/var-expectations/7/substitution.mod
new file mode 100644
index 000000000..43002de67
--- /dev/null
+++ b/tests/var-expectations/7/substitution.mod
@@ -0,0 +1,53 @@
+// --+ options: stochastic,transform_unary_ops,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;
+
+@#include "example/model/var-expectations/varexp-parameters.inc"
+
+beta = 1/(1+.02);
+
+// Define a VAR_EXPECTATION_MODEL
+var_model(model_name = toto, eqtags = [ 'X' 'Y' 'Z' ]);
+
+model;
+[ name = 'X' ]
+diff(x) = a*diff(x(-1)) + b*diff(x(-2)) + c*z(-2) + e_x;
+[ name = 'Z' ]
+z = f*z(-1) + e_z;
+[ name = 'Y' ]
+log(y) = d*log(y(-2)) + e*z(-1) + e_y;
+
+foo = .5*foo(-1) +
+@#include "example/model/var-expectations/varexp-expression.inc"
+;
+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');
+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
diff --git a/tests/var-expectations/8/example.mod b/tests/var-expectations/8/example.mod
new file mode 100644
index 000000000..44b81430f
--- /dev/null
+++ b/tests/var-expectations/8/example.mod
@@ -0,0 +1,51 @@
+// --+ 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
+var_expectation_model(model_name = varexp, variable = x, auxiliary_model_name = toto, horizon = 1, discount = beta)  ;
+
+model;
+[ name = 'X' ]
+diff(log(x)) = a*diff(log(x(-1))) + b*diff(log(x(-2))) + c*diff(z(-2)) + e_x;
+[ name = 'Z' ]
+diff(z) = f*diff(z(-1)) + e_z;
+[ name = 'Y' ]
+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');
+
+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');
+end;
\ No newline at end of file
diff --git a/tests/var-expectations/8/substitution.mod b/tests/var-expectations/8/substitution.mod
new file mode 100644
index 000000000..561d48034
--- /dev/null
+++ b/tests/var-expectations/8/substitution.mod
@@ -0,0 +1,50 @@
+// --+ options: stochastic,transform_unary_ops,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;
+
+@#include "example/model/var-expectations/varexp-parameters.inc"
+
+beta = 1/(1+.02);
+
+model;
+[ name = 'X' ]
+diff(log(x)) = a*diff(log(x(-1))) + b*diff(log(x(-2))) + c*diff(z(-2)) + e_x;
+[ name = 'Z' ]
+diff(z) = f*diff(z(-1)) + e_z;
+[ name = 'Y' ]
+log(y) = d*log(y(-2)) + e*diff(z(-1)) + e_y;
+
+foo = .5*foo(-1) +
+@#include "example/model/var-expectations/varexp-expression.inc"
+;
+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');
+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