# Correct way to generate large data sets (i.e.forward yield curve )

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:

Table[forwardYieldCurve[[1,j]]=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]
+((tenor0[[1,j+1]]-tenor0[[1,j]])/(dateArray[[1,j+1]]-dateArray[[1,j]]))tstep,{j,99}];//AbsoluteTiming


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

Table[forwardYieldCurve[[i+1,j]]=forwardYieldCurve[[i,j]]+driftM[[1,j]]tstep+
(volFit1[[1,j]]dX[[i+1,1]]+volFit2[[1,j]]dX[[i+1,2]]+volFit3[[1,j]]dX[[i+1,3]])Sqrt[tstep]+
((forwardYieldCurve[[i+1,j+1]]-forwardYieldCurve[[i+1,j]])/(dateArray[[1,j+1]]-dateArray[[1,j]]))tstep,{j,99},{i,**100**}];//AbsoluteTiming


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}
tenor0={{0.0050399,0.00537318,0.00578648,0.00614997,0.00633987,0.00637105,0.00632311,0.00625459,0.00622594,0.00631663,0.0065289,0.00679745,0.00706621,0.00731132,0.0075159,0.00766905,0.00778107,0.00786696,0.00793966,0.00800508,0.00806759,0.00813158,0.00820143,0.00828151,0.00837543,0.00848368,0.00860596,0.00874199,0.00889147,0.00905412,0.00922964,0.00941775,0.00961814,0.00983054,0.0100546,0.0102902,0.0105368,0.0107941,0.0110615,0.0113385,0.0116248,0.0119197,0.0122228,0.0125336,0.0128516,0.0131763,0.0135073,0.013844,0.0141859,0.0145327,0.0148838,0.0152389,0.0155975,0.0159592,0.0163236,0.0166903,0.0170588,0.0174287,0.0177995,0.0181709,0.0203931,0.0225666,0.0246436,0.0265946,0.0283977,0.0300428,0.0315247,0.0328461,0.0340124,0.035033,0.0359187,0.0366814,0.0373332,0.0378862,0.0383519,0.0387395,0.0390575,0.0393143,0.0395184,0.0396782,0.0398011,0.0398898,0.0399458,0.0399704,0.0399652,0.0399316,0.039871,0.0397848,0.0396746,0.0395418,0.0393879,0.0392142,0.0390222,0.0388134,0.0385892,0.0383511,0.0381006,0.037839,0.0375678,0.0372885}}
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
volFit1={{0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226,0.00718226}}
volFit2={{-5.97435*10^-6,-5.77847*10^-6,-5.58514*10^-6,-5.39435*10^-6,-5.20606*10^-6,-5.02029*10^-6,-4.83699*10^-6,-4.65618*10^-6,-4.47782*10^-6,-4.3019*10^-6,-4.12841*10^-6,-3.95734*10^-6,-3.78867*10^-6,-3.62238*10^-6,-3.45846*10^-6,-3.2969*10^-6,-3.13768*10^-6,-2.98078*10^-6,-2.8262*10^-6,-2.67392*10^-6,-2.52391*10^-6,-2.37618*10^-6,-2.23069*10^-6,-2.08745*10^-6,-1.94643*10^-6,-1.80762*10^-6,-1.671*10^-6,-1.53656*10^-6,-1.40429*10^-6,-1.27417*10^-6,-1.14618*10^-6,-1.02032*10^-6,-8.9656*10^-7,-7.74893*10^-7,-6.55303*10^-7,-5.37775*10^-7,-4.22296*10^-7,-3.08849*10^-7,-1.97421*10^-7,-8.79963*10^-8,1.94392*10^-8,1.249*10^-7,2.28402*10^-7,3.29958*10^-7,4.29584*10^-7,5.27294*10^-7,6.23103*10^-7,7.17026*10^-7,8.09077*10^-7,8.9927*10^-7,9.87622*10^-7,1.07415*10^-6,1.15886*10^-6,1.24177*10^-6,1.3229*10^-6,1.40226*10^-6,1.47986*10^-6,1.55573*10^-6,1.62987*10^-6,1.7023*10^-6,2.10178*10^-6,2.44338*10^-6,2.73027*10^-6,2.96562*10^-6,3.1526*10^-6,3.29436*10^-6,3.39409*10^-6,3.45494*10^-6,3.4801*10^-6,3.47271*10^-6,3.43596*10^-6,3.37301*10^-6,3.28703*10^-6,3.18118*10^-6,3.05863*10^-6,2.92256*10^-6,2.77613*10^-6,2.6225*10^-6,2.46485*10^-6,2.30635*10^-6,2.15015*10^-6,1.99943*10^-6,1.85736*10^-6,1.72711*10^-6,1.61183*10^-6,1.51471*10^-6,1.43891*10^-6,1.38759*10^-6,1.36392*10^-6,1.37108*10^-6,1.41223*10^-6,1.49054*10^-6,1.60917*10^-6,1.77129*10^-6,1.98008*10^-6,2.2387*10^-6,2.55031*10^-6,2.91809*10^-6,3.3452*10^-6,3.83481*10^-6}}
volFit3={{1.85601*10^-6,1.86013*10^-6,1.8634*10^-6,1.86583*10^-6,1.86742*10^-6,1.8682*10^-6,1.86815*10^-6,1.86729*10^-6,1.86563*10^-6,1.86316*10^-6,1.85991*10^-6,1.85587*10^-6,1.85105*10^-6,1.84546*10^-6,1.8391*10^-6,1.83199*10^-6,1.82412*10^-6,1.81551*10^-6,1.80616*10^-6,1.79608*10^-6,1.78528*10^-6,1.77375*10^-6,1.76152*10^-6,1.74858*10^-6,1.73494*10^-6,1.72062*10^-6,1.70561*10^-6,1.68992*10^-6,1.67356*10^-6,1.65653*10^-6,1.63885*10^-6,1.62052*10^-6,1.60154*10^-6,1.58193*10^-6,1.56169*10^-6,1.54082*10^-6,1.51934*10^-6,1.49724*10^-6,1.47454*10^-6,1.45125*10^-6,1.42737*10^-6,1.4029*10^-6,1.37786*10^-6,1.35224*10^-6,1.32607*10^-6,1.29933*10^-6,1.27205*10^-6,1.24422*10^-6,1.21586*10^-6,1.18697*10^-6,1.15756*10^-6,1.12763*10^-6,1.09719*10^-6,1.06625*10^-6,1.03481*10^-6,1.00289*10^-6,9.70481*10^-7,9.37598*10^-7,9.04246*10^-7,8.70433*10^-7,6.58225*10^-7,4.31082*10^-7,1.90459*10^-7,-6.21896*10^-8,-3.25408*10^-7,-5.97743*10^-7,-8.77739*10^-7,-1.16394*10^-6,-1.45489*10^-6,-1.74914*10^-6,-2.04524*10^-6,-2.34172*10^-6,-2.63713*10^-6,-2.93002*10^-6,-3.21893*10^-6,-3.50241*10^-6,-3.77901*10^-6,-4.04726*10^-6,-4.30572*10^-6,-4.55292*10^-6,-4.78742*10^-6,-5.00776*10^-6,-5.21248*10^-6,-5.40014*10^-6,-5.56927*10^-6,-5.71842*10^-6,-5.84613*10^-6,-5.95096*10^-6,-6.03144*10^-6,-6.08613*10^-6,-6.11356*10^-6,-6.11228*10^-6,-6.08084*10^-6,-6.01778*10^-6,-5.92165*10^-6,-5.791*10^-6,-5.62436*10^-6,-5.42028*10^-6,-5.17731*10^-6,-4.894*10^-6}}
randomWalkPCA[n_]:=  RandomVariate[NormalDistribution[0,1],n];
RandVarPCA[mcRun_]:=Table[randomWalkPCA[3],{mcRun}];
(dX:=RandVarPCA[1000])//MatrixForm;
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.}}

• +1 for a well written and formatted question. It took a while, but you did it, and that's all that matters! :)
– rm -rf
Jan 14, 2013 at 14:37
• 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, 2013 at 14:56
• Yeah +1. You're getting a grip on the site's etiquette. Jan 14, 2013 at 15:17
• 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. Jan 14, 2013 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. Jan 14, 2013 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

Differences[tenor0[[1]]]/Differences[dateArray[[1]]]


...and so on.

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

e.g.

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)

• Hi @Mike Honeychurch, could this be the reason why I am getting increased volatility with successive simulation runs as I noted: mathematica.stackexchange.com/questions/17781/… Jan 14, 2013 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? Jan 14, 2013 at 23:13

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]

Timing[
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 ---

• 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? Jan 14, 2013 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. Jan 14, 2013 at 23:44
• Hi Daniel, not sure why when trying the above I get: Dot::dotsh: "Tensors {0.00718226`,0.00718226...have incompatible shapes? Jan 15, 2013 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. Jan 15, 2013 at 0:15
• Yes Flatten it. Jan 15, 2013 at 0:23