Fast Simulations with Compile

Question

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

Answer 1

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 = Internal`Bag[{thc, pic}];
    Do[
     pic += -lambdaPi dt (pic - thc) + sigmaPi zpis[[k]];
     thc += -lambdaTheta dt (thc - thetaBar) + sigmaTheta zths[[k]];
     Internal`StuffBag[bag, {thc, pic}, 1];
     , {k, 1, steps}];
    Partition[Internal`BagPart[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.

Answer 2

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.

Given your definitions and parameters:

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.