surgibbs: pass prior as argument, allow limiting the number of parameters estimated, check inputs, change variable name
parent
b1d8e7aafa
commit
049be61ef5
|
@ -75,6 +75,18 @@ for i = 1:length(lhs)
|
|||
lastidx = lastidx + nobs + 1;
|
||||
end
|
||||
end
|
||||
|
||||
%% Return to surgibbs if called from there
|
||||
st = dbstack(1);
|
||||
if strcmp(st(1).name, 'surgibbs')
|
||||
varargout{1} = length(maxfp:minlp); %dof
|
||||
varargout{2} = pidxs;
|
||||
varargout{3} = newX;
|
||||
varargout{4} = newY;
|
||||
varargout{5} = length(lhs);
|
||||
return
|
||||
end
|
||||
|
||||
Y = newY;
|
||||
X = newX;
|
||||
oo_.sur.dof = length(maxfp:minlp);
|
||||
|
@ -91,19 +103,6 @@ kLeye = kron(chol(inv(M_.Sigma_e)), eye(oo_.sur.dof));
|
|||
[q, r] = qr(kLeye*X, 0);
|
||||
oo_.sur.beta = r\(q'*kLeye*Y);
|
||||
|
||||
%% Return to surgibbs if called from there
|
||||
st = dbstack(1);
|
||||
if strcmp(st(1).name, 'surgibbs')
|
||||
varargout{1} = oo_.sur.dof;
|
||||
varargout{2} = size(X, 2);
|
||||
varargout{3} = pidxs;
|
||||
varargout{4} = oo_.sur.beta;
|
||||
varargout{5} = X;
|
||||
varargout{6} = Y;
|
||||
varargout{7} = length(lhs);
|
||||
return
|
||||
end
|
||||
|
||||
M_.params(pidxs, 1) = oo_.sur.beta;
|
||||
|
||||
% Yhat
|
||||
|
|
|
@ -1,10 +1,13 @@
|
|||
function surgibbs(ds, A, ndraws, discarddraws)
|
||||
%function surgibbs(ds, ndraws, discarddraws)
|
||||
function surgibbs(ds, param_names, beta0, A, ndraws, discarddraws)
|
||||
%function surgibbs(ds, param_names, beta0, A, ndraws, discarddraws)
|
||||
% Implements Gibbs Samipling for SUR
|
||||
%
|
||||
% INPUTS
|
||||
% ds [dseries] data
|
||||
% A [matrix] prior distribution variance
|
||||
% param_names [cellstr] list of parameters to estimate
|
||||
% beta0 [vector] prior values (in order of param_names)
|
||||
% A [matrix] prior distribution variance (in order of
|
||||
% param_names)
|
||||
% ndraws [int] number of draws
|
||||
% discarddraws [int] number of draws to discard
|
||||
%
|
||||
|
@ -40,22 +43,35 @@ function surgibbs(ds, A, ndraws, discarddraws)
|
|||
global M_ oo_
|
||||
|
||||
%% Check input
|
||||
assert(nargin == 3 || nargin == 4, 'Incorrect number of arguments passed to surgibbs');
|
||||
if nargin == 3
|
||||
assert(nargin == 5 || nargin == 6, 'Incorrect number of arguments passed to surgibbs');
|
||||
assert(isdseries(ds), 'The 1st argument must be a dseries');
|
||||
assert(iscellstr(param_names), 'The 2nd argument must be a cellstr');
|
||||
assert(isvector(beta0) && length(beta0) == length(param_names), ...
|
||||
'The 3rd argument must be a vector with the same length as param_names and the same ');
|
||||
if isrow(beta0)
|
||||
beta0 = beta0';
|
||||
end
|
||||
assert(ismatrix(A) && all(all((A == A'))) && length(beta0) == size(A, 2), ...
|
||||
'The 4th argument must be a symmetric matrix with the same dimension as beta0');
|
||||
assert(isint(ndraws), 'The 5th argument must be an integer');
|
||||
if nargin == 5
|
||||
discarddraws = 0;
|
||||
else
|
||||
assert(isint(discarddraws), 'The 6th argument, if provided, must be an integer');
|
||||
end
|
||||
|
||||
%% Estimation
|
||||
beta0 = M_.params;
|
||||
[nobs, nvars, pidxs, beta, X, Y, m] = sur(ds);
|
||||
[nobs, pidxs, X, Y, m] = sur(ds);
|
||||
pnamesall = cellstr(M_.param_names(pidxs, :));
|
||||
if any(isnan(beta0))
|
||||
beta0 = beta;
|
||||
else
|
||||
beta = beta0;
|
||||
nparams = length(param_names);
|
||||
pidxs = zeros(nparams, 1);
|
||||
for i = 1:nparams
|
||||
pidxs(i) = find(strcmp(param_names{i}, pnamesall));
|
||||
end
|
||||
X = X(:, pidxs);
|
||||
beta = beta0;
|
||||
A = inv(A);
|
||||
oo_.surgibbs.betadraws = zeros(ndraws-discarddraws, nvars);
|
||||
oo_.surgibbs.betadraws = zeros(ndraws-discarddraws, nparams);
|
||||
for i = 1:ndraws
|
||||
% Draw Omega, given X, Y, Beta
|
||||
resid = reshape(Y - X*beta, nobs, m);
|
||||
|
@ -67,7 +83,7 @@ for i = 1:ndraws
|
|||
Omegabar = inv(tmp1 + A);
|
||||
betahat = tmp1\X'*tmp*Y;
|
||||
betabar = Omegabar*(tmp1*betahat+A*beta0);
|
||||
beta = rand_multivariate_normal(betabar', chol(Omegabar), nvars)';
|
||||
beta = rand_multivariate_normal(betabar', chol(Omegabar), nparams)';
|
||||
if i > discarddraws
|
||||
oo_.surgibbs.betadraws(i-discarddraws, :) = beta';
|
||||
end
|
||||
|
@ -78,21 +94,21 @@ oo_.surgibbs.beta = (sum(oo_.surgibbs.betadraws)/(ndraws-discarddraws))';
|
|||
M_.params(pidxs, 1) = oo_.surgibbs.beta;
|
||||
|
||||
%% Print Output
|
||||
dyn_table('Gibbs Sampling on SUR', {}, {}, pnamesall, ...
|
||||
dyn_table('Gibbs Sampling on SUR', {}, {}, param_names, ...
|
||||
{'Parameter Value'}, 4, oo_.surgibbs.beta);
|
||||
|
||||
%% Plot
|
||||
figure
|
||||
nrows = 5;
|
||||
ncols = floor(nvars/nrows);
|
||||
if mod(nvars, nrows) ~= 0
|
||||
ncols = floor(nparams/nrows);
|
||||
if mod(nparams, nrows) ~= 0
|
||||
ncols = ncols + 1;
|
||||
end
|
||||
for j = 1:length(pnamesall)
|
||||
M_.params(strmatch(pnamesall{j}, M_.param_names, 'exact')) = oo_.surgibbs.beta(j);
|
||||
for j = 1:length(param_names)
|
||||
M_.params(strmatch(param_names{j}, M_.param_names, 'exact')) = oo_.surgibbs.beta(j);
|
||||
subplot(nrows, ncols, j)
|
||||
histogram(oo_.surgibbs.betadraws(:, j))
|
||||
hc = histcounts(oo_.surgibbs.betadraws(:, j));
|
||||
line([oo_.surgibbs.beta(j) oo_.surgibbs.beta(j)], [min(hc) max(hc)], 'Color', 'red');
|
||||
title(pnamesall{j}, 'Interpreter', 'none')
|
||||
title(param_names{j}, 'Interpreter', 'none')
|
||||
end
|
||||
|
|
|
@ -17,8 +17,17 @@ bq0 = 6.7523;
|
|||
bq1 = -0.7969;
|
||||
|
||||
model(linear);
|
||||
price = bp0 + bp1*stormy + bp2*mixed + res_v;
|
||||
qty = bq0 + bq1*price + res_u;
|
||||
price = bp0 + bp1*stormy + bp2*mixed + res_v;
|
||||
end;
|
||||
|
||||
surgibbs(dseries('fishdata.csv'), 0.0005.*eye(M_.param_nbr), 10000, 5000);
|
||||
% Estimate all parameters
|
||||
%estparams = cellstr(M_.param_names);
|
||||
%estparamsval = M_.params;
|
||||
|
||||
% Estimate demand parameters
|
||||
estparams = {'bq1' 'bq0'};
|
||||
estparamsval = [bq1 bq0];
|
||||
|
||||
A = 0.0005.*eye(length(estparams));
|
||||
surgibbs(dseries('fishdata.csv'), estparams, estparamsval, A, 10000, 5000);
|
||||
|
|
Loading…
Reference in New Issue