Dear SPM list,
I'm trying to estimate the autocorrelation functions for my data using
MATLAB and the methods described in "To Smooth or Not to Smooth?" (2000).
I'm running into trouble defining the Va matrix. I'd like to use an
AR(16) or AR(32) model, but I'm unclear about the notation in Appendix A.
I assumed that rho(t) is the autocorrelation function, which is what I
thought would be given by xcorr(timeseries). If that were so, then Va =
toeplitz(xcorr(ts)), and there's no need for the AR model. So rho(t) must
be something more complex. I did:
if xc is xcorr(timeseries)...
xc = xc/max(xc);
xc = xc(200:end);
V = toeplitz(xc); % A.1
~ ~ ~ ~
Another option is to estimate the 'autoregression coefficients', a, using
eq. A.3. In this equation, I take the inverse of the autocorrelation
matrix, i think (?), and 'whiten' or deconvolve the autoregressive
function, it looks like. I'm assuming that rho(i) are the elements of
xcorr(timeseries), which gives me the autoregressive function.
I can compute this, if M is toeplitz(xc), with:
a = inv(M)*xc;
a = inv(M)*[xc(2:end); 0]; % A.3
...but I'm not sure which. is rho(1) the first element of xc? If so,
this should be 1, as the autocorrelation at time 0 = rho(1) = 1. Since
the diagonals of the V matrix are 1, and the first off-diagonal is rho(1),
I'm a little confused about exactly what to use. maybe rho(1) is the
autocorrelation function at timeshift = 1? If so, I need to add an
element to the xc vector to make it long enough to multiply; thus the
second of these two options. One clarifying question: what is matrix M
called in this case, if it's not the autocorrelation matrix V?
If I use the elements of a to make a triangular matrix, as in equation
A.2, and then compute the K and V matrices from that,
A = zeros(size(ts,1),size(ts,1));
for i = 1:size(ts,1)-1
A = A + diag(ones(200-i,1)*a(i),i);
end
K = inv(1 - A);
V = K * K';
In this case, the diagonals of V do not equal one, and the matrix doesn't
look much like what might be expected from an MRI autocorrelation
function, so I must be doing something wrong.
~ ~ ~ ~
In neither of these potential methods have I used the Yule-Walker method,
which I thought was necessary. I can estimate the AR model using:
coeff = aryule(timeseries,16) for an AR(16) model
but this gives me coefficients for the polynomial, and I'm not sure how to
convert those into an autocorrelation matrix V.
In summary, how can I compute Va?
I know this message may be very silly and confusing, but any help on this
would be gratefully appreciated!
Thanks,
Tor
_____________________________
Tor Wager
Department of Psychology
University of Michigan
Cognition and Perception Area
525 East University
Ann Arbor, MI 48109-1109
Office: 734-936-1295
Home: 734-995-8975
Email: [log in to unmask]
_____________________________
|