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)';

%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

d = sqrt(d);

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

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.*d));
 
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;

phi = norminv(x);

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

axis image

hold on;

plot(phi,upper,'r');

hold on;

plot(phi,lower,'r');

hold on;

plot(phi,phi);


%x = [ones(n,1),phi'];

%b = inv(x'*x)*x'*mid';

%MID = x*b;

%plot(phi',MID);