function [par,sig]=pacf(y,nr)

%PACF   Partial autocorrelation function
%       PACF computes partial autocorrelation function
%       of a time series and its standard deviation:
%       [C2 SIG] = PACF(Y,NR) -- C2 is a partial autocor.
%       function, SIG is its standard deviation.
%       ACFP(Y,NR) for column vector Y representing se-
%       ries, NR being the number of lags (cca 25)
%       (if not present, it's computed for all lags)

n=length(y);
if nargin==1
   nr=n;
end;

%Computing autocorrelation function
ar=acf(y,nr);
ar(nr+1)=NaN;

%Computing partial autocorrelation function
%Initial value of par an P matrix
par(1,1)=ar(1);
P=[1 ar(1); ar(1) 1];
for k=2:nr
   %Updating Ps matrix
   Ps=P;
   Ps(:,k)=ar(1:k);
   %Computing next par value
   par(k,1)=det(Ps)/det(P);
   %Updating P matrix
   P(k+1,1:k)=fliplr(ar(1:k)');
   P(1:k,k+1)=flipud(ar(1:k));
   P(k+1,k+1)=1;
end

%Computing its standard deviation
sig=ones(size(par))*sqrt(1/n);
