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