function [waldstat, waldprb] = waldtest(theta,R,q,vc)
% PURPOSE: computes wald test for two regressions
%---------------------------------------------------
% USAGE: [waldstat waldprob] = waldtest(theta,R,q,vc)
%    or: waldtest(theta,R,q,vc), which prints to the screen
%
% Where: Syntax 1 (user specifies the matrix of restrictions)
%        --------
%        theta = parameter estimate (kx1)
%        R = matrix of restrictions (dofxk)
%        q = R*theta (dofx1)
%        vc = var-cov matrix of theta (kxk)
%
%        Syntax 2 (when testing if some elements of theta = scalar)
%        --------
%        theta = parameter estimate (kx1)
%        R = vector (1xdof) identifying elements of theta
%        q = a scalar
%        vc = var-cov matrix of theta
%        Example: if you want to test the hypothesis
%                 that the 3rd and 5th elements of theta
%                 are equal to 0 waldtest(theta,[3,5],0,vc)
%
%        Syntax 3 (when testing if some elements of theta are =)
%        --------
%        theta = parameter estimate (kx1)
%        R = matrix identifying elements of theta
%            0's are simply ignored and used as place holders
%        q = []
%        vc = var-cov matrix of theta
%        Example: if you want to test the hypothesis
%                 that the 3rd and 5th elements of theta
%                 are equal and that the 1st, 6th, and 7th
%                 elements are also equal
%                 waldtest(theta,[3,5,0;1,6,7],[],vc)
%                 note that the 0 on the first line is simply ignored
%
%---------------------------------------------------
% RETURNS: waldstat = wald statistic 
%          waldprb  = marginal probability for waldstat
% NOTE:  large waldstat => reject the restrictions as inconsisent 
%                       with the data
% NOTE: requires chis_prb from Jim Lesage's Econometrics toolbox
%---------------------------------------------------

% written by:
% Guillaume Frechette
% Ohio State University
% Department of Economics
% 410 Arps Hall
% 1945 North High Street
% Columbus, OH 43210-1172
% frechette.6@osu.edu

pflag = 0;
if nargout == 0
   pflag = 1;
end
if nargin ~= 4 % flag incorrect arguments
   error('waldtest: Wrong # of input arguments');
end;
% flag if theta is 1xn
if length(theta)==1,error('waldtest: theta should be transposed');end;
[dof,k]=size(R);
if isempty(q)==0
   if k == length(theta) 
      % flag if q has the wrong dimensions
      if dof~=length(q), error('waldtest: Wrong # of restrictions in R or in q'); end;
   else
      % flag if R has the wrong dimensions
      if k>length(theta), error('waldtest: Wrong # of arguments in R'); end;
      temp=R;
      R=zeros(length(temp),length(theta));
      for i=1:k
         R(i,temp(1,i))=1;
      end
      [dof,k]=size(R);
      q=q*ones(dof,1);
   end
else
   temp=R; R=[];
   for i=1:dof
      rbr=temp(i,:);
      rbr=rbr(rbr~=0);
      [junkr,junkc]=size(rbr);
      [trashr,trashc]=size(R);
      R=[R;zeros(junkc-1,length(theta))];
      R(trashr+1:trashr+junkc-1,rbr(1,1))=1;
      for j=2:junkc
         R(trashr+j-1,rbr(1,j))=-1;
      end
   end
   [dof,k]=size(R);
   q=zeros(dof,1);
end

% compute waldstat
waldstat = (R*theta-q)'*inv(R*vc*R')*(R*theta-q);
waldprb = chis_prb(waldstat,dof);

if pflag == 1
   fprintf(1,'Wald stat        = %16.8f \n',waldstat);
   fprintf(1,'probability      =%16.4f\n',waldprb);
   fprintf(1,'num, dof         =         %4d\n',dof);
end;