Search code examples
c++eigen

Using Eigen class to sum certain numbers in a vector


I am new to C++ and I am using the Eigen library. I was wondering if there was a way to sum certain elements in a vector. For example, say I have a vector that is a 100 by 1 and I just want to sum the first 10 elements. Is there a way of doing that using the Eigen library?

What I am trying to do is this: say I have a vector that is 1000 by 1 and I want to take the mean of the first 10 elements, then the next 10 elements, and so on and store that in some vector. Hence I will have a vector of size 100 of the averages. Any thoughts or suggestions are greatly appreciated.

Here is the beginning steps I have in my code. I have a S_temp4vector that is 1000 by 1. Now I intialize a new vector S_A that I want to have as the vector of the means. Here is my messy sloppy code so far: (Note that my question resides in the crudeMonteCarlo function)

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

using namespace Eigen;
using namespace std;

void crudeMonteCarlo(int N,double K, double r, double S0, double sigma, double T, int n);
VectorXd time_vector(double min, double max, int n);
VectorXd call_payoff(VectorXd S, double K);

int main(){
    int N = 100;
    double K = 100;
    double r = 0.2;
    double S0 = 100;
    double sigma = 0.4;
    double T = 0.1;
    int n = 10;

    crudeMonteCarlo(N,K,r,S0,sigma,T,n);

    return 0;
}

VectorXd time_vector(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 generateGaussianNoise(int M, int N){
    MatrixXd Z(M,N);
    static random_device rd;
    static mt19937 e2(time(0));
    normal_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i < M; i++){
        for(int j = 0; j < N; j++){
            Z(i,j) = dist(e2);
        }
    }
    return Z;
}

VectorXd call_payoff(VectorXd S, double K){
    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;
}

void crudeMonteCarlo(int N,double K, double r, double S0, double sigma, double T, int n){
    // Create time vector
    VectorXd tt = time_vector(0.0,T,n);
    VectorXd t(n);
    double dt = T/n;
    for(int i = 0; i < n; i++){
            t(i) = tt(i+1);
    }

    // Generate standard normal Z matrix
    //MatrixXd Z = generateGaussianNoise(N,n);

    // Generate the log normal stock process N times to get S_A for crude Monte Carlo
    MatrixXd SS(N,n+1);
    MatrixXd Z = generateGaussianNoise(N,n);
    for(int i = 0; i < N; i++){
        SS(i,0) = S0;
        for(int j = 1; j <= n; j++){
                SS(i,j) = SS(i,j-1)*exp((double) (r - pow(sigma,2.0))*dt + sigma*sqrt(dt)*(double)Z(i,j-1));
        }
    }

    // This long bit of code gives me my S_A.....
    Map<RowVectorXd> S_temp1(SS.data(), SS.size());
    VectorXd S_temp2(S_temp1.size());
    for(int i = 0; i < S_temp2.size(); i++){
        S_temp2(i) = S_temp1(i);
    }
    VectorXd S_temp3(S_temp2.size() - N);
    int count = 0;
    for(int i = N; i < S_temp2.size(); i++){
        S_temp3(count) = S_temp2(i);
        count++;
    }
    VectorXd S_temp4(S_temp3.size());
    for(int i = 0; i < S_temp4.size(); i++){
        S_temp4(i) = S_temp3(i);
    }
    VectorXd S_A(N);
    S_A(0) = (S_temp4(0) + S_temp4(1) + S_temp4(2) + S_temp4(3) + S_temp4(4) + S_temp4(5) + S_temp4(6) + S_temp4(7) + S_temp4(8) + S_temp4(9))/(n);
    S_A(1) = (S_temp4(10) + S_temp4(11) + S_temp4(12) + S_temp4(13) + S_temp4(14) + S_temp4(15) + S_temp4(16) + S_temp4(17) + S_temp4(18) + S_temp4(19))/(n);
    int count1 = 0;
    for(int i = 0; i < S_temp4.size(); i++){
        S_A(count1) = 
    }



    // Calculate payoff of Asian option
    //VectorXd call_fun = call_payoff(S_A,K);




}

Solution

  • This question includes a lot of code, which makes it hard to understand the question you're trying to ask. Consider including only the code specific to your question.

    In any case, you can use Eigen directly to do all of these things quite simply. In Eigen, Vectors are just matrices with 1 column, so all of the reasoning here is directly applicable to what you've written.

    const Eigen::Matrix<double, 100, 1> v = Eigen::Matrix<double, 100, 1>::Random();
    
    const int num_rows = 10;
    const int num_cols = 1;
    
    const int starting_row = 0;
    const int starting_col = 0;
    
    const double sum_of_first_ten = v.block(starting_row, starting_col, num_rows, num_cols).sum();
    const double mean_of_first_ten = sum_of_first_ten / num_rows;
    

    In summary: You can use .block to get a block object, .sum() to sum that block, and then conventional division to get the mean.