function [FuncVal varargout] = OLS_b(par,x,y,varargin)
% [FuncVal Gr] = OLS_b(par,x,y)
% par: Beta's
% x: n-by-p matrix
% y: Observed data

n = length(y);
if nargin > 3 && varargin{1} == 1
    x = [ones(n,1) x];
end
eps = 0.0000001;
beta = par;
beta = beta(:);
e = y - x*beta;
FuncVal = sum(e.^2);

if nargout > 1
    Gr = [];
    for i = 1:length(beta)
        beta_a = beta;
        beta_a(i) = beta(i) + eps;
        Funk(i,1) = grad_OLS(beta_a,x,y);
        Gr(i,1) = (Funk(i) - FuncVal)/eps;
    end
    varargout{1} = Gr;
end


function Funk = grad_OLS(par,x,y)
beta = par;
beta = beta(:);
e = y - x*beta;
Funk = sum(e.^2);