function [siglev,res]=arres(w,A,v,k)
%ARRES 	Residuals of fitted AR process.	
%
%	[siglev,res]=ARRES(w,A,v) computes the time series of 
%	residuals
% 
%	  res(k,:)' = v(k+p,:)'- w - A1*v(k+p-1,:)' - ... - Ap*v(k,:)'
%
%	of an AR(p) model with A=[A1 ... Ap]. 
%       
%	Also returned is the significance level siglev of the modified 
%	Li-McLeod portmanteau (LMP) statistic. 
%
% 	Correlation matrices for the LMP statistic are computed 
% 	up to lag k=20, which can be changed to lag k by using
% 	[siglev,res]=ARRES(w,A,v,k).

%	Modified 2-Sep-97
%	Author: Tapio Schneider
%		AOS Program, Princeton University
%		Princeton, NJ 08544-0710
%		tapio@splash.Princeton.EDU

m    = size(v,2);              		% dimension of state vectors
p    = size(A,2)/m;       		% order of process
n    = length(v);       		% number of available states
nres = n-p; 				% number of residuals

% Default value for k 
if (nargin < 4) 
   k = 20;
end
if (k <= p) 				% check if k is in valid range 
   error('Maximum lag of residual correlation matrices too small.'); 
end
if (k >= nres) 
   error('Maximum lag of residual correlation matrices too big.'); 
end

w = w(:)';     				% force w to be row vector

% Get time series of residuals 
l = 1:nres; 				% vectorized loop l=1,...,nres 
    res(l,:) = v(l+p,:) - ones(nres,1)*w;
    for j=1:p
        res(l,:) = res(l,:)-v(l-j+p,:)*real(A(:,(j-1)*m+1:j*m)');
    end
% end of loop over l

% Correct residuals by subtraction of mean 
res = res-ones(nres,1)*mean(res);

% Compute lag zero correlation matrix
c0 = res'*res;
d  = diag(c0);
dd = sqrt(d*d');
c0 = c0./dd;

% Get "covariance matrix" in LMP statistic
c0_inv  = inv(c0); 			% inverse of lag 0 correlation matrix
rr = kron(c0_inv, c0_inv);       	% "covariance matrix" in LMP statistic

% Initialize LMP statistic and correlation matrices
lmp = 0;                         	% LMP statistic
cl  = zeros(m,m);                 	% correlation matrix
x   = zeros(m*m,1);                	% correlation matrix arranged as vector

% Compute modified Li-McLeod portmanteau statistic
for l=1:k
   cl  = (res(1:nres-l,:)'*res(l+1:nres,:))./dd;  % lag l correlation matrix
   x   = reshape(cl,m*m,1);        	% arrange cl as vector by stacking columns
   lmp = lmp+x'*rr*x;           	% sum up LMP statistic
end
lmp = n*lmp+m^2*k*(k+1)/2/n;     	% add remaining term and scale

dof_lmp = m^2*(k-p);                 	% degrees of freedom for LMP statistic

% Significance level with which hypothesis of uncorrelatedness is rejected
siglev = 1-gammainc(lmp/2,dof_lmp/2); 	








