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')