[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
RE: generate AR(1) data
From: |
Ted Harding |
Subject: |
RE: generate AR(1) data |
Date: |
Thu, 17 Jul 2003 08:37:04 +0100 (BST) |
On 17-Jul-03 Heber Farnsworth wrote:
> I need to generate a large number of random (Gaussian) vectors that
> have an AR(1) correlation structure. Obviously I could do this in a
> loop but this is slow. I would rather rotate a vector of IID normal
> draws so that the result will have the right covariance matrix. So
> what I need is the matrix square root of a symmetric toeplitz matrix
>
> 1 .5 .25 .125 ...
> .5 1 .5 .25 ...
> .25 .5 1 .5 ...
>
> etc (here the correlation is .5). It seems inefficient to try to
> compute this square root numerically. Isn't it known in closed form?
You could try the Choleski factorisation of T into T=V'*V where V is
upper triangular. In this case the top row of V is the same as that
of T; each of the remaining rows has zeros on the left of the main
diagonal, the on-diagonal term is sqrt(3/4), and subsequent terms are
all 1/2 of the one on their left.
In general, where your original T has top row 1, r, r^2, r^3, ... ,
the top row of V is as in T, and (reading rightwards from the
diagonal term) each subsequent row is of the form a*(1, r, r^2, r^3, ... )
where a^2 = 1 - r^2. Terms to the left of the diagonal are zero. If you
work it out (using a^2 = 1 - r^2) you can easily see that
for i <= j, [V'*V](i,j) = r^(j-i) and of course V'*V is symmetric.
Once you have V, let U be a ROW vector of independent N(0,1)s.
The covariance matrix of X = U*V is E(X'*X) = V'*E(U'*U)*V = V'*V = T.
Hope this helps!
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <address@hidden>
Fax-to-email: +44 (0)870 167 1972
Date: 17-Jul-03 Time: 08:37:04
------------------------------ XFMail ------------------------------
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------