function result = iv(y,x,w,noconst,display)
% IV estimation of multiple linear regression
%   result = iv(y,x,w)
% Inputs:
%   y: dependent variable, n by 1
%   x: independent variable, n by m
%   w: instrumental variable, n by m
%   noconst: logical variable, 1 for the case of no constant, 0 for the
%       case of having a constant term (default: 0).
% Output:
%   result: a structure containing
%       result.b : IV estimator
%       result.ssr : SSR
%       result.s2 : estimator for var(u)
%       result.bvar : Asymptotic covariance matrix of b
%       result.se : Asympototic standard error of b
%       result.bhrvar : Asymptotic heteroscedasticity-robust covariance matrix of b
%       result.hrse : Asympototic heteroscedasticity-robust standard error of b
%       result.R2 : R square (can be negative, meaningless for iv estimation)
%       result.Wald: Wald statistic for the significance test of the model
%       result.pvalue: p-value for the significance test of parameters

if nargin<4 || isempty(noconst)
    noconst = 0;
end

if nargin<5 || isempty(display)
    display = 0;
end

[n,m]=size(x);
b = (w'*x)\(w'*y);
e = y-x*b;
ssr = e'*e;
y1 = y-mean(y);
sst = y1'*y1;
R2 = 1-ssr/sst;
s2 = ssr/(n-m);
bvar =  s2/n*(w'*x)\(w'*w)/(x'*w)/n;
se = sqrt(diag(bvar));
bhrvar = (w'*x)\((repmat(e,1,m).*w)'*(repmat(e,1,m).*w))/(x'*w);
hrse = sqrt(diag(bhrvar));
t = b./hrse;
pvalue = 2*(1-tcdf(abs(t),n-m));
if noconst
    Wald = b'*(bhrvar\b);
else
    Wald = b(2:m)'*(bhrvar(2:m,2:m)\b(2:m));
end

result.b = b;
result.ssr = ssr;
result.s2 = s2;
result.bvar = bvar;
result.se = se;
result.bhrvar = bhrvar;
result.hrse = hrse;
result.pvalue = pvalue;
result.R2 = R2;
result.Wald = Wald;

if nargout < 1 || display
    disp('     Coef Est    S.E.    H.R.S.E.   p-value')
    disp([b se hrse pvalue])
    disp('N      SSR      Wald')
    disp([n ssr Wald])
end

