Update sur and surgibbs function outputs.
parent
7b1c61f63c
commit
75a929051f
|
@ -98,7 +98,7 @@ ast = handle_constant_eqs(get_ast(eqtags));
|
||||||
neqs = length(ast);
|
neqs = length(ast);
|
||||||
|
|
||||||
%% Find parameters and variable names in equations and setup estimation matrices
|
%% Find parameters and variable names in equations and setup estimation matrices
|
||||||
[Y, lhssub, X, ~, ~, residnames] = common_parsing(ds(ds_range), ast, true, param_names);
|
[Y, lhssub, X, fp, lp, residnames] = common_parsing(ds(ds_range), ast, true, param_names);
|
||||||
clear ast
|
clear ast
|
||||||
nobs = Y{1}.nobs;
|
nobs = Y{1}.nobs;
|
||||||
[Y, lhssub, X, constrained] = put_in_sur_form(Y, lhssub, X);
|
[Y, lhssub, X, constrained] = put_in_sur_form(Y, lhssub, X);
|
||||||
|
@ -114,6 +114,9 @@ if ~isempty(st) && strcmp(st(1).name, 'surgibbs')
|
||||||
varargout{2} = X{param_names{:}}.data;
|
varargout{2} = X{param_names{:}}.data;
|
||||||
varargout{3} = Y.data;
|
varargout{3} = Y.data;
|
||||||
varargout{4} = neqs;
|
varargout{4} = neqs;
|
||||||
|
varargout{5} = lhssub.data;
|
||||||
|
varargout{6} = fp{1};
|
||||||
|
varargout{7} = lp{1};
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -176,6 +179,9 @@ oo_.sur.(model_name).resid = Y.data - oo_.sur.(model_name).Yhat;
|
||||||
|
|
||||||
% Correct Yhat reported back to user
|
% Correct Yhat reported back to user
|
||||||
oo_.sur.(model_name).Yhat = oo_.sur.(model_name).Yhat + lhssub;
|
oo_.sur.(model_name).Yhat = oo_.sur.(model_name).Yhat + lhssub;
|
||||||
|
yhatname = [model_name '_FIT'];
|
||||||
|
ds.(yhatname) = dseries(oo_.sur.(model_name).Yhat.data, fp{1}, yhatname);
|
||||||
|
varargout{1} = ds;
|
||||||
|
|
||||||
%% Calculate statistics
|
%% Calculate statistics
|
||||||
% Estimate for sigma^2
|
% Estimate for sigma^2
|
||||||
|
|
|
@ -1,5 +1,5 @@
|
||||||
function surgibbs(ds, param_names, beta0, A, ndraws, discarddraws, thin, eqtags, model_name)
|
function ds = surgibbs(ds, param_names, beta0, A, ndraws, discarddraws, thin, eqtags, model_name)
|
||||||
%function surgibbs(ds, param_names, beta0, A, ndraws, discarddraws, thin, eqtags, model_name)
|
%function ds = surgibbs(ds, param_names, beta0, A, ndraws, discarddraws, thin, eqtags, model_name)
|
||||||
% Implements Gibbs Samipling for SUR
|
% Implements Gibbs Samipling for SUR
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
|
@ -88,11 +88,13 @@ end
|
||||||
% Using a Combination of Direct Monte Carlo and Importance Sampling
|
% Using a Combination of Direct Monte Carlo and Importance Sampling
|
||||||
% Techniques. Bayesian Analysis. 2010. pp 67-70.
|
% Techniques. Bayesian Analysis. 2010. pp 67-70.
|
||||||
if nargin == 8
|
if nargin == 8
|
||||||
[nobs, X, Y, m] = sur(ds, param_names, eqtags);
|
[nobs, X, Y, m, lhssub, fp, ~] = sur(ds, param_names, eqtags);
|
||||||
else
|
else
|
||||||
[nobs, X, Y, m] = sur(ds, param_names);
|
[nobs, X, Y, m, lhssub, fp, ~] = sur(ds, param_names);
|
||||||
end
|
end
|
||||||
|
|
||||||
|
oo_.surgibbs.(model_name).dof = nobs;
|
||||||
|
|
||||||
beta = beta0;
|
beta = beta0;
|
||||||
A = inv(A);
|
A = inv(A);
|
||||||
thinidx = 1;
|
thinidx = 1;
|
||||||
|
@ -116,7 +118,7 @@ for i = 1:ndraws
|
||||||
beta = rand_multivariate_normal(betabar', chol(Omegabar), nparams)';
|
beta = rand_multivariate_normal(betabar', chol(Omegabar), nparams)';
|
||||||
if i > discarddraws
|
if i > discarddraws
|
||||||
if thinidx == thin
|
if thinidx == thin
|
||||||
oo_.surgibbs.(model_name).betadraws(drawidx, :) = beta';
|
oo_.surgibbs.(model_name).betadraws(drawidx, 1:nparams) = beta';
|
||||||
thinidx = 1;
|
thinidx = 1;
|
||||||
drawidx = drawidx + 1;
|
drawidx = drawidx + 1;
|
||||||
else
|
else
|
||||||
|
@ -125,22 +127,65 @@ for i = 1:ndraws
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
% save parameter values
|
%% Save posterior moments.
|
||||||
oo_.surgibbs.(model_name).beta = (sum(oo_.surgibbs.(model_name).betadraws)/rows(oo_.surgibbs.(model_name).betadraws))';
|
oo_.surgibbs.(model_name).posterior.mean.beta = (sum(oo_.surgibbs.(model_name).betadraws)/rows(oo_.surgibbs.(model_name).betadraws))';
|
||||||
|
oo_.surgibbs.(model_name).posterior.variance.beta = cov(oo_.surgibbs.(model_name).betadraws);
|
||||||
|
|
||||||
incidxs = zeros(length(param_names), 1);
|
% Yhat
|
||||||
for i = 1:length(param_names)
|
oo_.surgibbs.(model_name).Yhat = X*oo_.surgibbs.(model_name).posterior.mean.beta;
|
||||||
incidxs(i) = strmatch(param_names{i}, M_.param_names, 'exact');
|
|
||||||
M_.params(incidxs(i)) = oo_.surgibbs.(model_name).beta(i);
|
% Residuals
|
||||||
|
oo_.surgibbs.(model_name).resid = Y - oo_.surgibbs.(model_name).Yhat;
|
||||||
|
|
||||||
|
% Correct Yhat reported back to user
|
||||||
|
oo_.surgibbs.(model_name).Yhat = oo_.surgibbs.(model_name).Yhat + lhssub;
|
||||||
|
yhatname = [model_name '_FIT'];
|
||||||
|
ds.(yhatname) = dseries(oo_.surgibbs.(model_name).Yhat, fp, yhatname);
|
||||||
|
|
||||||
|
% Compute and save posterior densities.
|
||||||
|
for i=1:nparams
|
||||||
|
xx = oo_.surgibbs.(model_name).betadraws(:,i);
|
||||||
|
nn = length(xx);
|
||||||
|
bandwidth = mh_optimal_bandwidth(xx, nn, 0, 'gaussian');
|
||||||
|
[x, f] = kernel_density_estimate(xx, 512, nn, bandwidth, 'gaussian');
|
||||||
|
oo_.surgibbs.(model_name).posterior.density.(param_names{i}) = [x, f];
|
||||||
end
|
end
|
||||||
|
|
||||||
|
% Update model1s parameters with posterior mean.
|
||||||
|
oo_.surgibbs.(model_name).param_idxs = zeros(length(param_names), 1);
|
||||||
|
for i = 1:length(param_names)
|
||||||
|
if ~strcmp(param_names{i}, 'intercept')
|
||||||
|
oo_.surgibbs.(model_name).param_idxs(i) = find(strcmp(M_.param_names, param_names{i}));
|
||||||
|
M_.params(oo_.surgibbs.(model_name).param_idxs(i)) = oo_.surgibbs.(model_name).posterior.mean.beta(i);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
oo_.surgibbs.(model_name).pnames = param_names;
|
||||||
|
oo_.surgibbs.(model_name).neqs = m;
|
||||||
|
|
||||||
|
% Estimate for sigma^2
|
||||||
|
SS_res = oo_.surgibbs.(model_name).resid'*oo_.surgibbs.(model_name).resid;
|
||||||
|
oo_.surgibbs.(model_name).s2 = SS_res/oo_.surgibbs.(model_name).dof;
|
||||||
|
|
||||||
|
% R^2
|
||||||
|
ym = Y - mean(Y);
|
||||||
|
SS_tot = ym'*ym;
|
||||||
|
oo_.surgibbs.(model_name).R2 = 1 - SS_res/SS_tot;
|
||||||
|
|
||||||
% Write .inc file
|
% Write .inc file
|
||||||
write_param_init_inc_file('surgibbs', model_name, incidxs, oo_.surgibbs.(model_name).beta);
|
write_param_init_inc_file('surgibbs', model_name, oo_.surgibbs.(model_name).param_idxs, oo_.surgibbs.(model_name).posterior.mean.beta);
|
||||||
|
|
||||||
%% Print Output
|
%% Print Output
|
||||||
if ~options_.noprint
|
if ~options_.noprint
|
||||||
dyn_table('Gibbs Sampling on SUR', {}, {}, param_names, ...
|
ttitle = ['Gibbs Sampling on SUR'];
|
||||||
{'Parameter Value'}, 4, oo_.surgibbs.(model_name).beta);
|
preamble = {['Model name: ' model_name], ...
|
||||||
|
sprintf('No. Equations: %d', oo_.surgibbs.(model_name).neqs), ...
|
||||||
|
sprintf('No. Independent Variables: %d', size(X, 2)), ...
|
||||||
|
sprintf('Observations: %d', oo_.surgibbs.(model_name).dof)};
|
||||||
|
|
||||||
|
afterward = {sprintf('s^2: %f', oo_.surgibbs.(model_name).s2), sprintf('R^2: %f', oo_.surgibbs.(model_name).R2)};
|
||||||
|
dyn_table(ttitle, preamble, afterward, param_names,...
|
||||||
|
{'Posterior mean', 'Posterior std.'}, 4,...
|
||||||
|
[oo_.surgibbs.(model_name).posterior.mean.beta, sqrt(diag(oo_.surgibbs.(model_name).posterior.variance.beta))]);
|
||||||
end
|
end
|
||||||
|
|
||||||
%% Plot
|
%% Plot
|
||||||
|
@ -155,7 +200,7 @@ if ~options_.nograph
|
||||||
subplot(nrows, ncols, j)
|
subplot(nrows, ncols, j)
|
||||||
histogram(oo_.surgibbs.(model_name).betadraws(:, j))
|
histogram(oo_.surgibbs.(model_name).betadraws(:, j))
|
||||||
hc = histcounts(oo_.surgibbs.(model_name).betadraws(:, j));
|
hc = histcounts(oo_.surgibbs.(model_name).betadraws(:, j));
|
||||||
line([oo_.surgibbs.(model_name).beta(j) oo_.surgibbs.(model_name).beta(j)], [min(hc) max(hc)], 'Color', 'red');
|
line([oo_.surgibbs.(model_name).posterior.mean.beta(j) oo_.surgibbs.(model_name).posterior.mean.beta(j)], [min(hc) max(hc)], 'Color', 'red');
|
||||||
title(param_names{j}, 'Interpreter', 'none')
|
title(param_names{j}, 'Interpreter', 'none')
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
|
@ -172,7 +172,7 @@ for i=1:NSIMS
|
||||||
idxs = [idxs j];
|
idxs = [idxs j];
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
sur(simdata{idxs});
|
simdata = sur(simdata{idxs});
|
||||||
BETA(i, :) = M_.params';
|
BETA(i, :) = M_.params';
|
||||||
oo_ = rmfield(oo_, 'sur');
|
oo_ = rmfield(oo_, 'sur');
|
||||||
end
|
end
|
||||||
|
|
|
@ -171,7 +171,7 @@ for i=1:NSIMS
|
||||||
idxs = [idxs j];
|
idxs = [idxs j];
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
sur(simdata{idxs});
|
simdata = sur(simdata{idxs});
|
||||||
BETA(i, :) = M_.params';
|
BETA(i, :) = M_.params';
|
||||||
oo_ = rmfield(oo_, 'sur');
|
oo_ = rmfield(oo_, 'sur');
|
||||||
end
|
end
|
||||||
|
|
|
@ -172,7 +172,7 @@ for i=1:NSIMS
|
||||||
idxs = [idxs j];
|
idxs = [idxs j];
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
sur(simdata{idxs}, {}, {}, 'mymodel', true);
|
simdata = sur(simdata{idxs}, {}, {}, 'mymodel', true);
|
||||||
BETA(i, :) = M_.params';
|
BETA(i, :) = M_.params';
|
||||||
oo_ = rmfield(oo_, 'sur');
|
oo_ = rmfield(oo_, 'sur');
|
||||||
end
|
end
|
||||||
|
|
|
@ -172,7 +172,7 @@ for i=1:NSIMS
|
||||||
idxs = [idxs j];
|
idxs = [idxs j];
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
sur(simdata{idxs}, {'u2_q_yed_u2_g_yer_L1', 'u2_estn_u2_estn_L1', 'u2_ehic_u2_ehic_L1', 'de_q_yed_de_g_yer_L1', 'u2_g_yer_u2_g_yer_L1', 'u2_stn_u2_q_yed_L1', 'u2_q_yed_ecm_u2_q_yed_L1', 'de_g_yer_u2_stn_L1'}, {}, 'mymodel', true);
|
simdata = sur(simdata{idxs}, {'u2_q_yed_u2_g_yer_L1', 'u2_estn_u2_estn_L1', 'u2_ehic_u2_ehic_L1', 'de_q_yed_de_g_yer_L1', 'u2_g_yer_u2_g_yer_L1', 'u2_stn_u2_q_yed_L1', 'u2_q_yed_ecm_u2_q_yed_L1', 'de_g_yer_u2_stn_L1'}, {}, 'mymodel', true);
|
||||||
BETA(i, :) = M_.params';
|
BETA(i, :) = M_.params';
|
||||||
oo_ = rmfield(oo_, 'sur');
|
oo_ = rmfield(oo_, 'sur');
|
||||||
end
|
end
|
||||||
|
|
|
@ -30,7 +30,7 @@ estparams = {'bq1' 'bq0'};
|
||||||
estparamsval = [bq1 bq0];
|
estparamsval = [bq1 bq0];
|
||||||
|
|
||||||
A = 0.0005.*eye(length(estparams));
|
A = 0.0005.*eye(length(estparams));
|
||||||
surgibbs(dseries('fishdata.csv'), estparams, estparamsval, A, 20000, 5000, 7);
|
simdata = surgibbs(dseries('fishdata.csv'), estparams, estparamsval, A, 20000, 5000, 7);
|
||||||
|
|
||||||
good = [6.791587808530124
|
good = [6.791587808530124
|
||||||
8.552700000000000
|
8.552700000000000
|
||||||
|
|
Loading…
Reference in New Issue