Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I would like to generate a set of forward yield curve matrix of size 1000 x 100. First I defined my SparseArray of 1000 x100:

(forwardYieldCurve=Normal[SparseArray[{{1,1}->0,{1000,100}-> 0}]])//MatrixForm;  

then initial first row of forwardYieldCurve using:


then for the second row and iterate with i (in BOLD) up to 100 rows of the forwardYieldCurve matrix:


takes around 4 minutes to do to obtain results of 100 x 100, which will be projected take it to around 40 minutes to run this single set of simulations. When I eventually set i to iterate up to 1000 x 100, and furthermore, I will repeat this many times to get a statistically monte-carlo simulation of distributions. How to optimize this to reduce run time.

My input data dimensions:

forwardYieldCurve -> {1000,100}  
driftM = {{4.29874*10^-6,8.59748*10^-6,0.0000128962,0.000017195,0.0000214937,0.0000257924,0.0000300912,0.0000343899,0.0000386887,0.0000429874,0.0000472861,0.0000515849,0.0000558836,0.0000601824,0.0000644811,0.0000687798,0.0000730786,0.0000773773,0.000081676,0.0000859748,0.0000902735,0.0000945723,0.000098871,0.00010317,0.000107468,0.000111767,0.000116066,0.000120365,0.000124663,0.000128962,0.000133261,0.00013756,0.000141858,0.000146157,0.000150456,0.000154755,0.000159053,0.000163352,0.000167651,0.00017195,0.000176248,0.000180547,0.000184846,0.000189144,0.000193443,0.000197742,0.000202041,0.000206339,0.000210638,0.000214937,0.000219236,0.000223534,0.000227833,0.000232132,0.000236431,0.000240729,0.000245028,0.000249327,0.000253626,0.000257924,0.000283717,0.000309509,0.000335302,0.000361094,0.000386886,0.000412679,0.000438471,0.000464264,0.000490056,0.000515849,0.000541641,0.000567433,0.000593226,0.000619018,0.000644811,0.000670603,0.000696396,0.000722188,0.000747981,0.000773773,0.000799565,0.000825358,0.00085115,0.000876943,0.000902735,0.000928528,0.00095432,0.000980113,0.0010059,0.0010317,0.00105749,0.00108328,0.00110907,0.00113487,0.00116066,0.00118645,0.00121224,0.00123804,0.00126383,0.00128962}}  
tstep = 0.01  
randomWalkPCA[n_]:=  RandomVariate[NormalDistribution[0,1],n];  
share|improve this question
+1 for a well written and formatted question. It took a while, but you did it, and that's all that matters! :) – R. M. Jan 14 '13 at 14:37
Thanks, it is making sense... – sebastian c. Jan 14 '13 at 14:40
Try to redefine your RandomWalk function to: randomWalk[x_] := Accumulate[Prepend[RandomVariate[NormalDistribution[0, 1], x], 0]]. So if you need to generate 1000 Random Walks of Length[] 100 try this: ListLinePlot[Table[randomWalk[100], {1000}]]. It takes only 1.8 seconds here... – Rod Jan 14 '13 at 14:56
I am about to go to work so short on time but I think you might have gone off the rails with the opening line. The point of using a sparse array is that it uses less memory for large matrices and runs faster for calculations. Wrapping it in Normal makes it a "normal" matrix which seems to defeat the purpose (...other than making it easier to create a big matrix). If your matrix is truly sparse then try and work up a method that takes advantage of sparse array calculations. – Mike Honeychurch Jan 14 '13 at 20:22
I notice that the lists tenor0, driftm, the three volfits and dateArray are all wrapped in an extra layer of List. IMO your code would be a lot easier to read (and easier to optimise) if you stored 1D lists as 1D lists. – Simon Woods Jan 14 '13 at 23:12

Not a complete solution but a few comments too detailed for a comment.

a) Firstly make use of listability. Coinicdentally mentioned this the other day as well. It is important because listable functions thread themselves athrough lists -- for want of a better description -- and as a rule perform their operations on lists much faster than comparable use of Map or Table

So for example this code fragment:

Table[tenor0[[1, j]] + 
  driftM[[1, j]]*
   tstep + (volFit1[[1, j]]*dX[[1, 1]] + volFit2[[1, j]]*dX[[1, 2]] + 
     volFit3[[1, j]]*dX[[1, 3]])*Sqrt[tstep], {j, 99}]

can be re-written as

tenor0[[1, 1 ;; 99]] + 
 driftM[[1, 1 ;; 99]]*
  tstep + (volFit1[[1, 1 ;; 99]]*dX[[1, 1]] + 
    volFit2[[1, 1 ;; 99]]*dX[[1, 2]] + 
    volFit3[[1, 1 ;; 99]]*dX[[1, 3]]) Sqrt[tstep]

You could also replace 99 in the index with -2.

Also consider this fragment:

Table[(tenor0[[1, j + 1]] - tenor0[[1, j]])/(
 dateArray[[1, j + 1]] - dateArray[[1, j]]), {j, 99}]

this is the same as


...and so on.

b) A major slowdown in the code fragment above seem to be your use of random numbers.


randomWalkPCA[n_] := RandomVariate[NormalDistribution[0, 1], n];
RandVarPCA[mcRun_] := Table[randomWalkPCA[3], {mcRun}];
(dX := RandVarPCA[1000]

So when you use e.g. dX[[1, 3]] you repeatedly regenerate these large number of random numbers only to take {1,3} from that large list.

The more efficient way of running Monte Carlo simulations is to create all your random numbers once only and sample all of them rather than generate a large amount of numbers, a lot of times and only take small sample from them each time. (FWIW I was asked to speed up some MC code a couple of years back with a brief for 10 times improvement and got 250 times primarily with proper handling of random number generation and listability.)

c) I don't think it is necessary or advisable to create a blank matrix which gets filled with values. This is basically procedural thinking and as above it is best to start thinking in terms of entire lists.

d) I am sure there are many other things that can be altered but these are intended to help point you in the right direction rather than being an exhaustive analysis (I am at work so do not have the time)

share|improve this answer
Hi @Mike Honeychurch, could this be the reason why I am getting increased volatility with successive simulation runs as I noted:… – sebastian c. Jan 14 '13 at 23:10
@sebastianc. What I have described should make your code run faster based on an alternative to your current code. By definition it should not change the results you obtain (once identical random numbers are used for comparison, see SeedRandom). If your results are not as you expected then you need to consider your implementation is correct or whether your underlying forumlas are correct. Are you a student or doing this commercially? – Mike Honeychurch Jan 14 '13 at 23:13
student, studying quant finance, but come from engineering background hoping to learn more about this for a successful career transition.. – sebastian c. Jan 14 '13 at 23:17

This overlaps with @Mike Honeychurch' reply. Define the 1xn matrices as simple vectors. For example:

dateArray = {0.0833333, 0.166667, 0.25, 0.333333, 0.416667, 0.5, 0.583333, 0.666667, 0.75, 0.833333, 0.916667, 1., 1.08333, 1.16667, 1.25, 1.33333, 1.41667, 1.5, 1.58333, 1.66667, 1.75, 1.83333, 1.91667, 2., 2.08333, 2.16667, 2.25, 2.33333, 2.41667, 2.5, 2.58333, 2.66667, 2.75, 2.83333, 2.91667, 3., 3.08333, 3.16667, 3.25, 3.33333, 3.41667, 3.5, 3.58333, 3.66667, 3.75, 3.83333, 3.91667, 4., 4.08333, 4.16667, 4.25, 4.33333, 4.41667, 4.5, 4.58333, 4.66667, 4.75, 4.83333, 4.91667, 5., 5.5, 6., 6.5, 7., 7.5, 8., 8.5, 9., 9.5, 10., 10.5, 11., 11.5, 12., 12.5, 13., 13.5, 14., 14.5, 15., 15.5, 16., 16.5, 17., 17.5, 18., 18.5, 19., 19.5, 20., 20.5, 21., 21.5, 22., 22.5, 23., 23.5, 24., 24.5, 25.};

Use Differences and some similar ideas to avoid recomputations.

tdiffs = Differences[tenor0];
ddiffs = Differences[dateArray];
qdiffs = tdiffs/ddiffs;

Also combine the volFitxxx stuff so we can use Dot instead of iterated multiply-and-add.

volFit = Transpose[{volFit1, volFit2, volFit3}];

Most importantly, define dX one time.

Here is a slight recoding of your example. It runs in a split second.

Timing[Module[{dx = dX, sqrtt = Sqrt[tstep]}, 
   Do[forwardYieldCurve[[1, j]] = 
     tenor0[[j]] + driftM[[j]] tstep + (volFit[[j]].dx[[1]]) sqrtt + 
      qdiffs[[j]] tstep, {j, 99}];
   Do[forwardYieldCurve[[i + 1, j]] = forwardYieldCurve[[i, j]] +
      driftM[[j]] tstep + volFit[[j]].dx[[i + 1]]*sqrtt +
      (forwardYieldCurve[[i + 1, j + 1]] - 
          forwardYieldCurve[[i + 1, j]])/ddiffs[[j]] tstep
    , {j, 99}, {i, 100}]

(* {0.100000, Null} *)

--- edit ---

To get different triples in every use of dX, can do as follows.

dX2 := randomWalkPCA[3]

 Module[{sqrtt = Sqrt[tstep]}, 
   Do[forwardYieldCurve[[1, j]] = 
     tenor0[[j]] + driftM[[j]] tstep + (volFit[[j]].dX2) sqrtt + 
      qdiffs[[j]] tstep, {j, 99}];
   Do[forwardYieldCurve[[i + 1, j]] = forwardYieldCurve[[i, j]] +
      driftM[[j]] tstep + volFit[[j]].dX2*sqrtt +
      (forwardYieldCurve[[i + 1, j + 1]] - 
          forwardYieldCurve[[i + 1, j]])/ddiffs[[j]] tstep
    , {j, 99}, {i, 100}]

(* Out[126]= {0.150000, Null} *)

--- end edit ---

share|improve this answer
Hi @Daniel, thanks I will try the above. But don't you normally have to regenerate random numbers each time when simulating a forward yield curve or else how would you get a Gaussian distribution? Or am I mistaken? – sebastian c. Jan 14 '13 at 23:38
You did not seem to be reusing them. More specifically, you were using one (new) row at a time. So you can generate all at once, or generate a row (3, that is) at a time. Now I'm not sure if you want to use the same set of values in each inner loop. If not, then best would be to generate only 3 at a time, not 3000. See edit for details. – Daniel Lichtblau Jan 14 '13 at 23:44
Hi Daniel, not sure why when trying the above I get: Dot::dotsh: "Tensors {0.00718226`,0.00718226...have incompatible shapes? – sebastian c. Jan 15 '13 at 0:03
Did you change all the 1 x n matrices into simple vectors? Also I forgot the defn of volFit, which I'll edit in right now. – Daniel Lichtblau Jan 15 '13 at 0:15
Yes Flatten it. – sebastian c. Jan 15 '13 at 0:23

Your Answer


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.