returns an (n x n) autocorrelation matrix for an AR(p) process FORMAT [Q] = pr_spm_q(A,n) A - vector pf p AR coeficients n - size of Q __________________________________________________________________________ spm_Q uses a Yule-Walker device to compute K where: y = K*z such that y is an AR(n) process generated from an i.i.d innovation z. This means cov(y) = <K*z*z'*K> = K*K' Critically, this is not the correlation because if cov(z) = eye(n) then trace(cov(y)) ~= n. This is why the normalization is required corr(y) = D*K*K'*D'; The reason the diagonals of corr(y) are not constant is that we are modeling finite length AR sequences, which incur boundary effects at the beginning and end of the sequence. __________________________________________________________________________ Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0001 function [Q] = pr_spm_q(A,n) 0002 % returns an (n x n) autocorrelation matrix for an AR(p) process 0003 % FORMAT [Q] = pr_spm_q(A,n) 0004 % 0005 % A - vector pf p AR coeficients 0006 % n - size of Q 0007 %__________________________________________________________________________ 0008 % spm_Q uses a Yule-Walker device to compute K where: 0009 % 0010 % y = K*z 0011 % 0012 % such that y is an AR(n) process generated from an i.i.d innovation 0013 % z. This means 0014 % 0015 % cov(y) = <K*z*z'*K> = K*K' 0016 % 0017 % Critically, this is not the correlation because if cov(z) = eye(n) 0018 % then trace(cov(y)) ~= n. This is why the normalization is required 0019 % 0020 % corr(y) = D*K*K'*D'; 0021 % 0022 % The reason the diagonals of corr(y) are not constant is that we 0023 % are modeling finite length AR sequences, which incur boundary effects 0024 % at the beginning and end of the sequence. 0025 %__________________________________________________________________________ 0026 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience 0027 0028 % Karl Friston 0029 % $Id: spm_Q.m 372 2005-12-08 17:12:13Z karl $ 0030 0031 0032 % compute Q 0033 %-------------------------------------------------------------------------- 0034 p = length(A); 0035 A = [1 -A(:)']; 0036 K = inv(spdiags(ones(n,1)*A,-[0:p],n,n)); 0037 K = K.*(abs(K) > 1e-4); 0038 Q = K*K'; 0039 D = spdiags(sqrt(1./diag(Q)),0,n,n); 0040 Q = D*Q*D; 0041 Q = Q.*(abs(Q) > 1e-4);