Allow exogenous variables in pac.estimation.iterative_ols routine.

The parameters associated to these additional variables can be estimated or calibrated.
time-shift
Stéphane Adjemian (Charybdis) 2019-03-01 23:32:47 +01:00
parent 21e6260011
commit d00b57541e
Signed by untrusted user who does not match committer: stepan
GPG Key ID: A6D44CB9C64CE77B
13 changed files with 818 additions and 32 deletions

View File

@ -1,4 +1,4 @@
function [pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, pnames_, ipnames_, params, data, islaggedvariables, eqtag] = ...
function [pacmodl, lhs, rhs, pnames, enames, xnames, rname, pid, eid, xid, pnames_, ipnames_, params, data, islaggedvariables, eqtag] = ...
init(M_, oo_, eqname, params, data, range)
% Copyright (C) 2018-2019 Dynare Team
@ -70,7 +70,10 @@ end
stack = dbstack;
ipnames__ = ipnames_; % The user provided order.
if isequal(stack(2).name, 'iterative_ols')
ipnames_ = [M_.pac.(pacmodl).equations.(eqtag).ec.params; M_.pac.(pacmodl).equations.(eqtag).ar.params']; % The correct order.
ipnames_ = [M_.pac.(pacmodl).equations.(eqtag).ec.params; M_.pac.(pacmodl).equations.(eqtag).ar.params'];
if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive')
ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).additive.params(~isnan(M_.pac.(pacmodl).equations.(eqtag).additive.params))'];
end
if isfield(M_.pac.(pacmodl).equations.(eqtag), 'share_of_optimizing_agents_index')
ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index];
end

View File

@ -43,7 +43,7 @@ global M_ oo_
debug = false; % If true, prints the value of SSR for the initial guess.
[pacmodl, ~, rhs, ~, ~, ~, ~, ~, ~, ~, ipnames_, params, data, ~, eqtag] = ...
[pacmodl, ~, rhs, ~, ~, ~, rname, ~, ~, ~, ~, ipnames_, params, data, ~, eqtag] = ...
pac.estimate.init(M_, oo_, eqname, params, data, range);
% Set initial condition.
@ -52,12 +52,19 @@ params0 = cell2mat(struct2cell(params));
% Set flag for models with non optimizing agents.
is_non_optimizing_agents = isfield(M_.pac.(pacmodl).equations.(eqtag), 'non_optimizing_behaviour');
% Set flag for models with exogenous variables (outside of non optimizing agents part)
if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive')
is_exogenous_variables = length(M_.pac.(pacmodl).equations.(eqtag).additive.vars)>1;
else
is_exogenous_variables = false;
end
if is_non_optimizing_agents
non_optimizing_behaviour = M_.pac.(pacmodl).equations.(eqtag).non_optimizing_behaviour;
non_optimizing_behaviour_params = NaN(length(non_optimizing_behaviour.params), 1);
noparams = isnan(non_optimizing_behaviour.params);
if ~all(noparams)
% Find estimated non optimizing behaviour parameters (if any.
% Find estimated non optimizing behaviour parameters (if any).
non_optimizing_behaviour_estimated_params = ismember(M_.param_names(non_optimizing_behaviour.params), fieldnames(params));
if any(non_optimizing_behaviour_estimated_params)
error('The estimation of non optimizing behaviour parameters is not yet allowed.')
@ -77,11 +84,32 @@ if is_non_optimizing_agents
error('The share of optimizing agents shoud be in (0,1).')
end
end
share_of_optimizing_agents_index = M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index;
else
share_of_optimizing_agents = 1.0;
share_of_optimizing_agents_index = [];
estimate_non_optimizing_agents_share = false;
end
if is_exogenous_variables
additive = M_.pac.(pacmodl).equations.(eqtag).additive;
residual_id = find(strcmp(rname, M_.exo_names));
residual_jd = find(additive.vars==residual_id & ~additive.isendo);
additive.params(residual_jd) = [];
additive.vars(residual_jd) = [];
additive.isendo(residual_jd) = [];
additive.lags(residual_jd) = [];
additive.scaling_factor(residual_jd) = [];
additive.estimation = ismember(additive.params, ipnames_);
else
additive.params = [];
additive.vars = [];
additive.isendo = [];
additive.lags = [];
additive.scaling_factor = [];
additive.estimation = [];
end
% Build PAC expectation matrix expression.
dataForPACExpectation0 = dseries();
listofvariables0 = {};
@ -144,6 +172,42 @@ else
dataForNonOptimizingBehaviour = dseries();
end
% Build data for exogenous variables (out of non optimizing behaviour term).
if is_exogenous_variables
listofvariables2 = {}; j = 0;
dataForExogenousVariables = dseries(); % Estimated parameters
dataForExogenousVariables_ = 0; % Calibrated parameters
is_any_calibrated_parameter_x = false;
is_any_estimated_parameter_x = false;
for i=1:length(additive.vars)
if additive.isendo(i)
variable = M_.endo_names{additive.vars(i)};
else
variable = M_.exo_names{additive.vars(i)};
end
if additive.estimation(i)
j = j+1;
is_any_estimated_parameter_x = true;
listofvariables2{j} = variable;
dataForExogenousVariables = [dataForExogenousVariables, additive.scaling_factor(i)*data{variable}.lag(additive.lags(i))];
else
is_any_calibrated_parameter_x = true;
tmp = data{variable}.lag(additive.lags(i)).data;
if ~isnan(additive.params(i))
tmp = M_.params(additive.params(i))*tmp;
end
tmp = additive.scaling_factor(i)*tmp;
dataForExogenousVariables_ = dataForExogenousVariables_+tmp;
end
end
if is_any_calibrated_parameter_x
dataForExogenousVariables_ = dseries(dataForExogenousVariables_, data.dates(1), 'exogenous_variables_associated_with_calibrated_parameters');
end
else
dataForExogenousVariables = dseries();
dataForExogenousVariables_ = dseries();
end
% Reorder ec.vars locally if necessary. Second variable must be the
% endogenous variable, while the first must be the associated trend.
if M_.pac.(pacmodl).equations.(eqtag).ec.isendo(2)
@ -152,7 +216,7 @@ else
ecvars = flip(M_.pac.(pacmodl).equations.(eqtag).ec.vars);
end
% Build matrix for EC and AR terms.
%% Build matrix for EC and AR terms.
DataForOLS = dseries();
% Error correction term is trend minus the level of the endogenous variable.
@ -169,16 +233,39 @@ end
XDATA = DataForOLS{listofvariables3{:}}(range).data;
if estimate_non_optimizing_agents_share
[~, ecm_params_id] = setdiff(ipnames_, M_.pac.(pacmodl).equations.(eqtag).share_of_optimizing_agents_index);
[~, share_param_id] = setdiff(1:length(ipnames_), ecm_params_id);
params0_ = params0(ecm_params_id);
share_of_optimizing_agents = params0(share_param_id);
if share_of_optimizing_agents>1 || share_of_optimizing_agents<0
error('Initial value for the share of optimizing agents shoud be in (0,1).')
end
% Get index in params0 for share of optimizing agents parameter (if
% not estimated, params_id_0 is empty).
if is_non_optimizing_agents
params_id_0 = find(ipnames_==share_of_optimizing_agents_index);
else
params0_ = params0;
params_id_0 = [];
end
% Get indices in params0 for EC and AR parameters
[~, params_id_1] = setdiff(ipnames_, [share_of_optimizing_agents_index, additive.params]);
% Get indices in params0 for other parameters (optimizing agents
% share plus parameters related to exogenous variables).
[~, params_id_2] = setdiff(1:length(ipnames_), params_id_1);
% Get indices in params0 for the parameters associated to the
% exogenous variables.
params_id_3 = setdiff(params_id_2, params_id_0);
% Get values for EC and AR parameters
params0_ = params0(params_id_1);
% Get values for parameters associated to the exogenous variables.
params0__ = params0(params_id_3);
% Get value of the share of optimizing agents.
if estimate_non_optimizing_agents_share
share_of_optimizing_agents = params0(params_id_0);
end
% Check that the share is in (0,1)
if share_of_optimizing_agents>1 || share_of_optimizing_agents<0
error('Initial value for the share of optimizing agents shoud be in (0,1).')
end
% Update the vector of parameters.
@ -187,19 +274,6 @@ M_.params(ipnames_) = params0;
% Update the reduced form PAC expectation parameters and compute the expectations.
[PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
if debug
YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
if is_non_optimizing_agents
YDATA = YDATA-share_of_optimizing_agents*PacExpectations;
YDATA = YDATA-(1-share_of_optimizing_agents)*(dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params);
else
YDATA = YDATA-PacExpectations;
end
r = YDATA-XDATA*(params0_*share_of_optimizing_agents);
ssr = r'*r;
fprintf('\nInitial value of the objective (SSR) is %s.\n\n', num2str(ssr));
end
noconvergence = true;
counter = 0;
@ -213,6 +287,14 @@ while noconvergence
else
YDATA = YDATA-PacExpectations;
end
if is_exogenous_variables
if is_any_calibrated_parameter_x
YDATA = YDATA-dataForExogenousVariables_(range).data;
end
if is_any_estimated_parameter_x
YDATA = YDATA-dataForExogenousVariables{listofvariables2{:}}(range).data*params0__;
end
end
% Run OLS to estimate PAC parameters (autoregressive parameters and error correction parameter).
params1_ = (XDATA\YDATA)/share_of_optimizing_agents;
% Compute residuals and sum of squareed residuals.
@ -226,18 +308,53 @@ while noconvergence
% Update the value of the share of non optimizing agents (if estimated)
if estimate_non_optimizing_agents_share
% First update the parameters and compute the PAC expectation reduced form parameters.
M_.params(ipnames_(ecm_params_id)) = params0_;
M_.params(ipnames_(params_id_1)) = params0_;
[PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
% Set vector for left handside variable.
YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
YDATA = YDATA-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params;
if is_exogenous_variables && is_any_calibrated_parameter_x
YDATA = YDATA-dataForExogenousVariables_(range).data;
end
% Set vector for regressor.
ZDATA = XDATA*params0_+PacExpectations-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params;
if is_exogenous_variables && is_any_estimated_parameter_x
ZDATA = [ZDATA, dataForExogenousVariables{listofvariables2{:}}(range).data];
end
% Update the (estimated) share of optimizing agents by running OLS
share_of_optimizing_agents = (ZDATA\YDATA);
beta = (ZDATA\YDATA);
share_of_optimizing_agents = beta(1);
if is_exogenous_variables && is_any_estimated_parameter_x
params0__ = beta(2:end);
end
% Force the share of optimizing agents to be in [0,1].
share_of_optimizing_agents = max(min(share_of_optimizing_agents, 1.0), 0.0);
M_.params(ipnames_(share_param_id)) = share_of_optimizing_agents;
M_.params(ipnames_(params_id_0)) = share_of_optimizing_agents;
if is_exogenous_variables && is_any_estimated_parameter_x
M_.params(ipnames_(params_id_3)) = params0__;
end
elseif is_exogenous_variables && is_any_estimated_parameter_x
% First update the parameters and compute the PAC expectation reduced form parameters.
M_.params(ipnames_(params_id_1)) = params0_;
[PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);
% Set vector for left handside variable.
YDATA = data{M_.endo_names{M_.pac.(pacmodl).equations.(eqtag).lhs_var}}(range).data;
if is_any_calibrated_parameter_x
YDATA = YDATA-dataForExogenousVariables_(range).data;
end
if is_non_optimizing_agents
YDATA = YDATA-share_of_optimizing_agents*(XDATA*params0_+PacExpectations) - ...
(1-share_of_optimizing_agents)*(dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params);
else
YDATA = YDATA-XDATA*params0_-PacExpectations;
end
% Set vector for regressor
ZDATA = dataForExogenousVariables{listofvariables2{:}}(range).data;
% Update the (estimated) parameters associated to the
% exogenous variables.
beta = (ZDATA\YDATA);
params0__ = beta;
M_.params(ipnames_(params_id_3)) = params0__;
else
M_.params(ipnames_) = params0_;
[PacExpectations, M_] = UpdatePacExpectationsData(dataPAC0, dataPAC1, data, range, pacmodl, eqtag, M_, oo_);

View File

@ -76,7 +76,7 @@ if nargin>4 && (isequal(optimizer, 'GaussNewton') || isequal(optimizer, 'lsqnonl
objective = 'r_';
end
[pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables, eqtag] = ...
[pacmodl, lhs, rhs, pnames, enames, xnames, ~, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables, eqtag] = ...
pac.estimate.init(M_, oo_, eqname, params, data, range);
% Check that the error correction term is correct.

@ -1 +1 @@
Subproject commit 2312ce13dc58d3645e0f37f1a757fbe50b6932fd
Subproject commit 5f013756f23fe89ec946fcff2cfdad7627693e12

View File

@ -426,6 +426,10 @@ MODFILES = \
pac/trend-component-20-2/example.mod \
pac/trend-component-20-3/example.mod \
pac/trend-component-20-4/example.mod \
pac/trend-component-21/example.mod \
pac/trend-component-22/example.mod \
pac/trend-component-23/example.mod \
pac/trend-component-24/example.mod \
ecb/backward-models/irf/solow_1.mod \
ecb/backward-models/irf/solow_2.mod \
dynare-command-options/ramst.mod

View File

@ -0,0 +1,8 @@
#!/bin/sh
rm -rf example
rm -rf +example
rm -f example.log
rm -f *.mat
rm -f *.m
rm -f *.dat

View File

@ -0,0 +1,160 @@
// --+ options: json=compute, stochastic +--
var x1 x2 x1bar x2bar z y x;
varexo ex1 ex2 ex1bar ex2bar ez ey ex;
parameters
rho_1 rho_2 rho_3 rho_4
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 c_z_dx2 beta
lambda;
rho_1 = .9;
rho_2 = -.2;
rho_3 = .4;
rho_4 = -.3;
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;
c_z_dx2 = .3;
lambda = 0.5; // Share of optimizing agents.
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:x']
x = rho_3*x(-1) + rho_4*x(-2) + ex;
[name='eq:x1']
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']
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']
x1bar = x1bar(-1) + ex1bar;
[name='eq:x2bar']
x2bar = x2bar(-1) + ex2bar;
[name='zpac']
diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1)) + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + 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;
var ex = 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
TrueData = simul_backward_model(initialconditions, 5000);
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.c_z_dx2 = .5;
eparams.lambda = .7;
// Define the dataset used for estimation
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols))
disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols))
skipline(2)
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.c_z_dx2 = .5;
eparams.lambda = .7;
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
c_z_dx2_nls = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
lambda_nls = M_.params(strmatch('lambda', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls))
disp(sprintf('Estimate of lambda: %f', lambda_nls))
if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
end
if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
end
if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
end
if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)')
end
if abs(lambda_nls-lambda_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (lambda)')
end

View File

@ -0,0 +1,8 @@
#!/bin/sh
rm -rf example
rm -rf +example
rm -f example.log
rm -f *.mat
rm -f *.m
rm -f *.dat

View File

@ -0,0 +1,150 @@
// --+ options: json=compute, stochastic +--
var x1 x2 x1bar x2bar z y x;
varexo ex1 ex2 ex1bar ex2bar ez ey ex;
parameters
rho_1 rho_2 rho_3 rho_4
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 c_z_dx2 beta
lambda;
rho_1 = .9;
rho_2 = -.2;
rho_3 = .4;
rho_4 = -.3;
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;
c_z_dx2 = .3;
lambda = 0.5; // Share of optimizing agents.
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:x']
x = rho_3*x(-1) + rho_4*x(-2) + ex;
[name='eq:x1']
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']
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']
x1bar = x1bar(-1) + ex1bar;
[name='eq:x2bar']
x2bar = x2bar(-1) + ex2bar;
[name='zpac']
diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1)) + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + 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;
var ex = 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
TrueData = simul_backward_model(initialconditions, 5000);
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.c_z_dx2 = .5;
// Define the dataset used for estimation
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols))
skipline(2)
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.c_z_dx2 = .5;
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
c_z_dx2_nls = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls))
if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
end
if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
end
if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
end
if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)')
end

View File

@ -0,0 +1,8 @@
#!/bin/sh
rm -rf example
rm -rf +example
rm -f example.log
rm -f *.mat
rm -f *.m
rm -f *.dat

View File

@ -0,0 +1,165 @@
// --+ options: json=compute, stochastic +--
var x1 x2 x1bar x2bar z y x u;
varexo ex1 ex2 ex1bar ex2bar ez ey ex eu;
parameters
rho_1 rho_2 rho_3 rho_4
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 c_z_dx2 c_z_u beta
lambda;
rho_1 = .9;
rho_2 = -.2;
rho_3 = .4;
rho_4 = -.3;
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;
c_z_dx2 = .3;
c_z_u = .3;
lambda = 0.5; // Share of optimizing agents.
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:u']
u = .5*u(-1) - .2*u(-2) + eu;
[name='eq:y']
y = rho_1*y(-1) + rho_2*y(-2) + ey;
[name='eq:x']
x = rho_3*x(-1) + rho_4*x(-2) + ex;
[name='eq:x1']
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']
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']
x1bar = x1bar(-1) + ex1bar;
[name='eq:x2bar']
x2bar = x2bar(-1) + ex2bar;
[name='zpac']
diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1)) + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + c_z_u*u + 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;
var ex = 0.1;
var eu = 0.05;
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
TrueData = simul_backward_model(initialconditions, 5000);
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.c_z_dx2 = .5;
eparams.lambda = .7;
// Define the dataset used for estimation
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols))
disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols))
skipline(2)
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.c_z_dx2 = .5;
eparams.lambda = .7;
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
c_z_dx2_nls = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact'));
lambda_nls = M_.params(strmatch('lambda', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls))
disp(sprintf('Estimate of lambda: %f', lambda_nls))
if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
end
if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
end
if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
end
if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)')
end
if abs(lambda_nls-lambda_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (lambda)')
end

View File

@ -0,0 +1,8 @@
#!/bin/sh
rm -rf example
rm -rf +example
rm -f example.log
rm -f *.mat
rm -f *.m
rm -f *.dat

View File

@ -0,0 +1,155 @@
// --+ options: json=compute, stochastic +--
var x1 x2 x1bar x2bar z y x u;
varexo ex1 ex2 ex1bar ex2bar ez ey ex eu;
parameters
rho_1 rho_2 rho_3 rho_4
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 c_z_dx2 c_z_u beta
lambda;
rho_1 = .9;
rho_2 = -.2;
rho_3 = .4;
rho_4 = -.3;
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;
c_z_dx2 = .3;
c_z_u = .3;
lambda = 0.5; // Share of optimizing agents.
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:u']
u = .5*u(-1) - .2*u(-2) + eu;
[name='eq:y']
y = rho_1*y(-1) + rho_2*y(-2) + ey;
[name='eq:x']
x = rho_3*x(-1) + rho_4*x(-2) + ex;
[name='eq:x1']
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']
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']
x1bar = x1bar(-1) + ex1bar;
[name='eq:x2bar']
x2bar = x2bar(-1) + ex2bar;
[name='zpac']
diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1)) + c_z_2*diff(z(-2)) + pac_expectation(pacman)) + (1-lambda)*( y + x) + c_z_dx2*diff(x2) + c_z_u*u + 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;
var ex = 0.1;
var eu = 0.05;
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
TrueData = simul_backward_model(initialconditions, 5000);
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.lambda = .7;
// Define the dataset used for estimation
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000);
e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols))
disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols))
skipline(2)
// Define a structure describing the parameters to be estimated (with initial conditions).
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
eparams.lambda = .7;
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
lambda_nls = M_.params(strmatch('lambda', M_.param_names, 'exact'));
disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
disp(sprintf('Estimate of lambda: %f', lambda_nls))
if abs(e_c_m_nls-e_c_m_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)')
end
if abs(c_z_1_nls-c_z_1_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)')
end
if abs(c_z_2_nls-c_z_2_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)')
end
if abs(lambda_nls-lambda_iterative_ols)>.01
error('Iterative OLS and NLS do not provide consistent estimates (lambda)')
end