Unknown errors within the C++ Eigen library for setting vectors/matrices

Question

I am trying to implement the Longstaff and Scwartz algorithm that essentially is a Least-Squares method to price American options. I am using Matlab code as a reference to program this in C++.

What I have done thus far is first I set a time step variable named

dt = T/N where T is the expiration time, and N is the number of time steps.

Then a time vector (in Matlab denoted t = 0:dt:T). I created a function for this in C++ that is posted below:

    VectorXd range(double min, double max, int N){
    VectorXd m(N + 1);
     double delta = (max-min)/N;
     for(int i = 0; i <= N; i++){
             m(i) = min + i*delta;
     }
    return m;
}

Then I have to generate a vector that generates a set of normally distributed random numbers. In Matlab the function is called z = randn(M/2,1) where M is the number of paths. In C++ I made a function for this which is listed below:

MatrixXd generateGaussianNoise(int n, int m){
    MatrixXd M(n,m);
    VectorXd V(n);
    normal_distribution<double> nd(0.0, 1.0);
    random_device rd;
    mt19937 gen(rd());
    for(int i = 0; i < n; i++){
        for(int j = 0; j < m; j++){
            M(i,j) = nd(gen);
        }
    }
    return M;
}

Now we have to create another vector called w = (r-sigma^2/2)*T + sigma*sqrt(T)*[z;z]. Before this I had to create another C++ function for handling the Matlab syntax [z;z] This is listed below:

MatrixXd generateMatrix(MatrixXd z){
    MatrixXd Dm(2*z.rows(),z.cols());
        Dm << z, -z;
        return Dm;
}

Then I can create the w variable with another C++ function I created listed below:

MatrixXd generateWMatrix(double r, double sigma, double T, MatrixXd Dm){
    MatrixXd w(Dm.rows(),Dm.cols());
    for(int i = 0; i < Dm.rows(); i++){
        for(int j = 0; j < Dm.cols(); j++){
            w(i,j) = (r - pow(sigma,2)/2)*T + sigma*sqrt(T)*Dm(i,j);
        }
    }
        return w;
}

Now we have to create a new variable S = S0*exp(w) where S0 is the initial asset price. We can do this in a for loop which is listed below:

MatrixXd S(w.rows(),w.cols());
        for(int i = 0; i < w.rows(); i++){
            for(int j = 0; j < w.cols(); j++){
                S(i,j) = exp(w(i,j));
            }
        }

Now this is where things for me get a bit ugly which is the Main Point of this post. We have to do a backward for loop for the actual algorithm. I will present the beginning part of the Matlab code for reference which is where things go wrong for me in C++:

for i = N:-1:2

    z = randn(M/2,1);
    w = t(i)*w/(t(i+1)) + sigma*sqrt(dt*t(i)/(t(i+1)))*[z;z];

Now in C++ what I try to do is this:

for(int i = N; i >= 2; i--){
      z(i) = generateGaussianNoise(M/2, 1);
      zz(i) = generateMatrix(z);
      w(i) = t(i)*w(i)/(t(i+1)) + sigma*sqrt(dt*t(i)/(t(i+1)))*zz(i);
}

Although I get the following errors for z(i):

cannot convert 'Eigen::MatrixXd {aka Eigen::Matrix<double, -1, -1>}' to 

'Eigen::DenseCoeffsBase<Eigen::Matrix<double, -1, -1>, 1>::Scalar {aka 
double}' in assignment

and for zz(i) I get the following error:

cannot convert 'Eigen::MatrixXd {aka Eigen::Matrix<double, -1, -1>}' to 

'Eigen::DenseCoeffsBase<Eigen::Matrix<double, -1, -1>, 1>::Scalar {aka double}' in assignment

and for w(i) I get the following error:

invalid arguments '
Candidates are:
double sqrt(double)
float sqrt(float)
long double sqrt(long double)
__gnu_cxx::__enable_if<74 0 __value 14 std::__is_integer 1 #074 0 __value 14 std::__is_integer 1 #0,double>::__type sqrt(#0)
const Eigen::CwiseUnaryOp<Eigen::internal::scalar_sqrt_op<#0::Scalar>,const #0> sqrt(const Eigen::ArrayBase<#0> &)
float sqrt(float)
long double sqrt(long double)
std::complex<#0> sqrt(const std::complex<#0> &)
__gnu_cxx::__enable_if<74 0 __value 14 std::__is_integer 1 #074 0 __value 14 std::__is_integer 1 #0,double>::__type 

I am not sure what I am doing wrong in the latter above or what these errors mean. Of course I have tried googling them to see what the problem is but I don't get any useful information. For completeness I will post my whole code below:

#include <iostream>
#include <cmath>
#include <math.h>
#include <limits>
#include <algorithm>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <random>

using namespace Eigen;
using namespace std;


double LaguerreExplicit(int R, double x); // Generates the (weighted) laguerre value
double payoff_Call(double S, double K); // Pay off of a call option
MatrixXd generateGaussianNoise(int n, int m); // Generates Normally distributed random numbers
double LSM(double T, double r, double sigma, double K, double S0, int N, int M, int R);
// T        Expiration time
// r        Riskless interest rate
// sigma    Volatility
// K        Strike price
// S0       Initial asset price
// N        Number of time steps
// M        Number of paths
// R        Number of basis functions

VectorXd range(double min, double max, int N);
MatrixXd generateMatrix(MatrixXd z);
MatrixXd generateWMatrix(double r, double sigma, double T, MatrixXd z);

int main(){
double r = 0.06;
double sigma = 0.25;
int T = 1.0;
int N = 2;
double dt = T/N;
    VectorXd t = range(0.0, T, N);

    MatrixXd result1 = generateGaussianNoise(2,1);


    MatrixXd result2 = generateMatrix(result1);

    MatrixXd w = generateWMatrix(r,sigma,T, result2);







}


double payoff_Call(double S, double K){
    double payoff;
    if((S - K) > 0)
    {
        payoff = S - K;
    }else
    {
        payoff = 0.0;
    }
    return payoff;
}

double LaguerreExplicit(int R, double x){
    double value;
    if(R==0)
    {
        value = 1;
    }
    else if(R==1)
    {
        value = 0.5*(pow(x,2) - 4.0*x + 2);
    }
    else if(R==3)
    {
        value = (1.0/6.0)*(-1*pow(x,3) + 9*pow(x,2) - 18*x + 6);
    }
    else if(R==4)
    {
        value = (1.0/24.0)*(pow(x,4) - 16*pow(x,3) + 72*pow(x,2) - 96*x + 24);
    }
    else if(R==5)
    {
        value = (1.0/120.0)*(-1*pow(x,5) + 25*pow(x,4) - 200*pow(x,3) + 600*pow(x,2) - 600*x + 120);
    }
    else if (R==6)
    {
        value = (1.0/720.0)*(pow(x,6) - 36*pow(x,5) + 450*pow(x,4) - 2400*pow(x,3) + 5400*pow(x,2) - 4320*x + 720);
    }
    else{
        cout << "Error!, R is out of range" << endl;
        value  = 0;
    }
    value = exp(-0.5*x)*value; // Weighted used in Longstaff-Scwartz
    return value;
}

MatrixXd generateGaussianNoise(int n, int m){
    MatrixXd M(n,m);
    normal_distribution<double> nd(0.0, 1.0);
    random_device rd;
    mt19937 gen(rd());
    for(int i = 0; i < n; i++){
        for(int j = 0; j < m; j++){
            M(i,j) = nd(gen);
        }
    }
    return M;
}

VectorXd range(double min, double max, int N){
    VectorXd m(N + 1);
     double delta = (max-min)/N;
     for(int i = 0; i <= N; i++){
             m(i) = min + i*delta;
     }
    return m;
}

MatrixXd generateMatrix(MatrixXd z){
    MatrixXd Dm(2*z.rows(),z.cols());
        Dm << z, -z;
        return Dm;
}


MatrixXd generateWMatrix(double r, double sigma, double T, MatrixXd Dm){
    MatrixXd w(Dm.rows(),Dm.cols());
    for(int i = 0; i < Dm.rows(); i++){
        for(int j = 0; j < Dm.cols(); j++){
            w(i,j) = (r - pow(sigma,2)/2)*T + sigma*sqrt(T)*Dm(i,j);
        }
    }
        return w;
}



double LSM(double T, double r, double sigma, double K, double S0, int N, int M, int R){
    double dt = T/N;
    VectorXd t = range(0,T,N);
    MatrixXd z = generateGaussianNoise(M/2, 1);
    MatrixXd zz = generateMatrix(z); // Need to do this for the [z;-z]
    MatrixXd w = generateWMatrix(r, sigma, T, zz);
    MatrixXd S(w.rows(),w.cols());
        for(int i = 0; i < w.rows(); i++){
            for(int j = 0; j < w.cols(); j++){
                S(i,j) = exp(w(i,j));
            }
        }
    MatrixXd P(S.rows(), S.cols());

    MatrixXd z_new(M,1);
    for(int i = N; i >= 2; i--){
        z_new.topRows(M/2) = generateGaussianNoise(M/2,1);
        z_new.bottomRows(M/2) = -z_new.topRows(M/2);
        w = t(i-1)*w/t(i) + sigma*sqrt((dt*t(i-1)/t(i)))*z_new;
    }







    return 0;

}

Update: I added the new full code but I am still getting some errors in the exp() and sqrt() in the function called double LSM(double T, double r, double sigma, double K, double S0, int N, int M, int R) I have tried to google what I could but I don't have any idea how to fix the following error messages:

For exp() I get this message:

Invalid arguments '
Candidates are:
double exp(double)
float exp(float)
long double exp(long double)
__gnu_cxx::__enable_if<74 0 __value 14 std::__is_integer 1 #074 0 __value 14 std::__is_integer 1 #0,double>::__type exp(#0)
const Eigen::CwiseUnaryOp<Eigen::internal::scalar_exp_op<#0::Scalar>,const #0> exp(const Eigen::ArrayBase<#0> &)
float exp(float)
long double exp(long double)
std::complex<#0> exp(const std::complex<#0> &)
__gnu_cxx::__enable_if<74 0 __value 14 std::__is_integer 1 #074 0 __value 14 std::__is_integer 1 #0,double>::__type 

For the sqrt() I get this message:

Invalid arguments '
Candidates are:
double sqrt(double)
float sqrt(float)
long double sqrt(long double)
__gnu_cxx::__enable_if<74 0 __value 14 std::__is_integer 1 #074 0 __value 14 std::__is_integer 1 #0,double>::__type sqrt(#0)
const Eigen::CwiseUnaryOp<Eigen::internal::scalar_sqrt_op<#0::Scalar>,const #0> sqrt(const Eigen::ArrayBase<#0> &)
float sqrt(float)
long double sqrt(long double)
std::complex<#0> sqrt(const std::complex<#0> &)
__gnu_cxx::__enable_if<74 0 __value 14 std::__is_integer 1 #074 0 __value 14 std::__is_integer 1 #0,double>::__type

Here is the Matlab code I am trying to implement in Eigen. I will just post the parts in the Matlab code that I have done successfully and the parts that I having a hard time implementing:

function u = LSM(T,r,sigma,K,S0,N,M,k)
% T Expiration time
% r Riskless interest rate
% sigma Volatility
% K Strike price
% S0 Initial asset price
% N Number of time steps
% M Number of paths
% k Number of basis functions
dt = T/N; % Time steps
t = 0:dt:T; % Time vector
z = randn(M/2,1);
w = (r-sigmaˆ2/2)*T + sigma*sqrt(T)*[z;-z];
S = S0*exp(w);
P = max(K-S,0); % Payoff at time T
for i = N:-1:2
z = randn(M/2,1);
w = t(i)*w/t(i+1) + sigma*sqrt(dt*t(i)/t(i+1))*[z;-z];

Everything after the for loop is where I am struggling...

Answer 1

You can almost literally translate this Matlab code

for i = N:-1:2
    z = randn(M/2,1);
    w = t(i)*w/t(i+1) + sigma*sqrt(dt*t(i)/t(i+1))*[z;-z];
end

to Eigen (assuming everything outside the loop works):

MatrixXd zz(M,1); // allocate space for [z; -z] only once
for(int i=N; i>=2; --i){
    zz.topRows(M/2) = generateGaussianNoise(M/2,1);
    zz.bottomRows(M/2) = -zz.topRows(M/2);

    w = t(i-1)*w/t(i) + sigma*std::sqrt(dt*t(i-1)/t(i))*zz;
}

Essentially, the only difference is that Matlab starts indexing at 1 and Eigen starts at 0.

Also concatenating [z;-z] is not possible "inlined" in Eigen, but I'm wondering why you do the same calculations for z and -z anyway ...

Answer 2

The following line does not make sense:

zz(i) = generateMatrix(z);

zz(i) is a double& whereas generateMatrix returns a MatrixXd. Same for w(i).