% Function to calculate the distribution of R-squared.
% Assumes regressors are non-stochastic (can be modified as in Reiss
% (2016)).
% Implemented as rsq(arg vals,X matrix,coef,error var)
% Requires:
%   rvec = points at which to evaluate the r-squared cdf
%   X = Regressor matrix (excluding constant term; (T-K-1)/2 integer)
%   coef = (Population) Coef vector (excluding constant term)
%   sigsq_e = Error variance
% Calculated
%   delt = noncentrality parameter (Here only numerator delt_1 > 0 assumed.)
% 
function rsqcdf=rsq(r_vec,X,coef,sigsq_e)
[T, K]=size(X); a=K/2; b=(T-K-1)/2; c=(T-1)/2; 
isb_int=b-round(b);
if isb_int~=0
    'b not integer'
    return;
end
S_X=(X-repmat(mean(X),T,1))'*(X-repmat(mean(X),T,1));
delt=coef'*S_X*coef/sigsq_e
  fac=exp(-delt*(1-r_vec)/2).*gamma(c).*(r_vec.^(a)).*(1-r_vec).^(b-1);
  % fac2 tracks double sum in eqn (27)
  fac2=0;
  % Double sum
  for l=0:b-1;
    for j=0:b-l-1;
        fac2=fac2+((r_vec./(1-r_vec)).^l).*(delt*r_vec./2).^j./...
                (gamma(b-l-j)*gamma(a+l+j+1)*gamma(j+1));
    end;
  end;
% Putting the sums together with the leading terms in (27) ...
rsqcdf=[r_vec(:) (fac(:).*fac2(:))];
