stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

PAC decomposition between stationary and non-stationary components. 

The taget in PAC equation can be decomposed into an arbitrary number of components (variables
in the VAR auxiliary model).

TODO Iterative OLS estimation (which is not the preferred estimation routine).
TODO Decomposition in the routine evaluating the forecasts for each component.
pac-components
Stéphane Adjemian (Ryûk) 2021-11-25 16:16:42 +01:00
parent 4db2899868
commit b297353b06
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
21 changed files with 1568 additions and 207 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
matlab/+pac/+estimate/init.m
View File

@ -30,14 +30,28 @@ pattern = '(\(model_name\s*=\s*)(?<name>\w+)\)';
pacmodl = regexp(RHS, pattern, 'names'); pacmodl = regexp(RHS, pattern, 'names');
pacmodl = pacmodl.name; pacmodl = pacmodl.name;
pacmodel = M_.pac.(pacmodl);
% Get the transformed equation to be estimated. % Get the transformed equation to be estimated.
[lhs, rhs, json] = get_lhs_and_rhs(eqname, M_); [lhs, rhs, json] = get_lhs_and_rhs(eqname, M_);
% Get definition of aux. variable for pac expectation... % Get definition of aux. variable for pac expectation...
auxrhs = json.model{strmatch(sprintf('pac_expectation_%s', pacmodl), M_.lhs, 'exact')}.rhs; if isfield(pacmodel, 'aux_id')
auxrhs = {M_.lhs{pacmodel.aux_id}, json.model{pacmodel.aux_id}.rhs};
elseif isfield(pacmodel, 'components')
auxrhs = cell(length(pacmodel.components), 2);
for i=1:length(pacmodel.components)
auxrhs{i,1} = M_.lhs{pacmodel.components(i).aux_id};
auxrhs{i,2} = sprintf('(%s)', json.model{pacmodel.components(i).aux_id}.rhs);
end
else
error('Cannot find auxiliary variables for PAC expectation.')
end
% ... and substitute in rhs. % ... and substitute in rhs.
rhs = strrep(rhs, sprintf('pac_expectation_%s', pacmodl), auxrhs); for i=1:rows(auxrhs)
rhs = strrep(rhs, auxrhs{i,1}, auxrhs{i,2});
end
% Get pacmodel properties % Get pacmodel properties
pacmodel = M_.pac.(pacmodl); pacmodel = M_.pac.(pacmodl);

296
matlab/+pac/+update/parameters.m
View File

@ -66,56 +66,118 @@ if ~isfield(varcalib, 'CompanionMatrix') || any(isnan(varcalib.CompanionMatrix(:
end end
% Show the equations where this PAC model is used. % Show the equations where this PAC model is used.
fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.eq_name); if verbose
skipline() fprintf('PAC model %s is used in equation %s.\n', pacname, pacmodel.eq_name);
skipline()
end
% Do we need to decompose the PAC expectation?
if isfield(pacmodel, 'components')
numberofcomponents = length(pacmodel.components);
else
numberofcomponents = 0;
end
% Build the vector of PAC parameters (ECM parameter + autoregressive parameters). % Build the vector of PAC parameters (ECM parameter + autoregressive parameters).
pacvalues = DynareModel.params([pacmodel.ec.params; pacmodel.ar.params(1:pacmodel.max_lag)']); pacvalues = DynareModel.params([pacmodel.ec.params; pacmodel.ar.params(1:pacmodel.max_lag)']);
% Get the indices for the stationary/nonstationary variables in the VAR system.
id = find(strcmp(DynareModel.endo_names{pacmodel.ec.vars(pacmodel.ec.istarget)}, varmodel.list_of_variables_in_companion_var));
if isempty(id)
% Find the auxiliary variables if any
ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
if isempty(ad)
error('Cannot find the trend variable in the Companion VAR/VECM model.')
else
for i=1:length(ad)
auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)}));
if isequal(auxinfo.endo_index, pacmodel.ec.vars(pacmodel.ec.istarget))
id = ad(i);
break
end
if isequal(auxinfo.type, 8) && isequal(auxinfo.orig_index, pacmodel.ec.vars(pacmodel.ec.istarget))
id = ad(i);
break
end
end
end
if isempty(id)
error('Cannot find the trend variable in the Companion VAR/VECM model.')
end
end
% Infer the kind of PAC exoectation % Get the indices for the stationary/nonstationary variables in the VAR system.
if isequal(pacmodel.auxiliary_model_type, 'var') if numberofcomponents
if varmodel.nonstationary(id) id = cell(numberofcomponents, 1);
kind = 'dd'; for i=1:numberofcomponents
if varmodel.isconstant id(i) = {find(strcmp(DynareModel.endo_names{pacmodel.components(i).endo_var}, varmodel.list_of_variables_in_companion_var))};
id = id+1; if isempty(id{i})
end % Find the auxiliary variables if any
else ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
kind = 'll' if isempty(ad)
if varmodel.isconstant error('Cannot find the trend variable in the Companion VAR/VECM model.')
id = id+1; else
for j=1:length(ad)
auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(j)}));
if isequal(auxinfo.endo_index, pacmodel.components(i).endo_var)
id(i) = {ad(j)};
break
end
if isequal(auxinfo.type, 8) && isequal(auxinfo.orig_index, pacmodel.components(i).endo_var)
id(i) = {ad(j)};
break
end
end
end
if isempty(id{i})
error('Cannot find the trend variable in the Companion VAR/VECM model.')
end
end end
end end
else else
% Trend component model is assumed. id = {find(strcmp(DynareModel.endo_names{pacmodel.ec.vars(pacmodel.ec.istarget)}, varmodel.list_of_variables_in_companion_var))};
kind = 'dd'; if isempty(id{1})
% Find the auxiliary variables if any
ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
if isempty(ad)
error('Cannot find the trend variable in the Companion VAR/VECM model.')
else
for i=1:length(ad)
auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)}));
if isequal(auxinfo.endo_index, pacmodel.ec.vars(pacmodel.ec.istarget))
id = {ad(i)};
break
end
if isequal(auxinfo.type, 8) && isequal(auxinfo.orig_index, pacmodel.ec.vars(pacmodel.ec.istarget))
id = {ad(i)};
break
end
end
end
if isempty(id{1})
error('Cannot find the trend variable in the Companion VAR/VECM model.')
end
end
end
if ~numberofcomponents
% Infer the kind of PAC exoectation
if isequal(pacmodel.auxiliary_model_type, 'var')
if varmodel.nonstationary(id{1})
kind = {'dd'};
if varmodel.isconstant
id{1} = id{1}+1;
end
else
kind = {'ll'};
if varmodel.isconstant
id{1} = id{1}+1;
end
end
else
% Trend component model is assumed.
kind = {'dd'};
end
else
if varmodel.isconstant
for i=1:numberofcomponents
id{i} = id{i}+1;
end
end
end end
% Override kind with the information provided by the user or update M_.pac % Override kind with the information provided by the user or update M_.pac
if ~numberofcomponents
if ~isempty(pacmodel.kind)
kind = {pacmodel.kind};
else
pacmodel.kind = kind{1};
end
else
kind = cell(numberofcomponents,1);
for i=1:numberofcomponents
if isempty(pacmodel.components(i).kind)
error('kind declaration is mandatory for each component in pac_target_info.')
else
kind{i} = pacmodel.components(i).kind;
end
end
end
% Get the value of the discount factor. % Get the value of the discount factor.
beta = DynareModel.params(pacmodel.discount_index); beta = DynareModel.params(pacmodel.discount_index);
@ -125,6 +187,12 @@ if isfield(pacmodel, 'growth_str')
growth_flag = true; growth_flag = true;
else else
growth_flag = false; growth_flag = false;
for i=1:numberofcomponents
if isfield(pacmodel.components(i), 'growth_str')
growth_flag = true;
break
end
end
end end
% Do we have rule of thumb agents? γ is the share of optimizing agents. % Do we have rule of thumb agents? γ is the share of optimizing agents.
@ -136,84 +204,114 @@ end
% Get h vector (plus the parameter for the growth neutrality correction). % Get h vector (plus the parameter for the growth neutrality correction).
if growth_flag if growth_flag
[h, growthneutrality] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind, id); h = cell(1,length(id));
growthneutrality = cell(1,length(id));
for i=1:length(id)
[h{i}, growthneutrality{i}] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind{i}, id{i});
end
else else
h = hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind, id); h = cell(1,length(id));
for i=1:length(id)
h(i) = {hVectors([pacvalues; beta], varcalib.CompanionMatrix, pacmodel.auxiliary_model_type, kind{i}, id{i})};
end
end end
% Update M_.params with h % Update M_.params with h
if isequal(pacmodel.auxiliary_model_type, 'var') if isequal(pacmodel.auxiliary_model_type, 'var')
if DynareModel.var.(pacmodel.auxiliary_model_name).isconstant if DynareModel.var.(pacmodel.auxiliary_model_name).isconstant
DynareModel.params(pacmodel.h_param_indices) = h; if isfield(pacmodel, 'h_param_indices')
% No decomposition
DynareModel.params(pacmodel.h_param_indices) = h{1};
else
for i=1:numberofcomponents
DynareModel.params(pacmodel.components(i).h_param_indices) = h{i};
end
end
else else
DynareModel.params(pacmodel.h_param_indices(1)) = .0; if isfield(pacmodel, 'h_param_indices')
DynareModel.params(pacmodel.h_param_indices(2:end)) = h; % No decomposition
end DynareModel.params(pacmodel.h_param_indices(1)) = .0;
DynareModel.params(pacmodel.h_param_indices(2:end)) = h{1};
else
for i=1:numberofcomponents
DynareModel.params(pacmodel.components(i).h_param_indices(1)) = .0;
DynareModel.params(pacmodel.components(i).h_param_indices(2:end)) = h{i};
end
end
end % If the auxiliary model (VAR) has no constant.
else else
DynareModel.params(pacmodel.h_param_indices) = h; DynareModel.params(pacmodel.h_param_indices) = h{1};
end end % if auxiliary model is a VAR
% Update the parameter related to the growth neutrality correction. % Update the parameter related to the growth neutrality correction.
if growth_flag if growth_flag
% Growth neutrality as returned by hVector is valid iff % Growth neutrality as returned by hVector is valid iff
% there is no exogenous variables in the model and in the % there is no exogenous variables in the model and in the
% absence of non optimizing agents. % absence of non optimizing agents.
gg = -(growthneutrality-1); % Finite sum of autoregressive parameters + infinite sum of the coefficients in the PAC expectation term. for j=1:length(id)
cc = 1.0-gg*gamma; % First adjustment of the growth neutrality correction (should also be divided by gamma, done below at the end of this section). if isnan(growthneutrality{j})
% We may have to further change the correction if we have nonzero mean exogenous variables. continue
ll = 0.0;
if isfield(pacmodel, 'optim_additive')
% Exogenous variables are present in the λ part (optimizing agents).
tmp0 = 0;
for i=1:length(pacmodel.optim_additive.params)
if isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
tmp0 = tmp0 + pacmodel.optim_additive.scaling_factor(i);
elseif ~isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
tmp0 = tmp0 + DynareModel.params(pacmodel.optim_additive.params(i))*pacmodel.optim_additive.scaling_factor(i);
elseif ~islogical(pacmodel.optim_additive.bgp{i})
error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
end
end end
cc = cc - tmp0*gamma; gg = -(growthneutrality{j}-1); % Finite sum of autoregressive parameters + infinite sum of the coefficients in the PAC expectation term.
end cc = 1.0-gg*gamma; % First adjustment of the growth neutrality correction (should also be divided by gamma, done below at the end of this section).
if gamma<1 % We may have to further change the correction if we have nonzero mean exogenous variables.
if isfield(pacmodel, 'non_optimizing_behaviour') && isfield(pacmodel.non_optimizing_behaviour, 'params') ll = 0.0;
% Exogenous variables are present in the 1-λ part (rule of thumb agents). if isfield(pacmodel, 'optim_additive')
% Exogenous variables are present in the λ part (optimizing agents).
tmp0 = 0; tmp0 = 0;
tmp1 = 0; for i=1:length(pacmodel.optim_additive.params)
for i=1:length(pacmodel.non_optimizing_behaviour.params) if isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
if isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i} tmp0 = tmp0 + pacmodel.optim_additive.scaling_factor(i);
tmp0 = tmp0 + pacmodel.non_optimizing_behaviour.scaling_factor(i); elseif ~isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
elseif ~isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i} tmp0 = tmp0 + DynareModel.params(pacmodel.optim_additive.params(i))*pacmodel.optim_additive.scaling_factor(i);
tmp0 = tmp0 + DynareModel.params(pacmodel.non_optimizing_behaviour.params(i))*pacmodel.non_optimizing_behaviour.scaling_factor(i); elseif ~islogical(pacmodel.optim_additive.bgp{i})
elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && isnan(pacmodel.non_optimizing_behaviour.params(i)) error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.bgp{i};
elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && ~isnan(pacmodel.non_optimizing_behaviour.params(i))
tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.params(i)*pacmodel.non_optimizing_behaviour.bgp{i};
end end
end end
cc = cc - (1.0-gamma)*tmp0; cc = cc - tmp0*gamma;
ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end end
end if gamma<1
if isfield(pacmodel, 'additive') if isfield(pacmodel, 'non_optimizing_behaviour') && isfield(pacmodel.non_optimizing_behaviour, 'params')
% Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation. % Exogenous variables are present in the 1-λ part (rule of thumb agents).
tmp0 = 0; tmp0 = 0;
tmp1 = 0; tmp1 = 0;
for i=1:length(pacmodel.additive.params) for i=1:length(pacmodel.non_optimizing_behaviour.params)
if isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i} if isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
tmp0 = tmp0 + pacmodel.additive.scaling_factor(i); tmp0 = tmp0 + pacmodel.non_optimizing_behaviour.scaling_factor(i);
elseif ~isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i} elseif ~isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
tmp0 = tmp0 + DynareModel.params(pacmodel.additive.params(i))*pacmodel.additive.scaling_factor(i); tmp0 = tmp0 + DynareModel.params(pacmodel.non_optimizing_behaviour.params(i))*pacmodel.non_optimizing_behaviour.scaling_factor(i);
elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && isnan(pacmodel.additive.params(i)) elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && isnan(pacmodel.non_optimizing_behaviour.params(i))
tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.bgp{i}; tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.bgp{i};
elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && ~isnan(pacmodel.additive.params(i)) elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && ~isnan(pacmodel.non_optimizing_behaviour.params(i))
tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.params(i)*pacmodel.additive.bgp{i}; tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.params(i)*pacmodel.non_optimizing_behaviour.bgp{i};
end
end
cc = cc - (1.0-gamma)*tmp0;
ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end end
end end
cc = cc - tmp0; if isfield(pacmodel, 'additive')
ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation. % Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
tmp0 = 0;
tmp1 = 0;
for i=1:length(pacmodel.additive.params)
if isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
tmp0 = tmp0 + pacmodel.additive.scaling_factor(i);
elseif ~isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
tmp0 = tmp0 + DynareModel.params(pacmodel.additive.params(i))*pacmodel.additive.scaling_factor(i);
elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && isnan(pacmodel.additive.params(i))
tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.bgp{i};
elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && ~isnan(pacmodel.additive.params(i))
tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.params(i)*pacmodel.additive.bgp{i};
end
end
cc = cc - tmp0;
ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end
if isfield(pacmodel, 'growth_neutrality_param_index')
DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
else
DynareModel.params(pacmodel.components(j).growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
end
end end
DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model. end
end

2
matlab/aggregate.m
View File

@ -90,7 +90,7 @@ for i=1:length(varargin)
model{j} = strrep(model{j}, eqtagname{1}, eqtagname_); model{j} = strrep(model{j}, eqtagname{1}, eqtagname_);
end end
else else
model{j} = eqtagname_; model{j} = strcat('[', eqtagname_, ']');
end end
% Add equation tag with block name. % Add equation tag with block name.
if ~isempty(rootfolder) if ~isempty(rootfolder)

95
matlab/cherrypick.m
View File

@ -82,10 +82,10 @@ try
% Get the original equation. % Get the original equation.
[LHS, RHS] = get_lhs_and_rhs(eqtags{i}, M_, true, json); [LHS, RHS] = get_lhs_and_rhs(eqtags{i}, M_, true, json);
% Get the parameters, endogenous and exogenous variables in the current equation. % Get the parameters, endogenous and exogenous variables in the current equation.
[pnames, ~, xnames] = get_variables_and_parameters_in_equation(LHS, RHS, M_); [pnames, enames, xnames] = get_variables_and_parameters_in_equation(LHS, RHS, M_);
lhs_expression = LHS; lhs_expression = LHS;
LHS = get_variables_and_parameters_in_expression(LHS); LHS = get_variables_and_parameters_in_expression(LHS);
enames = LHS; enames = union(enames, LHS);
if length(LHS)>1 if length(LHS)>1
error('Expressions with more than one variable on the LHS are not allowed.') error('Expressions with more than one variable on the LHS are not allowed.')
end end
@ -111,7 +111,11 @@ try
end end
% Remove residual from equation if required. % Remove residual from equation if required.
if noresids if noresids
exogenous_variables_to_be_removed = ~ismember(xnames, M_.simulation_exo_names); if isfield(M_, 'simulation_exo_names')
exogenous_variables_to_be_removed = ~ismember(xnames, M_.simulation_exo_names);
else
exogenous_variables_to_be_removed = false(size(xnames));
end
if any(exogenous_variables_to_be_removed) if any(exogenous_variables_to_be_removed)
switch sum(exogenous_variables_to_be_removed) switch sum(exogenous_variables_to_be_removed)
case 1 case 1
@ -128,31 +132,77 @@ try
% Unroll expectation terms if any. % Unroll expectation terms if any.
isvar = regexp(RHS, 'var_expectation\(model_name = (?<name>\w+)\)', 'names'); isvar = regexp(RHS, 'var_expectation\(model_name = (?<name>\w+)\)', 'names');
ispac = regexp(RHS, 'pac_expectation\(model_name = (?<name>\w+)\)', 'names'); ispac = regexp(RHS, 'pac_expectation\(model_name = (?<name>\w+)\)', 'names');
istar = regexp(RHS, 'pac_target_nonstationary\(model_name = (?<name>\w+)\)', 'names');
if ~isempty(isvar) if ~isempty(isvar)
rhs = write_expectations(isvar.name, 'var'); rhs = write_expectations(isvar.name, 'var');
lhs = sprintf('%s_VE', eqtags{i}); auxlhs = sprintf('%s_VE', eqtags{i});
RHS = strrep(RHS, sprintf('var_expectation(model_name = %s)', isvar.name), lhs); RHS = strrep(RHS, sprintf('var_expectation(model_name = %s)', isvar.name), auxlhs);
end end
if ~isempty(ispac) if ~isempty(ispac)
[rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac'); if isfield(M_.pac.(ispac.name), 'components')
[rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac', false, false);
else
[rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac');
end
if ~isempty(rhs) if ~isempty(rhs)
lhs = sprintf('%s_PE', eqtags{i}); if iscell(rhs) % PAC expectations are decomposed.
if isempty(growthneutralitycorrection) auxlhs = cell(size(rhs));
for k=1:length(M_.pac.(ispac.name).components)
if isequal(k, 1)
lhs = M_.lhs{M_.pac.(ispac.name).components(k).aux_id};
auxlhs{k} = lhs;
if ~isempty(growthneutralitycorrection{k})
lhs = sprintf('%s+%s', lhs, growthneutralitycorrection{k});
end
if ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
if isempty(growthneutralitycorrection{k})
lhs = sprintf('%s*%s', M_.pac.(ispac.name).components(k).coeff_str, lhs);
else
lhs = sprintf('%s*(%s)', M_.pac.(ispac.name).components(k).coeff_str, lhs);
end
end
else
lhs_ = M_.lhs{M_.pac.(ispac.name).components(k).aux_id};
auxlhs{k} = lhs_;
if ~isempty(growthneutralitycorrection{k})
lhs_ = sprintf('%s+%s', lhs_, growthneutralitycorrection{k});
end
if ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
if isempty(growthneutralitycorrection{k})
lhs_ = sprintf('%s*%s', M_.pac.(ispac.name).components(k).coeff_str, lhs_);
else
lhs_ = sprintf('%s*(%s)', M_.pac.(ispac.name).components(k).coeff_str, lhs_);
end
end
lhs = sprintf('%s+%s', lhs, lhs_);
end
end
RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), lhs); RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), lhs);
else else
RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', lhs, growthneutralitycorrection)); auxlhs = M_.lhs{M_.pac.(ispac.name).aux_id};
if isempty(growthneutralitycorrection)
RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), auxlhs);
else
RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', auxlhs, growthneutralitycorrection));
end
end end
else else
% MCE version of the PAC equation. % MCE version of the PAC equation.
[rhs, growthneutralitycorrection] = write_pac_mce_expectations(eqtags{i}, ispac.name); [rhs, growthneutralitycorrection] = write_pac_mce_expectations(eqtags{i}, ispac.name);
lhs = sprintf('%s_Z', eqtags{i}); auxlhs = sprintf('%s_Z', eqtags{i});
if isempty(growthneutralitycorrection) if isempty(growthneutralitycorrection)
RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), lhs); RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), auxlhs);
else else
RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', lhs, growthneutralitycorrection)); RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', auxlhs, growthneutralitycorrection));
end end
end end
end end
if ~isempty(istar)
% Note that auxlhs and rhs are already defined in the previous block, it is not possible to be here if isempty(ispac) is true.
auxlhs{end+1} = M_.endo_names{M_.pac.(ispac.name).ec.vars(M_.pac.(ispac.name).ec.istarget)};
rhs{end+1} = M_.aux_vars(strmatch(auxlhs{end}, M_.endo_names, 'exact')==[M_.aux_vars(:).endo_index]).orig_expr;
RHS = strrep(RHS, sprintf('pac_target_nonstationary(model_name = %s)', ispac.name), sprintf('%s(-1)', auxlhs{end}));
end
% Print equation for unrolled PAC/VAR-expectation and update % Print equation for unrolled PAC/VAR-expectation and update
% list of parameters and endogenous variables (if any). % list of parameters and endogenous variables (if any).
if ~isempty(rhs) if ~isempty(rhs)
@ -160,13 +210,28 @@ try
% will not return the lhs variable in expectation_enames since % will not return the lhs variable in expectation_enames since
% the name is created on the fly and is not a member of M_.endo_names. % the name is created on the fly and is not a member of M_.endo_names.
expectation_pnames = get_variables_and_parameters_in_equation('', rhs, M_); expectation_pnames = get_variables_and_parameters_in_equation('', rhs, M_);
expectation_enames = get_variables_and_parameters_in_expression(lhs); expectation_enames = get_variables_and_parameters_in_expression(auxlhs);
expectation_xnames = get_variables_and_parameters_in_expression(rhs); expectation_xnames = get_variables_and_parameters_in_expression(rhs);
pnames = union(pnames, expectation_pnames); pnames = union(pnames, expectation_pnames);
xnames = union(xnames, setdiff(expectation_xnames, expectation_pnames)); xnames = union(xnames, setdiff(expectation_xnames, expectation_pnames));
enames = union(enames, expectation_enames); enames = union(enames, expectation_enames);
fprintf(fid, '[name=''%s'']\n', lhs); if ischar(rhs)
fprintf(fid, '%s = %s;\n\n', lhs, rhs); fprintf(fid, '[name=''%s'']\n', auxlhs);
fprintf(fid, '%s = %s;\n\n', auxlhs, rhs);
else
for k=1:length(rhs)
fprintf(fid, '[name=''%s'']\n', auxlhs{k});
fprintf(fid, '%s = %s;\n\n', auxlhs{k}, rhs{k});
end
end
if isfield(M_.pac.(ispac.name), 'components')
for k=1:length(M_.pac.(ispac.name).components)
if ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
params = get_variables_and_parameters_in_expression(M_.pac.(ispac.name).components(k).coeff_str);
pnames = union(pnames, params);
end
end
end
else else
pRHS = get_variables_and_parameters_in_equation('', RHS, M_); pRHS = get_variables_and_parameters_in_equation('', RHS, M_);
xRHS = get_variables_and_parameters_in_expression(RHS); xRHS = get_variables_and_parameters_in_expression(RHS);

29
matlab/get_variables_and_parameters_in_expression.m
View File

@ -8,7 +8,7 @@ function objects = get_variables_and_parameters_in_expression(expr)
% OUTPUTS % OUTPUTS
% - objects [cell] cell of row char arrays, names of the variables and parameters in expr. % - objects [cell] cell of row char arrays, names of the variables and parameters in expr.
% Copyright (C) 2020 Dynare Team % Copyright © 2020-2021 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -25,16 +25,27 @@ function objects = get_variables_and_parameters_in_expression(expr)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
objects = strsplit(expr, {'+','-','*','/','^', ... if iscell(expr)
objects = splitexpr(expr{1});
for i=2:length(expr)
objects = [objects, splitexpr(expr{i})];
end
else
objects = splitexpr(expr);
end
% Filter out the numbers, punctuation.
objects(cellfun(@(x) all(isstrprop(x, 'digit')+isstrprop(x, 'punct')), objects)) = [];
% Filter out empty elements.
objects(cellfun(@(x) all(isempty(x)), objects)) = [];
function o = splitexpr(expr)
o = strsplit(expr, {'+','-','*','/','^', ...
'log(', 'log10(', 'ln(', 'exp(', ... 'log(', 'log10(', 'ln(', 'exp(', ...
'sqrt(', 'abs(', 'sign(', ... 'sqrt(', 'abs(', 'sign(', ...
'sin(', 'cos(', 'tan(', 'asin(', 'acos(', 'atan(', ... 'sin(', 'cos(', 'tan(', 'asin(', 'acos(', 'atan(', ...
'min(', 'max(', ... 'min(', 'max(', ...
'normcdf(', 'normpdf(', 'erf(', ... 'normcdf(', 'normpdf(', 'erf(', ...
'diff(', 'adl(', '(', ')', '\n', '\t', ' '}); 'diff(', 'adl(', '(', ')', '\n', '\t', ' '});
% Filter out the numbers, punctuation.
objects(cellfun(@(x) all(isstrprop(x, 'digit')+isstrprop(x, 'punct')), objects)) = [];
% Filter out empty elements.
objects(cellfun(@(x) all(isempty(x)), objects)) = [];

186
matlab/print_expectations.m
View File

@ -127,13 +127,29 @@ switch expectationmodelkind
end end
case 'pac' case 'pac'
parameter_declaration = 'parameters'; parameter_declaration = 'parameters';
for i=1:length(expectationmodel.h_param_indices) if isfield(expectationmodel, 'h_param_indices')
parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.h_param_indices(i)}); for i=1:length(expectationmodel.h_param_indices)
parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.h_param_indices(i)});
end
else
for j=1:length(expectationmodel.components)
for i=1:length(expectationmodel.components(j).h_param_indices)
parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.components(j).h_param_indices(i)});
end
end
end end
fprintf(fid, '%s;\n\n', parameter_declaration); fprintf(fid, '%s;\n\n', parameter_declaration);
if withcalibration if withcalibration
for i=1:length(expectationmodel.h_param_indices) if isfield(expectationmodel, 'h_param_indices')
fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.h_param_indices(i)}, M_.params(expectationmodel.h_param_indices(i))); for i=1:length(expectationmodel.h_param_indices)
fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.h_param_indices(i)}, M_.params(expectationmodel.h_param_indices(i)));
end
else
for j=1:length(expectationmodel.components)
for i=1:length(expectationmodel.components(j).h_param_indices)
fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.components(j).h_param_indices(i)}, M_.params(expectationmodel.components(j).h_param_indices(i)));
end
end
end end
end end
if isfield(expectationmodel, 'growth_neutrality_param_index') if isfield(expectationmodel, 'growth_neutrality_param_index')
@ -145,10 +161,20 @@ switch expectationmodelkind
growth_correction = true; growth_correction = true;
else else
growth_correction = false; growth_correction = false;
if isfield(expectationmodel, 'components')
for j=1:length(expectationmodel.components)
if isfield(expectationmodel.components(j), 'growth_neutrality_param_index') && ~isempty(expectationmodel.components(j).growth_neutrality_param_index)
fprintf(fid, '\n');
fprintf(fid, 'parameters %s;\n\n', M_.param_names{expectationmodel.components(j).growth_neutrality_param_index});
if withcalibration
fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.components(j).growth_neutrality_param_index}, M_.params(expectationmodel.components(j).growth_neutrality_param_index));
end
growth_correction = true;
end
end
end
end end
otherwise
end end
fclose(fid); fclose(fid);
skipline() skipline()
@ -211,6 +237,13 @@ fprintf(fid, 'ds = dseries();\n\n');
id = 0; id = 0;
clear('expression');
% Get coefficient values in the target (if any)
if exist(sprintf('+%s/pac_target_coefficients.m', M_.fname), 'file')
targetcoefficients = feval(sprintf('%s.pac_target_coefficients', M_.fname), expectationmodelname, M_.params);
end
maxlag = max(auxmodel.max_lag); maxlag = max(auxmodel.max_lag);
if isequal(expectationmodel.auxiliary_model_type, 'trend_component') if isequal(expectationmodel.auxiliary_model_type, 'trend_component')
% Need to add a lag since the error correction equations are rewritten in levels. % Need to add a lag since the error correction equations are rewritten in levels.
@ -222,13 +255,23 @@ if isequal(expectationmodelkind, 'var')
end end
if isequal(expectationmodelkind, 'var') && isequal(expectationmodel.auxiliary_model_type, 'var') if isequal(expectationmodelkind, 'var') && isequal(expectationmodel.auxiliary_model_type, 'var')
% Constant in the VAR auxiliary model
id = id+1; id = id+1;
expression = sprintf('%1.16f', M_.params(expectationmodel.param_indices(id))); expression = sprintf('%1.16f', M_.params(expectationmodel.param_indices(id)));
end end
if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var') if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
% Constant in the VAR auxiliary model
id = id+1; id = id+1;
expression = sprintf('%1.16f', M_.params(expectationmodel.h_param_indices(id))); if isfield(expectationmodel, 'h_param_indices')
constant = M_.params(expectationmodel.h_param_indices(id));
else
constant = 0;
for j=1:length(expectationmodel.components)
constant = constant + targetcoefficients(j)*M_.params(expectationmodel.components(j).h_param_indices(id));
end
end
expression = sprintf('%1.16f', constant);
end end
for i=1:maxlag for i=1:maxlag
@ -256,7 +299,14 @@ for i=1:maxlag
case 'var' case 'var'
parameter = M_.params(expectationmodel.param_indices(id)); parameter = M_.params(expectationmodel.param_indices(id));
case 'pac' case 'pac'
parameter = M_.params(expectationmodel.h_param_indices(id)); if isfield(expectationmodel, 'h_param_indices')
parameter = M_.params(expectationmodel.h_param_indices(id));
else
parameter = 0;
for k=1:length(expectationmodel.components)
parameter = parameter+targetcoefficients(k)*M_.params(expectationmodel.components(k).h_param_indices(id));
end
end
otherwise otherwise
end end
switch expectationmodelkind switch expectationmodelkind
@ -275,61 +325,115 @@ for i=1:maxlag
variable = sprintf('%s.%s()', variable, transformations{k}); variable = sprintf('%s.%s()', variable, transformations{k});
end end
end end
if isequal(id, 1) if exist('expression','var')
if isequal(expectationmodelkind, 'pac') && growth_correction if parameter>=0
pgrowth = M_.params(expectationmodel.growth_neutrality_param_index); expression = sprintf('%s+%1.16f*%s', expression, parameter, variable);
linearCombination = ''; elseif parameter<0
for iter = 1:numel(expectationmodel.growth_linear_comb) expression = sprintf('%s-%1.16f*%s', expression, -parameter, variable);
end
else
if parameter>=0
expression = sprintf('%1.16f*%s', parameter, variable);
elseif parameter<0
expression = sprintf('-%1.16f*%s', -parameter, variable);
end
end
end
end
if isequal(expectationmodelkind, 'pac') && growth_correction
if isfield(expectationmodel, 'growth_neutrality_param_index')
pgrowth = M_.params(expectationmodel.growth_neutrality_param_index);
for iter = 1:numel(expectationmodel.growth_linear_comb)
vgrowth='';
if expectationmodel.growth_linear_comb(iter).exo_id > 0
vgrowth = strcat('dbase.', M_.exo_names{expectationmodel.growth_linear_comb(iter).exo_id});
elseif expectationmodel.growth_linear_comb(iter).endo_id > 0
vgrowth = strcat('dbase.', M_.endo_names{expectationmodel.growth_linear_comb(iter).endo_id});
end
if expectationmodel.growth_linear_comb(iter).lag ~= 0
vgrowth = sprintf('%s(%d)', vgrowth, expectationmodel.growth_linear_comb(iter).lag);
end
if expectationmodel.growth_linear_comb(iter).param_id > 0
if ~isempty(vgrowth)
vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.growth_linear_comb(iter).param_id), vgrowth);
else
vgrowth = num2str(M_.params(expectationmodel.growth_linear_comb(iter).param_id), '%1.16f');
end
end
if abs(expectationmodel.growth_linear_comb(iter).constant) ~= 1
if ~isempty(vgrowth)
vgrowth = sprintf('%1.16f*%s', expectationmodel.growth_linear_comb(iter).constant, vgrowth);
else
vgrowth = num2str(expectationmodel.growth_linear_comb(iter).constant, '%1.16f');
end
end
if iter > 1
if expectationmodel.growth_linear_comb(iter).constant > 0
linearCombination = sprintf('%s+%s', linearCombination, vgrowth);
else
linearCombination = sprintf('%s-%s', linearCombination, vgrowth);
end
else
linearCombination = vgrowth;
end
end % loop over growth linear combination elements
growthcorrection = sprintf('%1.16f*(%s)', pgrowth, linearCombination);
else
first = true;
for i=1:length(expectationmodel.components)
if ~isequal(expectationmodel.components(i).kind, 'll') && isfield(expectationmodel.components(i), 'growth_neutrality_param_index') && isfield(expectationmodel.components(i), 'growth_linear_comb') && ~isempty(expectationmodel.components(i).growth_linear_comb)
pgrowth = targetcoefficients(i)*M_.params(expectationmodel.components(i).growth_neutrality_param_index);
for iter = 1:numel(expectationmodel.components(i).growth_linear_comb)
vgrowth=''; vgrowth='';
if expectationmodel.growth_linear_comb(iter).exo_id > 0 if expectationmodel.components(i).growth_linear_comb(iter).exo_id > 0
vgrowth = strcat('dbase.', M_.exo_names{expectationmodel.growth_linear_comb(iter).exo_id}); vgrowth = strcat('dbase.', M_.exo_names{expectationmodel.components(i).growth_linear_comb(iter).exo_id});
elseif expectationmodel.growth_linear_comb(iter).endo_id > 0 elseif expectationmodel.components(i).growth_linear_comb(iter).endo_id > 0
vgrowth = strcat('dbase.', M_.endo_names{expectationmodel.growth_linear_comb(iter).endo_id}); vgrowth = strcat('dbase.', M_.endo_names{expectationmodel.components(i).growth_linear_comb(iter).endo_id});
% TODO Check if we should not substitute auxiliary variables with original transformed variables here.
end end
if expectationmodel.growth_linear_comb(iter).lag ~= 0 if expectationmodel.components(i).growth_linear_comb(iter).lag ~= 0
vgrowth = sprintf('%s(%d)', vgrowth, expectationmodel.growth_linear_comb(iter).lag); vgrowth = sprintf('%s(%d)', vgrowth, expectationmodel.components(i).growth_linear_comb(iter).lag);
end end
if expectationmodel.growth_linear_comb(iter).param_id > 0 if expectationmodel.components(i).growth_linear_comb(iter).param_id > 0
if ~isempty(vgrowth) if ~isempty(vgrowth)
vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.growth_linear_comb(iter).param_id), vgrowth); vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.components(i).growth_linear_comb(iter).param_id), vgrowth);
else else
vgrowth = num2str(M_.params(expectationmodel.growth_linear_comb(iter).param_id), '%1.16f'); vgrowth = num2str(M_.params(expectationmodel.components(i).growth_linear_comb(iter).param_id), '%1.16f');
end end
end end
if abs(expectationmodel.growth_linear_comb(iter).constant) ~= 1 if abs(expectationmodel.components(i).growth_linear_comb(iter).constant) ~= 1
if ~isempty(vgrowth) if ~isempty(vgrowth)
vgrowth = sprintf('%1.16f*%s', expectationmodel.growth_linear_comb(iter).constant, vgrowth); vgrowth = sprintf('%1.16f*%s', expectationmodel.components(i).growth_linear_comb(iter).constant, vgrowth);
else else
vgrowth = num2str(expectationmodel.growth_linear_comb(iter).constant, '%1.16f'); vgrowth = num2str(expectationmodel.components(i).growth_linear_comb(iter).constant, '%1.16f');
end end
end end
if iter > 1 if iter > 1
if expectationmodel.growth_linear_comb(iter).constant > 0 if expectationmodel.components(i).growth_linear_comb(iter).constant > 0
linearCombination = sprintf('%s+%s', linearCombination, vgrowth); linearCombination = sprintf('%s+%s', linearCombination, vgrowth);
else else
linearCombination = sprintf('%s-%s', linearCombination, vgrowth); linearCombination = sprintf('%s-%s', linearCombination, vgrowth);
end end
else else
linearCombination=vgrowth; linearCombination = vgrowth;
end
end % loop over growth linear combination elements
if first
growthcorrection = sprintf('%1.16f*(%s)', pgrowth, linearCombination);
first = false;
else
if pgrowth>0
growthcorrection = sprintf('%s+%1.16f*(%s)', growthcorrection, pgrowth, linearCombination);
elseif pgrowth<0
growthcorrection = sprintf('%s-%1.16f*(%s)', growthcorrection, -pgrowth, linearCombination);
end end
end end
if parameter >= 0
expression = sprintf('%1.16f*(%s)+%1.16f*%s', pgrowth, linearCombination, parameter, variable);
else
expression = sprintf('%1.16f*(%s)-%1.16f*%s', pgrowth, linearCombination, -parameter, variable);
end
else
expression = sprintf('%1.16f*%s', parameter, variable);
end
else
if parameter>=0
expression = sprintf('%s+%1.16f*%s', expression, parameter, variable);
else
expression = sprintf('%s-%1.16f*%s', expression, -parameter, variable);
end end
end end
end end
end expression = sprintf('%s+%s', expression, growthcorrection);
end % growth_correction
fprintf(fid, 'ds.%s = %s;', expectationmodelname, expression); fprintf(fid, 'ds.%s = %s;', expectationmodelname, expression);
fclose(fid); fclose(fid);

160
matlab/write_expectations.m
View File

@ -1,4 +1,4 @@
function [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, iscrlf) function [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, iscrlf, aggregate)
% Prints the exansion of the VAR_EXPECTATION or PAC_EXPECTATION term in files. % Prints the exansion of the VAR_EXPECTATION or PAC_EXPECTATION term in files.
% %
@ -48,6 +48,16 @@ end
if nargin<3 if nargin<3
iscrlf = false; iscrlf = false;
aggregate = true;
end
if nargin<4
aggregate = true;
end
if isfield(expectationmodel, 'h_param_indices')
% Disaggregation requires components...
aggregate = true;
end end
% Get the name of the associated VAR model and test its existence. % Get the name of the associated VAR model and test its existence.
@ -82,7 +92,29 @@ end
if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var') if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
id = id+1; id = id+1;
expression = sprintf('%s', M_.param_names{expectationmodel.h_param_indices(id)}); if isfield(expectationmodel, 'h_param_indices')
expression = sprintf('%s', M_.param_names{expectationmodel.h_param_indices(id)});
else
if aggregate
if isequal(expectationmodel.components(1).coeff_str, '1')
expression = sprintf('%s', M_.param_names{expectationmodel.components(1).h_param_indices(id)});
else
expression = sprintf('%s*%s', expectationmodel.components(1).coeff_str, M_.param_names{expectationmodel.components(1).h_param_indices(id)});
end
for i=2:length(expectationmodel.components)
if isequal(expectationmodel.components(i).coeff_str, '1')
expression = sprintf('%s+%s', expression, M_.param_names{expectationmodel.components(i).h_param_indices(id)});
else
expression = sprintf('%s+%s*%s', expression, expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).h_param_indices(id)});
end
end
else
expression = cell(length(expectationmodel.components), 1);
for i=1:length(expectationmodel.components)
expression(i) = {M_.param_names{expectationmodel.components(i).h_param_indices(id)}};
end
end
end
end end
for i=1:maxlag for i=1:maxlag
@ -110,7 +142,31 @@ for i=1:maxlag
case 'var' case 'var'
parameter = M_.param_names{expectationmodel.param_indices(id)}; parameter = M_.param_names{expectationmodel.param_indices(id)};
case 'pac' case 'pac'
parameter = M_.param_names{expectationmodel.h_param_indices(id)}; if isfield(expectationmodel, 'h_param_indices')
parameter = M_.param_names{expectationmodel.h_param_indices(id)};
else
if aggregate
% TODO Check if we can have parameters entering with a minus sign in the linear combination defining the target.
if isequal(expectationmodel.components(1).coeff_str, '1')
parameter = M_.param_names{expectationmodel.components(1).h_param_indices(id)};
else
parameter = sprintf('%s*%s', expectationmodel.components(1).coeff_str, M_.param_names{expectationmodel.components(1).h_param_indices(id)});
end
for k=2:length(expectationmodel.components)
if isequal(expectationmodel.components(k).coeff_str, '1')
parameter = sprintf('%s+%s', parameter, M_.param_names{expectationmodel.components(k).h_param_indices(id)});
else
parameter = sprintf('%s+%s*%s', parameter, expectationmodel.components(k).coeff_str, M_.param_names{expectationmodel.components(k).h_param_indices(id)});
end
end
parameter = sprintf('(%s)', parameter);
else
parameter = cell(length(expectationmodel.components), 1);
for k=1:length(expectationmodel.components)
parameter(k) = {M_.param_names{expectationmodel.components(k).h_param_indices(id)}};
end
end
end
otherwise otherwise
end end
switch expectationmodelkind switch expectationmodelkind
@ -128,21 +184,47 @@ for i=1:maxlag
end end
end end
if isequal(id, 1) if isequal(id, 1)
if iscrlf if aggregate
expression = sprintf('%s*%s\n', parameter, variable); if iscrlf
expression = sprintf('%s*%s\n', parameter, variable);
else
expression = sprintf('%s*%s', parameter, variable);
end
else else
expression = sprintf('%s*%s', parameter, variable); for k=1:length(expectationmodel.components)
if iscrlf
expression(k) = {sprintf('%s*%s\n', parameter{k}, variable)};
else
expression(k) = {sprintf('%s*%s', parameter{k}, variable)};
end
end
end end
else else
if iscrlf if aggregate
expression = sprintf('%s + %s*%s\n', expression, parameter, variable); if iscrlf
expression = sprintf('%s + %s*%s\n', expression, parameter, variable);
else
expression = sprintf('%s + %s*%s', expression, parameter, variable);
end
else else
expression = sprintf('%s + %s*%s', expression, parameter, variable); for k=1:length(expectationmodel.components)
if iscrlf
expression(k) = {sprintf('%s + %s*%s\n', expression{k}, parameter{k}, variable)};
else
expression(k) = {sprintf('%s + %s*%s', expression{k}, parameter{k}, variable)};
end
end
end end
end end
end end
end end
if aggregate
growthneutralitycorrection = ''
else
growthneutralitycorrection = {};
end
if isfield(expectationmodel, 'growth_neutrality_param_index') if isfield(expectationmodel, 'growth_neutrality_param_index')
if numel(expectationmodel.growth_linear_comb) == 1 if numel(expectationmodel.growth_linear_comb) == 1
growthneutralitycorrection = sprintf('%s*%s', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str); growthneutralitycorrection = sprintf('%s*%s', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str);
@ -150,7 +232,65 @@ if isfield(expectationmodel, 'growth_neutrality_param_index')
growthneutralitycorrection = sprintf('%s*(%s)', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str); growthneutralitycorrection = sprintf('%s*(%s)', M_.param_names{expectationmodel.growth_neutrality_param_index}, expectationmodel.growth_str);
end end
else else
growthneutralitycorrection = ''; if isfield(expectationmodel, 'components')
if aggregate
growthneutralitycorrection = '';
for i=1:length(expectationmodel.components)
if ~isequal(expectationmodel.components(i).kind, 'll')
if isfield(expectationmodel.components(i), 'growth_neutrality_param_index')
if isempty(growthneutralitycorrection)
if ~isempty(expectationmodel.components(i).growth_str)
if isequal(expectationmodel.components(i).coeff_str, '1')
if numel(expectationmodel.components(i).growth_linear_comb) == 1
growthneutralitycorrection = sprintf('%s*%s', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
else
growthneutralitycorrection = sprintf('%s*(%s)', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
end
else
if numel(expectationmodel.components(i).growth_linear_comb) == 1
growthneutralitycorrection = sprintf('%s*%s*%s', expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
else
growthneutralitycorrection = sprintf('%s*%s*(%s)', expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
end
end
end
else
if ~isempty(expectationmodel.components(i).growth_str)
if isequal(expectationmodel.components(i).coeff_str, '1')
if numel(expectationmodel.components(i).growth_linear_comb) == 1
growthneutralitycorrection = sprintf('%s+%s*%s', growthneutralitycorrection, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
else
growthneutralitycorrection = sprintf('%s+%s*(%s)', growthneutralitycorrection, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
end
else
if numel(expectationmodel.components(i).growth_linear_comb) == 1
growthneutralitycorrection = sprintf('%s+%s*%s*%s', growthneutralitycorrection, expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
else
growthneutralitycorrection = sprintf('%s+%s*%s*(%s)', growthneutralitycorrection, expectationmodel.components(i).coeff_str, M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str);
end
end
end
end
end % if growth neutrality correction for this component
end % if non stationary component
end
else
growthneutralitycorrection = repmat({''}, length(expectationmodel.components), 1);
for i=1:length(growthneutralitycorrection)
if ~isequal(expectationmodel.components(i).kind, 'll')
if isfield(expectationmodel.components(i), 'growth_neutrality_param_index')
if ~isempty(expectationmodel.components(i).growth_str)
if numel(expectationmodel.components(i).growth_linear_comb) == 1
growthneutralitycorrection(i) = {sprintf('%s*%s', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str)};
else
growthneutralitycorrection(i) = {sprintf('%s*(%s)', M_.param_names{expectationmodel.components(i).growth_neutrality_param_index}, expectationmodel.components(i).growth_str)};
end
end
end % if growth neutrality correction for this component
end % if non stationary component
end
end % if aggregate
end
end end
if nargout==1 && ~isempty(growthneutralitycorrection) if nargout==1 && ~isempty(growthneutralitycorrection)

2
preprocessor

@ -1 +1 @@
Subproject commit 7dde09169ecd4b5e32b3cd35bef88958a262ec11 Subproject commit a8fce06dc46ca09329e6899e1fde47f2cd81809b

2
scripts/dynare.el
View File

@ -102,7 +102,7 @@
"osr_params_bounds" "observation_trends" "deterministic_trends" "optim_weights" "osr_params_bounds" "observation_trends" "deterministic_trends" "optim_weights"
"homotopy_setup" "conditional_forecast_paths" "svar_identification" "homotopy_setup" "conditional_forecast_paths" "svar_identification"
"moment_calibration" "irf_calibration" "ramsey_constraints" "generate_irfs" "moment_calibration" "irf_calibration" "ramsey_constraints" "generate_irfs"
"matched_moments" "occbin_constraints" "model_replace" "verbatim") "matched_moments" "occbin_constraints" "model_replace" "verbatim" "pac_target_info")
"Dynare block keywords.")) "Dynare block keywords."))
;; Mathematical functions and operators used in model equations (see "hand_side" in Bison file) ;; Mathematical functions and operators used in model equations (see "hand_side" in Bison file)

15
tests/Makefile.am
View File

@ -518,6 +518,15 @@ ECB_MODFILES = \
pac/var-6/substitution.mod \ pac/var-6/substitution.mod \
pac/var-7/example1.mod \ pac/var-7/example1.mod \
pac/var-7/substitution.mod \ pac/var-7/substitution.mod \
pac/var-8/example1.mod \
pac/var-8/substitution.mod \
pac/var-8e/example1.mod \
pac/var-9/example1.mod \
pac/var-9/substitution.mod \
pac/var-9e/example1.mod \
pac/var-10e/example1.mod \
pac/var-10e/example2.mod \
pac/var-11e/example1.mod \
pac/trend-component-1/example1.mod \ pac/trend-component-1/example1.mod \
pac/trend-component-2/example1.mod \ pac/trend-component-2/example1.mod \
pac/trend-component-3/example1.mod \ pac/trend-component-3/example1.mod \
@ -878,6 +887,12 @@ pac/var-6/substitution.m.trs: pac/var-6/example1.m.trs
pac/var-6/substitution.o.trs: pac/var-6/example1.o.trs pac/var-6/substitution.o.trs: pac/var-6/example1.o.trs
pac/var-7/substitution.m.trs: pac/var-7/example1.m.trs pac/var-7/substitution.m.trs: pac/var-7/example1.m.trs
pac/var-7/substitution.o.trs: pac/var-7/example1.o.trs pac/var-7/substitution.o.trs: pac/var-7/example1.o.trs
pac/var-8/substitution.m.trs: pac/var-8/example1.m.trs
pac/var-8/substitution.o.trs: pac/var-8/example1.o.trs
pac/var-9/substitution.m.trs: pac/var-9/example1.m.trs
pac/var-9/substitution.o.trs: pac/var-9/example1.o.trs
pac/var-10e/example2.m.trs: pac/var-10e/example1.m.trs
pac/var-10e/example2.o.trs: pac/var-10e/example1.o.trs
pac/trend-component-14/substitution.m.trs: pac/trend-component-14/example1.m.trs pac/trend-component-14/substitution.m.trs: pac/trend-component-14/example1.m.trs
pac/trend-component-14/substitution.o.trs: pac/trend-component-14/example1.o.trs pac/trend-component-14/substitution.o.trs: pac/trend-component-14/example1.o.trs

10
tests/ecb/cherrypick/test1.mod
View File

@ -45,7 +45,7 @@ c_z_2 = -.1;
c_z_dx2 = .3; c_z_dx2 = .3;
c_z_u = .3; c_z_u = .3;
c_z_dv = .4; c_z_dv = .4;
c_z_s = -.2; c_z_s = -.2;
cx = 1.0; cx = 1.0;
cy = 1.0; cy = 1.0;
@ -54,7 +54,7 @@ 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']); 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); pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman, auxname=rototo);
model; model;
@ -125,7 +125,7 @@ pac.update.expectation('pacman');
// variables and exogenous variables in .inc files under ./simulation-files folder. Note that // variables and exogenous variables in .inc files under ./simulation-files folder. Note that
// innovations ex1bar and ex2bar will not appear in the equations. // innovations ex1bar and ex2bar will not appear in the equations.
verbatim; verbatim;
cherrypick('test1', 'simulation-files1', {'z1', 'z2', 'z3'}, true); cherrypick('test1', 'simulation-files-1-a', {'z1', 'z2', 'z3'}, true);
cherrypick('test1', 'simulation-files2', {'zpac', 'eq:x1', 'eq:x', 'eq:y'}, true); cherrypick('test1', 'simulation-files-1-b', {'zpac', 'eq:x1', 'eq:x', 'eq:y'}, true);
aggregate('toto.mod', {}, '', 'simulation-files1', 'simulation-files2'); aggregate('toto1.mod', {}, '', 'simulation-files-1-a', 'simulation-files-1-b');
end; end;

110
tests/ecb/cherrypick/test2.mod
View File

@ -1,36 +1,100 @@
var x y z; // --+ options: json=compute, stochastic +--
model; var y x z v u w z1 z2 z3;
x = y(-1)+z(1); varexo ex ey ez eu ew;
y = x-z; parameters a_y_1 a_y_2 b_y_1 b_y_2 a_u_1 a_u_2 b_u_1 b_u_2 b_x_1 b_x_2 d_y;
diff(z) = 0; parameters 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;
a_u_1 = .2;
a_u_2 = .3;
b_u_1 = .1;
b_u_2 = .4;
b_x_1 = -.1;
b_x_2 = -.2;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
auxname pv_y_;
kind ll;
component x;
auxname pv_dx_;
growth diff(x(-1));
kind dl;
end; end;
model;
if ~isfield(M_, 'exo_names') [name='eq:u']
error('M_.exo_names not defined.') u = a_u_1*u(-1) + a_u_2*diff(x(-1)) + b_u_1*y(-2) + b_u_2*diff(x(-2)) + eu ;
end
if ~iscell(M_.exo_names) [name='eq:y']
error('M_.exo_names not a cell array.') y = a_y_1*y(-1) + a_y_2*diff(x(-1)) + b_y_1*y(-2) + b_y_2*diff(x(-2)) + ey ;
end
if ~isempty(M_.exo_names) [name='eq:x']
error('M_.exo_names not an empty cell array.') diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
end
if ~isfield(M_, 'param_names') [name='eq:v']
error('M_.param_names not defined.') v = x + d_y*y ;
end
if ~iscell(M_.param_names) [name='eq:pac']
error('M_.param_names not a cell array.') diff(z) = e_c_m*(pac_target_nonstationary(pacman)-z(-1)) + c_z_1*diff(z(-1)) + c_z_2*diff(z(-2)) + pac_expectation(pacman) + ez;
end
if ~isempty(M_.param_names) [name='eq:w']
error('M_.param_names not an empty cell array.') diff(w) = -.1*w + ew ;
end
[name='eq:z1']
z1 = z+y-x+u;
[name='eq:z2']
z2 = z-y+x-u;
[name='eq:z3']
z3 = u-diff(v);
end;
shocks;
var ex = 1.0;
var ey = 1.0;
var ez = 1.0;
var ew = 1.0;
var eu = 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');
// Select a subset of the equations and print the equations, the list of parameters, endogenous
// variables and exogenous variables in .inc files under ./simulation-files folder. Note that
// innovations ex1bar and ex2bar will not appear in the equations.
verbatim;
cherrypick('test2', 'simulation-files-2-a', {'eq:z1', 'eq:z2', 'eq:z3'}, true);
cherrypick('test2', 'simulation-files-2-b', {'eq:pac', 'eq:x', 'eq:y'}, true);
aggregate('toto2.mod', {}, '', 'simulation-files-2-a', 'simulation-files-2-b');
end;

103
tests/pac/var-10e/example1.mod Normal file
View File

@ -0,0 +1,103 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
growth diff(y(-1));
auxname pv_y_;
kind dd;
component x;
growth diff(x(-1));
auxname pv_dx_;
kind dd;
end;
model(use_dll);
[name='eq:y']
diff(y) = .01 + a_y_1*diff(y(-1)) + a_y_2*diff(x(-1)) + b_y_1*diff(y(-2)) + b_y_2*diff(x(-2)) + ey ;
[name='eq:x']
diff(x) = .05 + b_x_1*diff(y(-2)) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[name='zpac']
diff(z) = e_c_m*(pac_target_nonstationary(pacman)-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 500 periods
TrueData = simul_backward_model(initialconditions, 500);
TrueData.save('example1.data')
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
// Define the dataset used for estimation
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+400, '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'));
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))
save('example1.estimation.mat', 'e_c_m_nls', 'c_z_1_nls', 'c_z_2_nls')

104
tests/pac/var-10e/example2.mod Normal file
View File

@ -0,0 +1,104 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y', 'eq:v']);
pac_model(auxiliary_model_name=toto, discount=beta, kind=dd, growth=diff(v(-1)), model_name=pacman);
model(use_dll);
[name='eq:y']
diff(y) = .01 + a_y_1*diff(y(-1)) + a_y_2*diff(x(-1)) + b_y_1*diff(y(-2)) + b_y_2*diff(x(-2)) + ey ;
[name='eq:x']
diff(x) = .05 + b_x_1*diff(y(-2)) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
diff(v) = diff(x) + d_y*diff(y) ;
[name='zpac']
diff(z) = e_c_m*(v(-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 500 periods
TrueData = simul_backward_model(initialconditions, 500);
TrueData.save('example2.data')
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
// Define the dataset used for estimation
edata = TrueData;
edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+400, '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'));
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))
% Check consistency with disaggregated target
ts1 = dseries('example1.data.mat');
ts2 = dseries('example2.data.mat');
if max(abs(ts1.z.data-ts2.z.data))>1e-12
error('Simulations in example1 and example2 are not consistent.')
end
e1 = load('example1.estimation.mat');
if abs(e_c_m_nls-e1.e_c_m_nls)>1e-6 || abs(c_z_1_nls-e1.c_z_1_nls)>1e-6 || abs(c_z_2_nls-e1.c_z_2_nls)>1e-6
error('Estimations in example1 and example2 are not consistent.')
end
delete example1.data.mat
delete example2.data.mat
delete example1.estimation.mat

100
tests/pac/var-11e/example1.mod Normal file
View File

@ -0,0 +1,100 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y cst; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
cst = .1;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
auxname pv_y_;
kind ll;
component x;
growth diff(x(-1));
auxname pv_dx_;
kind dd;
end;
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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = .01 + x + d_y*y ;
[name='zpac']
diff(z) = e_c_m*(pac_target_nonstationary(pacman)-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 500 periods
TrueData = simul_backward_model(initialconditions, 5000);
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
// Define the dataset used for estimation
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'));
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))

90
tests/pac/var-8/example1.mod Normal file
View File

@ -0,0 +1,90 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
auxname pv_y_;
kind ll;
component x;
auxname pv_dx_;
kind dd;
end;
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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[name='eq:pac']
diff(z) = e_c_m*(pac_target_nonstationary(pacman)-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 500 periods
TrueData = simul_backward_model(initialconditions, 20);
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
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 = example1.set_auxiliary_variables(y, x, M_.params);
y = [y(find(M_.lead_lag_incidence(1,:))); y];
[residual, g1] = example1.dynamic(y, x', M_.params, oo_.steady_state, 1);
save('example1.mat', 'residual', 'g1', 'TrueData');
end;

83
tests/pac/var-8/substitution.mod Normal file
View File

@ -0,0 +1,83 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
@#include "example1/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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[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 "example1/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];
example1 = load('example1.mat');
[residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
if max(abs(example1.TrueData.data(:)-TrueData.data(:)))>1e-9
error('Simulations do not match.')
end
if ~isequal(length(residual), length(example1.residual)) || max(abs(example1.residual-residual))>1e-8
warning('Residuals do not match!')
end
if ~isequal(length(g1(:)), length(example1.g1(:))) || max(abs(example1.g1(:)-g1(:)))>1e-8
warning('Jacobian matrices do not match!')
end
end;

94
tests/pac/var-8e/example1.mod Normal file
View File

@ -0,0 +1,94 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
auxname pv_y_;
kind ll;
component x;
auxname pv_dx_;
kind dd;
end;
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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[name='zpac']
diff(z) = e_c_m*(pac_target_nonstationary(pacman)-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 500 periods
TrueData = simul_backward_model(initialconditions, 5000);
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
// Define the dataset used for estimation
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'));
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))

94
tests/pac/var-9/example1.mod Normal file
View File

@ -0,0 +1,94 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
auxname pv_y_;
kind ll;
component x;
growth diff(x(-1));
auxname pv_dx_;
kind dd;
end;
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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[name='eq:pac']
diff(z) = e_c_m*(pac_target_nonstationary(pacman)-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 500 periods
TrueData = simul_backward_model(initialconditions, 20);
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
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 = example1.set_auxiliary_variables(y, x, M_.params);
y = [y(find(M_.lead_lag_incidence(1,:))); y];
[residual, g1] = example1.dynamic(y, x', M_.params, oo_.steady_state, 1);
save('example1.mat', 'residual', 'g1', 'TrueData');
end;

84
tests/pac/var-9/substitution.mod Normal file
View File

@ -0,0 +1,84 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
@#include "example1/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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[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 "example1/model/pac-expectations/pacman-expression.inc"
+
@#include "example1/model/pac-expectations/pacman-growth-neutrality-correction.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];
example1 = load('example1.mat');
[residual, g1] = substitution.dynamic(y, x', M_.params, oo_.steady_state, 1);
if max(abs(example1.TrueData.data(:)-TrueData.data(:)))>1e-9
error('Simulations do not match.')
end
if ~isequal(length(residual), length(example1.residual)) || max(abs(example1.residual-residual))>1e-8
warning('Residuals do not match!')
end
if ~isequal(length(g1(:)), length(example1.g1(:))) || max(abs(example1.g1(:)-g1(:)))>1e-8
warning('Jacobian matrices do not match!')
end
end;

98
tests/pac/var-9e/example1.mod Normal file
View File

@ -0,0 +1,98 @@
// --+ options: json=compute, stochastic +--
var y x z v;
varexo ex ey ez ;
parameters a_y_1 a_y_2 b_y_1 b_y_2 b_x_1 b_x_2 d_y; // VAR parameters
parameters 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;
d_y = .5;
beta = .9;
e_c_m = .1;
c_z_1 = .7;
c_z_2 = -.3;
var_model(model_name=toto, structural, eqtags=['eq:x', 'eq:y']);
pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman);
pac_target_info(pacman);
target v;
auxname_target_nonstationary vns;
component y;
auxname pv_y_;
kind ll;
component x;
growth diff(x(-1));
auxname pv_dx_;
kind dd;
end;
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']
diff(x) = b_x_1*y(-2) + b_x_2*diff(x(-1)) + ex ;
[name='eq:v']
v = x + d_y*y ;
[name='zpac']
diff(z) = e_c_m*(pac_target_nonstationary(pacman)-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 500 periods
TrueData = simul_backward_model(initialconditions, 5000);
// Print expanded PAC_EXPECTATION term.
pac.print('pacman', 'eq:pac');
clear eparams
eparams.e_c_m = .9;
eparams.c_z_1 = .5;
eparams.c_z_2 = .2;
// Define the dataset used for estimation
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'));
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))