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.

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}]
share|improve this question
    
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
add comment

1 Answer 1

up vote 12 down vote accepted

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.

share|improve this answer
2  
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
add comment

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.