load WAGE1.mat;
dd = WAGE1;
n=size(dd,1);%number of rows of dd

  %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); %years of educ
x_2 = dd(:,3); % years of experience
x_3 = dd(:,4); % years with current employer
x_4 = dd(:,8);
one = ones(size(y,1),1);

% First Regression Model

[b, bint,r, rint, stats]=regress(y,[one x_1]);
sigma_hat2 =(1/(size(y,1)-2))*sum(r.^2);
b
bint
sigma_hat2 
stats

% defining a matrix of regressors
X=[one x_1];

% estimating the variance of the LS estimators
var_b = sigma_hat2*inv(X'*X);

% taking the square root of the diagonal elements
est_sdb = sqrt(diag(var_b));

% forming the t-statistic 
t = (b(2)-0.07)./est_sdb(2);

% obtaining the critical value for t

t_critical=tinv(.975,n-2);

t
t_critical


%Forecast
yhat = b(1)+b(2)*x_1;

% Graph
xf = [x_1 y yhat];
xf=sortrows(xf,1);

figure (1)

plot(xf(:,1),xf(:,2),'o',xf(:,1),xf(:,3))
title('Estimated regression of Log(wages) on Education')
xlabel('Educ')
ylabel('Log(wages)')

% % Second Regression Model
% [beta, betaint,r2, stats2]=regress(y,[one x_1 x_2 x_3]);
% nsigma_hat2 =(1/(size(y,1)-4))*sum(r2.^2);
% nsigma_hat2 
% beta
% betaint
% % Graph
% figure (2)
% scatter3(x_1,x_3,y,'filled')
% hold on
% [X1FIT,X2FIT] = meshgrid(x_1,x_3);
% YFIT = beta(1) + beta(2)*X1FIT + beta(3)*X2FIT;
% surf(X1FIT,X2FIT,YFIT)
% xlabel('Educ')
% ylabel('Experience')
% zlabel('Log(wages)')
% %view(50,10)
% hold off

fname = sprintf('results_OLSWage.txt');
save(fname,'b','bint','sigma_hat2','-ascii') 