2018-01-24 15:22:30 +01:00
|
|
|
function varargout = sur(ds, param_names, eqtags)
|
|
|
|
%function varargout = sur(ds, param_names, eqtags)
|
2017-11-17 14:41:27 +01:00
|
|
|
% Seemingly Unrelated Regressions
|
|
|
|
%
|
|
|
|
% INPUTS
|
2018-01-16 18:42:15 +01:00
|
|
|
% ds [dseries] data to use in estimation
|
2018-01-24 15:22:30 +01:00
|
|
|
% param_names [cellstr] list of parameters to estimate
|
2018-01-16 18:42:15 +01:00
|
|
|
% eqtags [cellstr] names of equation tags to estimate. If empty,
|
|
|
|
% estimate all equations
|
2017-11-17 14:41:27 +01:00
|
|
|
%
|
|
|
|
% OUTPUTS
|
|
|
|
% none
|
|
|
|
%
|
|
|
|
% SPECIAL REQUIREMENTS
|
2017-12-12 12:13:14 +01:00
|
|
|
% dynare must be run with the option: json=compute
|
2017-11-17 14:41:27 +01:00
|
|
|
|
2018-01-05 17:13:46 +01:00
|
|
|
% Copyright (C) 2017-2018 Dynare Team
|
2017-11-17 14:41:27 +01:00
|
|
|
%
|
|
|
|
% 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_ options_
|
|
|
|
|
|
|
|
%% Check input argument
|
2018-01-24 15:22:30 +01:00
|
|
|
assert(nargin >= 1 && nargin <= 3, 'You must provide one, two, or three arguments');
|
2017-11-17 14:41:27 +01:00
|
|
|
assert(~isempty(ds) && isdseries(ds), 'The first argument must be a dseries');
|
2018-01-24 15:22:30 +01:00
|
|
|
if nargin >= 2
|
|
|
|
assert(iscellstr(param_names), 'The 2nd argument must be a cellstr');
|
|
|
|
else
|
|
|
|
param_names = {};
|
|
|
|
end
|
2017-11-17 14:41:27 +01:00
|
|
|
|
|
|
|
%% Read JSON
|
2018-07-10 20:15:53 +02:00
|
|
|
jsonfile = [M_.fname filesep() 'model' filesep() 'json' filesep() 'modfile-original.json'];
|
2017-11-17 14:41:27 +01:00
|
|
|
if exist(jsonfile, 'file') ~= 2
|
2017-12-12 12:13:14 +01:00
|
|
|
error('Could not find %s! Please use the json=compute option (See the Dynare invocation section in the reference manual).', jsonfile);
|
2017-11-17 14:41:27 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
jsonmodel = loadjson(jsonfile);
|
|
|
|
jsonmodel = jsonmodel.model;
|
2018-01-24 15:22:30 +01:00
|
|
|
if nargin == 3
|
2018-01-19 12:51:51 +01:00
|
|
|
jsonmodel = getEquationsByTags(jsonmodel, 'name', eqtags);
|
2018-01-16 18:42:15 +01:00
|
|
|
end
|
2017-11-17 14:41:27 +01:00
|
|
|
|
|
|
|
%% Find parameters and variable names in equations and setup estimation matrices
|
2018-01-24 15:22:30 +01:00
|
|
|
[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, opidxs, surconstrainedparams] = ...
|
2018-01-19 12:51:51 +01:00
|
|
|
pooled_sur_common(ds, jsonmodel);
|
2017-11-17 14:41:27 +01:00
|
|
|
|
2018-01-16 18:42:15 +01:00
|
|
|
if nargin == 1 && size(X, 2) ~= M_.param_nbr
|
2017-12-12 12:13:14 +01:00
|
|
|
warning(['Not all parameters were used in model: ' ...
|
2018-01-11 17:10:12 +01:00
|
|
|
sprintf('%s', strjoin(setdiff(M_.param_names, pbeta), ', '))]);
|
2017-11-17 14:41:27 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
%% Force equations to have the same sample range
|
|
|
|
maxfp = max([startdates{:}]);
|
|
|
|
minlp = min([enddates{:}]);
|
|
|
|
nobs = minlp - maxfp;
|
2018-01-19 12:51:51 +01:00
|
|
|
newY = zeros(nobs*length(jsonmodel), 1);
|
|
|
|
newX = zeros(nobs*length(jsonmodel), columns(X));
|
2017-11-17 14:41:27 +01:00
|
|
|
lastidx = 1;
|
2018-01-19 12:51:51 +01:00
|
|
|
for i = 1:length(jsonmodel)
|
|
|
|
if i == length(jsonmodel)
|
2017-11-17 14:41:27 +01:00
|
|
|
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;
|
2018-01-19 12:51:51 +01:00
|
|
|
if i ~= length(jsonmodel)
|
2017-11-17 14:41:27 +01:00
|
|
|
lastidx = lastidx + nobs + 1;
|
|
|
|
end
|
|
|
|
end
|
2018-01-11 11:18:47 +01:00
|
|
|
|
2018-01-24 15:22:30 +01:00
|
|
|
if ~isempty(param_names)
|
|
|
|
pnamesall = M_.param_names(opidxs);
|
|
|
|
nparams = length(param_names);
|
|
|
|
pidxs = zeros(nparams, 1);
|
|
|
|
for i = 1:nparams
|
|
|
|
idxs = find(strcmp(param_names{i}, pnamesall));
|
|
|
|
if isempty(idxs)
|
|
|
|
if ~isempty(eqtags)
|
|
|
|
error(['Could not find ' param_names{i} ...
|
|
|
|
' in the provided equations specified by ' strjoin(eqtags, ',')]);
|
|
|
|
end
|
|
|
|
error('Unspecified error. Please report');
|
|
|
|
end
|
|
|
|
pidxs(i) = idxs;
|
|
|
|
end
|
|
|
|
vars = [vars{:}];
|
|
|
|
vars = {vars(pidxs)};
|
|
|
|
newY = newY - newX(:, setdiff(1:size(newX, 2), pidxs)) * M_.params(setdiff(opidxs, opidxs(pidxs), 'stable'));
|
|
|
|
newX = newX(:, pidxs);
|
2018-01-29 15:26:28 +01:00
|
|
|
opidxs = opidxs(pidxs);
|
2018-01-24 15:22:30 +01:00
|
|
|
end
|
|
|
|
|
2018-01-11 11:18:47 +01:00
|
|
|
%% Return to surgibbs if called from there
|
|
|
|
st = dbstack(1);
|
|
|
|
if strcmp(st(1).name, 'surgibbs')
|
|
|
|
varargout{1} = length(maxfp:minlp); %dof
|
2018-01-29 15:26:28 +01:00
|
|
|
varargout{2} = opidxs;
|
2018-01-11 11:18:47 +01:00
|
|
|
varargout{3} = newX;
|
|
|
|
varargout{4} = newY;
|
2018-01-19 12:51:51 +01:00
|
|
|
varargout{5} = length(jsonmodel);
|
2018-01-11 11:18:47 +01:00
|
|
|
return
|
|
|
|
end
|
|
|
|
|
2017-11-17 14:41:27 +01:00
|
|
|
Y = newY;
|
|
|
|
X = newX;
|
2018-01-05 17:13:46 +01:00
|
|
|
oo_.sur.dof = length(maxfp:minlp);
|
2017-11-17 14:41:27 +01:00
|
|
|
|
|
|
|
%% Estimation
|
|
|
|
% Estimated Parameters
|
|
|
|
[q, r] = qr(X, 0);
|
2017-12-12 13:07:25 +01:00
|
|
|
xpxi = (r'*r)\eye(size(X, 2));
|
2017-11-17 14:41:27 +01:00
|
|
|
resid = Y - X * (r\(q'*Y));
|
2018-01-19 12:51:51 +01:00
|
|
|
resid = reshape(resid, oo_.sur.dof, length(jsonmodel));
|
2017-11-17 14:41:27 +01:00
|
|
|
|
|
|
|
M_.Sigma_e = resid'*resid/oo_.sur.dof;
|
2017-12-18 15:01:09 +01:00
|
|
|
kLeye = kron(chol(inv(M_.Sigma_e)), eye(oo_.sur.dof));
|
2017-11-17 14:41:27 +01:00
|
|
|
[q, r] = qr(kLeye*X, 0);
|
|
|
|
oo_.sur.beta = r\(q'*kLeye*Y);
|
|
|
|
|
2018-01-29 15:26:28 +01:00
|
|
|
M_.params(opidxs) = oo_.sur.beta;
|
2018-01-05 17:13:46 +01:00
|
|
|
|
|
|
|
% Yhat
|
|
|
|
oo_.sur.Yhat = X * oo_.sur.beta;
|
|
|
|
|
|
|
|
% Residuals
|
|
|
|
oo_.sur.resid = Y - oo_.sur.Yhat;
|
|
|
|
|
2017-11-17 14:41:27 +01:00
|
|
|
%% Calculate statistics
|
|
|
|
% Estimate for sigma^2
|
|
|
|
SS_res = oo_.sur.resid'*oo_.sur.resid;
|
|
|
|
oo_.sur.s2 = SS_res/oo_.sur.dof;
|
|
|
|
|
|
|
|
% R^2
|
|
|
|
ym = Y - mean(Y);
|
|
|
|
SS_tot = ym'*ym;
|
|
|
|
oo_.sur.R2 = 1 - SS_res/SS_tot;
|
|
|
|
|
|
|
|
% Adjusted R^2
|
|
|
|
oo_.sur.adjR2 = oo_.sur.R2 - (1 - oo_.sur.R2)*M_.param_nbr/(oo_.sur.dof - 1);
|
|
|
|
|
|
|
|
% Durbin-Watson
|
|
|
|
ediff = oo_.sur.resid(2:oo_.sur.dof) - oo_.sur.resid(1:oo_.sur.dof - 1);
|
|
|
|
oo_.sur.dw = (ediff'*ediff)/SS_res;
|
|
|
|
|
|
|
|
% Standard Error
|
|
|
|
oo_.sur.stderr = sqrt(oo_.sur.s2*diag(xpxi));
|
|
|
|
|
|
|
|
% T-Stat
|
|
|
|
oo_.sur.tstat = oo_.sur.beta./oo_.sur.stderr;
|
|
|
|
|
|
|
|
%% Print Output
|
|
|
|
if ~options_.noprint
|
2018-01-29 15:36:56 +01:00
|
|
|
preamble = {sprintf('No. Equations: %d', length(jsonmodel)), ...
|
|
|
|
sprintf('No. Independent Variables: %d', size(X, 2)), ...
|
2017-11-17 14:41:27 +01:00
|
|
|
sprintf('Observations: %d', oo_.sur.dof)};
|
|
|
|
|
|
|
|
afterward = {sprintf('R^2: %f', oo_.sur.R2), ...
|
|
|
|
sprintf('R^2 Adjusted: %f', oo_.sur.adjR2), ...
|
|
|
|
sprintf('s^2: %f', oo_.sur.s2), ...
|
|
|
|
sprintf('Durbin-Watson: %f', oo_.sur.dw)};
|
|
|
|
|
2017-12-12 13:07:25 +01:00
|
|
|
if ~isempty(surconstrainedparams)
|
|
|
|
afterward = [afterward, ...
|
|
|
|
sprintf('Constrained parameters: %s', ...
|
|
|
|
strjoin(pbeta(surconstrainedparams), ', '))];
|
|
|
|
end
|
|
|
|
|
2017-11-17 14:41:27 +01:00
|
|
|
dyn_table('SUR Estimation', preamble, afterward, [vars{:}], ...
|
|
|
|
{'Coefficients','t-statistic','Std. Error'}, 4, ...
|
|
|
|
[oo_.sur.beta oo_.sur.tstat oo_.sur.stderr]);
|
|
|
|
end
|
|
|
|
end
|