help-octave
[Top][All Lists]
Advanced

[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
-------------------------------------------------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]