% import data from the file runner.mat

load runner;

% compute sample size

n=length(Y);

%%%%%%%%%%%%%%%%% First linear regession model, only using %%%%%%%%%%%%%%%%
%                  the variable X5, aerobic capacity, V_O2.        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Contruct matrix X for the independent variable X5, aerobic capacity, V_O2.

X=[ones(n,1) X5];

p=1; % number of independent variables

% the least square estimate of beta
beta=inv(X'*X)*X'*Y;

% Compute totol variation of responsible variable of Y, SSTO
SSTO=sum((Y-mean(Y))'*(Y-mean(Y))); 

% Compute the variation of Regression, SSR
Y_hat=X*beta;
SSR=sum((Y_hat-mean(Y_hat))'*(Y_hat-mean(Y_hat)));

%Compute the variation of residuals
Error=Y-Y_hat; % compute residuals
SSE=Error'*Error;

% Compute Mean square of SSR and Mean square of SSE
MSR=SSR/p;
MSE=SSE/(n-p-1);

% Compute R^2 for this model
R_square=SSR/SSTO;

% Compute F-statistics
F=MSR/MSE;



% compute the estimate of sigma, standard deviation of white noise

s_square=Error'*Error/(n-2);
s=sqrt(s_square);



% compute estimated covariance matrix of the least-square estimate beta_hat
beta_var=s_square*inv(X'*X);

% compute estimate standard error of the least-square estimate beta_hat
beta_error=sqrt(diag(beta_var));

% compute t-statistics for beta

t_beta=beta./sqrt(diag(beta_var));
   
t_alpha=-tinv(0.95,n-p-1);


% 95 percent prediction interval for the running time when x_k=55.  
x_k=55;
y_pred=[1 x_k]*beta;

% Compute forecast error
forecast_error=s*(1+[1 x_k]*inv(X'*X)*[1 x_k]')^0.5;

% Compute confidence interval [CI_down, CI_up]
CI_up=y_pred+tinv(0.975,n-p-1)*forecast_error;
CI_down=y_pred-tinv(0.975,n-p-1)*forecast_error;


%%%%%%%%%%%%%%Second Model using all independent variables %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Contruct matrix X_all for all independent variables from X1 to X5.

X_all=[ones(n,1) X1 X2 X3 X4 X5];

% the least square estimate of beta for the parameter from beta_0 to beta_5

beta_all=inv(X_all'*X_all)*X_all'*Y;
p_all=5; % number of parameters

% compute fitted value y_hat
Y_hat_all=X_all*beta_all;

% compute variation of regession 
SSR_all=sum((Y_hat_all-mean(Y_hat_all))'*(Y_hat_all-mean(Y_hat_all)));

% compute mean square of SSR
MSR_all=SSR_all/p_all;

% compute residuals 
Error_all=Y-Y_hat_all;

% compute variation of residuals
SSE_all=Error_all'*Error_all;

% compute mean square of SSE
MSE_all=SSE_all/(n-p_all-1);

% compute the estimate standar deviation, sigma, of white noise, 
s_square_all=Error_all'*Error_all/(n-p_all-1);
s_all=sqrt(s_square_all);

% compute the estimate covariance matrix of the least estimate of beta_all 
beta_var_all=s_square_all*inv(X_all'*X_all);

% compute estimate standard error of the least-square estimate beta_hat
beta_error_all=sqrt(diag(beta_var_all));

% compute the t-statistics for beta0 to beta5;
t_beta_all=beta_all./sqrt(diag(beta_var_all));

% compute R^2 and F-statistics for this model with all independent
% variables
R_square_all=SSR_all/SSTO;
F_all=MSR_all/MSE_all;

