Performance issues with Trinomial Tree to calculate price of option

Question

A program I wrote isn't performing quite to where I had hoped so I went looking to find where the majority of the time of the program was being spent and it was in the function below so I would like to get some eyes on it to help tune it up and make it run a little quicker. Please note I realize it isn't the most readable but I heavily commented it to help prevent this and wrote it in the manner I did since my sole goal of the function is to get the result as fast as possible.

The function uses an inlined function of the "Option" class called getPayoffType() which simply returns a value of the class as well as this function:

inline double Option::getPayoff(double ulPrice)
{
return ((cp == Call) ? ((ulPrice + 0.000005 > strike) ? ulPrice - strike : 0.0) 
        : ((strike + 0.000005 > ulPrice) ? strike - ulPrice : 0.0));
}

Here's the code:

void TrinomialEngine(Option* option, Underlying* underlying)
{
int timeSteps = 20; //set the number of "time steps" aka the number of steps in the tree, this will later be dynamic

double vol = option->vol; //get a parameter of our model from the Option class
double rate = option->rate; //"" ""
PayoffType payoff = option->getPayoffType(); //" "  " "

double dt = option->getTimeToExpiration() / timeSteps; //the length of time (in years) each time step represents
double v = rate - vol*vol*0.5; //working in log-units
double x = vol*sqrt(2.0*dt); //stock-price jump-size
double dis = exp(-rate*dt); //timestep discounting
double edx = exp(x); //precomputing this constant to save time

double pu = 0.5*(dt*(vol*vol + v*v*dt)/x/x + (v*dt/x));  //up probability
double pm = 1.0 - dt*(vol*vol + v*v*dt)/x/x;  //middle probability
double pd = 0.5*(dt*(vol*vol + v*v*dt)/x/x - (v*dt/x));  //down probability


int nodes = timeSteps * timeSteps; //how many "nodes" there are in the tree
double* tree = new double[nodes * 3];

/*how the tree is modeled in the array (each asterisk is a node, number is array position)

             11 Option
             10 Div    *
             9 Stock

  2 Option   8 Option
  1 Div   *  7 Div     *
  0 Stock    6 Stock

             5 Option       
             4 Div     *
             3 Stock    

timesteps   1            2

        0  3  12  28
        1  2   3   4
        1  4   9   16
*/

tree[0] = underlying->theo; //get a parameter from our "Underlying" class

for(int i = 1; i < timeSteps; ++i) //set up the "Stock" prices at each node
{
    tree[i*i*3] = tree[0]*exp(-i*x); //bottom node of each time step
    tree[i*i*3 + 1] = 0.0; //this will be changed later to actually hold a value

    for(int j = i*i*3 + 3; j < (i+1)*(i+1)*3; j += 3)
    {
        //working up the nodes for the current "time step"
        tree[j] = tree[j - 3] * edx;
        tree[j + 1] = 0.0; //this will be changed later to actually hold a value
    }
}

//value option at expiry
for(int i = timeSteps * timeSteps * 3 - 3; i >= (timeSteps - 1) * (timeSteps - 1) * 3; i -= 3)
{
    //calulating the "Option" value for each node on the last time step of the tree
    tree[i + 2] = option->getPayoff( tree[i] + tree[i + 1] ); //inlined function of "Option" class
}

int j( (timeSteps - 2) * 6 + 3 );
int nodeCount( (timeSteps - 1) * 2 - 1 );
int currStep( timeSteps - 1 );
double exVal(0.0);

if (payoff == American) //assume this is always true for now
{   
    for(int i = (timeSteps - 1) * (timeSteps - 1) * 3 - 3; i >= 0; i -= 3)
    {
        //calculate the value of the option at every other node, working backwards through time steps

        tree[i + 2] = dis * (tree[i + j + 2] * pd + tree[i + j + 5] * pm + tree[i + j + 8] * pu); //discounted value of option

        exVal = option->getPayoff( tree[i] + tree[i + 1] ); //value of option if exercised

        //the value of the option at these nodes is the maximum of the discounted value and the exercised value
        if (exVal > tree[i + 2] + 0.000005)
            tree[i + 2] = exVal;

        --nodeCount; //tick back how many nodes we ahve evaluated
        if (nodeCount < 1)
        {
            //this is if we've evaluated all the nodes in the timestep so reset our indices
            j -= 6;
            --currStep;
            nodeCount = currStep * 2 - 1;
        }
    }
} else
{
    for(int i = (timeSteps - 1) * (timeSteps - 1) * 3 - 3; i >= 0; i -= 3)
    {
        tree[i + 2] = dis * (tree[i + j + 2] * pd + tree[i + j + 5] * pm + tree[i + j + 8] * pu);

        --nodeCount;
        if (nodeCount < 1)
        {
            j -= 6;
            --currStep;
            nodeCount = (currStep - 1) * 2 - 1;
        }
    }
}

//grab our results!
option->theo = tree[2]; 
option->delta = ( tree[11] - tree[5] ) / ( tree[9] - tree[3] );

delete[] tree; //don't forget to free the memory
}

Answer 1

My advice:

1. Precompute whatever is possible to precompute. I would say: all those vol*vol, x*x (especially before the double division /x/x). Your compiler may be smart enough to make some optimization itself, but it's better not to risk and to check if some performance improvement can be achieved. It think, however, that you will spend most of the time computing what is inside the loops.

2. sqrt() is often a time expensive function. If you can get rid of it, that's good. But I don't know if it makes sense for you to compute sqrt(2.0*dt) outside TrinomialEngine. If getTimeToExpiration() is almost constant (as is timeSteps) in your application, you could try to precompute that sqrt somewhere else, in order to have it ready. A lot depends on how many times do you call TrinomialEngine and how often are "options" changed.

3. Try to avoid operations inside the loop condition, index declaration, and keep increment "easy".

   This is good:

   $for(int i = 1; i < timeSteps; ++i)
   

   This is bad:

   $for(int j = i*i*3 + 3; j < (i+1)*(i+1)*3; j += 3)
   

   The risk in having all those "complex" for loops is that your compiler will fail to optimize them properly, since it is (most of the time) designed to be good at optimizing loops from 0 to N.

   This should give you most of the performance improvement that you can get. I cannot see any other big issue.

4. It has nothing to do with performances, but I would put the brackets {} around tree[i + 2] = exVal; after if (exVal > tree[i + 2] + 0.000005). That's good practice.