sur: fix up and use common code to create matrices
parent
d8f26525b0
commit
f8c0282b01
|
@ -1,5 +1,5 @@
|
||||||
function [Y, lhssub, X, startdates, enddates, startidxs, residnames, pbeta, vars, surpidxs, surconestrainedparams] = common_parsing(ds, ast, jsonmodel, overlapping_dates)
|
function [Y, lhssub, X, startdates, enddates, startidxs, residnames, pbeta, vars] = common_parsing(ds, ast, jsonmodel, overlapping_dates)
|
||||||
%function [Y, lhssub, X, startdates, enddates, startidxs, residnames, pbeta, vars, surpidxs, surconstrainedparams] = common_parsing(ds, ast, jsonmodel, overlapping_dates)
|
%function [Y, lhssub, X, startdates, enddates, startidxs, residnames, pbeta, vars] = common_parsing(ds, ast, jsonmodel, overlapping_dates)
|
||||||
%
|
%
|
||||||
% Code common to sur.m and pooled_ols.m
|
% Code common to sur.m and pooled_ols.m
|
||||||
%
|
%
|
||||||
|
@ -84,6 +84,7 @@ end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
return
|
return
|
||||||
%%
|
%%
|
||||||
|
|
||||||
|
|
|
@ -9,7 +9,7 @@ function [Y, lhssub, X] = parse_ols_style_equation(ds, ast, line)
|
||||||
% line [int] equation line number
|
% line [int] equation line number
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% Y [dseries] LHS of the equation
|
% Y [dseries] LHS of the equation (with lhssub subtracted)
|
||||||
% lhssub [dseries] RHS subtracted from LHS
|
% lhssub [dseries] RHS subtracted from LHS
|
||||||
% X [dseries] RHS of the equation
|
% X [dseries] RHS of the equation
|
||||||
%
|
%
|
||||||
|
@ -203,7 +203,7 @@ end
|
||||||
function [param, X] = parseTimesNodeHelper(ds, node, line, param, X)
|
function [param, X] = parseTimesNodeHelper(ds, node, line, param, X)
|
||||||
if isOlsParamExpr(node, line)
|
if isOlsParamExpr(node, line)
|
||||||
param = assignParam(param, node, line);
|
param = assignParam(param, node, line);
|
||||||
elseif isOlsVarExpr(node, line)
|
elseif isOlsVarExpr(ds, node, line)
|
||||||
if isempty(X)
|
if isempty(X)
|
||||||
X = evalNode(ds, node, line, X);
|
X = evalNode(ds, node, line, X);
|
||||||
else
|
else
|
||||||
|
@ -237,7 +237,7 @@ else
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
function tf = isOlsVar(node)
|
function tf = isOlsVar(ds, node)
|
||||||
if strcmp(node.node_type, 'VariableNode') ...
|
if strcmp(node.node_type, 'VariableNode') ...
|
||||||
&& (strcmp(node.type, 'endogenous') ...
|
&& (strcmp(node.type, 'endogenous') ...
|
||||||
|| (strcmp(node.type, 'exogenous') && any(strcmp(ds.name, node.name))))
|
|| (strcmp(node.type, 'exogenous') && any(strcmp(ds.name, node.name))))
|
||||||
|
@ -251,11 +251,11 @@ else
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
function tf = isOlsVarExpr(node, line)
|
function tf = isOlsVarExpr(ds, node, line)
|
||||||
if strcmp(node.node_type, 'VariableNode') || strcmp(node.node_type, 'UnaryOpNode')
|
if strcmp(node.node_type, 'VariableNode') || strcmp(node.node_type, 'UnaryOpNode')
|
||||||
tf = isOlsVar(node);
|
tf = isOlsVar(ds, node);
|
||||||
elseif strcmp(node.node_type, 'BinaryOpNode')
|
elseif strcmp(node.node_type, 'BinaryOpNode')
|
||||||
tf = isOlsVarExpr(node.arg1, line) || isOlsVarExpr(node.arg2, line);
|
tf = isOlsVarExpr(ds, node.arg1, line) || isOlsVarExpr(ds, node.arg2, line);
|
||||||
else
|
else
|
||||||
ols_error(['got unexpected type ' node.node_type], line);
|
ols_error(['got unexpected type ' node.node_type], line);
|
||||||
end
|
end
|
||||||
|
|
|
@ -0,0 +1,67 @@
|
||||||
|
function [Yvec, Xmat, constrained] = put_in_sur_form(Y, X)
|
||||||
|
%function [Yvec, Xmat, constrained] = put_in_sur_form(Y, X)
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% Y [cell array] dependent variables
|
||||||
|
% X [cell array] regressors
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% Yvec [vector] dependent variables
|
||||||
|
% Xmat [matrix] regressors
|
||||||
|
% constrained [cellstr] names of parameters that were constrained
|
||||||
|
%
|
||||||
|
% SPECIAL REQUIREMENTS
|
||||||
|
% none
|
||||||
|
|
||||||
|
% Copyright (C) 2019 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/>.
|
||||||
|
|
||||||
|
%% Check inputs
|
||||||
|
if nargin ~= 2
|
||||||
|
error('put_in_sur_form expects 2 arguments');
|
||||||
|
end
|
||||||
|
|
||||||
|
if isempty(Y) || ~iscell(Y) || isempty(X) || ~iscell(X) || length(Y) ~= length(X)
|
||||||
|
error('put_in_sur_form arguments should be cells of the same size');
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Organize output
|
||||||
|
nobs = size(X{1}, 1);
|
||||||
|
neqs = length(X);
|
||||||
|
Xmat = dseries([X{1}.data; zeros(nobs*(neqs-1), size(X{1}, 2))], X{1}.firstdate, X{1}.name);
|
||||||
|
Yvec = Y{1};
|
||||||
|
constrained = {};
|
||||||
|
for i = 2:neqs
|
||||||
|
to_remove = [];
|
||||||
|
Xtmp = dseries([zeros(nobs*(i-1), size(X{i}, 2)); X{i}.data; zeros(nobs*(neqs-i), size(X{i}, 2))], X{i}.firstdate, X{i}.name);
|
||||||
|
for j = 1:length(X{i}.name)
|
||||||
|
idx = find(strcmp(Xmat.name, X{i}.name{j}));
|
||||||
|
if ~isempty(idx)
|
||||||
|
Xmat.(Xmat.name{idx}) = Xmat{idx} + Xtmp{j};
|
||||||
|
to_remove = [to_remove j];
|
||||||
|
constrained{end+1} = Xmat.name{idx};
|
||||||
|
end
|
||||||
|
end
|
||||||
|
for j = length(to_remove):-1:1
|
||||||
|
Xtmp = Xtmp.remove(Xtmp.name{j});
|
||||||
|
end
|
||||||
|
if ~isempty(Xtmp)
|
||||||
|
Xmat = [Xmat Xtmp];
|
||||||
|
end
|
||||||
|
Yvec = dseries([Yvec.data; Y{i}.data], Yvec.firstdate);
|
||||||
|
end
|
||||||
|
end
|
133
matlab/ols/sur.m
133
matlab/ols/sur.m
|
@ -14,7 +14,7 @@ function varargout = sur(ds, param_names, eqtags)
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% dynare must be run with the option: json=compute
|
% dynare must be run with the option: json=compute
|
||||||
|
|
||||||
% Copyright (C) 2017-2018 Dynare Team
|
% Copyright (C) 2017-2019 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -34,113 +34,100 @@ function varargout = sur(ds, param_names, eqtags)
|
||||||
global M_ oo_ options_
|
global M_ oo_ options_
|
||||||
|
|
||||||
%% Check input argument
|
%% Check input argument
|
||||||
assert(nargin >= 1 && nargin <= 3, 'You must provide one, two, or three arguments');
|
assert(nargin >= 1 && nargin <= 3, 'sur() takes between 1 and 3 arguments');
|
||||||
assert(~isempty(ds) && isdseries(ds), 'The first argument must be a dseries');
|
|
||||||
if nargin >= 2
|
if nargin < 3
|
||||||
assert(iscellstr(param_names), 'The 2nd argument must be a cellstr');
|
eqtags = {};
|
||||||
else
|
end
|
||||||
|
|
||||||
|
if nargin < 2
|
||||||
param_names = {};
|
param_names = {};
|
||||||
|
else
|
||||||
|
assert(iscellstr(param_names), 'sur: the 2nd argument must be a cellstr');
|
||||||
end
|
end
|
||||||
|
|
||||||
%% Read JSON
|
%% Get Equation(s)
|
||||||
jsonfile = [M_.fname filesep() 'model' filesep() 'json' filesep() 'modfile-original.json'];
|
[ast, jsonmodel] = get_ast_jsonmodel(eqtags);
|
||||||
if exist(jsonfile, 'file') ~= 2
|
neqs = length(jsonmodel);
|
||||||
error('Could not find %s! Please use the json=compute option (See the Dynare invocation section in the reference manual).', jsonfile);
|
|
||||||
end
|
|
||||||
|
|
||||||
jsonmodel = loadjson(jsonfile);
|
|
||||||
jsonmodel = jsonmodel.model;
|
|
||||||
if nargin == 3
|
|
||||||
jsonmodel = getEquationsByTags(jsonmodel, 'name', eqtags);
|
|
||||||
end
|
|
||||||
|
|
||||||
%% Find parameters and variable names in equations and setup estimation matrices
|
%% Find parameters and variable names in equations and setup estimation matrices
|
||||||
[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, opidxs, surconstrainedparams] = ...
|
[Y, ~, X] = common_parsing(ds, ast, jsonmodel, true);
|
||||||
pooled_sur_common(ds, jsonmodel);
|
clear ast jsonmodel;
|
||||||
|
nobs = Y{1}.nobs;
|
||||||
|
[Y, X, constrained] = put_in_sur_form(Y, X);
|
||||||
|
|
||||||
if nargin == 1 && size(X, 2) ~= M_.param_nbr
|
if nargin == 1 && size(X, 2) ~= M_.param_nbr
|
||||||
warning(['Not all parameters were used in model: ' ...
|
warning(['Not all parameters were used in model: ' strjoin(setdiff(M_.param_names, X.name), ', ')]);
|
||||||
sprintf('%s', strjoin(setdiff(M_.param_names, pbeta), ', '))]);
|
|
||||||
end
|
end
|
||||||
|
|
||||||
%% Force equations to have the same sample range
|
% constrained_param_idxs: indexes in X.name of parameters that were constrained
|
||||||
maxfp = max([startdates{:}]);
|
constrained_param_idxs = zeros(length(constrained), 1);
|
||||||
minlp = min([enddates{:}]);
|
for i = 1:length(constrained)
|
||||||
nobs = minlp - maxfp;
|
constrained_param_idxs(i, 1) = find(strcmp(X.name, constrained{i}));
|
||||||
newY = zeros(nobs*length(jsonmodel), 1);
|
|
||||||
newX = zeros(nobs*length(jsonmodel), columns(X));
|
|
||||||
lastidx = 1;
|
|
||||||
for i = 1:length(jsonmodel)
|
|
||||||
if i == length(jsonmodel)
|
|
||||||
yds = dseries(Y(startidxs(i):end), startdates{i});
|
|
||||||
xds = dseries(X(startidxs(i):end, :), startdates{i});
|
|
||||||
else
|
|
||||||
yds = dseries(Y(startidxs(i):startidxs(i+1)-1), startdates{i});
|
|
||||||
xds = dseries(X(startidxs(i):startidxs(i+1)-1, :), startdates{i});
|
|
||||||
end
|
|
||||||
newY(lastidx:lastidx + nobs, 1) = yds(maxfp:minlp).data;
|
|
||||||
newX(lastidx:lastidx + nobs, :) = xds(maxfp:minlp, :).data;
|
|
||||||
if i ~= length(jsonmodel)
|
|
||||||
lastidx = lastidx + nobs + 1;
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
constrained_params_str = strjoin(X.name(constrained_param_idxs), ', ');
|
||||||
|
|
||||||
if ~isempty(param_names)
|
if ~isempty(param_names)
|
||||||
pnamesall = M_.param_names(opidxs);
|
newX = dseries();
|
||||||
nparams = length(param_names);
|
nparams = length(param_names);
|
||||||
pidxs = zeros(nparams, 1);
|
pidxs = zeros(nparams, 1);
|
||||||
for i = 1:nparams
|
for i = 1:nparams
|
||||||
idxs = find(strcmp(param_names{i}, pnamesall));
|
idx = find(strcmp(param_names{i}, X.name));
|
||||||
if isempty(idxs)
|
if isempty(idx)
|
||||||
if ~isempty(eqtags)
|
if ~isempty(eqtags)
|
||||||
error(['Could not find ' param_names{i} ...
|
error(['Could not find ' param_names{i} ...
|
||||||
' in the provided equations specified by ' strjoin(eqtags, ',')]);
|
' in the provided equations specified by ' strjoin(eqtags, ',')]);
|
||||||
end
|
end
|
||||||
error('Unspecified error. Please report');
|
error('Unspecified error. Please report');
|
||||||
end
|
end
|
||||||
pidxs(i) = idxs;
|
pidxs(i) = idx;
|
||||||
|
newX = [newX X.(X.name{idx})];
|
||||||
end
|
end
|
||||||
vars = [vars{:}];
|
subcols = setdiff(1:length(X.name), pidxs);
|
||||||
vars = {vars(pidxs)};
|
for i = length(subcols):-1:1
|
||||||
newY = newY - newX(:, setdiff(1:size(newX, 2), pidxs)) * M_.params(setdiff(opidxs, opidxs(pidxs), 'stable'));
|
Y = Y - M_.params(strcmp(X.name{subcols(i)}, M_.param_names))*X.(X.name{subcols(i)});
|
||||||
newX = newX(:, pidxs);
|
end
|
||||||
opidxs = opidxs(pidxs);
|
X = newX;
|
||||||
|
end
|
||||||
|
|
||||||
|
% opidxs: indexes in M_.params associated with columns of X
|
||||||
|
opidxs = zeros(length(X.name), 1);
|
||||||
|
for i = 1:length(X.name)
|
||||||
|
opidxs(i, 1) = find(strcmp(X.name{i}, M_.param_names));
|
||||||
end
|
end
|
||||||
|
|
||||||
%% Return to surgibbs if called from there
|
%% Return to surgibbs if called from there
|
||||||
st = dbstack(1);
|
st = dbstack(1);
|
||||||
if strcmp(st(1).name, 'surgibbs')
|
if strcmp(st(1).name, 'surgibbs')
|
||||||
varargout{1} = length(maxfp:minlp); %dof
|
varargout{1} = nobs; %dof
|
||||||
varargout{2} = opidxs;
|
varargout{2} = opidxs;
|
||||||
varargout{3} = newX;
|
varargout{3} = X.data;
|
||||||
varargout{4} = newY;
|
varargout{4} = Y.data;
|
||||||
varargout{5} = length(jsonmodel);
|
varargout{5} = neqs;
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
Y = newY;
|
|
||||||
X = newX;
|
|
||||||
oo_.sur.dof = length(maxfp:minlp);
|
|
||||||
|
|
||||||
%% Estimation
|
%% Estimation
|
||||||
|
oo_.sur.dof = nobs;
|
||||||
|
|
||||||
% Estimated Parameters
|
% Estimated Parameters
|
||||||
[q, r] = qr(X, 0);
|
[q, r] = qr(X.data, 0);
|
||||||
xpxi = (r'*r)\eye(size(X, 2));
|
xpxi = (r'*r)\eye(size(X.data, 2));
|
||||||
resid = Y - X * (r\(q'*Y));
|
resid = Y.data - X.data * (r\(q'*Y.data));
|
||||||
resid = reshape(resid, oo_.sur.dof, length(jsonmodel));
|
resid = reshape(resid, oo_.sur.dof, neqs);
|
||||||
|
|
||||||
M_.Sigma_e = resid'*resid/oo_.sur.dof;
|
M_.Sigma_e = resid'*resid/oo_.sur.dof;
|
||||||
kLeye = kron(chol(inv(M_.Sigma_e)), eye(oo_.sur.dof));
|
kLeye = kron(chol(inv(M_.Sigma_e)), eye(oo_.sur.dof));
|
||||||
[q, r] = qr(kLeye*X, 0);
|
[q, r] = qr(kLeye*X.data, 0);
|
||||||
oo_.sur.beta = r\(q'*kLeye*Y);
|
oo_.sur.beta = r\(q'*kLeye*Y.data);
|
||||||
|
|
||||||
M_.params(opidxs) = oo_.sur.beta;
|
M_.params(opidxs) = oo_.sur.beta;
|
||||||
|
|
||||||
% Yhat
|
% Yhat
|
||||||
oo_.sur.Yhat = X * oo_.sur.beta;
|
oo_.sur.Yhat = X.data * oo_.sur.beta;
|
||||||
|
|
||||||
% Residuals
|
% Residuals
|
||||||
oo_.sur.resid = Y - oo_.sur.Yhat;
|
oo_.sur.resid = Y.data - oo_.sur.Yhat;
|
||||||
|
|
||||||
%% Calculate statistics
|
%% Calculate statistics
|
||||||
% Estimate for sigma^2
|
% Estimate for sigma^2
|
||||||
|
@ -148,7 +135,7 @@ SS_res = oo_.sur.resid'*oo_.sur.resid;
|
||||||
oo_.sur.s2 = SS_res/oo_.sur.dof;
|
oo_.sur.s2 = SS_res/oo_.sur.dof;
|
||||||
|
|
||||||
% R^2
|
% R^2
|
||||||
ym = Y - mean(Y);
|
ym = Y.data - mean(Y.data);
|
||||||
SS_tot = ym'*ym;
|
SS_tot = ym'*ym;
|
||||||
oo_.sur.R2 = 1 - SS_res/SS_tot;
|
oo_.sur.R2 = 1 - SS_res/SS_tot;
|
||||||
|
|
||||||
|
@ -167,7 +154,7 @@ oo_.sur.tstat = oo_.sur.beta./oo_.sur.stderr;
|
||||||
|
|
||||||
%% Print Output
|
%% Print Output
|
||||||
if ~options_.noprint
|
if ~options_.noprint
|
||||||
preamble = {sprintf('No. Equations: %d', length(jsonmodel)), ...
|
preamble = {sprintf('No. Equations: %d', neqs), ...
|
||||||
sprintf('No. Independent Variables: %d', size(X, 2)), ...
|
sprintf('No. Independent Variables: %d', size(X, 2)), ...
|
||||||
sprintf('Observations: %d', oo_.sur.dof)};
|
sprintf('Observations: %d', oo_.sur.dof)};
|
||||||
|
|
||||||
|
@ -176,14 +163,12 @@ if ~options_.noprint
|
||||||
sprintf('s^2: %f', oo_.sur.s2), ...
|
sprintf('s^2: %f', oo_.sur.s2), ...
|
||||||
sprintf('Durbin-Watson: %f', oo_.sur.dw)};
|
sprintf('Durbin-Watson: %f', oo_.sur.dw)};
|
||||||
|
|
||||||
if ~isempty(surconstrainedparams)
|
if ~isempty(constrained_param_idxs)
|
||||||
afterward = [afterward, ...
|
afterward = [afterward, ['Constrained parameters: ' constrained_params_str]];
|
||||||
sprintf('Constrained parameters: %s', ...
|
|
||||||
strjoin(pbeta(surconstrainedparams), ', '))];
|
|
||||||
end
|
end
|
||||||
|
|
||||||
dyn_table('SUR Estimation', preamble, afterward, [vars{:}], ...
|
dyn_table('SUR Estimation', preamble, afterward, X.name, ...
|
||||||
{'Coefficients','t-statistic','Std. Error'}, 4, ...
|
{'Estimates','t-statistic','Std. Error'}, 4, ...
|
||||||
[oo_.sur.beta oo_.sur.tstat oo_.sur.stderr]);
|
[oo_.sur.beta oo_.sur.tstat oo_.sur.stderr]);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
Loading…
Reference in New Issue