This question is really a specific problem and a methodological one concerning MMA best practices. I want to simulate a system of stochastic processes. If this were a geometric Brownian motion or random walk then the recursive (univariate) nature of the problem means FoldList or Nestlist is easy and neat to implement.

However, what if I wanted to simulate an mean-reverting process which itself had a mean-reverting long run mean? Now, I cannot work out the MMA efficient way of doing it other than within a For loop.

For example, the following code simulates an OU process for inflation with a stochastic central tendency :


norθ[mu_, sigma_] := Random[NormalDistribution[mu, sigma]]; 
norπ[mu_, sigma_] := Random[NormalDistribution[mu, sigma]]; 

deltaθt[θnow_] :=-λθ*(θnow-θbar)*deltaT + σθ*norθ[0, 1]*Sqrt[deltaT]    
deltaπt[πnow_] := -λπ*(πnow - θnow)*deltaT + σπ*norπ[0, 1]*Sqrt[deltaT]

λθ     = 0.07;
σθ     = 1.2;
θbar   = 2;
θnow   = 2;
λπ     = 1.;
σπ     = 1.25;
πnow   = 2;
deltaT = 1/12;
noYear = 100*12;

Process = Reap[For[i = 1, i < noYear, i++, 
    Sow[{i, πnow, θnow}]; 
    πnow = deltaπt[πnow] + πnow; 
    θnow = deltaθt[θnow] + θnow; 
 ]][[2, 1]]; 

inflation = Process[[All, {1, 2}]]; 
target    = Process[[All, {1, 3}]];

ListLinePlot[{inflation, target}]
  • $\begingroup$ Welcome to Mathematica.SE, Luap! (Nicely obsfucated username there ;-) ) Have you looked at FoldList? Also note that you can use RandomVariate to get a list of numbers distributed according to a desired distribution. $\endgroup$
    – Verbeia
    Commented Oct 3, 2012 at 22:42
  • $\begingroup$ I wonder if you could do this with a MixtureDistribution[] see: reference.wolfram.com/mathematica/ref/MixtureDistribution.html $\endgroup$
    – Jagra
    Commented Oct 3, 2012 at 23:13
  • $\begingroup$ Happy to be hear. I divorced Matlab and Gauss recently through boredom. They'll get over it. RandomVariate is certainty useful if you know / can work out the distribution. But in more complicated examples it might not be so easy. $\endgroup$ Commented Oct 3, 2012 at 23:45

1 Answer 1


For the example you give there is no reason you can't use NestList, you just need to make two simple changes:

  1. Don't use the side effect in deltaπt to get the value for θnow, give it as an explicit second argument
  2. Then you just do:

    NestList[{deltaπt[#[[1]], #[[2]]] + #[[1]], deltaθt[#[[2]]] + #[[2]]} &, {2, 2}, noYear]

Similar changes would allow for use of FoldList to pass in the random values instead of using the function calls.

Not being an expert in Brownian motion, I can't tell you if this covers all examples, but I can't see why there is any recursive stochastic system you can't model using NestList/FoldList.

Good luck.

  • 2
    $\begingroup$ Also if you are worried about performance be careful with using fractions of integers as parameters, since mathematica will deal with these symbolically, ie make deltaT=1.0/12.0 and you should notice a speed bump $\endgroup$
    – Gabriel
    Commented Oct 3, 2012 at 23:42
  • $\begingroup$ Thanks a lot, great fractions of integers tip also. $\endgroup$ Commented Oct 4, 2012 at 8:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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