[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
RE: generate AR(1) data
From: |
Mike Miller |
Subject: |
RE: generate AR(1) data |
Date: |
Thu, 17 Jul 2003 09:38:24 -0500 (CDT) |
If it wasn't clear -- Ted and I are saying the same thing about how to use
Cholesky factorization but we're using slightly different notation. All
I'm adding is that there is a chol() function in Octave. --Mike
On Thu, 17 Jul 2003 address@hidden wrote:
> 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.
-------------------------------------------------------------
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
-------------------------------------------------------------