// Related Software: Scilab
// License     : BSD (4.4): https://fileexchange.scilab.org/licenses/7
// Functions   : OLS()
// Author (C)  : Mario MAGGI - Italy
// Release     : 1.0
// Date        : 2019-08-30
// Distribution: https://fileexchange.scilab.org/toolboxes/561000
// History     :
//   1.0 - 2019-08-30: first publication on FileExchange
//
function [betas,std_err,t_stats,p_values,sigma,resid,R2,R2_adj,F_stat,F_prob]=OLS(Y,X)
//-------------------------------------------------------------
// linear regression:
// OLS estimate of the model
//        Y = X*betas + epsilon,
// where epsilon is an homoskedastic Gaussian white noise
// INPUTS
// Y the (nx1) vector of the realizations of the endogenous variable
// X the (nxp) matrix of the realizations of p exogenous variables
// OUTPUTS
// betas, the OLS parameters estimates
// std_err, standard errors of the parameter estimates
// T_stats, the t statistics of the parameters estimate
// p_values, the p-values of the parameters estimate
// sigma, the standard deviation of the estimated epsilon
// resid, the residuals
// R2, the R squared
// R2_adj, adjusted R squared
// F_stat, the F statistic of the regression
// F_prob, the p-value of the F statistic
//-------------------------------------------------------------

[n,p] = size(X); //n number of observations, p number of exogenous variables

betas = pinv(X)*Y; //ols parameter estimate
Y_hat = X*betas; //endogenous fitted values
resid = Y-Y_hat; //fitted residuals
RSS = sum((Y_hat-mean(Y)).^2); //regression sum of squares
ESS = sum(resid.^2); //error sum of squares
TSS = sum((Y-mean(Y)).^2); //total sum of squares
sigma = sqrt(ESS/(n-p)); //error standard error estimate
covbetas = sigma^2*inv(X'*X); //covariance matrix of the parameter estimates
std_err = sqrt(diag(covbetas)); //standard errors of the parameter estimates
t_stats = betas./std_err; //t statistic
[altro, p_values] = cdft("PQ",abs(t_stats),(n-p)*ones(t_stats));
p_values = 2*p_values; //p values
R2 = RSS/TSS; //R squared
R2_adj = 1-(1-R2)*(n-1)/(n-p); //adjusted R squared
if p>1
  F_stat = R2*(n-p)/((1-R2)*(p-1)); //F statistic
  [altro, F_prob] = cdff("PQ",F_stat,p-1,n-p); //p-value of F statistic
end //if

rango = rank(X,n*norm(X)*1e-6);
if rango<p
  mprintf('         WARNING!\n the exogenous matrix is rank deficient\n its robust rank is %i < %i',rango,p)
end //if

disp('------------------------------------------------------')
disp('   parameter    std. err.     t stat.      p-value')
disp([betas std_err t_stats p_values])
disp('------------------------------------------------------')
mprintf('noise std. err.         %f\n',sigma)
mprintf('R squaded (adjusted)    %f (%f)\n',R2,R2_adj)
if p>1
  mprintf('F statistic (p-value)   %f (%f)\n',F_stat,F_prob)
end//if

endfunction
