function [beta,se,pv,S,ssr,R2,R2a,F,pvF] = myols(Y,X,hr)
% [beta,se,pv,S,ssr,R2,R2a,F,pvF] = myols(Y,X,hr)
%   Estimate a linear multiple regression of the following form,
%       y = beta_0 + beta_1 x_1 + ... + beta_k x_k + u
%   Inputs: 
%       Y:  an n-by-1 vector of dependent variable 
%       X:  an n-by-m matrix of independent variables. The first column is
%               a vector of 1s. 
%       hr: a logical variable, if hr==1 then heteroscedasticity-robust
%               covariance/standard errors are calculated. (default: hr==1)

if nargin < 3 || isempty(hr)
    hr = 0;
end

[n,m] = size(X);
beta = (X'*X)\X'*Y;
e = Y - X*beta;
ssr = e'*e;
s2 = ssr/(n - m);

if hr==1
    S = inv(X'*X) * (X'*(X.*repmat(e.^2,1,m))) * inv(X'*X);
else
    S = inv(X'*X)*s2;
end
se = sqrt(diag(S));

Y1 = Y-mean(Y);
sst = Y1'*Y1;
R2 = 1 - ssr / sst;
R2a = 1 - (ssr/(n-m)) / (sst/(n-1));

pv = 2*(1-tcdf(abs(beta./se),n-m));
F = ((sst-ssr)/(m-1)) / (ssr/(n-m));
pvF = 1-fcdf(F,m-1,n-m);

if nargout == 0
    disp('Estimation Results:');
    disp(['N = ' num2str(n)]);
    disp('    Parameter    S.E.     P-value');
    disp([beta se pv]);
    disp(['R2 = ' num2str(R2)]);
    disp(['Adjusted R2 = ' num2str(R2a)]);
    disp(['SSR = ' num2str(ssr)]);
    disp(['F statistic of the model (P-value): ' num2str(F) ' (' num2str(pvF) ')']);
end
    



