Verified Commit ac2b3f2b authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Add integration tests for model inversion.

parent 0df38369
......@@ -463,6 +463,31 @@ MODFILES = \
bgp/solow-1/solow.mod \
bgp/nk-1/nk.mod \
bgp/ramsey-1/ramsey.mod \
model-inversion/bk-1/solow_ces.mod \
model-inversion/bk-2/solow_ces.mod \
model-inversion/bk-3a/varmodel.mod \
model-inversion/bk-3a-nlsolve/varmodel.mod \
model-inversion/bk-3b/varmodel.mod \
model-inversion/bk-3b-nlsolve/varmodel.mod \
model-inversion/bk-4/varmodel.mod \
model-inversion/bk-4-nlsolve/varmodel.mod \
model-inversion/bk-5/varmodel.mod \
model-inversion/bk-5-nlsolve/varmodel.mod \
model-inversion/bk-6/varmodel.mod \
model-inversion/bk-7/varmodel.mod \
model-inversion/bk-8/varmodel.mod \
model-inversion/bk-9/msm1.mod \
model-inversion/bk-9/msm2.mod \
model-inversion/bk-diff-1/varmodel.mod \
model-inversion/bk-diff-2/varmodel.mod \
model-inversion/bk-diff-3/varmodel.mod \
model-inversion/bk-diff-4/varmodel.mod \
model-inversion/bk-diff-5/varmodel.mod \
model-inversion/nk-1/simulate.mod \
model-inversion/nk-1/invert.mod \
model-inversion/nk-2/simulate.mod \
model-inversion/nk-2/invert.mod \
model-inversion/nk-2/z_check_inversion.mod \
dynare-command-options/ramst.mod \
log_transform/example1.mod \
log_transform/ramst.mod \
......@@ -947,6 +972,14 @@ ramst_model_edit.o.trs: ramst.o.trs
log_transform/ramst.m.trs: ramst.m.trs
log_transform/ramst.o.trs: ramst.o.trs
model-inversion/nk-1/invert.m.trs: model-inversion/nk-1/simulate.m.trs
model-inversion/nk-1/invert.o.trs: model-inversion/nk-1/simulate.o.trs
model-inversion/nk-2/invert.m.trs: model-inversion/nk-2/simulate.m.trs
model-inversion/nk-2/invert.o.trs: model-inversion/nk-2/simulate.o.trs
model-inversion/nk-2/z_check_inversion.m.trs: model-inversion/nk-2/invert.m.trs
model-inversion/nk-2/z_check_inversion.o.trs: model-inversion/nk-2/invert.o.trs
observation_trends_and_prefiltering/MCMC: m/observation_trends_and_prefiltering/MCMC o/observation_trends_and_prefiltering/MCMC
m/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.m.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
......@@ -1160,6 +1193,10 @@ estimation/univariate: m/estimation/univariate o/estimation/univariate
m/estimation/univariate: $(patsubst %.mod, %.m.trs, $(filter estimation/univariate/%.mod, $(MODFILES)))
o/estimation/univariate: $(patsubst %.mod, %.o.trs, $(filter estimation/univariate/%.mod, $(MODFILES)))
particle: m/model-inversion o/model-inversion
m/model-inversion: $(patsubst %.mod, %.m.trs, $(filter model-inversion/%.mod, $(MODFILES)))
o/model-inversion: $(patsubst %.mod, %.o.trs, $(filter model-inversion/%.mod, $(MODFILES)))
# ECB files
M_ECB_TRS_FILES = $(patsubst %.mod, %.m.trs, $(ECB_MODFILES))
......
// --+ options: stochastic +--
/* © 2022 Dynare Team
*
* This file is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* It is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with the file. If not, see <http://www.gnu.org/licenses/>.
*/
var Efficiency // $A$
EfficiencyGrowth // $X$
Population // $L$
PopulationGrowth // $N$
Output // $Y$
PhysicalCapitalStock ; // $K$
varexo e_x // $\varepsilon_x$
e_n ; // $\varepsilon_n$
parameters alpha // $\alpha$
epsilon // $\varepsilon$
delta // $\delta$
s // $s$
rho_x // $\rho_x$
rho_n // $\rho_n$
EfficiencyGrowth_ss // $X^{\star}$
PopulationGrowth_ss ; // $N^{\star}$
alpha = .33;
epsilon = .70;
delta = .02;
s = .20;
rho_x = .90;
rho_n = .95;
EfficiencyGrowth_ss = 1.02;
PopulationGrowth_ss = 1.02;
if s>delta*alpha^(-epsilon/(epsilon-1))
disp('The model admits a unique positive steady state.')
end
model;
Efficiency = EfficiencyGrowth*Efficiency(-1);
EfficiencyGrowth/EfficiencyGrowth_ss = (EfficiencyGrowth(-1)/EfficiencyGrowth_ss)^(rho_x)*exp(e_x);
Population = PopulationGrowth*Population(-1);
PopulationGrowth/PopulationGrowth_ss = (PopulationGrowth(-1)/PopulationGrowth_ss)^(rho_n)*exp(e_n);
Output = (alpha*PhysicalCapitalStock(-1)^((epsilon-1)/epsilon)+(1-alpha)*(Efficiency*Population)^((epsilon-1)/epsilon))^(epsilon/(epsilon-1));
PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
end;
histval;
Efficiency(0) = 1;
EfficiencyGrowth(0) = 1.02;
Population(0) = 1;
PopulationGrowth(0) = 1.02;
PhysicalCapitalStock(0) = 1;
end;
shocks;
var e_x = 0.005;
var e_n = 0.001;
end;
TrueData = simul_backward_nonlinear_model([], 100, options_, M_, oo_);
// Set the periods where some of the endogenous variables will be constrained.
subsample = 2Y:100Y;
// Copy the generated data
SimulatedData = copy(TrueData);
// Set the constrained paths for the endogenous variables (Output and PhysicalCapitalStock).
constrainedpaths = SimulatedData{'Output','PhysicalCapitalStock'}(subsample);
// Set the instruments (innovations used to control the paths for the endogenous variables).
exogenousvariables = dseries(NaN(100, 2), 1Y, {'e_x';'e_n'});
// Invert the model by calling the model_inversion routine.
[endogenousvariables, exogenousvariables] = model_inversion(constrainedpaths, exogenousvariables, SimulatedData, M_, options_, oo_);
// Check that all the constraints are satisfied.
if max(abs(constrainedpaths.Output(subsample).data-endogenousvariables.Output(subsample).data))>1e-12
error('Constraint on Output path is not satisfied!')
end
if max(abs(constrainedpaths.PhysicalCapitalStock(subsample).data-endogenousvariables.PhysicalCapitalStock(subsample).data))>1e-12
error('Constraint on PhysicalCapitalStock path is not satisfied!')
end
// Check the consistency of the deduced innovations with the true innovations.
/* REMARK
**
** In this model the two innovatons, on population growth and efficiency growth, cannot be identified simulatneously because
** what matters for the dynamic of the physical capital stock and the output is the same non linear function of the two innovations.
** This explains why we only check that the product of Efficiency and Population are the same in the simulated (true) data and in the
** data returned by the inversion routine.
*/
TrueEfficiency = TrueData.Efficiency;
TruePopulation = TrueData.Population;
TrueEfficiencyTimesPopulation = TrueEfficiency*TruePopulation;
EfficiencyTimesPopulation = endogenousvariables.Efficiency*endogenousvariables.Population;
if max(abs(TrueEfficiencyTimesPopulation(2Y:100Y).data-EfficiencyTimesPopulation(2Y:100Y).data))>1e-5
error('Model inversion is not consitent with true innovations.')
end
// --+ options: stochastic +--
/* © 2022 Dynare Team
*
* This file is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* It is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with the file. If not, see <http://www.gnu.org/licenses/>.
*/
var Efficiency // $A$
EfficiencyGrowth // $X$
Population // $L$
PopulationGrowth // $N$
Output // $Y$
PhysicalCapitalStock ; // $K$
varexo e_x // $\varepsilon_x$
e_n ; // $\varepsilon_n$
parameters alpha // $\alpha$
epsilon // $\varepsilon$
delta // $\delta$
s // $s$
rho_x // $\rho_x$
rho_n // $\rho_n$
EfficiencyGrowth_ss // $X^{\star}$
PopulationGrowth_ss ; // $N^{\star}$
alpha = .33;
epsilon = .70;
delta = .02;
s = .20;
rho_x = .90;
rho_n = .95;
EfficiencyGrowth_ss = 1.02;
PopulationGrowth_ss = 1.02;
if s>delta*alpha^(-epsilon/(epsilon-1))
disp('The model admits a unique positive steady state.')
end
model;
Efficiency = EfficiencyGrowth*Efficiency(-1);
EfficiencyGrowth/EfficiencyGrowth_ss = (EfficiencyGrowth(-1)/EfficiencyGrowth_ss)^(rho_x)*exp(e_x);
Population = PopulationGrowth*Population(-1);
PopulationGrowth/PopulationGrowth_ss = (PopulationGrowth(-1)/PopulationGrowth_ss)^(rho_n)*exp(e_n);
Output = (alpha*PhysicalCapitalStock(-1)^((epsilon-1)/epsilon)+(1-alpha)*(Efficiency*Population)^((epsilon-1)/epsilon))^(epsilon/(epsilon-1));
PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
end;
histval;
Efficiency(0) = 1;
EfficiencyGrowth(0) = 1.02;
Population(0) = 1;
PopulationGrowth(0) = 1.02;
PhysicalCapitalStock(0) = 1;
end;
shocks;
var e_x = 0.005;
var e_n = 0.001;
end;
TrueData = simul_backward_nonlinear_model([], 200, options_, M_, oo_);
// Set the periods where some of the endogenous variables will be constrained.
subsample = 2Y:100Y;
// Copy the generated data
SimulatedData = copy(TrueData);
// Set the constrained paths for the endogenous variables (Output and PhysicalCapitalStock).
constrainedpaths = SimulatedData{'Output'}(subsample);
// Set the instruments (innovations used to control the paths for the endogenous variables).
exogenousvariables = dseries([NaN(100, 1), TrueData.e_n.data(1:100)], 1Y, M_.exo_names);
// Invert the model by calling the model_inversion routine.
[endogenousvariables, exogenousvariables] = model_inversion(constrainedpaths, exogenousvariables, SimulatedData, M_, options_, oo_);
// Check that all the constraints are satisfied.
if max(abs(constrainedpaths.Output(subsample).data-endogenousvariables.Output(subsample).data))>1e-12
error('Constraint on Output path is not satisfied!')
end
// Check that the solution for the PhysicalCapitalStock is consistent with the simulated data
if max(abs(SimulatedData.PhysicalCapitalStock(2Y:100Y).data-endogenousvariables.PhysicalCapitalStock(2Y:100Y).data))>1e-9
error('Results are not consistent with true data (Physical capital stock).')
end
// Check the consistency of the deduced innovations with the true innovations.
/* REMARK
**
** In this model the two innovatons, on population growth and efficiency growth, cannot be identified simulatneously because
** what matters for the dynamic of the physical capital stock and the output is the same non linear function of the two innovations.
** This explains why we only check that the product of Efficiency and Population are the same in the simulated (true) data and in the
** data returned by the inversion routine.
*/
if max(abs(SimulatedData.e_x(2Y:100Y).data-exogenousvariables.e_x(2Y:100Y).data))>1e-5
error('Model inversion is not consitent with true innovations (e_x).')
end
if max(abs(SimulatedData.e_n(2Y:100Y).data-exogenousvariables.e_n(2Y:100Y).data))>1e-16
error('Model inversion changed a calibrated innovation (e_n).')
end
// --+ options: stochastic +--
/* © 2022 Dynare Team
*
* This file is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* It is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with the file. If not, see <http://www.gnu.org/licenses/>.
*/
var y1 y2 y3 ;
varexo e1 e2 e3 ;
parameters a11 a12 a13 a21 a22 a23 a31 a32 a33 ;
/*
** Simulate the elements of the first order autoregressive matrix (we impose stability of the model, note that
** inversion fails if the model is explosive, ie the autoregressive matrix has at least one root greater than
** one in modulus) probably because of the propagation of roundoff errors.
*/
D = diag([.9 .7 .8]);
P = randn(3,3);
A = P*D*inv(P);
a11 = A(1,1);
a12 = A(1,2);
a13 = A(1,3);
a21 = A(2,1);
a22 = A(2,2);
a23 = A(2,3);
a31 = A(3,1);
a32 = A(3,2);
a33 = A(3,3);
model;
y1 = a11*y1(-1) + a12*y2(-1) + a13*y3(-1) + e1 ;
y2 = a21*y1(-1) + a22*y2(-1) + a23*y3(-1) + e2 ;
y3 = a31*y1(-1) + a32*y2(-1) + a33*y3(-1) + e3 ;
end;
histval;
y1(0) = 0;
y2(0) = 0;
y3(0) = 0;
end;
shocks;
var e1 = 1.0;
var e2 = 1.0;
var e3 = 1.0;
end;
steady;
check;
TrueData = simul_backward_model([], 200);
// Set the periods where some of the endogenous variables will be constrained.
subsample = 3Y:100Y;
// Load the generated data
SimulatedData = copy(TrueData);
// Set the constrained paths for the endogenous variables (Output and PhysicalCapitalStock).
constrainedpaths = SimulatedData{'y1'}(subsample);
// Set the instruments (innovations used to control the paths for the endogenous variables).
exogenousvariables = dseries([NaN(100, 1) TrueData{'e2','e3'}.data(1:100,:)], '1Y', M_.exo_names);
// Invert the model by calling the model_inversion routine.
[endogenousvariables, exogenousvariables] = model_inversion(constrainedpaths, exogenousvariables, SimulatedData, M_, options_, oo_);
// Check that all the constraints are satisfied.
if max(abs(constrainedpaths(subsample).y1.data-endogenousvariables(subsample).y1.data))>1e-12
error('Constraint on y1 path is not satisfied!')
end
if max(abs(exogenousvariables(subsample).e2.data-SimulatedData(subsample).e2.data))>1e-12
error('Constraint on e1 path is not satisfied!')
end
if max(abs(exogenousvariables(subsample).e3.data-SimulatedData(subsample).e3.data))>1e-12
error('Constraint on e2 path is not satisfied!')
end
// Check consistency of the results.
if max(abs(SimulatedData(subsample).y2.data-endogenousvariables(subsample).y2.data))>1e-12
error('Model inversion is not consistent with respect to y2')
end
if max(abs(SimulatedData(subsample).y3.data-endogenousvariables(subsample).y3.data))>1e-12
error('Model inversion is not consistent with respect to y3')
end
if max(abs(exogenousvariables(subsample).e1.data-SimulatedData(subsample).e1.data))>1e-12
error('Model inversion is not consistent with true innovations (e3)')
end
// --+ options: stochastic +--
/* © 2022 Dynare Team
*
* This file is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* It is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with the file. If not, see <http://www.gnu.org/licenses/>.
*/
var y1 y2 y3 ;
varexo e1 e2 e3 ;
parameters a11 a12 a13 a21 a22 a23 a31 a32 a33 ;
/*
** Simulate the elements of the first order autoregressive matrix (we impose stability of the model, note that
** inversion fails if the model is explosive, ie the autoregressive matrix has at least one root greater than
** one in modulus) probably because of the propagation of roundoff errors.
*/
D = diag([.9 .7 .8]);
P = randn(3,3);
A = P*D*inv(P);
a11 = A(1,1);
a12 = A(1,2);
a13 = A(1,3);
a21 = A(2,1);
a22 = A(2,2);
a23 = A(2,3);
a31 = A(3,1);
a32 = A(3,2);
a33 = A(3,3);
model(linear);
y1 = a11*y1(-1) + a12*y2(-1) + a13*y3(-1) + e1 ;
y2 = a21*y1(-1) + a22*y2(-1) + a23*y3(-1) + e2 ;
y3 = a31*y1(-1) + a32*y2(-1) + a33*y3(-1) + e3 ;
end;
histval;
y1(0) = 0;
y2(0) = 0;
y3(0) = 0;
end;
shocks;
var e1 = 1.0;
var e2 = 1.0;
var e3 = 1.0;
end;
steady;
check;
TrueData = simul_backward_model([], 200);//, options_, M_, oo_);
// Set the periods where some of the endogenous variables will be constrained.
subsample = 3Y:100Y;
// Load the generated data
SimulatedData = copy(TrueData);
// Set the constrained paths for the endogenous variables (Output and PhysicalCapitalStock).
constrainedpaths = SimulatedData{'y1'}(subsample);
// Set the instruments (innovations used to control the paths for the endogenous variables).
exogenousvariables = dseries([NaN(100, 1) TrueData{'e2','e3'}.data(1:100,:)], '1Y', M_.exo_names);
// Invert the model by calling the model_inversion routine.
[endogenousvariables, exogenousvariables] = model_inversion(constrainedpaths, exogenousvariables, SimulatedData, M_, options_, oo_);
// Check that all the constraints are satisfied.
if max(abs(constrainedpaths(subsample).y1.data-endogenousvariables(subsample).y1.data))>1e-12
error('Constraint on y1 path is not satisfied!')
end
if max(abs(exogenousvariables(subsample).e2.data-SimulatedData(subsample).e2.data))>1e-12
error('Constraint on e1 path is not satisfied!')
end
if max(abs(exogenousvariables(subsample).e3.data-SimulatedData(subsample).e3.data))>1e-12
error('Constraint on e2 path is not satisfied!')
end
// Check consistency of the results.
if max(abs(SimulatedData(subsample).y2.data-endogenousvariables(subsample).y2.data))>1e-12
error('Model inversion is not consistent with respect to y2')
end
if max(abs(SimulatedData(subsample).y3.data-endogenousvariables(subsample).y3.data))>1e-12
error('Model inversion is not consistent with respect to y3')
end
if max(abs(exogenousvariables(subsample).e1.data-SimulatedData(subsample).e1.data))>1e-12
error('Model inversion is not consistent with true innovations (e3)')
end
// --+ options: stochastic +--
/* © 2022 Dynare Team
*
* This file is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* It is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with the file. If not, see <http://www.gnu.org/licenses/>.
*/
var y1 y2 y3 ;
varexo e1 e2 e3 ;
parameters a11 a12 a13 a21 a22 a23 a31 a32 a33 b11 b12 b13 b22 b23 b33 ;
/*
** Simulate the elements of the first order autoregressive matrix (we impose stability of the model, note that
** inversion fails if the model is explosive, ie the autoregressive matrix has at least one root greater than
** one in modulus) probably because of the propagation of roundoff errors.
*/
D = diag([.9 .7 .8]);
P = randn(3,3);
A = P*D*inv(P);
a11 = A(1,1);
a12 = A(1,2);
a13 = A(1,3);
a21 = A(2,1);
a22 = A(2,2);
a23 = A(2,3);
a31 = A(3,1);
a32 = A(3,2);
a33 = A(3,3);
b11 = .10;
b12 = -.30;
b13 = .05;
b22 = .20;
b23 = -.05;
b33 = .10;
model;
y1 = a11*y1(-1) + a12*y2(-1) + a13*y3(-1) + b11*e1 + b12*e2 + b13*e3 ;
y2 = a21*y1(-1) + a22*y2(-1) + a23*y3(-1) + b22*e2 + b23*e3 ;
y3 = a31*y1(-1) + a32*y2(-1) + a33*y3(-1) + b33*e3 ;
end;
histval;
y1(0) = 0;
y2(0) = 0;
y3(0) = 0;
end;
shocks;
var e1 = 1.0;
var e2 = 1.0;
var e3 = 1.0;
end;
steady;
check;
TrueData = simul_backward_model([], 200);
// Set the periods where some of the endogenous variables will be constrained.
subsample = 3Y:100Y;
// Load the generated data
SimulatedData = copy(TrueData);
// Set the constrained paths for the endogenous variables (Output and PhysicalCapitalStock).
constrainedpaths = SimulatedData{'y1'}(subsample);
// Set the instruments (innovations used to control the paths for the endogenous variables).
exogenousvariables = dseries([NaN(100, 1) TrueData{'e2','e3'}.data(1:100,:)], '1Y', M_.exo_names);
// Invert the model by calling the model_inversion routine.
[endogenousvariables, exogenousvariables] = model_inversion(constrainedpaths, exogenousvariables, SimulatedData, M_, options_, oo_);
// Check that all the constraints are satisfied.
if max(abs(constrainedpaths(subsample).y1.data-endogenousvariables(subsample).y1.data))>1e-12
error('Constraint on y1 path is not satisfied!')
end
if max(abs(exogenousvariables(subsample).e2.data-SimulatedData(subsample).e2.data))>1e-12
error('Constraint on e1 path is not satisfied!')
end
if max(abs(exogenousvariables(subsample).e3.data-SimulatedData(subsample).e3.data))>1e-12
error('Constraint on e2 path is not satisfied!')
end
// Check consistency of the results.
if max(abs(SimulatedData(subsample).y2.data-endogenousvariables(subsample).y2.data))>1e-12
error('Model inversion is not consistent with respect to y2')
end
if max(abs(SimulatedData(subsample).y3.data-endogenousvariables(subsample).y3.data))>1e-12
error('Model inversion is not consistent with respect to y3')
end
if max(abs(exogenousvariables(subsample).e1.data-SimulatedData(subsample).e1.data))>1e-12
error('Model inversion is not consistent with true innovations (e3)')
end
// --+ options: stochastic +--
/* © 2022 Dynare Team