% start with a vector y of dependent variable observations, and a matrix x of independent variable observations
clearvars -except x y

% establish the number of observations and independent variables
[ n K ] = size(x);

% create the K by 1 vector b, which estimates the beta values
b = inv( x' * x ) * x' * y;  

% using b, generate the predicted y values for each observation
yhat = x * b;

% e is the vector of residuals for each observation
e = y - yhat;

% SSR is the sum of squared residuals
e2 = zeros(n,1);
for i = 1:n
    e2(i) = e(i) ^ 2;
end
SSR = sum(e2);

% s2 is the estimated variance of the error term
s2 = SSR / ( n - K );

% varbk is a K by 1 vector that gives the estimated variance of each b value
varb = s2 * inv( x' * x ); 
varbk = zeros(K,1);
for k = 1:K
    varbk(k) = varb(k,k);
end

% SE is the standard error, i.e. the square root of varbk
SE = zeros(K,1);
for k = 1:K
    SE(k) = sqrt( varbk(k) );
end

% the t-statistic for each b(k) is the point estimate divided by the standard error
t = zeros(K,1);
for k = 1:K
    t(k) = b(k) / SE(k);
end

% p values generated using the assumption that the b values are normally distributed
% note that the normcdf( ) function here relies on the Statistics Toolbox; if you don't have that, you should delete this section of the code
pnorm = zeros(K,1);
for k = 1:K
    pnorm(k) = 2 * normcdf( -abs( t(k) ) );
end

% p values generated using the assumption that the b values are distributed according to a Student's t distribution. (this is much more standard, but rarely very different from the normal approximation)
% note that the tcdf( ) function here relies on the Statistics Toolbox; if you don't have that, you should delete this section of the code
pstud = zeros(K,1);
for k = 1:K
    pstud(k) = 2 * tcdf( -abs( t(k) ) , n-K );
end

% calculate R squared
ybar = mean(y);
yminusybar2 = zeros(n,1);   yhatminusybar2 = zeros(n,1);
for i = 1:n
    yminusybar2(i) = ( y(i) - ybar )^2;
    yhatminusybar2(i) = ( yhat(i) - ybar )^2;
end
R2 = sum(yhatminusybar2) / sum(yminusybar2);

display(b)
%#ok<*MINV>