%% Regression Analysis with Matlab
% Michael Bar, San Francisco State University
% Table of Contents
% Introduction
% Statistics and Machine Learning Toolbox
%     Reading the data from the web
%     Estimating the model
% Matrix Algebra
%     OLS problem in matrix notation
%     Solving the OLS problem with matrix algebra
%% Introduction
% In this note we demonstrate the use of matrices for regression analysis. The multivariate regression model is:
% Y_i =\beta {\;}_1 +\beta_2 X_{2i} +\ldotp \ldotp \ldotp +\beta_k X_{\textrm{ki}} +u_i ,\;\;\;i=1,\ldotp \ldotp \ldotp ,n
% Thus, the outcome of observation i, Y_i, is a linear function of regressors \left(X_2 ,\ldotp \ldotp \ldotp ,X_k \right), plus an error term u_i. Let \left(b_1 ,\ldotp \ldotp \ldotp ,b_k \right) be some estimators of the unknown coefficients \left(\beta_1 ,\ldotp \ldotp \ldotp ,\beta_k \right). We define the fitted eqaution, which gives the predicted values of the outcome for given estimates:
%  \hat{\;Y_i } =b{\;}_1 +b_2 X_{2i} +\ldotp \ldotp \ldotp +b_k X_{\textrm{ki}}
% The residual of observation i is e_i =Y_i -\hat{\;Y_i }. The OLS (Ordinary Least Square) estimator of \left(\beta_1 ,\ldotp \ldotp \ldotp ,\beta_k \right) is the vector \left(b_1 ,\ldotp \ldotp \ldotp ,b_k \right) which minimizes the sum of squared residuals, i.e. solves:
% \min_{b_1 ,\ldotp \ldotp \ldotp ,b_k } \textrm{RSS}=\sum_{i=1}^n e_i^2 =\sum_{i=1}^n {\left(Y_i -b{\;}_1 -b_2 X_{2i} -\ldotp \ldotp \ldotp -b_k X_{\textrm{ki}} \right)}^2
% The first order necessary conditions for this this quadratic optimization problem is a system of k linear equations with k unknowns, b_1 ,\ldotp \ldotp \ldotp ,b_k.

clear
close all

%
%% Statistics and Machine Learning Toolbox
% Matlab's Statistics and Machine Learning Toolbox can estimate the multivariate linear regression models using the command fitlm, which is similar to the lm() in R.
%% Reading the data from the web

wage = webread("https://dl.dropboxusercontent.com/s/32m1t3iyno319j6/wage21.csv?dl=0");

%
%% Estimating the model
% The next command generates output very similar to R.

m = fitlm(wage,'EARNINGS ~ S + EXP');
m

%
%% Matrix Algebra
% In this section we demonstrate what statistical packages are doing when they are estimating the multivariate regression model.
%% OLS problem in matrix notation
% The model in the introduction can be presented in matrix form:
% Y={\left\lbrack \begin{array}{c}
% Y_1 \\
% \vdots \\
% Y_n
% \end{array}\right\rbrack }_{n\times 1} ,X={\left\lbrack \begin{array}{cccc}
% 1 & X_{21}  & \cdots  & X_{\textrm{k1}} \\
% \vdots  & \vdots  & \ddots  & \vdots \\
% 1 & X_{2n}  & \cdots  & X_{\textrm{kn}}
% \end{array}\right\rbrack }_{n\times k} ,\beta \;={\left\lbrack \begin{array}{c}
% \beta_1 \\
% \vdots \\
% \beta_k
% \end{array}\right\rbrack }_{k\times 1} ,\;u={\left\lbrack \begin{array}{c}
% u_1 \\
% \vdots \\
% u_n
% \end{array}\right\rbrack }_{n\times 1}
% the model becomes
% Y=X\beta +u
% The fitted equation and residuals of vecotr of estimates b={\left\lbrack \begin{array}{ccc}
% b_1  & \cdots  & b_k
% \end{array}\right\rbrack }^{\prime } are
% \hat{\;Y} =\textrm{Xb},\;\;e=Y-\hat{\;Y}
% The OLS problem in matrix form is therefore:
% \min_b \textrm{RSS}=e^{\prime } e={\left(Y-\textrm{Xb}\right)}^{\prime } \left(Y-\textrm{Xb}\right)
% The solution is
% b={\left(X^{\prime } X\right)}^{-1} X^{\prime } Y
%
%% Solving the OLS problem with matrix algebra
% We prepare the vector of dependent variable Y and the design matrix X, containing column of 1s and the regressors X_2 ,\ldotp \ldotp \ldotp ,X_k.

Y = wage(:,{'EARNINGS'}); %Extracting the dependent variable
X = wage(:,{'S','EXP'}); %Extracting the needed regressors
Y = table2array(Y); %Converting to matrix
X = table2array(X); %Converting to matrix
X = [ones(length(X),1), X]; %Adding vector of 1s for X

%
% Next we use the analytical solution for OLS coefficients, b={\left(X^{\prime } X\right)}^{-1} X^{\prime } Y.

b = (X'*X)\X'*Y %Gaussian elimination (faster, more accurate)
b_inv = inv(X'*X)*X'*Y %Matrix inversion

%
% Matrix algebra is used to obtain elegan formulas for decomposition of sum of squares:
% \begin{array}{l}
% \sum_{i=1}^n {\left(Y_i -\bar{\;Y} \right)}^2 =\sum_{i=1}^n {\left({\hat{\;Y} }_i -\bar{\;Y} \right)}^2 +\sum_{i=1}^n {\left(Y_i -\hat{\;Y_i } \right)}^2 \\
% \textrm{TSS}=\textrm{ESS}+\textrm{RSS}
% \end{array}
% In matrix form, these become:
% \begin{array}{l}
% \textrm{TSS}=Y^{\prime } Y-n\bar{\;Y^2 } \\
% \textrm{ESS}=b^{\prime } X^{\prime } \textrm{Xb}-n\bar{\;Y^2 } \\
% \textrm{RSS}=e^{\prime } e\\
% R^2 =\frac{\textrm{ESS}}{\textrm{TSS}}
% \end{array}
% The estimated covariance matrix of OLS estimators:
% \textrm{cov}\left(b\right)={\left(X^{\prime } X\right)}^{-1} \frac{e^{\prime } e}{n-k}
% where k is the number of number of estimated coefficients.
% Thus, the standard errors of OLS estimators is the the square root of the diagonal of the covariance matrix.
% Implementing the above in Matlab is straightforward.

Y_hat = X*b; %Fitted (predicted) values
e = Y - Y_hat; %Residuals
n = size(X,1); %Number of observations
k = size(X,2); %Number of regressors

C = (X'*X)\eye(k)*(e'*e)/(n-k); % cov(b)
se = sqrt(diag(C)) %Standard errors of b

TSS = Y'*Y - n*mean(Y)^2; %Total Sum of Squares
ESS = b'*(X'*X)*b - n*mean(Y)^2; %Explained Sum of Squares
RSS = e'*e; %Residuals Sum of Squares
R2 = ESS/TSS; %Ordinary R-squared
R2_adj = 1 - (1-R2)*(n-1)/(n-k); %Adjusted R-squared

% Presenting the results
disp(['TSS = ', num2str(TSS)]) 
disp(['ESS = ', num2str(ESS)]) 
disp(['RSS = ', num2str(RSS)]) 
disp(['R^2 = ', num2str(R2)]) 
disp(['R^2_adj = ', num2str(R2_adj)]) 

%
% Comparing our results to the ones obtained by Matlab's fitlm() function:

disp(['TSS = ', num2str(m.SST)]) %Total Sum of Squares
disp(['ESS = ', num2str(m.SSR)]) %Explained Sum of Squares
disp(['RSS = ', num2str(m.SSE)]) %Residuals Sum of Squares
disp('R^2 = ') %R-squared (ordinary and adjusted)
disp(m.Rsquared)