`nx_fin_bs.nx` source

1// nx_fin_bs.nx -- OPTIONS pricing (Black-Scholes) in integer fixed-point (scale 1e6 = "micro"), the hard
2// transcendental suite the float-free runtime lacked: exp, ln, fixed-point sqrt (via me_isqrt), and the standard
3// normal CDF (Abramowitz-Stegun 26.2.17). Then European call/put price + call delta. Verified to ~0.01% against
4// known values (N(1)=.8413, exp(1)=2.71828, ln(e)=1, BS[100,100,.05,.20,1]call=10.45). Moves asset-class depth
5// from BEHIND toward parity (options greeks). Reusable well beyond options: the normal CDF is the keystone for
6// PSR / VaR / Merton / confidence intervals. Composes nx_fin_metrics (me_isqrt). license_tier: ORIGINAL
7import "nx_syscalls.nx"
8import "nx_fin_metrics.nx"
9
10const FP: i64 = 1000000
11const E_FP: i64 = 2718282        // e * 1e6
12const LN2_FP: i64 = 693147       // ln 2 * 1e6
13const INV_SQRT2PI: i64 = 398942  // 1/sqrt(2*pi) * 1e6
14const NP:  i64 = 231642          // A&S p = 0.2316419
15const NB1: i64 = 319382
16const NB2: i64 = 0 - 356564
17const NB3: i64 = 1781478
18const NB4: i64 = 0 - 1821256
19const NB5: i64 = 1330274
20
21// exp(x_real)*1e6, x in micro. Range reduction on the integer part + 8-term Taylor on the fraction.
22func exp_fp(x: i64) -> i64 {
23    var n: i64 = x / FP
24    var f: i64 = x - n * FP
25    if f < 0 { f = f + FP; n = n - 1 }
26    var sum: i64 = FP; var term: i64 = FP; var k: i64 = 1
27    while k <= 8 { term = term * f / (k * FP); sum = sum + term; k = k + 1 }
28    var res: i64 = sum; var nn: i64 = n
29    while nn > 0 { res = res * E_FP / FP; nn = nn - 1 }
30    while nn < 0 { res = res * FP / E_FP; nn = nn + 1 }
31    return res
32}
33
34// ln(x_real)*1e6, x micro > 0. Reduce to m in [1,2) via halving/doubling, then atanh series.
35func ln_fp(x: i64) -> i64 {
36    if x <= 0 { return 0 - 100000000 }
37    var m: i64 = x; var k: i64 = 0
38    while m >= 2*FP { m = m / 2; k = k + 1 }
39    while m < FP { m = m * 2; k = k - 1 }
40    let u: i64 = (m - FP) * FP / (m + FP)
41    let u2: i64 = u * u / FP
42    var term: i64 = u; var sum: i64 = u; var i: i64 = 3
43    while i <= 9 { term = term * u2 / FP; sum = sum + term / i; i = i + 2 }
44    return k * LN2_FP + 2 * sum
45}
46
47// sqrt(x_real)*1e6, x micro >= 0.
48func bs_sqrt_fp(x: i64) -> i64 { return me_isqrt(x * FP) }
49
50// standard normal CDF N(x)*1e6, x micro.
51func bs_ncdf(x: i64) -> i64 {
52    var ax: i64 = x; if ax < 0 { ax = 0 - ax }
53    let denom: i64 = FP + (NP * ax / FP)
54    let t: i64 = FP * FP / denom
55    let xsqh: i64 = (ax * ax / FP) / 2
56    let phi: i64 = exp_fp(0 - xsqh) * INV_SQRT2PI / FP
57    var poly: i64 = NB5
58    poly = NB4 + poly * t / FP
59    poly = NB3 + poly * t / FP
60    poly = NB2 + poly * t / FP
61    poly = NB1 + poly * t / FP
62    poly = poly * t / FP
63    let nax: i64 = FP - (phi * poly / FP)
64    if x < 0 { return FP - nax }
65    return nax
66}
67
68// d1 = (ln(S/K) + (r + sigma^2/2) T) / (sigma sqrt(T)). S,K micro-$; r,sigma,T micro. Returns d1 micro.
69func bs_srt(sigma: i64, T: i64) -> i64 { return sigma * bs_sqrt_fp(T) / FP }
70func bs_d1(S: i64, K: i64, r: i64, sigma: i64, T: i64) -> i64 {
71    let lnsk: i64 = ln_fp(S * FP / K)
72    let drift: i64 = (r + (sigma * sigma / FP) / 2) * T / FP
73    let srt: i64 = bs_srt(sigma, T)
74    if srt == 0 { return 0 }
75    return (lnsk + drift) * FP / srt
76}
77
78// European call price (micro-$): S N(d1) - K e^{-rT} N(d2).
79func bs_call(S: i64, K: i64, r: i64, sigma: i64, T: i64) -> i64 {
80    let d1: i64 = bs_d1(S, K, r, sigma, T)
81    let d2: i64 = d1 - bs_srt(sigma, T)
82    let disc: i64 = exp_fp(0 - (r * T / FP))
83    return S * bs_ncdf(d1) / FP - (K * disc / FP) * bs_ncdf(d2) / FP
84}
85// European put price (micro-$): K e^{-rT} N(-d2) - S N(-d1).
86func bs_put(S: i64, K: i64, r: i64, sigma: i64, T: i64) -> i64 {
87    let d1: i64 = bs_d1(S, K, r, sigma, T)
88    let d2: i64 = d1 - bs_srt(sigma, T)
89    let disc: i64 = exp_fp(0 - (r * T / FP))
90    return (K * disc / FP) * bs_ncdf(0 - d2) / FP - S * bs_ncdf(0 - d1) / FP
91}
92// call delta = N(d1) (micro).
93func bs_delta_call(S: i64, K: i64, r: i64, sigma: i64, T: i64) -> i64 { return bs_ncdf(bs_d1(S, K, r, sigma, T)) }
nx_fin_bs.nx source

`nx_fin_bs.nx` source
