# The model

## Black Scholes model

Under the risk neutral probability measure $\mathbb{Q}$, the underlying is defined as a geometric Brownian motion solving the following SDE:

\begin{equation}
\begin{cases}
\frac{d S_t}{S_t} = r dt + \sigma d W_t \\
S_0 = x
\end{cases}
\end{equation}

From that SDE, one can explcitiely compute the first two moments of $S$:

\begin{equation}
\begin{cases}
\mathbb{E}_t^{\mathbb{Q}} (S_{t+h}) = S_t e^{r h} \\
\mathbb{E}_t^{\mathbb{Q}} (S_{t+h}^2) = S_t^2 e^{2 r h + \sigma^2 h}
\end{cases}
\end{equation}

## The trinomial tree model matching the Black-Scholes model

In the sequel, we will use a recombining tree to model the underlying moves.
We define by $u, u > 1$, some multiplicative factor of the underlying from time $t$ to time $t+h$.
In the trinomial model, the underlying at time $t+h$ can take 3 values knowing the $S_t$:

\begin{equation}
\begin{cases}
S_{t+h} = u \times S_t > S_t \text { with probability } p_u \\
S_{t+h} = d \times S_t < S_t \text { with probability } p_d \\
S_{t+h} = S_t \text { with probability } p_m = 1 - p_u - p_d \\
\end{cases}
\end{equation}

with by vertue of our recombining assumption, $d = \frac{1}{u}$.

That is, given $u$, we must find the matching probabilities $p_u, p_d, p_m$ matching the Black-Scholes model.
Thus we have 3 unknown variables. We must at least give 3 independent equations to solve the system.

The probabilities should be non negative and summing to one, which is already written for the definition of $p_m$. We are left with 2 unknown variables $p_u$ and $p_d$. 
Matching the first two moments given above we have:

\begin{equation}
\begin{cases}
\mathbb{E}_t^{\mathbb{Q}} (S_{t+h}) = S_t e^{r h} = S_t \left(u p_u + (1 - p_u - p_m) + d p_d \right) \\
\mathbb{E}_t^{\mathbb{Q}} (S_{t+h}^2) = S_t^2 e^{2 r h + \sigma^2 h} = S_t^2 \left(u^2 p_u + (1 - p_u - p_m) + d^2 p_d \right)
\end{cases}
\end{equation}

So we have:

$$
\begin{pmatrix}
e^{r h} - 1 \\
e^{2rh + \sigma^2 h} - 1
\end{pmatrix}
= 
\begin{pmatrix}
u - 1 & d - 1 \\
u^2 - 1 & d^2 - 1
\end{pmatrix}
\begin{pmatrix}
p_u \\
p_d
\end{pmatrix}
$$

or equivalently, solving the linear system:

$$
\begin{pmatrix}
u - 1 & d - 1 \\
u^2 - 1 & d^2 - 1
\end{pmatrix}^{-1}
\begin{pmatrix}
e^{r h} - 1 \\
e^{2rh + \sigma^2 h} - 1
\end{pmatrix}
= 
\begin{pmatrix}
p_u \\
p_d
\end{pmatrix}
$$

And using the classical formula of the inverse of a 2x2 matrix, we have:
$$
\begin{pmatrix}
u - 1 & d - 1 \\
u^2 - 1 & d^2 - 1
\end{pmatrix}^{-1}
= \frac{1}{(u - 1) (d^2 - 1) - (u^2 - 1) (d - 1)}
\begin{pmatrix}
d^2 - 1 & 1 - d \\
1 - u^2 & u - 1
\end{pmatrix}
$$

and $p_m := 1 - p_u - p_d$.
Moreover, once the computation is done, one should check that all the probabilities are non negative.


```python
import numpy as np

from scipy.stats import norm
```


```python
class TrinomialBSModel(object):
    def __init__(self, S0=100., r=0.02, sigma=0.2, mat=1.):
        self.__s0 = S0
        self.__r = r
        self.__sigma = sigma
        self.__T = mat
        
    def __compute_probs(self):
        B = np.array([-1. + np.exp(self.__r * self.__h), 
                      -1. + np.exp(2. * self.__r * self.__h + self.__sigma**2 * self.__h)])
        
        d = self.__down
        u = self.__up
        
        A = np.array([[u - 1., d - 1.], 
                      [u**2 - 1., d**2 - 1.]])
        
        det = (u - 1.) * (d**2 - 1.) - (u**2 - 1.) * (d - 1.)
        invA = 1. / det * np.array([[d**2 - 1., 1. - d], 
                                    [1. - u**2, u - 1.]])
        
        res = invA.dot(B)
        
        self.__pu = res[0]
        self.__pd = res[1]
        self.__pm = 1. - self.__pu - self.__pd
        
        assert 0 <= self.__pu <= 1., 'p_u should lie in [0, 1] given %s' % self.__pu
        assert 0 <= self.__pd <= 1., 'p_d should lie in [0, 1] given %s' % self.__pd
        assert 0 <= self.__pm <= 1., 'p_m should lie in [0, 1] given %s' % self.__pm
        
    def __check_up_value(self, up):
        if up is None:
            lbda = np.sqrt(0.5 * np.pi)
            up = np.exp(lbda * self.__sigma * np.sqrt(self.__h))
            
        assert up > 0., 'up should be non negative'
        
        down = 1. / up
        assert down < up, 'up <= 1. / up = down'
                
        self.__up = up
        self.__down = down
        
    def __gen_stock_vec(self, nb):
        vec_u = self.__up * np.ones(nb)
        np.cumprod(vec_u, out=vec_u)

        vec_d = self.__down * np.ones(nb)
        np.cumprod(vec_d, out=vec_d)
        
        res = np.concatenate((vec_d[::-1], [1.], vec_u))
        res *= self.__s0
        
        return res
    
    def payoff(self, stock_vec):
        raise NotImplementedError()
        
    def compute_current_price(self, crt_vec_stock, nxt_vec_prices):
        expectation = np.zeros(crt_vec_stock.size)
        for i in range(expectation.size):
            tmp = nxt_vec_prices[i] * self.__pd
            tmp += nxt_vec_prices[i + 1] * self.__pm
            tmp += nxt_vec_prices[i + 2] * self.__pu
            
            expectation[i] = tmp
            
        return self.__discount * expectation
        
    def price(self, nb_steps, up=None):
        assert nb_steps > 0, 'nb_steps shoud be > 0'
        
        nb_steps = int(nb_steps)
        
        self.__h = self.__T / nb_steps
        self.__check_up_value(up)
        self.__compute_probs()
        
        self.__discount = np.exp(-self.__r * self.__h)
       
        final_vec_stock = self.__gen_stock_vec(nb_steps)
        final_payoff = self.payoff(final_vec_stock)
        nxt_vec_prices = final_payoff
        
        for i in range(1, nb_steps + 1):
            vec_stock = self.__gen_stock_vec(nb_steps - i)
            nxt_vec_prices = self.compute_current_price(vec_stock, nxt_vec_prices)
            
        return nxt_vec_prices[0]        
```


```python
class TrinomialBSCall(TrinomialBSModel):
    def __init__(self, S0=100., r=0.02, sigma=0.2, mat=1., K=100.):
        super(TrinomialBSCall, self).__init__(S0, r, sigma, mat)
        self.__K = K
    
    def payoff(self, s):
        return np.maximum(s - self.__K, 0.)
        
class TrinomialBSAmericanCall(TrinomialBSCall):
    def compute_current_price(self, crt_vec_stock, nxt_vec_prices):
        crt_payoff = self.payoff(crt_vec_stock)
        crt_prices = super(TrinomialBSAmericanCall, self).compute_current_price(crt_vec_stock, nxt_vec_prices)
        
        return np.maximum(crt_payoff, crt_prices)
```


```python
class TrinomialBSPut(TrinomialBSModel):
    def __init__(self, S0=100., r=0.02, sigma=0.2, mat=1., K=100.):
        super(TrinomialBSPut, self).__init__(S0, r, sigma, mat)
        self.__K = K
    
    def payoff(self, s):
        return np.maximum(self.__K - s, 0.)
    
class TrinomialBSAmericanPut(TrinomialBSPut):
    def compute_current_price(self, crt_vec_stock, nxt_vec_prices):
        crt_payoff = self.payoff(crt_vec_stock)
        crt_prices = super(TrinomialBSAmericanPut, self).compute_current_price(crt_vec_stock, nxt_vec_prices)
        
        return np.maximum(crt_payoff, crt_prices)
```


```python
def bs_call_price(S=100., r=0.02, sigma=0.2, t=0., T=1., K=100.):
    ttm = T - t
    
    if ttm < 0:
        return 0.
    elif ttm == 0.:
        return np.maximum(S - K, 0.)

    vol = sigma * np.sqrt(ttm)

    d_minus = np.log(S / K) + (r - 0.5 * sigma**2) * ttm
    d_minus /= vol

    d_plus = d_minus + vol

    res = S * norm.cdf(d_plus)
    res -= K * np.exp(-r * ttm) * norm.cdf(d_minus)

    return res
```


```python
def bs_put_price(S=100., r=0.02, sigma=0.2, t=0., T=1., K=100.):
    # Using call-put parity :)
    ttm = T - t
    
    if ttm < 0:
        return 0.
    elif ttm == 0.:
        return np.maximum(K - S, 0.)

    dsct_strike = K * np.exp(-r * ttm)

    cap_S = S * np.exp(r * t)

    call = bs_call_price(S, r, sigma, t, T, K)
    return call - cap_S + dsct_strike
```


```python
tree = TrinomialBSCall()

print(tree.price(1000))
print(bs_call_price())
```

    8.91559126943
    8.91603727857



```python
tree = TrinomialBSPut()

print(tree.price(1000))
print(bs_put_price())
```

    6.93545860012
    6.93590460925



```python
tree = TrinomialBSAmericanCall()

print(tree.price(1000))
```

    8.91559126943



```python
tree = TrinomialBSAmericanPut()

print(tree.price(1000))
```

    7.11070940695


2efae404e62d703d64685b7090d7518db8458bb5
ipynb
Jupyter Notebook
Black Scholes trinomial tree.ipynb
yarmenti/bs_trinomial_tree
bad9c22b44e23881be15de4000345cd44b381edb
MIT
Black Scholes trinomial tree.ipynb
yarmenti/bs_trinomial_tree
bad9c22b44e23881be15de4000345cd44b381edb
MIT
Black Scholes trinomial tree.ipynb
yarmenti/bs_trinomial_tree
bad9c22b44e23881be15de4000345cd44b381edb
MIT
Qwen/Qwen-72B
1. YES
2. YES