function [beta sigma2] = bmareg(y,x,bhatprior,g)
% Estimates a Bayesian regression model with a g-prior. Formulas are based
% on Koop (2003), p. 271 and p. 37
% INPUTS: y - dependent variable, of size Tx1
%         x - independent variable, of size Txk
%         bhatprior - prior on the regression coefficients, of size kx1,
%         g - the value of g in the g prior on the regression coefficients, scalar
%         the prior on var-covar is assumed to be uninformative
% OUTPUTS: beta - posterior (conditional) mean of the regression coefficients, of size kx1
%          sigma2 - posterior shape parameter on the Gamma distribution for
%          the var-covar
[n,k] = size(x);

bhat = pinv(x'*x)*x'*y;
dy = y - (x*bhatprior);
px = eye(n) - x*pinv(x'*x)*x';

beta =  bhat/(1+g) + bhatprior*g/(1+g);
sigma2 = (1/(1+g)*y'*px*y + (g/(1+g))*(dy)'*dy)/n;