# Simulating Stochastic Processes the MMA way

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 :

Clear["Global*"];

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}]

-
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. –  Verbeia Oct 3 '12 at 22:42
I wonder if you could do this with a MixtureDistribution[] see: reference.wolfram.com/mathematica/ref/MixtureDistribution.html –  Jagra Oct 3 '12 at 23:13
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. –  Luap Nalehw Oct 3 '12 at 23:45

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.

-
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 –  Gabriel Oct 3 '12 at 23:42
Thanks a lot, great fractions of integers tip also. –  Luap Nalehw Oct 4 '12 at 8:11