clear all;  close all;  clc
%read series
[LFP, WHRS, KL6, K618, WA, WE, WW, RPWG, HHRS, HA, HE, HW, FAMINC, MTR, WMED, WFED, UN, CIT, EX] = textread('mroz.asc','%d %d %d %d %d %d %f %f %d %d %d %f %d %f %d %d %f %d %d');
%generate x's
ON=ones(size(LFP)); PRIN=FAMINC-WHRS.*WW;   EX2=EX.*EX;
%estimate LWW
LW=log(WW(1:428));  Y=LW;   X=[ON WA WE CIT EX EX2];    X=X(1:428,:);
%OLS on subsample
beta0=inv(X'*X)*(X'*Y);     FLW=[ON WA WE CIT EX EX2]*beta0;
%bombine fitted and actuals
FLW=FLW(429:753);       LWW=[LW;FLW];
%construct main variables
X=[ON KL6 K618 WA WE UN CIT PRIN/100000 LWW];      Y=LFP;
%LPM model
disp('LPM model');      beta1=inv(X'*X)*(X'*Y);
V=(((Y-X*beta1)'*(Y-X*beta1))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));      ste1=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t1=GE*beta1./ste1;      pv1=2-2*tcdf(abs(t1),size(X,1)-size(X,2));
disp('    c        KL6      K618       WA       WE        UN        CIT       PRIN       LWW ');
disp([beta1 ste1 t1 pv1]');     disp('R-squared');      R_2=zeros(1,2);     
R_2(1,1)=1-((Y-X*beta1)'*(Y-X*beta1))/((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);

%marginal effect of LWW
Xav=X'*ON/size(ON,1);
mf1=beta1;  se1=ste1;
disp(['marginal effect of LWW    ' num2str([mf1(9) se1(9)])]);

%Probit model
disp('Probit model');   
beta2=fminsearch(@myf2,beta1,... 
                    optimset('TolX',1e-4),[Y X]);
Om=zeros(size(X,2));    H=zeros(size(X,2));
for i=1:size(X,1)
G=normcdf(X(i,:)*beta2);
C=normpdf(X(i,:)*beta2);
Om=Om+X(i,:)'*X(i,:)*((Y(i,:)-G)*C*(Y(i,:)-G)*C)/size(X,1);
H=H+X(i,:)'*(C*C+(Y(i,:)-G)*(X(i,:)*beta2)*C)*X(i,:)/size(X,1);
end
V=inv(H)*Om*inv(H)/size(X,1);
GE=eye(size(X,2));      ste2=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t2=GE*beta2./ste2;      pv2=2-2*normcdf(abs(t2));
disp('    c        KL6      K618       WA       WE        UN        CIT       PRIN       LWW ');
disp([beta2 ste2 t2 pv2]');
%marginal effect of LWW
Xav=X'*ON/size(ON,1);
a=Xav'*beta2;   f=normpdf(a);   mf2=f*beta2;  se2=f*(a*a+1)*ste2;
disp(['marginal effect of LWW    ' num2str([mf2(9) se2(9)])]);

%Logit model
disp('Logit model');
beta3=fminsearch(@myf3,beta1,... 
                    optimset('TolX',1e-4),[Y X]);
Om=zeros(size(X,2));    H=zeros(size(X,2));
for i=1:size(X,1)
G=exp(X(i,:)*beta3)./(1+exp(X(i,:)*beta3));
C=exp(X(i,:)*beta3)./((1+exp(X(i,:)*beta3)).^2);
Om=Om+X(i,:)'*X(i,:)*((Y(i,:)-G)*C*(Y(i,:)-G)*C)/size(X,1);
H=H+X(i,:)'*(C*C+(Y(i,:)-G)*(X(i,:)*beta2)*C)*X(i,:)/size(X,1);
end
V=inv(H)*Om*inv(H)/size(X,1);
GE=eye(size(X,2));      ste3=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t3=GE*beta3./ste3;      pv3=2-2*normcdf(abs(t3));
disp('    c        KL6      K618       WA       WE        UN        CIT       PRIN       LWW ');
disp([beta3 ste3 t3 pv3]');
disp('Logit model corrected by 0.5513');
disp('    c        KL6      K618       WA       WE        UN        CIT       PRIN       LWW ');
disp([beta3*0.5513 ste3*0.5513 t3 pv3]');

%marginal effect of LWW
Xav=X'*ON/size(ON,1);
a=Xav'*beta3;   ea=exp(a);   f=ea/(1+ea);  mf3=f*beta3;  
fi=(ea*(1-ea)/((ea+1)^3))*a+ea/((1+ea)^2);  se3=fi*ste3;
disp(['marginal effect of LWW    ' num2str([mf3(9) se3(9)])]);

%rescaled marginal effects
Xav=X'*ON/size(ON,1);
a=Xav'*beta3*0.5513;   ea=exp(a);   f=ea/(1+ea);  mf3=f*beta3*0.5513;  
fi=(ea*(1-ea)/((ea+1)^3))*a+ea/((1+ea)^2);  se3=fi*ste3;
disp(['marginal effect of LWW    ' num2str([mf3(9) se3(9)])]);


er1=Y-X*beta1;
er2=Y-normcdf(X*beta2);
er3=Y-exp(X*beta3)./(1+exp(X*beta3));
[t1,x1]=ecdf(er1);
[t2,x2]=ecdf(er2);
[t3,x3]=ecdf(er3);
plot(x1,t1,'-',x2,t2,'.-',x3,t3,':',...
                    'LineWidth',2)
legend('LPM','Probit','Logit')