Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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... –  rm -rf 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}]]
 ListLinePlot[%]

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, 
  RandomVariate[NormalDistribution[0, 1], {100}]]

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

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

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

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

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

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

ListLinePlot[x /@ Range[300]]

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.