clear;
%define variables
[EDU,SOUTH,NONWH,HISP,FE,MAR,MARRFE,EX,EXSQ,UNIO,LNWAGE,AGE,D1,D2,D3,D4,D5,D6,D7] = textread('CPS85.txt','%d %d %d %d %d %d %d %d %d %d %f %d %d %d %d %d %d %d %d');
ON=ones(size(EDU));

%check scan of data
size(EDU);
ON'*(MAR.*FE-MARRFE);
ON'*(EX.*EX-EXSQ);

% (a) average hourly wage: algebraic and geometric
disp('(a)');
WAGE=exp(LNWAGE);
wage_arith=(WAGE'*ON)/(ON'*ON);
wage_geom=exp((LNWAGE'*ON)/(ON'*ON));
disp('               arith        geom');
disp(['     wage     ' num2str([wage_arith wage_geom])]);

% (b) sample means: 
disp('(b)');
south_fr=(SOUTH'*ON)/(ON'*ON);
fe_fr=(FE'*ON)/(ON'*ON);
unio_fr=(UNIO'*ON)/(ON'*ON);
nonwh_fr=(NONWH'*ON)/(ON'*ON);
hisp_fr=(HISP'*ON)/(ON'*ON);
disp('               males   females   south     non-south   union  non-union    white   hispanic    black');
disp([1 1-fe_fr    fe_fr   1-south_fr  south_fr    unio_fr 1-unio_fr   1-nonwh_fr-hisp_fr  hisp_fr  nonwh_fr]);
% (c) means and stdevs
disp('(c)');
%LOG WAGES
lnwage_av=(LNWAGE'*ON)/(ON'*ON);
lnwage_sd=sqrt((((LNWAGE-lnwage_av).^2)'*ON)/(ON'*ON));
%female
lnwage_f_av=(LNWAGE'*FE)/(FE'*FE);
lnwage_f_sd=sqrt((((LNWAGE-lnwage_f_av).^2)'*FE)/(FE'*FE));
%male
lnwage_m_av=(LNWAGE'*(1-FE))/((1-FE)'*(1-FE));
lnwage_m_sd=sqrt((((LNWAGE-lnwage_m_av).^2)'*(1-FE))/((1-FE)'*(1-FE)));
%white
lnwage_w_av=(LNWAGE'*(1-NONWH-HISP))/((1-NONWH-HISP)'*(1-NONWH-HISP));
lnwage_w_sd=sqrt((((LNWAGE-lnwage_w_av).^2)'*(1-NONWH-HISP))/((1-NONWH-HISP)'*(1-NONWH-HISP)));
%hispanic
lnwage_h_av=(LNWAGE'*HISP)/(HISP'*HISP);
lnwage_h_sd=sqrt((((LNWAGE-lnwage_h_av).^2)'*HISP)/(HISP'*HISP));
%black
lnwage_b_av=(LNWAGE'*NONWH)/(NONWH'*NONWH);
lnwage_b_sd=sqrt((((LNWAGE-lnwage_b_av).^2)'*NONWH)/(NONWH'*NONWH));
%union
lnwage_u_av=(LNWAGE'*UNIO)/(UNIO'*UNIO);
lnwage_u_sd=sqrt((((LNWAGE-lnwage_u_av).^2)'*UNIO)/(UNIO'*UNIO));
%non-union
lnwage_n_av=(LNWAGE'*(1-UNIO))/((1-UNIO)'*(1-UNIO));
lnwage_n_sd=sqrt((((LNWAGE-lnwage_n_av).^2)'*(1-UNIO))/((1-UNIO)'*(1-UNIO)));

% EDUCATION
edu_av=(EDU'*ON)/(ON'*ON);
edu_sd=sqrt((((EDU-edu_av).^2)'*ON)/(ON'*ON));
%female
edu_f_av=(EDU'*FE)/(FE'*FE);
edu_f_sd=sqrt((((EDU-edu_f_av).^2)'*FE)/(FE'*FE));
%male
edu_m_av=(EDU'*(1-FE))/((1-FE)'*(1-FE));
edu_m_sd=sqrt((((EDU-edu_m_av).^2)'*(1-FE))/((1-FE)'*(1-FE)));
%white
edu_w_av=(EDU'*(1-NONWH-HISP))/((1-NONWH-HISP)'*(1-NONWH-HISP));
edu_w_sd=sqrt((((EDU-edu_w_av).^2)'*(1-NONWH-HISP))/((1-NONWH-HISP)'*(1-NONWH-HISP)));
%hispanic
edu_h_av=(EDU'*HISP)/(HISP'*HISP);
edu_h_sd=sqrt((((EDU-edu_h_av).^2)'*HISP)/(HISP'*HISP));
%black
edu_b_av=(EDU'*NONWH)/(NONWH'*NONWH);
edu_b_sd=sqrt((((EDU-edu_b_av).^2)'*NONWH)/(NONWH'*NONWH));
%union
edu_u_av=(EDU'*UNIO)/(UNIO'*UNIO);
edu_u_sd=sqrt((((EDU-edu_u_av).^2)'*UNIO)/(UNIO'*UNIO));
%non-union
edu_n_av=(EDU'*(1-UNIO))/((1-UNIO)'*(1-UNIO));
edu_n_sd=sqrt((((EDU-edu_n_av).^2)'*(1-UNIO))/((1-UNIO)'*(1-UNIO)));

%EXPERIENCE
ex_av=(EX'*ON)/(ON'*ON);
ex_sd=sqrt((((EX-ex_av).^2)'*ON)/(ON'*ON));
%female
ex_f_av=(EX'*FE)/(FE'*FE);
ex_f_sd=sqrt((((EX-ex_f_av).^2)'*FE)/(FE'*FE));
%male
ex_m_av=(EX'*(1-FE))/((1-FE)'*(1-FE));
ex_m_sd=sqrt((((EX-ex_m_av).^2)'*(1-FE))/((1-FE)'*(1-FE)));
%white
ex_w_av=(EX'*(1-NONWH-HISP))/((1-NONWH-HISP)'*(1-NONWH-HISP));
ex_w_sd=sqrt((((EX-ex_w_av).^2)'*(1-NONWH-HISP))/((1-NONWH-HISP)'*(1-NONWH-HISP)));
%hispanic
ex_h_av=(EX'*HISP)/(HISP'*HISP);
ex_h_sd=sqrt((((EX-ex_h_av).^2)'*HISP)/(HISP'*HISP));
%black
ex_b_av=(EX'*NONWH)/(NONWH'*NONWH);
ex_b_sd=sqrt((((EX-ex_b_av).^2)'*NONWH)/(NONWH'*NONWH));
%union
ex_u_av=(EX'*UNIO)/(UNIO'*UNIO);
ex_u_sd=sqrt((((EX-ex_u_av).^2)'*UNIO)/(UNIO'*UNIO));
%non-union
ex_n_av=(EX'*(1-UNIO))/((1-UNIO)'*(1-UNIO));
ex_n_sd=sqrt((((EX-ex_n_av).^2)'*(1-UNIO))/((1-UNIO)'*(1-UNIO)));

%           total   male    female  white   hisp    black   unio   non-unio
% lnwage
% edu
% ex

AV=[lnwage_av lnwage_m_av lnwage_f_av lnwage_w_av lnwage_h_av lnwage_b_av lnwage_u_av lnwage_n_av;
edu_av edu_m_av edu_f_av edu_w_av edu_h_av edu_b_av edu_u_av edu_n_av;
ex_av ex_m_av ex_f_av ex_w_av ex_h_av ex_b_av ex_u_av ex_n_av
];
disp('     Total   males   females     white     hispanic    black   union   non-union');
disp(AV);

SD=[lnwage_sd lnwage_m_sd lnwage_f_sd lnwage_w_sd lnwage_h_sd lnwage_b_sd lnwage_u_sd lnwage_n_sd;
edu_sd edu_m_sd edu_f_sd edu_w_sd edu_h_sd edu_b_sd edu_u_sd edu_n_sd;
ex_sd ex_m_sd ex_f_sd ex_w_sd ex_h_sd ex_b_sd ex_u_sd ex_n_sd
];
disp('     Total   males   females     white     hispanic    black   union   non-union');
disp(SD);


%(d) 
disp('(d)');
Y=LNWAGE;
X=[ON EDU EX];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(3)));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const    EDU        EX');
disp([beta ste t pv]');

%(e)
disp('(e)   R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

% (f)
disp('(f) P-value of coef on EX');
disp(2-2*tcdf(abs(t(3)),size(X,1)-size(X,2)))

% (g) 
disp('(g) ');
Y=LNWAGE;
X=[ON EDU EX EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const    EDU        EX       EX^2');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);
disp('Age of maximum');
disp(-0.5*beta(3)/beta(4));


% (h) 
disp('(h) ');
Y=LNWAGE;
X=[ON FE EDU EX EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const    FE       EDU        EX       EX^2');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

GE=[0 1 0 0 0];
t=GE*beta./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: no gender differences    ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);

% (i)
disp('(i) ');
Y=LNWAGE;
X=[(1-FE) FE EDU EX EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     MA        FE       EDU        EX       EX^2');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

GE=[1 -1 0 0 0];
t=GE*beta./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: no gender differences    ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);

% (j)
disp('(j) ');
Y=LNWAGE;
X=[(1-FE) FE EDU.*(1-FE) EDU.*FE EX EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     MA        FE       EDU*MA    EDU*FE      EX       EX^2');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

GE=[0 0 1 -1 0 0];
t=GE*beta./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: no change   ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);

% (k)
disp('(k) ');
Y=LNWAGE;
X=[ON EDU EX EX.*EDU EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const     EDU      EX     EDU*EX       EX^2 ');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

GE=[0 0 0 1 0];
t=GE*beta./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: no change   ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);

Y=LNWAGE;
X=[ON EDU EX EXSQ EDU.*EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const     EDU      EX     EX^2       EDU*EX^2 ');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);
beta(5)
ste(5)

GE=[0 0 0 0 1];
t=GE*beta./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: no change   ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);

Y=LNWAGE;
X=[ON EDU EX EXSQ EDU.*EX EDU.*EXSQ];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const     EDU      EX     EX^2       EDU*EX      EDU*EX^2 ');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

GE=[0 0 0 0 1 0
    0 0 0 0 0 1];
F=(((GE*beta)')*inv(GE*V*GE')*(GE*beta))/size(GE,1);
disp(['P-value of hypothesis: no change   ' num2str(fpdf(F,size(GE,1),size(X,1)-size(X,2)))]);
disp(' ');
disp('corr EDU      EX        EX^2      EDU*EX    EDU*EX^2 ');
disp(corrcoef([EDU EX EXSQ EDU.*EX EDU.*EXSQ]));

% (l)
disp('(k) ');
Y=LNWAGE;
X=[ON EDU EX EXSQ UNIO];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));
ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     const     EDU      EX        EX^2        UNIO');
disp([beta ste t pv]');
disp('R-squared');
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

GE=[0 0 0 0 1];
t=GE*beta./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: no change   ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);
