This was a funny challenge: putting the Black-Scholes equation into HP-12C. It calculates the theoretical premium for calls and puts.

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 shown value is the call value. Press x⇄y to see the put value, and again to go back to 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
01	RCL n	t
02	g √x	√t
03	RCL PMT	100.σ
04	%	√t.σ
05	STO 2	Let mem2 = √t.σ
06	RCL n	t
07	RCL i	100.r
08	%	r.t
09	RCL 2	√t.σ
10	2
11	Yx	t.σ2
12	2
13	÷	t.σ2/2
14	+	(r + σ2/2).t
15	RCL PV	S
16	RCL FV	K
17	÷	S/K
18	g LN	ln(S/K)
19	+	ln(S/K) + (r + σ2/2).t
20	RCL 2	√t.σ
21	÷	d1
22	STO 1	Let mem1 = d1
23	STO 0	Let mem0 = d1
24	RCL 2	√t.σ
25	STO - 0	Let mem0 = d1 - √t.σ = d2
26	RCL 0	swap mem0 and mem1
27	RCL 1
28	STO 0
29	R↓
30	STO 1	swap done
31	RCL 0	Get d1 or d2 (whichever is at mem0) and
32	4	calculate cummulative normal distribution
33	Yx	d4
34	RCL 6	c4
35	x⇄y
36	×	c4.d4
37	g LSTx	d4
38	g √x	d2
39	RCL 4	c2
40	x⇄y
41	×	c2.d2
42	g LSTx	d2
43	g √x	d (guaranteed to be positive)
44	RCL 3	c1
45	×	c1.d
46	+	c1.d+c2.d2
47	+	c1.d+c2.d2+c4.d4
48	RCL 0	d
49	6
50	Yx	d6
51	g √x	d3 (guaranteed to be positive)
52	RCL 5	c3
53	×	c3.d3
54	+	c1.d+c2.d2+c3.d3+c4.d4
55	1
56	+	1+c1.d+c2.d2+c3.d3+c4.d4
57	4
58	CHS
59	Yx	(1+c1.d+c2.d2+c3.d3+c4.d4)-4
60	2
61	÷	(1+c1.d+c2.d2+c3.d3+c4.d4)-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↓
74	STO 0	Let mem0 = N(d)
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	0
79	STO 2	Let mem2 = 0
80	GTO 26	Go back to 2nd round of N(d) calculation
81	RCL PV	S
82	STO × 1	Let mem1 = S.N(d1)
83	RCL n	t
84	RCL i	r
85	%	r.t
86	CHS	-r.t
87	g ex	e-r.t
88	RCL FV	K
89	×	K.e-rt
90	STO × 0	Let mem0 = K.e-rt.N(d2)
91	RCL PV	S
92	-	K.e-rt - S
93	RCL 1	S.N(d1)
94	RCL 0	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 make the equation to fit into HP-12C, since a program can have at most 99 steps, the "language" is poor, and using all 99 steps reduces the number of available memory positions, which in turn makes the program more complicated and forces usage of more opcodes.

The most complicated function is the cummulative normal distribution function, which we implement using a rational approximation with precision of 4 decimal digits. The coefficients of this function are entered as part of the program but they lie on memory positions instead of inside the program. A better, seven-coefficient function could not be used because of opcode and memory constraints.