Black-Scholes program for the HP-12C

That was a funny challenge: cramming the Black-Scholes equation in an HP-12C.

Usage

Fill "i" register with the interest, as usual (e.g. 12 for 12% per year).
Fill "PV" with the spot value of the underlying asset.
Fill "PMT" with the volatility, as a percentage (e.g. 25 for 25% per year).
Fill "n" with the time to expiration. Use the same time unit as interest and volatility. If you enter interest and volatility as % per year, you must supply the time to expiration expressed in years (e.g. 2 months = 0.1666 years).
Fill "FV" with the strike of the option.
Press R/S to get the call value. It takes several seconds.
The value on display is the call value. Press x⇄y to see the put value, and again to go back to the call value.

The option parameters are not changed during calculation, so you can fiddle with them and press R/S as many times as you want.

The program

Opcode #	Operation	Remarks
	0.196854
	STO 3	N(x) 1st coefficient
	0.115194
	STO 4	N(x) 2nd coefficient
	0.000344
	STO 5	N(x) 3rd coefficient
	0.019527
	STO 6	N(x) 4th coefficient
	f P/R	Goes into programming mode
	f CLEAR PRGM	Clears programming memory
06	RCL n	t
07	RCL i	100.r
08	%	rt
08	CHS	-rt
08	e^x	e^-rt
16	RCL FV	K
13	×	K.e^-rt
02	STO EPX	√t
15	RCL PV	S
15	÷	K.e^-rt/S
15	1/x	S.e^rt/K
18	g LN	ln(S/K) + rt
02	STO 0	√t
01	RCL n	t
02	g √x	√t
03	RCL PMT	100.σ
04	%	√t.σ
17	STO ÷ 0	S/K
05	STO 2	Let mem2 = √t.σ
10	2
17	÷	S/K
05	STO - 0	Let mem2 = √t.σ
05	RCL 0	Let mem2 = √t.σ
32	4	calculate cummulative normal distribution
33	Y^x	d⁴
34	RCL 6	c4
35	x⇄y
36	×	c4.d⁴
37	g LSTx	d⁴
38	g √x	d²
39	RCL 4	c2
40	x⇄y
41	×	c2.d²
42	g LSTx	d²
43	g √x	d (guaranteed to be positive)
44	RCL 3	c1
45	×	c1.d
46	+	c1.d+c2.d²
47	+	c1.d+c2.d²+c4.d⁴
48	RCL 0	d
49	6
50	Y^x	d⁶
51	g √x	d³ (guaranteed to be positive)
52	RCL 5	c3
53	×	c3.d³
54	+	c1.d+c2.d² + c3.d³ + c4.d⁴
55	1
56	+	1+c1.d+c2.d² + c3.d³ + c4.d⁴
57	4
58	CHS
59	Y^x	(1+c1.d+c2.d² + c3.d³ + c4.d⁴)^{- 4}
60	2
61	÷	(1+c1.d+c2.d² + c3.d³ + c4.d⁴)^{- 4}/2
62	0	Current result is 1-N(abs(d))
63	RCL 0
64	g x≤y	d < 0?
65	GTO 72	d<0, N(d)=1-N(-d), we are done, jump ahead
66	R↓	d>0, need to do 1-(1-(N(d))
67	R↓	Restore 1-N(d) into register X
68	CHS
69	1
70	+	N(d) corrected, we are done
71	GTO 74	Jump ahead
72	R↓	Restore N(d) into register X
73	R↓
75	RCL 2	mem2 is √t.σ at 1st round, zero otherwise
76	g x=0	mem2 is zero?
77	GTO 81	Jump to calculate call and put
78	RCL + 0	Let mem0 = d2 + √t.σ = d1
79	STO - 2	Let mem2 = √t.σ - √t.σ = zero
73	R↓	Recalls N(d2)
74	STO 1	Let mem1 = N(d2)
80	GTO 26	Go back to calculate N(d1)
81	R↓	Recalls N(d1)
81	RCL PV	S
36	×	S.N(d1)
82	STO 0	Let mem0 = S.N(d1)
89	RCL EPX	K.e^-rt
90	STO × 1	Let mem1 = K.e^-rt.N(d2)
91	RCL PV	S
92	-	K.e^-rt - S
93	RCL 0	S.N(d1)
94	RCL 1	K.e^-rt . N(d2)
95	-	S.N(d1) - K.e^-rt . N(d2) = Call
96	+	Call + K.e^-rt - S = Put
97	g LSTx	Call at X, put at Y
	f P/R	Back to normal mode

Implementation remarks

It was difficult to fit the equation in 99 programming steps, and using all steps reduces the available memory, which in turn makes the program more complicated.

The most critical function is the cummulative normal distribution, that we implement using a rational approximation with precision of 4 decimal digits. The coefficients are entered as part of the program but they are stored on memory positions (the STO 3, STO 4... at the beginning of the table). There are better approximations but they could not be fitted in such a small "computer".