need a big pascal's table

APL-related discussions - a stream of APL consciousness.
Not sure where to start a discussion ? Here's the place to be
Forum rules
This forum is for discussing APL-related issues. If you think that the subject is off-topic, then the Chat forum is probably a better place for your thoughts !

need a big pascal's table

Postby tclviii-dyalog on Fri Aug 19, 2022 3:52 pm

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] }
tclviii-dyalog
 
Posts: 28
Joined: Tue Mar 02, 2010 6:04 pm

Re: need a big pascal's table

Postby petermsiegel on Sat Aug 20, 2022 3:31 am

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

Code: Select all
 ∇p←pascal2 num
 ;c;k;n;q
 ;⎕IO;⎕FR

 ⍝ Pascal's triangle via an O(N2) alg 
 ⍝ Non-vectorized "naive" algorithm based on a Python pgm below.
 ⍝ Equiv. dfn trivial to write and similar performance...
 ⍝
 ⍝ There are ways to improve the algorithm, but watch out for WS full.
 ⍝ ⍝⍝⍝ ORIGINAL PYTHON...
 ⍝   def pascalTriangle(num):
 ⍝     for n in range(1, num + 1):
 ⍝        C = 1
 ⍝        for k in range(1, n + 1):
 ⍝          # value of first column is always 1
 ⍝          print(C, end = " ")
 ⍝          C = int(C * (n - k) / k)
 ⍝        print("")

 ⎕IO←1
 ⎕FR←1287             ⍝ Anticipate larger Pascal triangles

 ⍝ Set up array p in advance to limit copying
 p←num num⍴0
 :For n :In ⍳num
     q←,c←1
     :For k :In ⍳¯1+n
         c←⌊c×(n-k)÷k
         q,←c
     :EndFor
     p[n;⍳≢q]←q
 :EndFor
⍝∇⍣§./pascal2.aplf§0§ 2022 8 19 20 21 14 951 §bUUrM§0

Results, where original is your original function...
Code: Select all
       (original 500)≡pascal2 500
1
      ⎕TS ⋄ sink←⍬ ⋄ sink←pascal2 2048 ⋄ ⎕TS
2022 8 19 20 21 56 279
2022 8 19 20 21 58 728

and with WSMAX:"4G" (this size selected out of a hat)
Code: Select all
 ⎕TS ⋄ sink←⍬ ⋄ sink←pascal2 10000 ⋄ ⎕TS
2022 8 19 20 46 27 237
2022 8 19 20 47 29 727
Last edited by petermsiegel on Mon Aug 22, 2022 3:22 am, edited 4 times in total.
petermsiegel
 
Posts: 127
Joined: Thu Nov 11, 2010 11:04 pm

Re: need a big pascal's table

Postby tclviii-dyalog on Sat Aug 20, 2022 7:11 am

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
tclviii-dyalog
 
Posts: 28
Joined: Tue Mar 02, 2010 6:04 pm

Re: need a big pascal's table

Postby petermsiegel on Mon Aug 22, 2022 11:38 pm

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.
petermsiegel
 
Posts: 127
Joined: Thu Nov 11, 2010 11:04 pm

Re: need a big pascal's table

Postby Veli-Matti on Tue Aug 23, 2022 7:30 am

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
Veli-Matti
 
Posts: 92
Joined: Sat Nov 28, 2009 3:12 pm

Re: need a big pascal's table

Postby petermsiegel on Tue Aug 23, 2022 5:29 pm

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:
Code: Select all
      (pascalCD 10000)≡pascal3 10000    ⍝ ⎕FR←1287 required
1

      ]runtime "pascal3 10000"
* Benchmarking "pascal3 10000"
              (ms)
 CPU (avg):  15048
 Elapsed:    15341

        ]runtime "pascalCD 10000"
* Benchmarking "pascalCD 10000"
             (ms)
 CPU (avg):  9455
 Elapsed:    9509

And now, the function pascalCD:
Code: Select all
    ∇ pascalCD←{                                                 
[1]   ⍝ Pascal's Triangle via convolution (Dfn version)           
[2]   ⍝ For info/code on convolve, see notes.xtimes in ws "dfns".
[3]   ⍝ triangle(∆) ← pascalCD ndim(⍵)                           
[4]   ⍝    ndim: a pos. integer                                   
[5]   ⍝    ∆:    array of shape (ndim ndim)                       
[6]   ⍝ pascalCD 10000: under 10 sec on my  my Macbook Air 2020 
[7]                                                               
[8]        ⎕IO←0 ⋄ ⎕FR←1287                                       
[9]   ⍝ convolve:  Requires ⎕IO←0                                 
[10]       convolve←{+⌿(-⍳⍴⍺)⌽⍺∘.×⍵,0×1↓⍺}                       
[11]                                                             
[12]  ⍝ Set up the triangle ∆ and the power (⍣) "loop"           
[13]       ∆←⍵ ⍵⍴0                                               
[14]       r←0                                                   
[15]       ∆[r;0]←1   ⍝ row number (r) incremented in CNext       
[16]       CNext←{∆[r;⍳1+r⊣r+←1]←⌊1 1 convolve ⍵}                 
[17]                                                             
[18]     ⍝ Iterate: Each row based on convolution of the prior row
[19]       ∆⊣CNext⍣(⍵-1)⊣1                                       
[20]                                                             
[21]   }                                                         
     ∇   
Algorithmic fun!
petermsiegel
 
Posts: 127
Joined: Thu Nov 11, 2010 11:04 pm

Re: need a big pascal's table

Postby Veli-Matti on Tue Aug 23, 2022 7:32 pm

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 :)
Veli-Matti
 
Posts: 92
Joined: Sat Nov 28, 2009 3:12 pm

Re: need a big pascal's table

Postby tclviii-dyalog on Tue Aug 23, 2022 7:45 pm

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.
tclviii-dyalog
 
Posts: 28
Joined: Tue Mar 02, 2010 6:04 pm

Re: need a big pascal's table

Postby tclviii-dyalog on Tue Aug 23, 2022 7:50 pm

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.
tclviii-dyalog
 
Posts: 28
Joined: Tue Mar 02, 2010 6:04 pm

Re: need a big pascal's table

Postby Veli-Matti on Wed Aug 24, 2022 6:03 am

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
Veli-Matti
 
Posts: 92
Joined: Sat Nov 28, 2009 3:12 pm

Next

Return to APL Chat

Who is online

Users browsing this forum: No registered users and 1 guest