Black Scholes Merton option with only Python standard library (Work in progress)

option.py

from math import sqrt, erf, pi, exp, log

N = lambda x: (1.0 + erf(x / sqrt(2.0))) / 2.0
phi = lambda x: exp(-0.5*x*x)/sqrt(2.0*pi)

class Option:
    def __init__(self, St, K, t, s, r=0, q=0):
        self.St = St # 0.90
        self.K = K # 1.00
        self.s = s # 0.20
        self.r = r # 0.0
        self.q = q # 0.0
        self.t = t # 1.0

    @property
    def d1(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        d1 = (log(St/K) + (r - q + s*s/2) * t) / s / sqrt(t)
        return d1

    @property
    def d2(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        d1 = self.d1
        d2 = d1 - s * sqrt(t)
        return d2

    @property
    def C(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        d1 = self.d1; d2 = self.d2
        C = St * exp(-q*t) * N(d1) - K * exp(-r*t) * N(d2)
        return C
    
    @property
    def P(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        d1 = self.d1; d2 = self.d2
        P = -St * exp(-q*t) * N(-d1) + K * exp(-r*t) * N(-d2)
        return P

    @property
    def dCdK(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        d2 = self.d2
        return -exp(-r*t)*N(d2)

    @property
    def dPdK(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        return self.dCdK + exp(-r*t)

    @property
    def z(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        mu = (r - q - s*s/2.0) * t
        return ((log(K) - log(St)) - mu) / s / sqrt(t)

    @property
    def probability(self):
        return N(self.z)

    @property
    def dCdt(self):
        St = self.St; K = self.K; s = self.s; r = self.r; q = self.q; t = self.t
        d1 = self.d1; d2 = self.d2
        return (
            St * (-q) * exp(-q*t) * N(d1) 
            - K * (-r) * exp(-r*t) * N(d2) 
            + St * exp(-q*t) * phi(d1) * s / 2 / sqrt(t)
            )