# 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. – PlatoManiac Nov 14 '12 at 11:18
sorry, yes, I hadn't actually run the code. I've fixed and edited this above. – Luap Nalehw Nov 14 '12 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. – Mr.Wizard Nov 14 '12 at 12:45
ok - thanks, that is possible. – Luap Nalehw Nov 14 '12 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. – Mr.Wizard Nov 14 '12 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. :^) – Mr.Wizard Nov 14 '12 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 – Luap Nalehw Nov 14 '12 at 14:25
So for those who couldn't buy the book, what does the Ozaki method achieve? Thanks :) – AimForClarity Jan 5 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. – Sasha Jan 5 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 – Luap Nalehw Nov 14 '12 at 13:52
@LuapNalehw (Paul Whelan?) I simplified my method somewhat. – Mr.Wizard Nov 14 '12 at 14:25
friends call me Luap – Luap Nalehw Nov 14 '12 at 15:23
@Luap -- I'm Paul too, which is why I recognized that. :-) – Mr.Wizard Nov 15 '12 at 1:48
@Mr.Wizard you can define rands in one step using RandomReal[NormalDistribution[0, 1], {Maturity, 2}] – faysou Nov 15 '12 at 13:28