Simulating Stochastic Processes the MMA way

Question

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

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.