clear all; close all; clc; %parameters Omega=[ 4.456, -0.274, 0.227; -0.274, 5.323, 0.017; 0.227, 0.017, 5.247]; be=[1; 0.5; -0.5; 0.25]; mu=[2;2;3]; P=Omega^(1/2); N=100; J=500; for j=1:J %generate samples for j=1:J j ksi=randn(3,N); X=mu*ones(1,N)+P*ksi; X=[ones(N,1) X']; Y=X*be+trnd(5,size(X*be)); Z=[X (X(:,2:4).^2)]; %OLS beta0=inv(X'*X)*(X'*Y); %GMM, Identity matrix V3=eye(size(Z,2)); beta3=fminsearch(@gmm,beta0,[],X,Y,Z,V3); A=zeros(size(X,2),size(Z,2)); V3=zeros(size(Z,2)); for i=1:N A=A+X(i,:)'*Z(i,:)/N; V3=V3+Z(i,:)'*Z(i,:)*((Y(i)-X(i,:)*beta3)^2)/N; end %GMM optimal matrix 1 [beta3,fval3]=fminsearch(@gmm,beta3,[],X,Y,Z,V3); V3=zeros(size(Z,2)); M3=zeros(size(Z(1,:)')); for i=1:N V3=V3+Z(i,:)'*Z(i,:)*((Y(i)-X(i,:)*beta3)^2)/N; end %GMM optimal matrix 2 [beta3,fval3]=fminsearch(@gmm,beta0,[],X,Y,Z,V3); for i=1:N M3=M3+Z(i,:)'*(Y(i)-X(i,:)*beta3)/N; end VAR3=inv(A*inv(V3)*A')/N; %Wald test b(2)*b(3)+b(4)=0 beta=beta3; r=beta(2)*beta(3)+beta(4); R=[0 beta(3) beta(2) 1]; W(j)=r'*inv(R*VAR3*R')*r; %Restricted estimate options = optimset('LargeScale','off','Display','off'); [beta5,fval5]=fmincon(@gmm,beta3,[],[],[],[],[],[],@restr,options,X,Y,Z,V3); V5=zeros(size(Z,2)); M5=zeros(size(Z(1,:)')); for i=1:N V5=V5+Z(i,:)'*Z(i,:)*((Y(i)-X(i,:)*beta5)^2)/N; M5=M5+Z(i,:)'*(Y(i)-X(i,:)*beta5)/N; end %Lagrange Multiplier test LM(j)=N*M5'*inv(V3)*A'*inv(A*inv(V5)*A')*A*inv(V3)*M5; %Likelihood Ration test LR(j)=N*(M5'*inv(V3)*M5-M3'*inv(V3)*M3); end %plot results figure(1); hist(W,50) figure(2); hist(LM,50) figure(3); hist(LR,50) mean([W;LM;LR]') median([W;LM;LR]') [fw,xw]=ecdf(W); [fm,xm]=ecdf(LM); [fr,xr]=ecdf(LR); p1W=1-interp1(xw(2:J),fw(2:J),3.84) p1LM=1-interp1(xm(2:J),fm(2:J),3.84) p1LR=1-interp1(xr(2:J),fr(2:J),3.84) [xr,fr]=stairs(xr,fr); [xw,fw]=stairs(xw,fw); [xm,fm]=stairs(xm,fm); figure(4) plot(xw,fw,xm,fm,xr,fr) legend('Wald','LM','LR')