.py

def BinomHeston(self, n_steps, opt_type, American=False, comp_fun=None, kappa=None, var_mean=None, var_var=None, return_Hs=False):
    '''
    Calculates option value with stochastic volatility according to a
    discretization of the Heston model presented on p. 754
                                                                                                                                   
    ===========
    n_steps: number of steps in binomial calculation
    opt_type: 'call' or 'put'
    American: Boolean; if True then value of a `type` American option is calculated
                                                                                                                                   
    returns float, the option value
    '''
                                                                                                                                   
    '''
    Parameters that are fixed
    '''
    S, K, vol, r, T, div = self.S, self.K, self.vol, self.r, self.T, self.div
                                                                                                                                   
    # kappa = .01
    kappa = kappa if kappa is not None else .1
    var_mean = var_mean if var_mean is not None else self.vol**2
    var_variance = var_var if var_var is not None else .1**2
                                                                                                                                   
    h0 = T / n_steps
    # v0 = vol ** 2
    v0 = vol ** 2
    # print(v0)
                                                                                                                                   
    def dZ(dt):
        return np.random.choice([-1, 1]) * np.sqrt(dt)
                                                                                                                                   
    def dv(vt, h):
        # import ipdb; ipdb.set_trace()
        # return (kappa * (var_mean - vt)*h) + (var_variance * np.sqrt(vt) * np.sqrt(dZ(h)))
        return kappa*(var_mean - vt)*h + var_variance * np.sqrt(vt) * np.sqrt(h) * np.random.normal()
        # return (kappa * (var_mean - vt)*h) + (var_variance * np.sqrt(vt) * dZ(h))
                                                                                                                                   
    def simulate_vols(v0, h, var_var, n_steps):
        '''
        Simulate volatilities according to that kappa model.
        returns the volatilities
        '''
        # v0 = var_mean
        L = [v0]
        for t in range(n_steps):
            d = dv(v0, h)
            v0 = v0
            if v0 + d > 0:
                v0 = v0 + d
            else:
                while True:
                    N = simulate_vols(v0, h**2/100, var_var, n_steps=100)[-1]
                    if not np.isnan(N) and v0 + N > 0:
                        v0 = N
                        break
                                                                                                                                   
            L.append(v0)
        return np.sqrt(np.array(L))
                                                                                                                                   
    def get_up_down(vol, h):
        return (np.exp((r-div)*h + vol*np.sqrt(h)),
                np.exp((r-div)*h - vol*np.sqrt(h)))
                                                                                                                                   
    def replicating_portfolio_value(Cu, Cd, vol, h):
        u, d = get_up_down(vol, h)
        delta = np.exp(-div*h) * (Cu-Cd) / (S*(u-d))
        B = np.exp(-r*h) * (u*Cd - d*Cu) / (u - d)
        return delta*S + B
                                                                                                                                   
    def step_backward(C_T, vol, h):
        '''
        V_T: a list of values for the option at step T. Two values next
        to each other represent the binomial values that they might have.
                                                                                                                                   
        returns: Option values for the underlying given prices at step T-1
        '''
        S_t = []
        for Cu, Cd in zip(C_T, C_T[1:]):
            S_t.append(replicating_portfolio_value(Cu, Cd, vol, h))
        return S_t
                                                                                                                                   
    Hs, Vs = [], []  # h-steps and variances
    ###for step in range(n_steps):
    ###    Vs.append(np.sqrt(v0))
    ###    v1 = max(0.01, v0 + dv(v0, h0))  # ensure that it cannot turn negative.
    ###    h1 = v0/v1 * h0
    ###    Hs.append(h1)
    ###    h0 = Hs[-1]
    ###    v0 = v1
                                                                                                                                   
    ### NEW
    Hs.append(h0)
    Vs = simulate_vols(v0, h0, var_variance, n_steps)
    # for step in range(n_steps):
    #     Vs.append(np.sqrt(v0))
    #     v0 = max(.01, v0 + dv(v0, h0))
    for step in range(n_steps-1):
        h0 = Hs[-1] * (Vs[step]/Vs[step+1])**2  # TODO: shouldnt this be squared?
        Hs.append(h0)
    ### NEW END
    # import ipdb; ipdb.set_trace()
                                                                                                                                   
    Hs = np.array(Hs)
    Hs = Hs / Hs.sum() * self.T
    # import ipdb; ipdb.set_trace()
                                                                                                                                   
    u, d = get_up_down(Vs[-1], Hs[-1])
    S_T = [S * u**i * d**(n_steps-i) for i in np.arange(n_steps+1)[::-1]]
    if comp_fun is not None:
        V_T = [comp_fun(self, s, 0) for s in S_T]
    else:
        V_T = [max(0, x - K) if opt_type == 'call' else max(0, K - x) for x in S_T]
    V_t = V_T
                                                                                                                                   
    for h, vol, t in zip(Hs[::-1], Vs[::-1], np.arange(n_steps)[::-1]):
        u, d = get_up_down(vol, h)
                                                                                                                                   
        V_t = step_backward(V_t, vol, h)
        S_t = [S * u**i * d**(t-i) for i in np.arange(t+1)[::-1]]
        if comp_fun is not None:
            V_t = [comp_fun(self, s, v) for s, v in zip(S_t, V_t)]
            continue
                                                                                                                                   
        if American and opt_type == 'call':
            V_t = [max(c_t, s_t-K) for (c_t, s_t) in zip(V_t, S_t)]
        if American and opt_type == 'put':
            V_t = [max(p_t, K-s_t) for (p_t, s_t) in zip(V_t, S_t)]
                                                                                                                                   
    if return_Hs:
        return V_t[0], Hs
    # import ipdb; ipdb.set_trace()
    return V_t[0]