function pacorr(x,n,name);
% pacorr(x, n, cor, name);
% Calculates and plots the serial auocorrelation correlation function
% and the partial autocorrelation function of x, for
% up to n lags. n should be < N/4, where N is the length
% of x.  The unbiased estimate of r and  the progressive 
% Bonferroni correction of alpha are used.
% Written by Gerry Middleton, November 1996 & modified by
% J. Bradbury Sept 2000 and S. Vehrencamp Nov 2000.
P=[0.5 0.2 0.1 0.05 0.02 0.01 0.005 0.002 0.001]';
T=[0.6745 1.2816 1.6449 1.96 2.3263 2.5758 2.807 3.0902 3.2905]';
NN = length(x);
xbar = mean(x);
xd = x - xbar;
ss = xd'*xd;
for i = 1:n+1
   xdlag = x(1:NN+1-i) - xbar;
   r(i) = xd(i:NN)'*xdlag / ss;
end
rr=r(2:n+1);
se=sqrt(1+2*rr*rr')/NN;
CI05=1.96*se;
CI01=2.58*se;
for i=1:(n+1)
	d(i,:)=[i r(i) 1 0];
end
d=sortrows(d,2);
d=flipud(d);
for i=2:(n+1)
	if (abs(d(i,2))>CI05)
			d(i,4)=1;
	end
end
d=sortrows(d,1);
for i=1:(n+1)
	d(i,:)=[i r(i) 1 0];
end
d=sortrows(d,2);
d=flipud(d);
for i=2:(n+1)
	pr=0.05/(i-1);
	tt=interp1(P,T,pr);
	d(i,3)=tt*se;
	if (abs(d(i,2))>d(i,3))
		d(i,4)=1;
	end
end
d=sortrows(d,1);
lag = [0:n];
figure
plot(lag, r,'linewidth',1);
hold on
plot([lag(1) lag(n+1)],[0 0],'k:');
dx=0.02*(max(lag)-min(lag));
dy=0.02*(max(r)-min(r));
text((lag(1)+dx),(1+(2*dy)),name);
xlabel('Lag size in samples');
title('Serial correlation function with progressive Bonferroni');
ylabel('Autocorrelation');
for i=1:n
	if (d(i,4)==1)
		plot(lag(i),r(i),'k*');
	end
end
hold off

rg=cumsum(ones(n,1));
r=[rg r(2:(n+1))'];

rp = zeros(1,n);
rp(1) = 1.0;    % rp(1) is the (partial) correlation at lag 0   
rp(2) = r(2);   % rp(2) is the (partial) corelation at lag 1
for j = 2:n     % these are the true partial correlations
   R = zeros(j); 
   ri = r(1:j);
   R(1,:) = ri;
   for i = 2:j
      ri = [r(i) ri(1:j-1)];
      R(i,:) = ri;
   end
   phi = R\r(2:j+1)';
   rp(j+1) = phi(j);
end
lag = [0:n];
figure
plot(lag, rp,'linewidth',1);
hold on
dx=0.02*(max(lag)-min(lag));
dy=0.02*(max(rp)-min(rp));
y = 2*sqrt(1/NN);
text((lag(1)+dx),(max(rp)+(2*dy)),name);
plot([0,n],[-y,-y],':',[0,n],[y,y],':');
title('Partial autocorrelation');
xlabel('Lag size in samples');
ylabel('Partial r');
hold off
return