Added Iterative OLS routine for PAC equation.

time-shift
Stéphane Adjemian(Charybdis) 2018-07-23 16:12:46 +02:00
parent f03a9b98e1
commit 0edebbe5d7
3 changed files with 255 additions and 87 deletions

View File

@ -0,0 +1,105 @@
function [pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, pnames_, ipnames_, data, islaggedvariables] = ...
init(M_, oo_, eqname, params, data, range)
% Copyright (C) 2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Get the original equation to be estimated
[LHS, RHS] = get_lhs_and_rhs(eqname, M_, true);
% Check that the equation has a PAC expectation term.
if ~contains(RHS, 'pac_expectation', 'IgnoreCase', true)
error('This is not a PAC equation.')
end
% Get the name of the PAC model.
pattern = '(\(model_name\s*=\s*)(?<name>\w+)\)';
pacmodl = regexp(RHS, pattern, 'names');
pacmodl = pacmodl.name;
% Get the transformed equation to be estimated.
[lhs, rhs] = get_lhs_and_rhs(eqname, M_);
% Get the parameters and variables in the PAC equation.
[pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_);
% Get the list of estimated parameters
pnames_ = fieldnames(params);
% Check that the estimated parameters are used in the PAC equation.
ParametersNotInPAC = setdiff(pnames_, pnames);
if ~isempty(ParametersNotInPAC)
skipline()
if length(ParametersNotInPAC)>1
list = sprintf(' %s\n', ParametersNotInPAC{:});
remk = sprintf(' The following parameters:\n\n%s\n do not appear in the PAC equation.', list);
else
remk = sprintf(' Parameter %s does not appear in the PAC equation.', ParametersNotInPAC{1});
end
disp(remk)
skipline()
error('The estimated parameters must be used in equation %s.', eqname)
end
% Get indices of estimated parameters.
ipnames_ = zeros(size(pnames_));
for i=1:length(ipnames_)
ipnames_(i) = strmatch(pnames_{i}, M_.param_names, 'exact');
end
% Add the auxiliary variables in the dataset.
data = feval([M_.fname '.dynamic_set_auxiliary_series'], data, M_.params);
% Check that the data for endogenous variables have values.
if any(isnan(data{enames{:}}(range).data(:)))
error('Some variable values are missing in the database.')
end
% Set the number of exogenous variables.
xnbr = length(xnames);
% Test if we have a residual and get its name (-> rname).
if isequal(xnbr, 1)
rname = M_.exo_names{strmatch(xnames{1}, M_.exo_names, 'exact')};
irname = 1;
if ~all(isnan(data{xnames{1}}.data))
error('The residual (%s) must have NaN values in the provided database.', xnames{1})
end
else
% We have observed exogenous variables in the PAC equation.
tmp = data{xnames{:}}(range).data;
idx = find(all(~isnan(tmp))); % Indices for the observed exogenous variables.
if isequal(length(idx), length(xnames))
error('There is no residual in this equation, all the exogenous variables are observed.')
else
if length(idx)<length(xnames)-1
error('It is not allowed to have more than one residual in a PAC equation')
end
irname = setdiff(1:length(xnames), idx);
rname = xnames{irname};
end
end
% Remove residuals from the equation.
%
% Note that a plus or minus will remain in the equation, but this seems to
% be without consequence.
rhs = regexprep(rhs, rname, '');
% Create a dummy variable
islaggedvariables = contains(rhs, '(-1)');

View File

@ -0,0 +1,146 @@
function iterative_ols(eqname, params, data, range)
% Estimates the parameters of a PAC equation by Iterative Ordinary Least Squares.
%
% INPUTS
% - eqname [string] Name of the pac equation.
% - params [struct] Describes the parameters to be estimated.
% - data [dseries] Database for the estimation
% - range [dates] Range of dates for the estimation.
%
% OUTPUTS
% - none
%
% REMARKS
% [1] The estimation results are printed in the command line window, and the
% parameters are updated accordingly in M_.params.
% [2] The second input is a structure. Each fieldname corresponds to the
% name of an estimated parameter, the value of the field is the initial
% guess used for the estimation (by NLS).
% [3] The third input is a dseries object which must at least contain all
% the variables appearing in the estimated equation. The residual of the
% equation must have NaN values in the object.
% [4] It is assumed that the residual is additive.
% Copyright (C) 2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_
[pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, pnames_, ipnames_, data, islaggedvariables] = ...
pac.estimate.init(M_, oo_, eqname, params, data, range);
% Set initial condition.
params0 = cell2mat(struct2cell(params));
% Build PAC expectation matrix expression.
dataForPACExpectation0 = dseries();
listofvariables0 = {};
if ~isempty(M_.pac.(pacmodl).h0_param_indices)
for i=1:length(M_.pac.(pacmodl).h0_param_indices)
match = regexp(rhs, sprintf('(?<var>((\\w*)|\\w*\\(-1\\)))\\*%s', M_.param_names{M_.pac.(pacmodl).h0_param_indices(i)}), 'names');
if isempty(match)
match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*)|\\w*\\(-1\\)))', M_.param_names{M_.pac.(pacmodl).h0_param_indices(i)}), 'names');
end
if isempty(strfind(match.var, '(-1)'))
listofvariables0{i} = match.var;
dataForPACExpectation0 = [dataForPACExpectation0, data{listofvariables0{i}}];
else
listofvariables0{i} = match.var(1:end-4);
dataForPACExpectation0 = [dataForPACExpectation0, data{match.var(1:end-4)}.lag(1)];
end
end
end
dataForPACExpectation1 = dseries();
listofvariables1 = {};
if ~isempty(M_.pac.(pacmodl).h1_param_indices)
for i=1:length(M_.pac.(pacmodl).h1_param_indices)
match = regexp(rhs, sprintf('(?<var>((\\w*)|\\w*\\(-1\\)))\\*%s', M_.param_names{M_.pac.(pacmodl).h1_param_indices(i)}), 'names');
if isempty(match)
match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*)|\\w*\\(-1\\)))', M_.param_names{M_.pac.(pacmodl).h1_param_indices(i)}), 'names');
end
if isempty(strfind(match.var, '(-1)'))
listofvariables1{i} = match.var;
dataForPACExpectation1 = [dataForPACExpectation1, data{listofvariables1{i}}];
else
listofvariables1{i} = match.var(1:end-4);
dataForPACExpectation1 = [dataForPACExpectation1, data{match.var(1:end-4)}.lag(1)];
end
end
end
% Set correction for growth neutrality
correction = 0;
if isfield(M_.pac.(pacmodl), 'growth_type')
switch M_.pac.(pacmodl).growth_type
case 'parameter'
correction = M_.params(M_.pac.(pacmodl).growth_index)*M_.params(M_.pac.(pacmodl).growth_neutrality_param_index);
otherwise
error('Not yet implemented.')
end
end
% Build matrix for EC and AR terms.
DataForOLS = dseries();
DataForOLS{'ec-term'} = data{M_.endo_names{M_.pac.(pacmodl).ec.vars(1)}}.lag(1)-data{M_.endo_names{M_.pac.(pacmodl).ec.vars(2)}}.lag(1);
listofvariables3 = {'ec-term'};
xparm = { M_.param_names(M_.pac.(pacmodl).ec.params(1))};
for i = 1:length(M_.pac.(pacmodl).ar.params)
if islagof(M_.pac.(pacmodl).ar.vars(i), M_.pac.(pacmodl).lhs_var)
DataForOLS = [DataForOLS, data{M_.endo_names{M_.pac.(pacmodl).ar.vars(i)}}];
listofvariables3{i+1} = M_.endo_names{M_.pac.(pacmodl).ar.vars(i)};
xparm{i+1} = M_.param_names(M_.pac.(pacmodl).ar.params(i));
end
end
XDATA = DataForOLS{listofvariables3{:}}(range).data;
params0
% Iterative OLS
noconvergence = true;
counter = 0;
while noconvergence
counter = counter+1;
% Update the PAC parameter,s.
for i=1:length(ipnames_)
M_.params(ipnames_(i)) = params0(i);
end
M_ = pac.update.parameters(pacmodl, M_, oo_);
% Compute lhs conditional on params0
PACExpectations = 0;
if ~isempty(listofvariables0)
PACExpectations = dataForPACExpectation0{listofvariables0{:}}(range).data*M_.params(M_.pac.pacman.h0_param_indices);
end
if ~isempty(listofvariables1)
PACExpectations = dataForPACExpectation1{listofvariables1{:}}(range).data*M_.params(M_.pac.pacman.h1_param_indices);
end
YDATA = data{M_.endo_names{M_.pac.(pacmodl).lhs_var}}(range).data-correction-PACExpectations;
% Do OLS
params1 = XDATA\YDATA;
noconvergence = max(abs(params0-params1))>1e-6;
disp(sprintf('Iter. %u,\t crit: %10.5f', counter, max(abs(params0-params1))))
params0 = params1;
end
% Update M_.params
for i=1:length(params0)
M_.params(ipnames_(i)) = params0(i);
end
M_ = pac.update.parameters(pacmodl, M_, oo_);

View File

@ -41,24 +41,8 @@ function nls(eqname, params, data, range)
global M_ oo_
% Get the original equation to be estimated
[LHS, RHS] = get_lhs_and_rhs(eqname, M_, true);
% Check that the equation has a PAC expectation term.
if ~contains(RHS, 'pac_expectation', 'IgnoreCase', true)
error('This is not a PAC equation.')
end
% Get the name of the PAC model.
pattern = '(\(model_name\s*=\s*)(?<name>\w+)\)';
pacmodl = regexp(RHS, pattern, 'names');
pacmodl = pacmodl.name;
% Get the transformed equation to be estimated.
[lhs, rhs] = get_lhs_and_rhs(eqname, M_);
% Get the parameters and variables in the PAC equation.
[pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_);
[pacmodl, lhs, rhs, pnames, enames, xnames, pid, eid, xid, pnames_, ipnames_, data, islaggedvariables] = ...
pac.estimate.init(M_, oo_, eqname, params, data, range);
% List of objects to be replaced
objNames = [pnames; enames; xnames];
@ -70,75 +54,6 @@ objNames = objNames(I);
objIndex = objIndex(I);
objTypes = objTypes(I);
% Get the list of estimated parameters
pnames_ = fieldnames(params);
% Check that the estimated parameters are used in the PAC equation.
ParametersNotInPAC = setdiff(pnames_, pnames);
if ~isempty(ParametersNotInPAC)
skipline()
if length(ParametersNotInPAC)>1
list = sprintf(' %s\n', ParametersNotInPAC{:});
remk = sprintf(' The following parameters:\n\n%s\n do not appear in the PAC equation.', list);
else
remk = sprintf(' Parameter %s does not appear in the PAC equation.', ParametersNotInPAC{1});
end
disp(remk)
skipline()
error('The estimated parameters must be used in equation %s.', eqname)
end
% Get indices of estimated parameters.
ipnames_ = zeros(size(pnames_));
for i=1:length(ipnames_)
ipnames_(i) = strmatch(pnames_{i}, M_.param_names, 'exact');
end
% Add the auxiliary variables in the dataset.
data = feval([M_.fname '.dynamic_set_auxiliary_series'], data, M_.params);
% Check that the data for endogenous variables have values. Note that we
% also need data in range(1)-1 for the lagged variables, and we do not test
% for them here.
if any(isnan(data{enames{:}}(range).data(:)))
error('Some variable values are missing in the database.')
end
% Set the number of exogenous variables.
xnbr = length(xnames);
% Test if we have a residual and get its name (-> rname).
if isequal(xnbr, 1)
rname = M_.exo_names{strmatch(xnames{1}, M_.exo_names, 'exact')};
irname = 1;
if ~all(isnan(data{xnames{1}}.data))
error('The residual (%s) must have NaN values in the provided database.', xnames{1})
end
else
% We have observed exogenous variables in the PAC equation.
tmp = data{xnames{:}}(range).data;
idx = find(all(~isnan(tmp))); % Indices for the observed exogenous variables.
if isequal(length(idx), length(xnames))
error('There is no residual in this equation, all the exogenous variables are observed.')
else
if length(idx)<length(xnames)-1
error('It is not allowed to have more than one residual in a PAC equation')
end
irname = setdiff(1:length(xnames), idx);
rname = xnames{irname};
end
end
% Remove residuals from the equation.
%
% Note that a plus or minus will remain in the equation, but this seems to
% be without consequence.
rhs = regexprep(rhs, rname, '');
% Create a dummy variable
islaggedvariables = contains(rhs, '(-1)');
% Substitute parameters and variables.
for i=1:length(objNames)
switch objTypes(i)
@ -212,3 +127,5 @@ pparams1 = fminunc(ssr, params0, options);
for i=1:length(pparams1)
M_.params(ipnames_(i)) = pparams1(i);
end
M_ = pac.update.parameters(pacmodl, M_, oo_);