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

X=[X1 X2 X3 X4 X5];
Z=[X X(:,2:5).*X(:,2:5)];% X1.*X2 X1.*X3 X1.*X4 X1.*X5];
N=size(Y,1);

%OLS
beta0=inv(X'*X)*(X'*Y) %OLS
VAR0=(((Y-X*beta0)'*(Y-X*beta0))/(size(X,1)-size(X,2)))*inv(X'*X);  err0=Y-X*beta0;
GE=eye(size(X,2));  ste0=sqrt(spdiags((GE*VAR0*GE').*eye(size(X,2))))

%GMM Z type 1, I    %a1
V1=eye(size(X,2));
beta1=fminsearch(@gmm,beta0,[],X,Y,X,V1)
A=zeros(size(X,2),size(X,2));
V1=zeros(size(X,2));
for i=1:N
    A=A+X(i,:)'*X(i,:)/N;
    V1=V1+X(i,:)'*X(i,:)*((Y(i)-X(i,:)*beta1)^2)/N;
end
VAR1=inv(A*A')*A*inv(V1)*A'*inv(A*A')/N;    err1=Y-X*beta1;
GE=eye(size(X,2));  ste1=sqrt(spdiags((GE*VAR1*GE').*eye(size(X,2)))) 
%GMM Z type 1, Vn   %a2
beta2=fminsearch(@gmm,beta0,[],X,Y,X,V1)
VAR2=inv(A*inv(V1)*A')/N;    err2=Y-X*beta2;
GE=eye(size(X,2));  ste2=sqrt(spdiags((GE*VAR2*GE').*eye(size(X,2)))) 

%GMM Z type 2, I    %b1
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
VAR3=inv(A*A')*A*inv(V3)*A'*inv(A*A')/N;    err3=Y-X*beta3;
GE=eye(size(X,2));  ste3=sqrt(spdiags((GE*VAR3*GE').*eye(size(X,2))))
%GMM Z type 2, Vn   %b2
beta4=fminsearch(@gmm,beta0,[],X,Y,Z,V3)
VAR4=inv(A*inv(V3)*A')/N;    err4=Y-X*beta4;
GE=eye(size(X,2));  ste4=sqrt(spdiags((GE*VAR4*GE').*eye(size(X,2))))

%GMM Z type 2, Vn optimal   %b2
V5=zeros(size(Z,2));
for i=1:N
    V5=V5+Z(i,:)'*Z(i,:)*((Y(i)-X(i,:)*beta4)^2)/N;
end
beta5=fminsearch(@gmm,beta0,[],X,Y,Z,V5)
VAR5=inv(A*inv(V5)*A')/N;    err5=Y-X*beta5;
GE=eye(size(X,2));  ste5=sqrt(spdiags((GE*VAR5*GE').*eye(size(X,2))))

%disp('    c        X2       X3         X4');
%disp([beta2 ste2 t2 pv2]');

%Wald test
beta=beta5;
delta=beta(2)^2+beta(4)^2-beta(3)^2-beta(5)^2;
GE=[0 2*beta(2) -2*beta1(3) 2*beta(4) -2*beta(5)];
xi2=delta'*inv(GE*VAR5*GE')*delta;
disp(['P-value of hypothesis: ' num2str(1-chi2cdf(xi2,1))]);