load WAGE1.mat;
dd = WAGE1;
n=size(dd,1);

  %Obs:   526

  %1. wage                     average hourly earnings
  %2. educ                     years of education
  %3. exper                    years potential experience
  %4. tenure                   years with current employer
  %5. nonwhite                 =1 if nonwhite
  %6. female                   =1 if female
  %7. married                  =1 if married
  %8. numdep                   number of dependents
  %9. smsa                     =1 if live in SMSA
 %10. northcen                 =1 if live in north central U.S
 %11. south                    =1 if live in southern region
 %12. west                     =1 if live in western region
 %13. construc                 =1 if work in construc. indus.
 %14. ndurman                  =1 if in nondur. manuf. indus.
 %15. trcommpu                 =1 if in trans, commun, pub ut
 %16. trade                    =1 if in wholesale or retail
 %17. services                 =1 if in services indus.
 %18. profserv                 =1 if in prof. serv. indus.
 %19. profocc                  =1 if in profess. occupation
 %20. clerocc                  =1 if in clerical occupation
 %21. servocc                  =1 if in service occupation
 %22. lwage                    log(wage)
 %23. expersq                  exper^2
 %24. tenursq                  tenure^2
 
y = dd(:,22);
x_1 = dd(:,2);
x_2 = dd(:,3);
x_3 = dd(:,4);
%x_4 = dd(:,8);

% X is the matrix of regressors

bfone = ones(size(y,1),1);

fprintf('Number of observations is: %d', size(y,1));

X = [ bfone x_1 x_2 x_3 ];

% phi_hat is the least squares estimate

beta_hat = (X'*X)\(X'*y);
beta_hat

% sigma_hat2

sigma_hat2 = (1/(size(y,1)-size(X,2)))*(y-X*beta_hat)'*(y-X*beta_hat);
sigma_hat2

% Variance of beta_hat

format long

var_beta = sigma_hat2*(X'*X)^(-1);
var_beta

diag(var_beta,0)

% All t-statistics

t_vec=beta_hat.*(sqrt(diag(var_beta,0))).^(-1);

t_vec

% Degrees off reedom

df = size(y,1)-size(X,2);
df

% Critical value

alpha = 0.05;
tc=tinv(1-(alpha/2),df);
tc

%Confidence interval

u_limit = beta_hat + tc*sqrt(diag(var_beta,0));
l_limit = beta_hat - tc*sqrt(diag(var_beta,0));

%R-Squared

rss = (y-X*beta_hat)'*(y-X*beta_hat);
tss = (y-mean(y))'*(y-mean(y));
R_squared = 1- (rss/tss);
R_squared


% Multiple Hypothesis Testing
R = [0 1 0 0  ;
     0 0 1 0  ];
r = [0;0];
 
m= size(R,1);

dm = R*beta_hat-r;

f = ((dm'*((R*((X'*X)\R'))\dm))/m)/sigma_hat2;

f

fc=finv(1-alpha,m,df);
fc

otp = [beta_hat t_vec l_limit u_limit];
otp

%Forecast

y_hat = X*beta_hat;

xf = [y y_hat];
corr(xf)
xf=sortrows(xf,1);
plot(xf(:,1),xf(:,2),'+')

fname = sprintf('results.txt');
save(fname,'otp','-ascii') 