function [beta, std, res, r2] = ols(y,x,const)
% usage: [beta, std, res, r2] = ols(y,x,const)
% ols regression procedure
%   y = n*1 vector, dependent
%   x = n*m matrix, independent
%   const = 1: add constant, 0: no constant
% returns
%   beta = estimated linear params m*1 (or (m+1)*1 with constant)
%   std  = standard error for beta
%   res  = estimated residuals n*1
%   r2   = r-squared


% augment x if needed
if(const==1)
	xi = [ones(size(x,1),1), x];
else
	xi = x;
end

% store x'x matrix
q = xi'*xi;
% now invert it
qi = inv(q);
% OLS formula for beta
beta = qi*xi'*y;
% now get the estimated residuals
res = y-xi*beta;
% sum of squared errors  (size(x,1) = length of sample)
sq = sum(res.^2)/size(x,1);
% std errors for beta 
std = sqrt(sq * diag(qi));
% r-squared
r2 = 1-sum(res.^2)/sum((y-mean(y)).^2);