diff --git a/matlab/ols/sur.m b/matlab/ols/sur.m
index 08783aeda..7b34b0b95 100644
--- a/matlab/ols/sur.m
+++ b/matlab/ols/sur.m
@@ -98,7 +98,7 @@ ast = handle_constant_eqs(get_ast(eqtags));
 neqs = length(ast);
 
 %% 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
 nobs = Y{1}.nobs;
 [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{3} = Y.data;
     varargout{4} = neqs;
+    varargout{5} = lhssub.data;
+    varargout{6} = fp{1};
+    varargout{7} = lp{1};
     return
 end
 
@@ -176,6 +179,9 @@ oo_.sur.(model_name).resid = Y.data - oo_.sur.(model_name).Yhat;
 
 % Correct Yhat reported back to user
 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
 % Estimate for sigma^2
diff --git a/matlab/surgibbs.m b/matlab/surgibbs.m
index a770c5910..91d0c6a88 100644
--- a/matlab/surgibbs.m
+++ b/matlab/surgibbs.m
@@ -1,5 +1,5 @@
-function 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)
+%function ds = surgibbs(ds, param_names, beta0, A, ndraws, discarddraws, thin, eqtags, model_name)
 % Implements Gibbs Samipling for SUR
 %
 % INPUTS
@@ -88,11 +88,13 @@ end
 %  Using a Combination of Direct Monte Carlo and Importance Sampling
 %  Techniques. Bayesian Analysis. 2010. pp 67-70.
 if nargin == 8
-    [nobs, X, Y, m] = sur(ds, param_names, eqtags);
+    [nobs, X, Y, m, lhssub, fp, ~] = sur(ds, param_names, eqtags);
 else
-    [nobs, X, Y, m] = sur(ds, param_names);
+    [nobs, X, Y, m, lhssub, fp, ~] = sur(ds, param_names);
 end
 
+oo_.surgibbs.(model_name).dof = nobs;
+
 beta = beta0;
 A = inv(A);
 thinidx = 1;
@@ -116,7 +118,7 @@ for i = 1:ndraws
     beta = rand_multivariate_normal(betabar', chol(Omegabar), nparams)';
     if i > discarddraws
         if thinidx == thin
-            oo_.surgibbs.(model_name).betadraws(drawidx, :) = beta';
+            oo_.surgibbs.(model_name).betadraws(drawidx, 1:nparams) = beta';
             thinidx = 1;
             drawidx = drawidx + 1;
         else
@@ -125,22 +127,65 @@ for i = 1:ndraws
     end
 end
 
-% save parameter values
-oo_.surgibbs.(model_name).beta = (sum(oo_.surgibbs.(model_name).betadraws)/rows(oo_.surgibbs.(model_name).betadraws))';
+%% Save posterior moments.
+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);
-for i = 1:length(param_names)
-    incidxs(i) = strmatch(param_names{i}, M_.param_names, 'exact');
-    M_.params(incidxs(i)) = oo_.surgibbs.(model_name).beta(i);
+% Yhat
+oo_.surgibbs.(model_name).Yhat = X*oo_.surgibbs.(model_name).posterior.mean.beta;  
+
+% 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
 
+% 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_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
 if ~options_.noprint
-    dyn_table('Gibbs Sampling on SUR', {}, {}, param_names, ...
-        {'Parameter Value'}, 4, oo_.surgibbs.(model_name).beta);
+    ttitle = ['Gibbs Sampling on SUR'];
+    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
 
 %% Plot
@@ -155,7 +200,7 @@ if ~options_.nograph
         subplot(nrows, ncols, j)
         histogram(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')
     end
 end
diff --git a/tests/ecb/SUR/panel_var_diff_NB_simulation_test.mod b/tests/ecb/SUR/panel_var_diff_NB_simulation_test.mod
index f30206691..40995fb59 100644
--- a/tests/ecb/SUR/panel_var_diff_NB_simulation_test.mod
+++ b/tests/ecb/SUR/panel_var_diff_NB_simulation_test.mod
@@ -172,7 +172,7 @@ for i=1:NSIMS
             idxs = [idxs j];
         end
     end
-    sur(simdata{idxs});
+    simdata = sur(simdata{idxs});
     BETA(i, :) = M_.params';
     oo_ = rmfield(oo_, 'sur');
 end
diff --git a/tests/ecb/SUR/panel_var_diff_NB_simulation_zero_eq.mod b/tests/ecb/SUR/panel_var_diff_NB_simulation_zero_eq.mod
index 6cd386e5b..39bc2d146 100644
--- a/tests/ecb/SUR/panel_var_diff_NB_simulation_zero_eq.mod
+++ b/tests/ecb/SUR/panel_var_diff_NB_simulation_zero_eq.mod
@@ -171,7 +171,7 @@ for i=1:NSIMS
             idxs = [idxs j];
         end
     end
-    sur(simdata{idxs});
+    simdata = sur(simdata{idxs});
     BETA(i, :) = M_.params';
     oo_ = rmfield(oo_, 'sur');
 end
diff --git a/tests/ecb/SUR/sur_noniterative.mod b/tests/ecb/SUR/sur_noniterative.mod
index f374bd178..ef29f3ddb 100644
--- a/tests/ecb/SUR/sur_noniterative.mod
+++ b/tests/ecb/SUR/sur_noniterative.mod
@@ -172,7 +172,7 @@ for i=1:NSIMS
             idxs = [idxs j];
         end
     end
-    sur(simdata{idxs}, {}, {}, 'mymodel', true);
+    simdata = sur(simdata{idxs}, {}, {}, 'mymodel', true);
     BETA(i, :) = M_.params';
     oo_ = rmfield(oo_, 'sur');
 end
diff --git a/tests/ecb/SUR/sur_params_noniterative.mod b/tests/ecb/SUR/sur_params_noniterative.mod
index 884fe75f1..b56f65584 100644
--- a/tests/ecb/SUR/sur_params_noniterative.mod
+++ b/tests/ecb/SUR/sur_params_noniterative.mod
@@ -172,7 +172,7 @@ for i=1:NSIMS
             idxs = [idxs j];
         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';
     oo_ = rmfield(oo_, 'sur');
 end
diff --git a/tests/ecb/SURGibbs/fulton_fish.mod b/tests/ecb/SURGibbs/fulton_fish.mod
index d8de5e58c..23e8c8f41 100644
--- a/tests/ecb/SURGibbs/fulton_fish.mod
+++ b/tests/ecb/SURGibbs/fulton_fish.mod
@@ -30,7 +30,7 @@ estparams = {'bq1' 'bq0'};
 estparamsval = [bq1 bq0];
 
 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
    8.552700000000000