diff --git a/matlab/ols/sur.m b/matlab/ols/sur.m
index 2d9a466b6..6e45de722 100644
--- a/matlab/ols/sur.m
+++ b/matlab/ols/sur.m
@@ -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
diff --git a/matlab/surgibbs.m b/matlab/surgibbs.m
index c7b8398cb..684808ce5 100644
--- a/matlab/surgibbs.m
+++ b/matlab/surgibbs.m
@@ -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
diff --git a/tests/ecb/SURGibbs/fulton_fish.mod b/tests/ecb/SURGibbs/fulton_fish.mod
index 0466d5faa..c65fbbc84 100644
--- a/tests/ecb/SURGibbs/fulton_fish.mod
+++ b/tests/ecb/SURGibbs/fulton_fish.mod
@@ -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);