function quantile_plot_resids(y,X,sims);

% QUANTILE_PLOT_RESIDS does a quantile plot of regression residuals with 
% simulated 95% confidence bounds on the plot.

% INPUTS:      y = vector of dependent variables
%              X = explanatory data for the regressions
%           sims = (optional) number of simulations to be used in confidence bounds
%                  Default: = 1000
%
% Written by: T.E.SMITH, 3/16/03              
%


if nargin < 3
   
   sims = 1000;
   
end

n = length(y);

X = [ones(n,1),X]; %add intercept

k = size(X,2);

H = X*inv(X'*X)*X';

M = eye(n) - H;

res = M*y; 

s = sqrt((1/(n - k))*y'*M*y);

d = diag(M)';

dd = sqrt(d);

%calculate studentized residuals (Sen, p.157)

ss = zeros(1,n);

for i = 1:n
   
   ss(i) = sqrt( ( (n - k)*(s^2) - (res(i)^2)/d(i) )/(n - k - 1) );
   
end


res = res./(ss.*dd)';

res = sort(res);

%Now simulated studentized residuals

RES = M*normrnd(0,1,n,sims);

RES = RES';

for i = 1:sims
   
   s = sqrt((1/(n - k))*RES(i,:)*M*RES(i,:)');
   
   for j = 1:n
   
   ss(j) = sqrt( ( (n - k)*(s^2) - (RES(i,j)^2)/d(j) )/(n - k - 1) );
   
   end

   RES(i,:) = sort(RES(i,:)./(ss.*dd));
 
end

lower = prctile(RES,2.5);

mid = prctile(RES,50);

upper = prctile(RES,97.5);


x = ([1:n] - (3/8))/(n + .25);

%x = ([1:n] - (1/2))/n;

z = norminv(x);

for i = 1:n
   
   SE(i) = sqrt(x(i)*(1 - x(i))/n)/normpdf(z(i));
   
end

UPPER = z + 1.96*SE;

LOWER = z - 1.96*SE;

plot(z,res,'k.');

axis image

hold on;

plot(z,upper,'r');

hold on;

plot(z,lower,'r');

hold on;

plot(z,UPPER,'k');

hold on;

plot(z,LOWER,'k');

hold on;

plot(z,z);