Writting "A Simple Least-Squares Approach" by Longstaff and Schwartz; pricing American option contracts

Question

Background:

I have made a post about this on another account here. I am writing the least squares algorithm into a class in C++ and I want to make sure that what I am doing is the most efficient and hopefully fast. I used the Eigen library to write all the sub-routines to price the American option contracts. I have not completed the algorithm yet but I have a majority of the sub-routines done, and tested them to make sure they are working correctly.

Question:

I want to know if there is anything I can do to improve my code as it stands now or if there is anything I am doing wrong in terms of the syntax of writing the class. Here is the code:

Here is the header file:

/*
 * LSM.h
 *
 *  Created on: Oct 8, 2017
 *      Author: 
 */

#include <vector>
#include <Eigen/Dense>
#include <Eigen/Geometry>




#ifndef LSM_H
#define LSM_H




class LSM {
public:
    // Overload Constructor
    LSM(const double, const double, const double, const int, const int, const double, const double, const int, const int);

    // Destructor
    ~LSM();

    // Generate the Laguerre Polynomials
    Eigen::MatrixXd Laguerre(Eigen::VectorXd, const int);

    // Generate M paths of stock prices (Geometric Brownian Motion)
    Eigen::VectorXd GBM(const int, const int, const double, const double, const double, const double, const double);

    // Payoff of call option
    Eigen::VectorXd callPayoff(Eigen::VectorXd, const double);

    // Payoff of put option
    Eigen::VectorXd putPayoff(Eigen::VectorXd, const double);

    // Find function for finding the paths that are in the money (call option)
    Eigen::VectorXd Findcallpath(Eigen::VectorXd, const double);

    // Find function for finding the paths that are in the money (put option)
    Eigen::VectorXd Findputpath(Eigen::VectorXd, const double);

    // Find price of call given path
    Eigen::VectorXd Findcallprices(Eigen::VectorXd, Eigen::VectorXd);

    // Find price of put given path
    Eigen::VectorXd Findputprices(Eigen::VectorXd, Eigen::VectorXd);

    // Find return of call (stock price - strike price)
    Eigen::VectorXd Findcallreturn(Eigen::VectorXd, const double);

    // Find return of put (strike price - stock price)
    Eigen::VectorXd Findputreturn(Eigen::VectorXd, const double);

    // Using Two-sided Jacobi SVD decomposition of a rectangular matrix
    Eigen::VectorXd Jacobi(Eigen::MatrixXd, Eigen::VectorXd);






private:
    // Member variables
    double new_r;
    double new_q;
    double new_sigma;
    int new_T;
    int new_N;
    double new_K;
    double new_S0;
    int new_M;
    int new_R;

};






#endif

Here is the .cpp file associated with the header file:

#include <iostream>
#include <vector>
#include <random>
#include <time.h>
#include <math.h>
#include "LSM.h"
#include <Eigen/Dense>
#include <Eigen/Geometry>


LSM::LSM( const double r, const double q, const double sigma, const int T, const int N, const double K, const double S0, const int M, const int R){
    new_r = r;
    new_q = q;
    new_sigma = sigma;
    new_T = T;
    new_N = N;
    new_K = K;
    new_S0 = S0;
    new_M = M;
    new_R = R;


/*  Eigen::VectorXd V(4);
    V(0) = 100;
    V(1) = 102;
    V(2) = 103;
    V(3) = 104;

    Eigen::MatrixXd A = Laguerre(2,V);
    std::cout << A << std::endl;*/

/*  Eigen::VectorXd v;
    v = GBM(new_M, new_N, new_T, new_r, new_q, new_sigma, new_S0);
    std::cout << v << std::endl;*/


/*  Eigen::VectorXd S(3);
    S(0) = 101;
    S(1) = 102;
    S(2) = 105;
    S = Findcallpath(S,102);
    std::cout << S << std::endl;*/




}

LSM::~LSM(){

}

Eigen::MatrixXd LSM::Laguerre(Eigen::VectorXd X, const int R){
    int n = X.rows();
        int m = R + 1;
        Eigen::MatrixXd value(n, m);
        for(int i = 0; i < n; i++){
            for(int j = 0; j < m; j++){
                if(R == 1){
                    value(i,0) = 1.0;
                    value(i,1) = -X(i) + 1.0;
                }
                else if(R == 2){
                    value(i,0) = 1.0;
                    value(i,1) = -X(i) + 1.0;
                    value(i,2) = 1.0/2.0*(2 - 4*X(i) + X(i)*X(i));
                }
                else if(R == 3){
                    value(i,0) = 1.0;
                    value(i,1) = -X(i) + 1.0;
                    value(i,2) = 1.0/2.0*(2 - 4*X(i) + X(i)*X(i));
                    value(i,3) = 1.0/6.0*(6.0 - 18.0*X(i,0) + 9.0*X(i)*X(i) - pow((double)X(i,0),3.0));
                }
            }
        }
        return value;
}

Eigen::VectorXd LSM::GBM(const int M, const int N, const double T, const double r, const double q, const double sigma, const double S0){
    double dt = T/N;
    Eigen::VectorXd Z(M);
    Eigen::VectorXd S(M);
    S(0) = S0;
    std::mt19937 e2(time(0));
    std::normal_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i < M; i++){
        Z(i) = dist(e2);
    }
    double drift  = exp(dt*((r - q)-0.5*sigma*sigma));
    double vol = sqrt(sigma*sigma*dt);
    for(int i = 1; i < M; i++){
        S(i) = S(i-1) * drift * exp(vol * Z(i));
    }
    return S;
}

Eigen::VectorXd LSM::callPayoff(Eigen::VectorXd S, const double K){
    Eigen::VectorXd C(S.size());
        for(int i = 0; i < S.size(); i++){
            if(S(i) - K > 0){
                C(i) = S(i) - K;
            }else{
                C(i) = 0.0;
            }
        }
        return C;
}

Eigen::VectorXd LSM::putPayoff(Eigen::VectorXd S, const double K){
    Eigen::VectorXd P(S.size());
        for(int i = 0; i < S.size(); i++){
            if(K - S(i) > 0){
                P(i) = K - S(i);
            }else{
                P(i) = 0.0;
            }
        }
        return P;
}


Eigen::VectorXd LSM::Findcallpath(Eigen::VectorXd S, const double K){
    Eigen::VectorXd path(S.size());
    int count = 0;
    for(int i = 0; i < S.size(); i++){
        if(S(i) - K > 0){
            path(count) = i;
            count++;
        }
    }
    path.conservativeResize(count);
    return path;
}
Eigen::VectorXd LSM::Findputpath(Eigen::VectorXd S, const double K){
    Eigen::VectorXd path(S.size());
    int count = 0;
    for(int i = 0; i < S.size(); i++){
        if(K - S(i) > 0){
            path(count) = i;
            count++;
        }
    }
    path.conservativeResize(count);
    return path;
}

Eigen::VectorXd Findcallprices(Eigen::VectorXd path, Eigen::VectorXd S){
    Eigen::VectorXd C(path.size());
    for(int i = 0; i < path.size(); i++){
        C(i) = S(path(i));
    }
    return C;
}

Eigen::VectorXd Findputprices(Eigen::VectorXd path, Eigen::VectorXd S){
    Eigen::VectorXd P(path.size());
    for(int i = 0; i < path.size(); i++){
        P(i) = S(path(i));
    }
    return P;
}

Eigen::VectorXd LSM::Jacobi(Eigen::MatrixXd L, Eigen::VectorXd Y){
    Eigen::VectorXd reg(L.rows());
    return reg = L.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
}

Eigen::VectorXd LSM::Findcallreturn(Eigen::VectorXd S, const double K){
    Eigen::VectorXd C_return(S.size());
    for(int i = 0; i < S.size; i++){
        C_return(i) = (S(i) - K);
    }
    return C_return;
}

Eigen::VectorXd LSM::Findputreturn(Eigen::VectorXd S, const double K){
    Eigen::VectorXd P_return(S.size());
    for(int i = 0; i < S.size; i++){
        P_return(i) = (K - S(i));
    }
    return P_return;
}

Answer 1

Here are some general coding style tips and improvements you could (and should!) make. I did not actually review the inner workings of your code since I am neither familiar with Eigen nor with the algorithm you are trying to implement.

1. LSM(const double, const double, const double, const int, const int, const double, const double, const int, const int); give names to your parameters. As is, this line is unreadable because one cannot know which parameter has which meaning.
2. Do not define a destructor if it does not actually do anything. The compiler generates the empty constructor for you anyway.
3. Put your #includes inside your header guard. There is no sense in making compiler and preprocessor do additional work if you do not actually use these headers.
4. Order your #includes. As a rule of thumb, begin with the header this source file implements (if in a .cpp file), then all headers from the same project you are working on, then headers from other libraries and finally headers from the STL. This ensures that all your headers are self-contained, i.e. include all the necessary headers themselves. Also, you should sort your #includes (in their respective groups, of course) alphabetically to make it easier to check whether some header is actually included or not.
5. Do not include C standard library headers directly. All of those headers have a C++ version with a prefixed c that brings the definitions into the std as opposed to the global namespace (i.e., #include <cmath> instead of #include <math.h> and prefix everything from that header with std::).
6. Why do you prefix all your variable names with new_? If this has to do something with name clashes with parameters to your constructor, you should instead use the this->variable = variable; variant or, even better, if using C++11 or later, a member initialization list instead.
7. Choose better variable names. One letter variable names are high up on the list of things that hurt readability the most, because knowing what a variable contains and what it stands for is often crucial to understanding code. I realize that some of your variables represent mathematical symbols, but those who do not should definitely have a more descriptive name.
8. Make your naming scheme consistent. For example, why do you have some variables with a capital letter and some without? Why do some methods follow CamelCase and others a plain capitalize-the-first-letter-and-nothing-else-scheme? There is no general consensus as to which naming style is best, so you are free to pick one, but please stay consistent (and do not pick a style that is hard to read, e. g. all lowercase).

9. Adopt a sensible indentation scheme. Again, this is something where a lot of programmers do something slightly different, but there are some general rules that most people follow. Most importantly, you should always indent when creating a new block, which you already obey. However, when you are not creating a new block, there are seldom reasons to increase indentation. In particular, lines such as

   Eigen::VectorXd P(S.size());
        for(int i = 0; i < S.size(); i++){
   

   are confusing because of the wrong indent.

10. Leave some horizontal space. Generally, it is a good idea to leave a space around binary operators (which you already do most of the time, but not always). Also, Most people prefer to have a space between a control structure keyword and its condition block and between a closing parens and a following opening curly brace as well (e.g. if (R == 1) { instead of if(R == 1){).

11. In Jacobi, why do you have a variable reg that you only ever write to? You could just write return L.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y); and be done.