clear all;
close all;
clc;
[Y	X1	X2	X3 X4] = textread('hw4p4.txt','%f %f %f %f %f');

X=[X1 X2 X3 X4];

%ML estimates
beta0=inv(X'*X)*(X'*Y);
beta1=fminsearch(@myf2,beta0,[],X,Y);
options = optimset('Display','iter','LargeScale','on','Jacobian','off');
beta2=fsolve(@myf,beta0,options,X,Y);

%WAY1
G=normcdf(X*beta1);
F=normpdf(X*beta1);
E=(((Y.*F)./G)-(((1-Y).*F)./(1-G)));
Om=zeros(size(X,2));
for i=1:size(X,1)
    Om=Om+X(i,:)'*X(i,:)*E(i,1)*E(i,1);
end
V=inv(Om);
GE=eye(size(X,2));      ste1=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t1=GE*beta1./ste1;      pv1=2-2*normcdf(abs(t1));
disp('    c        X2       X3         X4');
disp([beta1 ste1 t1 pv1]');

%WAY2
G=normcdf(X*beta2);
F=normpdf(X*beta2);
E=(((Y.*F)./G)-(((1-Y).*F)./(1-G)));
Om=zeros(size(X,2));
for i=1:size(X,1)
    Om=Om+X(i,:)'*X(i,:)*E(i,1)*E(i,1);
end
V=inv(Om);
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        X2       X3         X4');
disp([beta2 ste2 t2 pv2]');

%Wald test
GE=[0 1 1 0];
t=(GE*beta1)./sqrt(spdiags((GE*V*GE')));
%disp(['P-value of hypothesis: b2+b3=0  ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);
disp(['P-value of hypothesis: b2+b3=0  ' num2str(1-chi2cdf(t^2,1))]);

%marginal effect of LWW
ON=X1;
Xav=X'*ON/size(ON,1);
Yav=Y'*ON/size(ON,1);
a=Xav'*beta1;   f=normpdf(a);   mf2=f*beta1.*Xav/Yav;  se2=f*(a*a+1)*ste2.*Xav/Yav;
disp('elasticities of Y w.r.t. X and their standard errors');
disp(num2str([mf2 se2]));

%Wald test
delta=beta1(2)*beta1(1)/(beta1(3))^2
GE=[beta1(2)/(beta1(3))^2 beta1(1)/(beta1(3))^2 -2*beta1(2)*beta1(1)/(beta1(3))^3 0];
D=sqrt(GE*V*GE')

