function result = tsls(y,x,w,noconst,display)
% TSLS (Two-Stage Least Square) estimation of multiple linear regression
%   result = iv(y,x,w)
% Inputs:
%   y: dependent variable, n by 1
%   x: independent variable, n by p
%   w: instrumental variable, n by q, q>=p
%   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 : TSLS 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
%       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,p]=size(x);
P = w*inv(w'*w)*(w');
b = (x'*P*x)\(x'*P*y);
e = y-x*b;
ssr = e'*e;
y1 = y-mean(y);
sst = y1'*y1;
R2 = 1-ssr/sst;
s2 = ssr/(n-p);
bvar =  s2/n*inv(x'*P*x);
se = sqrt(diag(bvar));
bhrvar = (x'*P*x)\((repmat(P*e,1,p).*x)'*(repmat(P*e,1,p).*x))/(x'*P*x);
hrse = sqrt(diag(bhrvar));
t = b./hrse;
pvalue = 2*(1-tcdf(abs(t),n-p));
if noconst
    Wald = b'*(bhrvar\b);
else
    Wald = b(2:p)'*(bhrvar(2:p,2:p)\b(2:p));
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

