Dyalog Forums

i need a big pascal's triangle (say, 10,000 rows)

tried with normal numbers
⎕FR←645 ⋄⎕ts⋄ sink←{⍉(⍳⍵)∘.!⍳⍵} 1029 ⋄ ⎕ts
2022 8 19 10 52 37 354
2022 8 19 10 52 37 571
⎕FR←645 ⋄⎕ts⋄ sink←{⍉(⍳⍵)∘.!⍳⍵} 1030 ⋄ ⎕ts
2022 8 19 10 52 50 483
DOMAIN ERROR
⎕FR←645 ⋄ ⎕TS ⋄ sink←{⍉(⍳⍵)∘.!⍳⍵}1030 ⋄ ⎕TS

bignums takes about 5 minutes for 2060 rows, about 19 minutes for 4096 rows, still running after an hour for 10240 rows

∧
⎕FR←1287 ⋄ ⎕ts⋄ sink←{⍉(⍳⍵)∘.!⍳⍵}2060 ⋄ ⎕ts
2022 8 19 10 53 4 848
2022 8 19 10 58 23 138

anybody have any clever tricks to speed it up ? [i suppose i could wait several hours to compute it once and store the result in a memory mapped file. but having to multiply a bignum pascal triangle by another matrix would take a long time ... but i was hoping someone whose comp sci is better than mine would say "here is a way to generate the logs of pascals triangle,and instead of multiplying a matrix by pascal triangle you could just add logs, and taking the antilogs is easy enough, and 64 bit logs will give you enough precision for 10,000 rows ")

the use case is the function below, in particular lines 14, 32, and 33

∇ crrtree2 [⎕]∇
[0] crrtree2←{
[1] ⍝ why separate crr tree generation from crr pricing?
[2] ⍝ cuz sometimes you just want to pathwise monte carlo through a crr tree not backwards loop over it
[3] ⍝ so break tree generation out from traditional crr pricing
[4]
[5] ⍝ ⍵[1]=strike, ⍵[2]=spotprice,
[6] ⍝ ⍵[3]=expiry in years, 1 month=12÷365.25 or 12÷tradingdays←12×21
[7] ⍝ ⍵[4]=annualized volatility, 21%/yr would = 0.21
[8] ⍝ ⍵[5]=NUMBER of timesteps N.B. root=timestep zero, so tree with 6 end leaves has 5 timesteps
[9] ⍝ ⍵[6]=irate, 5%/yr would = 0.05
[10] ⍝ ⍵[7]=optional (for stock) annual dividend rate
[11] ⍝ MANDITORY THAT ⍵[6]=⍵[7] FOR OPTIONS ON FUTURES
[12]
[13] ⍝ result[1;;] is stock price at node, result[2;;] prob ANY SINGLE PATH reaching node
[14] ⍝ result[3;;] total prob reaching node = (prob SINGLE PATH)×(NUMB PATHS=Pascals triangle)
[15]
[16] spot expiry vol timesteps irate divrate←1↓7↑⍵ ⍝ 1↓ as trees dont need strike arg of option pricer ⍵
[17] ⎕IO←1 ⍝ I've 5 fingers, not 4 fingers and zeroeth thumb
[18]
[19] mask←(⍳1+timesteps)∘.≤⍳1+timesteps⍝ "mask×" will be used to zero out tree below diagonal
[20]
[21] u←*vol×(∆t←expiry÷timesteps)*0.5 ⍝ multiplicative size up jump
[22] d←*-vol×∆t*0.5 ⍝ multiplicative size down jump, also = 1÷u
[23]
[24] uppowers←maskׯ1++\mask ⍝ numb up jumps to reach a node, ¯1+ sets tree root to zero jumps
[25] downpowers←mask×|uppowers-[2]¯1+⍳1+timesteps ⍝ numb down jumps to get to node
[26]
[27] prices←mask×spot×(u*uppowers)×d*downpowers ⍝ prices at node through time.
[28] ⍝ ⍝ yep, above prices table couldve been a tacit fork with some ⊢s, ⊣s, and ⍨s
[29] ⍝ ⍝ I DONT CARE, CLARITY OVER CODE GOLF
[30]
[31] Pd←1-Pu←((*∆t×irate-divrate)-d)÷u-d ⍝ Pu=up probability, Pd=down probability
[32] P1Path←mask×(Pu*uppowers)×Pd*downpowers ⍝ probability ANY ONE SINGLE path hits node
[33] PAnyPath←mask×P1Path×⍉pascal timesteps+1 ⍝ probability ANYofALL paths hits node
[34] ↑prices P1Path PAnyPath ⍝ result = rank 3 flat array ⍴=(3,(1+timesstep),1+timestep)
[35] }

crrtree2 uses this subfunction
∇ pascal [⎕]∇
[0] pascal←{⍝ returns ⍵ lines pof pascals triangle
[1]
[2] ⎕IO←0
[3] ⍉(⍳⍵)∘.!⍳⍵
[4] }

above two functions are subfunctions of a Cox Ross Rubenstein option pricer

∇ crrpricer2 [⎕]∇
[0] crrpricer2←{
[1] ⍝ return CRR value of an American or European option
[2] ⍝ Euro redundant to closed form Euro Black Scholes, but useful for futures & calibrating volatility smiles tables
[3]
[4] ⍝ ⍺[1]←1 for American exercise, 0 for European Exercise
[5] ⍝ ⍺[2] put/call ¯1=put, +1=call
[6]
[7] ⍝ ⍵[1]= strike, ⍵[2]=spotprice
[8] ⍝ ⍵[3]=expiry in years, 1 month=12÷year←365.25 or 12÷year←12×21 (21=tradingdays per month)
[9] ⍝ ⍵[4]=annualized volatility, 35 percent =0.35
[10] ⍝ ⍵[5]=NUMBER of timesteps N.B. root=timestep zero, so tree with 6 end leaves has 5 timesteps
[11] ⍝ ⍵[6]=irate, 5%/yr=0.05
[12] ⍝ ⍵[7]=optional (for stock) annual dividend rate
[13] ⍝ MANDITORY THAT ⍵[6]=⍵[7] FOR FUTURES
[14]
[15] isAmer PutCall←⍺
[16] strike spot expiry vol timesteps irate divrate←7↑⍵
[17]
[18] ⎕IO←1 ⍝ I've 5 fingers, not 4 fingers and zeroeth thumb
[19]
[20] mask←(⍳1+timesteps)∘.≤⍳1+timesteps ⍝ mask has 1s for valid nodes, 0s for invalid below diagonal
[21]
[22] trees←crrtree2 ⍵ ⍝ treeS[1;;]=prices, [2;;]=prob 1 single path hits node, [3;;]=sum of prob each individual path h
its node
[23]
[24] ⍝ disc∆t←÷*irate×expiry÷timesteps ⍝ present val of $1 1∆t in future
[25] disc∆t←÷1+irate×expiry÷timesteps ⍝ should be EXACTLY same as above, but small diff
[26]
[27] Pu Pd←2↑trees[3;;2] ⍝ prob up and down jumps
[28] prctree←(↓[1]mask)/¨↓[1]trees[1;;] ⍝ trimmed nested VectorOfVectors of price tree
[29] intrinsic←0⌈PutCall×prctree-|strike ⍝ in the money amount at each node
[30]
[31] {(⍺×isAmer)⌈disc∆t×(Pu,Pd)+.×0 ¯1↓0 1⌽⍵,[0.5]⍵}/intrinsic ⍝ THE BIG FOLD FROM THE RIGHT
[32] ⍝ 0k, wrote the above as a joke to prove i can write as cryptic a line of APL as anyone else.
[33] ⍝ IN PRODUCTION CODE I'd write "foo/intrinsics" while defining "foo" as below
[34]
[35] ⍝ foo←{
[36] ⍝ ⍝ whats going on with that reduction above
[37] ⍝ ⍝ ⍺="current i.e " now's " intrinsics, ⍵=intrinscs 1∆t forward in time
[38]
[39] ⍝ intrinsics←0 ¯1↓ 0 1 ⌽ ⍵,[0.5]⍵ ⍝ intrinsic 1∆t ahead [1;]=intrinsic if upjump, [2;]=intrinsic if downjump
[40] ⍝ expval←(Pu,Pd)+.×intrinsics ⍝ expval=(ProbUpJump×upIntrinsics)+(ProbDownJump×downIntrinsics)
[41] ⍝ discval←disc∆t×expval ⍝ discval=discounted expval to "now" of intrinsics[1∆t ahead]
[42]
[43] ⍝ result← (⍺×isAmer)⌈discval ⍝ IsAmer×MAX(intrinsic "now" vs expval intrinsic 1∆t ahead)
[44] ⍝ ⍝ if isAmer=0, makes now's intrinsics=0, equivalent to European exercise
[45]
[46] ⍝ ⎕←result
[47] ⍝ }
[48] }

Does this help? Tried a naive, scalar N*2 time-complexity algorithm...


Results, where original is your original function...
and with WSMAX:"4G" (this size selected out of a hat)

Squinting at the teeny tiny type, on the teeny tiny screen, of my teeny tiny phone, in the dimly lit back seat of the teeny tiny Toyota Prius uber, on my way from one smoky jazz club to another ... All i can say is thank you, and i'll have to contemplate this when i get up at the crack of noon tomorrow

If you're interested, it appears practical to calculate a 10000-row triangle in about 10 seconds in Dyalog APL using a simple convolution algorithm, as long as you have enough virtual memory.

I couldn't help but try to fine-tune Peter's original algorithm:

m←pascal3 n;r;⎕IO;⎕FR

 (⎕IO ⎕FR)←1 1287
 m←n n⍴n↑1

 :For r :In ⍳n
     m[r;1+⍳r-1]←⌊×\(⌽÷⊢)⍳r-1
 :EndFor

Benchmarking gives

]runtime "pascal2 10000"

* Benchmarking "pascal2 10000"
               (ms) 
 CPU (avg):   86063 
 Elapsed:    123997 
     ]runtime "pascal3 10000"

* Benchmarking "pascal3 10000"
              (ms) 
 CPU (avg):  39296 
 Elapsed:    39341

-Veli-Matti

Side Note: As a side note, my first attempt surprised me. I just assumed the naive Python version would have to run slowly in APL, but the algorithmic improvements dominated over the interpreter overhead.

After looking at random math[s] sites, I decided to use a convolution to replace the element-by-element solution with a row-at-a-time solution, similar to yours (V-M's). As a random site points out, "The simplest and most well-known convolution triangle is Pascal's triangle, which is formed by convolving the sequence {l, 1, 1, • • •} with itself repeatedly." The code pascalCD (convolution dfn version) is shown below.

First, a comparison of your (elegant and concise) pascal3 and pascalCD: And now, the function pascalCD: Algorithmic fun!

Oh, that's nice!

I just realised that the min function was unnecessary in pascal3 making the calculation
slightly quicker: m[r;1+i]←×\(⌽i)÷i←⍳r-1
..but it is still half as quick as your pascalCD.
Hats off!

Nevertheless, you might be amused with this different kind of approach

pascal4←{
     (⎕IO ⎕DIV)←1 ⋄ ⎕FR←(⍵>999)(↑⌽)645 1287
     1,0⍪×\(∘.≥⍨v)×((⌽v-1)⌽r⍴⌽v)÷(r←2⍴⍵-1)⍴v←⍳⍵-1
 }

-wm :)

petermsiegel wrote:If you're interested, it appears practical to calculate a 10000-row triangle in about 10 seconds in Dyalog APL using a simple convolution algorithm, as long as you have enough virtual memory.

Peter,

First, let me apologize for not getting back to you sooner. I spent a few days running around New York like a chicken without a head taking photographs of meeting rooms on behalf of somebody else on this forum (you know who you are).

Funny you should mention convolution. The whole point of my wanting such a big grid is to try to price variance gamma options without resorting to the Fourier transform convolution that Madan & Carr use to solve the problem of pricing options on an asset that undergoes a variance gamma process. Mostly because I can't remember how to do a a fast fourier transform. (There's a lot of things I used to be able to do that I can't do anymore. For instance, I used to play three sports back in school, now I run out of breath tying my shoes)

My idea about the options pricing is that just as Mark Rubenstein demonstrated that a trinomial options tree is just a binomial options tree with twice as many timesteps, (i.e delta t time step of the trinomial is 1/2 delta t time step of the binomial) you could evaluate a binomial options tree only at random columns where the space between each random column was a gamma variate, and that ought to yield you a variance gamma options price (since a variance gamma process is just a long normal process observed on a "gamma clock")

So yes, I would be very interested in a convolution solution. Not only would it make my code run faster,seeing convolution done in APL might spark some dim forgotten memory of how to do fourier transforms.

Veli-Matti wrote:I couldn't help but try to fine-tune Peter's original algorithm:

m←pascal3 n;r;⎕IO;⎕FR

 (⎕IO ⎕FR)←1 1287
 m←n n⍴n↑1

 :For r :In ⍳n
     m[r;1+⍳r-1]←⌊×\(⌽÷⊢)⍳r-1
 :EndFor

Benchmarking gives

]runtime "pascal2 10000"

* Benchmarking "pascal2 10000"
               (ms) 
 CPU (avg):   86063 
 Elapsed:    123997 
     ]runtime "pascal3 10000"

* Benchmarking "pascal3 10000"
              (ms) 
 CPU (avg):  39296 
 Elapsed:    39341

-Veli-Matti

A train. A train with a tack in it ... That's going to take me about an hour to figure out.

Being thrilled by pascalCD I noticed that you can pretty easily make that even quicker by just simplifying the code.

]runtime "pascal3 10000"   

* Benchmarking "pascal3 10000"
              (ms) 
 CPU (avg):  15875 
 Elapsed:    15877 
      ]runtime "pascalCD 10000"

* Benchmarking "pascalCD 10000"
             (ms) 
 CPU (avg):  9282 
 Elapsed:    9277 
      ]runtime "pascalCD2 10000"

* Benchmarking "pascalCD2 10000"
             (ms) 
 CPU (avg):  3609 
 Elapsed:    3608

Leaving the comments out:

pascalCD2←{
     ⎕IO←0 ⋄ ⎕FR←(⍵>999)(↑⌽)645 1287
     ∆←⍵ ⍵⍴⍵↑1 ⋄ r←0
     CNext←{∆[r;⍳1+r⊣r+←1]←(⍵,0)+0,⍵}
     ∆⊣CNext⍣(⍵-1)⊣1
 }

-wm