Fast Simulations with Compile

this post relates to another post that I didn't follow up propely. If I wanted to simulate a system of stochastic proesses like the following, and loop over this run many many times would writing the processes as 'compiled' pure functions speed up the run time? Or, is Nestlist already trying to do this for me?

norTheta[mu_, sigma_] := Random[NormalDistribution[mu, sigma]];
norPi[mu_, sigma_] := Random[NormalDistribution[mu, sigma]];

thetaNext[thetaNow_] :=
thetaNow + (-lambdaTheta*(thetaNow - thetaBar)*deltaT
+ sigmaTheta*norTheta[0, 1]*Sqrt[deltaT]);
piNext[piNow_, thetaNow_] := piNow + (-lambdaPi*(piNow - thetaNow)*deltaT +
sigmaPi*norPi[0, 1]*Sqrt[deltaT]);

lambdaTheta = 0.07; sigmaTheta = 1.2; thetaBar = 2; lambdaPi = 1.0;
sigmaPi = 1.25; deltaT = 1/12;
steps = 252;
T = 5;
deltaT = 1/steps; // N
Maturity = T*steps;

simulateRun = Transpose[NestList[{piNext[#[[1]], #[[2]]],
thetaNext[#[[2]]]} &, {2, 2}, Maturity]];

• Check simulateRun[[1, 3]] after evaluating the above code. There is something wrong in your function call thetaNext[#[[1]], #[[2]]]. You have defined it as a single argument function before hand. Commented Nov 14, 2012 at 11:18
• sorry, yes, I hadn't actually run the code. I've fixed and edited this above. Commented Nov 14, 2012 at 12:26
• If it is possible to generate all your random data in advance and then access it with an incrementing index (list[[i++]]) that should be somewhat faster, for what it's worth. Commented Nov 14, 2012 at 12:45
• ok - thanks, that is possible. Commented Nov 14, 2012 at 13:18
• I added an answer that as well as another optimization. Please let me know if these are applicable to your real problem, or if you have trouble understanding or implementing them. Commented Nov 14, 2012 at 13:26

You are trying to implement Euler-Maruyama simulation method for a 2-stage short-term interest rate model which is given by the following system of SDEs: $$\begin{eqnarray} \mathrm{d} \theta_t &=& -\lambda_\theta \left( \theta_t - \bar\theta\right) \mathrm{d}t + \sigma_\theta \mathrm{d}W_{\theta,t} \\ \mathrm{d} \pi_t &=& -\lambda_\pi\left( \theta_t - \pi_t \right) \mathrm{d}t + \sigma_\pi \mathrm{d} W_{\pi,t} \end{eqnarray}$$ where $W_\theta$ and $W_\pi$ are independent standard Wiener processes.

Here is the compiled code implementing the above.

cfEM = Compile[{{lambdaTheta, _Real}, {thetaBar, _Real}, {sigmaTheta, \
_Real}, {lambdaPi, _Real}, {sigmaPi, _Real}, {th0, _Real}, {pi0, \
_Real}, {dt, _Real}, {steps, _Integer}},
Module[{zs, bag, thc, pic, zths, zpis},
zths = RandomReal[NormalDistribution[0, Sqrt[dt]], steps];
zpis = RandomReal[NormalDistribution[0, Sqrt[dt]], steps];
thc = th0; pic = pi0;
bag = InternalBag[{thc, pic}];
Do[
pic += -lambdaPi dt (pic - thc) + sigmaPi zpis[[k]];
thc += -lambdaTheta dt (thc - thetaBar) + sigmaTheta zths[[k]];
InternalStuffBag[bag, {thc, pic}, 1];
, {k, 1, steps}];
Partition[InternalBagPart[bag, All], 2]
]
];


Call example:

In[14]:= AbsoluteTiming[
Length[data = cfEM[0.07, 2., 1.2, 1, 1.25, 2., 2., 1/252., 5 252]]]

Out[14]= {0., 1261}


Note, however, that Euler-Maruyama method is only approximate, while your system of SDEs is Gaussian. An exact simulation method, known as Ozaki method, is available and is not hard to code. See this book for the details of the 1D case. You would need to generalize it to the 2D case, but it is not very hard.

• ...Or I could just wait for Sasha to come along and make my answer look silly. :^) Commented Nov 14, 2012 at 13:44
• Wow, great answer. Well, actually its an inflation process with a stochastic target, but yes it is also a commonly deployed short rate model. I 'm trying to migrate to MMA --- this gives me good encouragement. My actual system is a bit more funky than this example but I shall check out the Ozaki method - thank a lot Commented Nov 14, 2012 at 14:25
• So for those who couldn't buy the book, what does the Ozaki method achieve? Thanks :) Commented Jan 5, 2016 at 2:43
• It's exact on SDEs with linear drift and constant diffusion, where Euler-Maruyama or Milstein are only approximate in the strong sense. Commented Jan 5, 2016 at 2:45

There are two areas for optimization that I see here.

The first, if possible, is to generate all your random data in advance and then access it with an incrementing index, e.g. list[[i++]].

The second is to partially evaluate the definitions of thetaNext and piNext for a given set of parameters.

A note: Random has been deprecated for some time now and may produce inferior results. You should be using RandomReal/RandomInteger in version 7 or RandomVariate in version 8.

Update

Here is a cleaner implementation of my recommendations. Should this be inapplicable my original code is visible in the edit history of this post.

thetaNext[thetaNow_] :=
thetaNow + (-lambdaTheta*(thetaNow - thetaBar)*deltaT +
sigmaTheta*norTheta[0, 1]*Sqrt[deltaT]);

piNext[piNow_, thetaNow_] :=
piNow + (-lambdaPi*(piNow - thetaNow)*deltaT +
sigmaPi*norPi[0, 1]*Sqrt[deltaT]);

lambdaTheta = 0.07; sigmaTheta = 1.2; thetaBar = 2; lambdaPi = 1.0;
sigmaPi = 1.25; deltaT = 1/12;
steps = 15000;
T = 5;
deltaT = 1/steps; // N
Maturity = T*steps;


We can reduce the core of your NestList function (for these specific parameters) as follows:

func =
Block[{norTheta, norPi, Part},
Function @@ FullSimplify @ {
{piNext[#[[1]], #[[2]]], thetaNext[#[[2]]]} /.
{norPi[0, 1] -> #2[[1]], norTheta[0, 1] -> #2[[2]]}
}
]

{0.916667 #1[[1]] + 0.0833333 #1[[2]] + 0.360844 #2[[1]], 0.0116667 + 0.994167 #1[[2]] + 0.34641 #2[[2]]} &


We will use the second argument of this function to insert the random data, which we create with:

piList = RandomReal[NormalDistribution[0, 1], Maturity];
thetaList = RandomReal[NormalDistribution[0, 1], Maturity];
rands = {piList, thetaList}\[Transpose];


(A shorter form exists for this specific data but I am trying to retain some generality.)

And which we use FoldList to provide:

FoldList[func, {2, 2}, rands]\[Transpose] // Timing // First
`

0.219

By comparison your code takes 4.368 seconds on my machine with the same parameters.

• Thanks a lot. Both answers provide excellent solutions. Thanks Mr Wizard Commented Nov 14, 2012 at 13:52
• @LuapNalehw (Paul Whelan?) I simplified my method somewhat. Commented Nov 14, 2012 at 14:25
• friends call me Luap Commented Nov 14, 2012 at 15:23
• @Luap -- I'm Paul too, which is why I recognized that. :-) Commented Nov 15, 2012 at 1:48
• @Mr.Wizard you can define rands in one step using RandomReal[NormalDistribution[0, 1], {Maturity, 2}] Commented Nov 15, 2012 at 13:28