Trinomial tree option pricing python

I've just started picking up Python and have built binomial and trinomial models just to test my understanding, especially about arrays.

For terminal stock price array, I have stepped down from InitialStock * u**[iSteps - i] for i going from zero to 2*iSteps+ 1. ie. I start at top-most terminal stock price InitialStock * u**[iSteps] and travel down until I get to InitalStock * u**[-iSteps] and populate the TerminalStock array from 0 to 2*iSteps+1 [ie. last array element is InitialStock * u**[-iSteps]. Because I come from the land of VBA, I haven't yet picked up the eloquent [and hence, brief and sometimes difficult-to-read Python style], my loops look like this:

for i in range[0, 2*iSteps+1] :

    # Terminal stock is initial with up staps and down steps
    dblStockTerminal = dblStock * u ** float[iSteps - i]

    # compute intrinsic values at terminal stock value
    dblTrinOpt[i][iSteps] = max[ dblSwitch * [dblStockTerminal - dblStrike], 0 ]

I have initialised dblTrinOpt array as dblTrinOpt = np.ndarray[ [ 2*iSteps+1, iSteps+1], float] so that it has 2*iSteps price elements at terminal stock price and iStep time steps, where iSteps = OptionTerm / NumberOfSteps. Yes, apology for my clunky python code! [but I have to read it].

after that, my option calc is as follows [yes, again the code is inelegant]:

#------------------------------------------
# steps in time from Terminal to Initial Stock price
#------------------------------------------
for i in range[iSteps-1, -1, -1] :

    #------------------------------------------
    # steps in price range from highest to lowest
    #------------------------------------------
    for j in range[0, 2*i+1] :

        #------------------------------------------
        # discount average of future option value
        # as usual, trinomial averages three values
        #------------------------------------------
        dblTrinOpt[j][i] = dblDisc * [pu * dblTrinOpt[j][i+1] + pm * \
            dblTrinOpt[j+1][i+1] + pd * dblTrinOpt[j+2][i+1]]

after that, my Trinomial function passes dblTrinOpt[0][0] as final discounted option price.

I wanted to write this in a more readable form so that it became easier to add two components: [a] an American option facility, hence tracking stock prices at each tree node [and so then things like barriers and Bermudans would be an easy addition, and [b] add two additional terminal prices one at each end and so that at time zero, I'd have three stock and option values so that I could compute delta and gamma without needing to recompute the option tree [i.e.. avoid the bump-and-recompute approach]. I hope that helps, even though I didn't fix your specific code, but explained the logic.

Trinomial trees in options pricing

In the binomial tree, each node leads to two other nodes in the next time step. Similarly in a trinomial tree, each node leads to three other nodes in the next time step. Besides having up and down states, the middle node of the trinomial tree indicates no change in state. When extended over more than two time steps, the trinomial tree can be thought of as a recombining tree, where the middle nodes always retain the same values as the previous time step.

Let's consider the Boyle trinomial tree, where the tree is calibrated such that the probability of up, down, and flat movements,

def setup_parameters[self]:    """ Required calculations for the model """    self.u = math.exp[self.sigma*math.sqrt[2.*self.dt]]    self.d = 1/self.u    self.m = 1    self.qu = [[math.exp[[self.r-self.div] *                         self.dt/2.] -                math.exp[-self.sigma *                         math.sqrt[self.dt/2.]]] /               [math.exp[self.sigma *                         math.sqrt[self.dt/2.]] -                math.exp[-self.sigma *                         math.sqrt[self.dt/2.]]]]**2    self.qd = [[math.exp[self.sigma *                         math.sqrt[self.dt/2.]] - math.exp[[self.r-self.div] ...

Trinomial Model: the Auld Triangle

A two-jump process for the asset price over each discrete time step was developed in the binomial lattice. Boyle expanded this frame of reference and explored the feasibility of option valuation by allowing for an extra jump in the stochastic process. In keeping with Black Scholes, Boyle examined an asset [S] with a lognormal distribution of returns. Over a small time interval, this distribution can be approximated by a three-point jump process in such a way that the expected return on the asset is the riskless rate, and the variance of the discrete distribution is equal to the variance of the corresponding lognormal distribution. The three point jump process was introduced by Phelim Boyle[1986] as a trinomial tree to price options and the effect has been momentous in the finance literature. Perhaps shamrock mythology or the well-known ballad associated with Brendan Behan inspired the Boyle insight to include a third jump in lattice valuation. His trinomial paper has spawned a huge amount of ground breaking research. In the trinomial model, the asset price S is assumed to jump uS or mS or dS after one time period [dt = T/n], where u > m > d. Joshi [2008] point out that the trinomial model is characterized by the following five parameters: [1] the probability of an up move pu, [2] the probability of an down move pd, [3] the multiplier on the stock price for an up move u, [4] the multiplier on the stock price for a middle move m, [5] the multiplier on the stock price for a down move d. A recombining tree is computationally more efficient so we require

ud = m*m

M = exp [r∆t],

V = exp [σ 2∆t],

dt or ∆t = T/N

where where N is the total number of steps of a trinomial tree. For a tree to be risk-neutral, the mean and variance across each time steps must be asymptotically correct. Boyle [1986] chose the parameters to be:

m = 1, u = exp[λσ√ ∆t], d = 1/u

pu =[ md − M[m + d] + [M^2]*V ]/ [u − d][u − m] ,

pd =[ um − M[u + m] + [M^2]*V ]/ [u − d][m − d]

Boyle suggested that the choice of value for λ should exceed 1 and the best results were obtained when λ is approximately 1.20. One approach to constructing trinomial trees is to develop two steps of a binomial in combination as a single step of a trinomial tree. This can be engineered with many binomials CRR[1979], JR[1979] and Tian [1993] where the volatility is constant. For example, a two-step presentation of the CRR[1979] is presented below in C++:

b = r - q;

dt = T/n;

u = exp[v*sqrt[2.0*dt]];

d = 1.0/u;

pu = pow[[exp[0.5*b*dt] - exp[-v*sqrt[0.5*dt]]] / [exp[v*sqrt[0.5*dt]] - exp[-v*sqrt[0.5*dt]]], 2];

pd = pow[[exp[v*sqrt[0.5*dt]] - exp[0.5*b*dt]] / [exp[v*sqrt[0.5*dt]] - exp[-v*sqrt[0.5*dt]]], 2];

pm = 1.0 - pu - pd;

or in R using Google Colab

# Time Interval:

dt = Time/n

# Up-and-down jump sizes:

u = exp[+sigma * sqrt[2*dt]]

d = exp[-sigma * sqrt[2*dt]]

# Probabilities of going up and down:

pu = [[exp[b * dt/2] - exp[ -sigma * sqrt[dt/2]]] / [exp[sigma * sqrt[dt/2]] - exp[-sigma * sqrt[dt/2]]]] ^ 2

pd = [[ exp[sigma * sqrt[dt/2]] - exp[ b * dt/2]] / [exp[sigma * sqrt[dt/2]] - exp[-sigma * sqrt[dt/2]]]] ^ 2

# Probability of staying at the same asset price level:

pm = 1 - pu - pd

We initially present below two snippets of C++ code which are investigated for efficiency. Given the numerically intensive nature of lattices we propose making use of dynamic memory techniques in the lattice estimation. We compare the performance of a static vs dynamic memory set up in pricing models not unlike the approach adapted for ESOs in the previous page and consistent with techniques previously outlined for the binomial estimation where memory usage was optimized. The static code was extracted from Volopta.com developed by Fabrice Rouah. The Dynamic representation was transposed from from Espen Haug's Excel VBA trinomial tree. Espen's dynamic tree design appears to be more efficient and can reach higher step size on the online C++ compilers. We also present how to manually set up the trinomial tree in an excel workbook so readers/viewers can compare the basic visual representation against the actual C++ code presented below.

Static Trinomial Tree in C++

// Static Trinomial Tree C++

// by Fabrice Douglas Rouah, FRouah.com and Volopta.com

//#include "StdAfx.h"

#include

#include

#include

#include

#include

#include

#include

#include

using namespace std;

double Trinomial[double Spot,double K,double r,double q,double v,double T,char PutCall,char EuroAmer,int n]

{

int i,j;

double b = r - q;

double dt = T/n;

double u = exp[v*sqrt[2.0*dt]];

double d = 1.0/u;

double pu = pow[[exp[0.5*b*dt] - exp[-v*sqrt[0.5*dt]]] / [exp[v*sqrt[0.5*dt]] - exp[-v*sqrt[0.5*dt]]], 2];

double pd = pow[[exp[v*sqrt[0.5*dt]] - exp[0.5*b*dt]] / [exp[v*sqrt[0.5*dt]] - exp[-v*sqrt[0.5*dt]]], 2];

double pm = 1.0 - pu - pd;

// Build the trinomial tree for the stock price evolution

vector S[2*n+1,vector [n+1]];

for [j=0; j