Take the 2-minute tour ×

Find $X_t = c_1 X_{t-1} + \dots + c_n X_{t-n} + n_t$ for given $c_i$, initial conditions $(X_1, \dots, X_n)$, and distribution for i.i.d. $n_t$.

I would like to know if there is a more efficient or elegant way to code this than using a for loop.

share|improve this question
I think something like my answer here is what you want. This is the same (or nearly the same) as that question, but yours is much clearer. I guess I'll hold off on voting to close... –  The Toad May 17 '12 at 1:09
All the answers were good. I learned something from you all; I'm glad I asked. Yes, I meant "generate" when I said "find". –  Emre May 17 '12 at 5:20

3 Answers 3

up vote 7 down vote accepted

Generating an autoregressive time series is not the same as "finding" it, since you have the parameters. If I understand you correctly, this is an application made for FoldList:

Here is a simple AR(1) process to illustrate the technique:

FoldList[c1 #1 + #2&, x0, RandomVariate[NormalDistribution[0,1],{100}]

For example

 FoldList[0.9 #1 + #2 &, 0.2, RandomVariate[NormalDistribution[0, 1], {100}]]

A simulated AR(1) process

For an AR($p$) process with $p>1$, something like this would work:

FoldList[Join[{params.#1+#2},Most[#1]]&, startingvalues, 

Where params is your list of c parameters and startingvalues is your list of starting values.

Multivariate processes (i.e. simulation of data from a VAR to get a vector of length $p$ at each time period) can be handled similarly except of course that params will be a matrix and you'll be saving a matrix consisting of the vector of data for the current period, plus the $p$ lags, at each step. Then you need to extract the matrix of data from the repetitions representing the lags at each step using Part, i.e.

FoldList[Join[{matrixofparams.#1+#2},Most[#1]]&, matrixofstartingvalues, 
   RandomVariate[NormalDistribution[0, 1], {100,p}]][[All, 1, All ]]
share|improve this answer
Do you think "finding" means something else or that he just expressed himself wrong? If it means something else, what would it be? –  Rojo May 17 '12 at 4:29
I would have assumed "finding" meant "estimating the parameters", but he has the parameters and the title says "generating", so I interpret him as meaning "simulate". I use the FoldList approach to generating an AR process all the time, to show people how functional programming in Mathematica works, and to generate fake economic-looking data for my test graphs. –  Verbeia May 17 '12 at 4:34
Yeah, it's a practical one-liner, +1 –  Rojo May 17 '12 at 4:42

A straightforward way would be using recursion and memoization. An example:

n = 5;
c = RandomReal[NormalDistribution[], n]/100;

Array[(x[#] = RandomReal[NormalDistribution[]]) &, n]; (* Initial conditions *)
x[t_Integer] /; t > n := x[t] = c.Table[x[i], {i, t - 1, t - n, -1}] +

ListLinePlot[x /@ Range[300]]

enter image description here

share|improve this answer

I like @RM 's approach. Another one, with a more "mathematical" notation (but I must say I'm not certain of the random behaviour of the way I'm getting the numbers) could be the following

First we create the discrete noise

i : noise["Seed"] := i = RandomInteger[{-# , #}] &[Developer`$MaxMachineInteger];
noise[n_Integer] := 
 BlockRandom[SeedRandom[# + n]; 
    RandomVariate[NormalDistribution[]]] &[noise["Seed"]]

Now, noise is a realisation of the process. You can realise another one by resetting the seed with noise["Seed"]=.;

If you do DiscretePlot[noise[n], {n, 0, 10}] several times, you'll see what I mean.

Now, just use RecurrenceTable

    x[0] == 0.5, 
    x[1] == 0.54, 
    x[2] == -2.3,
    x[n] == noise[n] + 0.75 x[n - 1] - 0.23 x[n - 2] + 0.2 x[n - 3]}, 
        x[n], {n, 0, 10}]
share|improve this answer

Your Answer


