main.cpp

enum Style{
  EU = 1,
  AM = 2
};

enum Payoff{
  PUT = 1,
  CALL = 2
};

/**
 * Gauss Siedel tridiagonal solver
 * @param x_end
 * @param x_start
 * @param bound
 * @param a
 * @param b
 * @param c
 * @param d
 * @param n
 * @param w
 * @return
 */
template <class V1, class V2, class V3>
static inline double gs_tridiag_step(
        V1& x_end,
        const V2& bound,
        const xt::xtensor<double,1>& a,
        const xt::xtensor<double,1>& b,
        const xt::xtensor<double,1>& c,
        const V3& d,
        unsigned long n,
        double w)
{
    double err = 0.0;
    auto init = x_end(0);
    auto next_gs = (d(0)-b(0)*x_end(1))/a(0);
    auto next_bound = std::max((1-w)*init+w*next_gs, bound(0));
    x_end(0) = next_bound;
    err = std::max(err, std::abs(next_bound - init));
    for (auto i = 1; i<n; i++) {
        init = x_end(i);
        next_gs = (d(i)-c(i-1)*next_gs-b(i)*x_end(i+1))/a(i);
        next_bound = std::max((1-w)*init+w*next_gs, bound(i));
        x_end(i) = next_bound;
        err = std::max(err, std::abs(next_bound - init));
    }
    init = x_end(n);
    next_gs = (d(n)-c(n-1)*next_gs)/a(n);
    next_bound = std::max((1-w)*init+w*next_gs, bound(n));
    x_end(n) = next_bound;
    err = std::max(err, std::abs(next_bound - init));
    return err; //inf norm
}

BSFDResult fd(unsigned long M, double K, double r, double sigma, double div, double T, Payoff payoff, Style style,
        double rec_dx,
        double rec_dt)
{
    // Step sizes:
    double dt = 0.; // size of time steps
    double dx = 0.; // size of space steps that guarantee explicit stability
    if (rec_dt!=0 && rec_dx!=0) {
        dt = rec_dt;
        dx = rec_dx;
    }
    else {
        dt = T/((double) M-1);
        dx = sigma*sqrt(dt);
    }

    const double sigma_sq = gcem::sq(sigma);
    const double alpha = 0.5*sigma_sq;  // diffusivity in heat equation
    const double dn = alpha*dt/(dx*dx);  // diffusion number

    // size of the grid
    //TODO: heuristic for defining those
    const double L1 = -5;
    const double L2 = 1.5;
    const unsigned long N = ceil((L2-L1)/dx+1);

    xt::xtensor<double, 2> x = xt::linspace<double>(L1, L2, N).reshape({1, N});
    xt::xtensor<double, 2> tau = xt::linspace<double>(0, T, M).reshape({M, 1});


    xt::xtensor<double, 2> S = (K*xt::exp(x-(r-div-0.5*sigma_sq)*tau));
    auto S_init = xt::view(S, xt::all(), N-1, xt::newaxis());
    // Solution matrix, solve backwards where idx = 0 is exp, idx = N is today
    xt::xtensor<double, 2> U ({ M, N});
    auto U_first_col_v = xt::view(U, xt::all(), 0);
    auto S_first_col_v = xt::view(S, xt::all(), 0);
    auto S_first_row_v = xt::view(S, 0, xt::all());
    auto U_first_row_v = xt::view(U, 0, xt::all());

    // Boundary condition and exercice region
    xt::xtensor<double, 2> bound;
    if (payoff==CALL) {
        auto U_last_col_v = xt::view(U, xt::all(), N-1);
        if (style==EU) {
            U_last_col_v = xt::flatten(xt::exp(r*tau)*xt::maximum(xt::exp(-div*tau)*S_init-K*xt::exp(-r*tau), 0));
        }
        else {
            bound = exp(r*tau)*(xt::view(S, xt::all(), xt::range(1,N-1))-K);
            U_last_col_v = xt::flatten(xt::maximum(xt::maximum(S_init-K, 0)*xt::exp(r*tau),
                    xt::exp(r*tau)*xt::maximum(xt::exp(-div*tau)*S_init-K*xt::exp(-r*tau), 0)));
        }
        U_first_col_v = 0.0;
        U_first_row_v = xt::maximum(S_first_row_v-K, 0);
    }
    else {
        if (style==EU) {
            U_first_col_v = K-S_first_col_v*xt::exp(r*tau); // Merton Bound
        }
        else {
            bound = exp(r*tau)*(K-xt::view(S, xt::all(), xt::range(1,N-1)));
            U_first_col_v = (K-S_first_col_v)*xt::exp(r*tau);  // early exercise-need to double check for Calls!
        }
        auto U_last_col_v = xt::view(U, xt::all(), N-1);
        U_last_col_v = 0;
        U_first_row_v = xt::maximum(K-S_first_row_v, 0);
    }

    // Tridiagonal scheme equations
    const double theta_ = 1./2.;  //TODO: may want to alter this through time, in which case just place in loop below
  //TODO:  convert  to views lazy eval
    xt::xtensor<double, 1> main_diag({N-2}, 1. + 2. * theta_*dn);
    xt::xtensor<double, 1> off_diag({N-2}, -theta_*dn);

    xt::xtensor<double, 2> D = xt::xtensor<double, 2>::from_shape({{M, N-2}});
    const double cutoff = 0.01;
    const double w = 1.5; // step size

    for (long i = 1; i<M; i++) {
        auto d = xt::view(D, i-1, xt::all());
        for (long j = 1; j<N-1; j++) {
            d(j-1) = (1-theta_)*dn*U(i-1, j+1)+(1-2*(1-theta_)*dn)*U(i-1, j)+(
                    1-theta_)*dn*U(i-1, j-1);
        }
        d(0) = d(0)+dn*theta_*U(i, 0); // from boundary conditions
        d(N-3) = d(N-3)+dn*theta_*U(i, N-1);

        if (style==EU) {
            xt::xtensor<double, 1> U_tmp_v = xt::view(U, i, xt::range(1, N));
            //TODO: verify and test EU case
            algo::thomas_tridiag(main_diag, off_diag, off_diag, d, U_tmp_v);
        }
        else {
            // Solve for american in explicit SOR
            double err = std::numeric_limits<double>::infinity();
            auto step_result_v = xt::view(U, i, xt::range(1, N-1));
            step_result_v = d;
            auto current_bound = xt::view(bound, i, xt::all());
            while (err>cutoff) {
                err = gs_tridiag_step(step_result_v, current_bound, main_diag, off_diag, off_diag, d, N-3, w);
                //TODO: add convergence test
            }

        }
    }
    // put back into C out of U
    xt::xtensor<double, 2> C = xt::exp(-r*tau)*U;

    // Greeks: chain rule with (x, tau) cordinates, central diff for Delta and Gamma
    // delta
    auto dc_dx_v1 = xt::view(C, xt::all(), xt::range(2, _));
    auto dc_dx_v2 = xt::view(C, xt::all(), xt::range(0, N-2));
    xt::xtensor<double, 2> dc_dx = (dc_dx_v1-dc_dx_v2)/(2*dx);
    //TODO: avoid aliasing extra copy by folding it above

    auto dc_dx_v = xt::concatenate(
            xtuple(xt::view(dc_dx, xt::all(), 0, xt::newaxis()), dc_dx, xt::view(dc_dx, xt::all(), N-3, xt::newaxis())),
            1);
    xt::xtensor<double, 2> delta = dc_dx_v/S;
    // gamma
    xt::xtensor<double, 2> dc2_dx2 = xt::diff(C, 2, 1)/(dx*dx);
    auto dc2_dx2_v = xt::concatenate(xtuple(xt::view(dc2_dx2, xt::all(), 0, xt::newaxis()), dc2_dx2,
            xt::view(dc2_dx2, xt::all(), N-3, xt::newaxis())), 1);
    auto chain = dc2_dx2_v-dc_dx_v;
    xt::xtensor<double, 2> gamma = chain/(S*S);
    //TODO: why hardcode the threshold to 2
    filtration(gamma, S<2.0) = 0.;  // fix boundary errors

    // Numerical theta, Forward difference for Thetas
    xt::xtensor<double, 2> zeros = xt::zeros<double>({1L, (long) N});

    xt::xtensor<double, 2> dc_dtau =
            xt::concatenate(xtuple(zeros, xt::diff(C, 1, 0)), 0)/dt;  //TODO: stability check here for dt ~ 0 ?
    xt::xtensor<double, 2> theta = -dc_dx_v*(r-div-0.5*sigma_sq)-dc_dtau;
    xt::view(theta, 0, xt::all()) = 0;
    // Analytic theta
    //theta = -0.5 * (sigma_sq) * (S * S) * gamma  // look into effect of IR/Div

    // Numeric vega
    xt::xtensor<double, 2> vega;
    if (style==EU && (rec_dt==0 && rec_dx==0)) {
        vega = (sigma*tau)*(S*S*gamma);
    }
    else if (rec_dt==0 && rec_dx==0) {
        auto dvol = 0.05; //TODO: what's the correct value here should be % rather than absolute
        auto res1 = bs::fd(M, K, r, sigma+dvol, div, T, payoff, style, dx, dt);
        auto res2 = bs::fd(M, K, r, sigma-dvol, div, T, payoff, style, dx, dt);
        vega = (res1.C-res2.C-delta*(res1.S-res2.S)-0.5*gamma*(gcem::sq(res1.S-S)-gcem::sq(res2.S-S)))/(2*dvol);

        //faster(only one more loop) but less accurate
        //auto res1 = fd(M, K, r, div, T, sigma + dvol, style, payoff, (dx, dt))
        //vega = (res1.C - C - delta * (res1.S - S) - 0.5 * gamma * ((res1.S - S) ** 2)) / dvol
    }
    return BSFDResult{
            .S=std::move(S),
            .C=std::move(C),
            .delta=std::move(delta),
            .gamma=std::move(gamma),
            .theta=std::move(theta),
            .vega=std::move(vega),

    };
}

test.cpp

TEST(BS_FD, CallAmWithDivAndRf)
{

    double K = 100.0; // strike
    double r = 0.0194;  // LIBOR IR
    double d = 0.017;  // SPY div
    double T = 1;  // 1 year
    double sigma = 0.2;  // 20% annual vol
    auto M = 3 * 365; // grid steps
    Payoff payoff = CALL;
    Style style = AM;

    auto res = fd(M, K, r, sigma, d, T, payoff, style, 0, 0);

    auto shape = res.S.shape();
    auto rix = shape[0];
    auto cix = shape[1];
    EXPECT_DOUBLE_EQ(0.67379469990854668, res.S(0,0));
    EXPECT_DOUBLE_EQ(0.715907850860794, res.S(10,10));
    EXPECT_DOUBLE_EQ(0.685758458883008, res.S(rix-1,0));
    EXPECT_DOUBLE_EQ(4.481689070338065e+02, res.S(0,cix-1));
    EXPECT_DOUBLE_EQ(4.561265012154291e+02, res.S(rix-1,cix-1));
    EXPECT_DOUBLE_EQ(4.319054901783607e+02, res.S(rix-10,cix-10));

    EXPECT_DOUBLE_EQ(0., res.C(0,0));
    EXPECT_DOUBLE_EQ(0., res.C(10,10));
    EXPECT_DOUBLE_EQ(0., res.C(rix-1,0));
    EXPECT_DOUBLE_EQ(3.481689070338065e+02, res.C(0,cix-1));
    EXPECT_DOUBLE_EQ(3.561265012154291e+02, res.C(rix-1,cix-1));
    EXPECT_DOUBLE_EQ(3.319054901783607e+02, res.C(rix-10,cix-10));

    EXPECT_DOUBLE_EQ(0., res.delta(0,0));
    EXPECT_DOUBLE_EQ(0., res.delta(10,10));
    EXPECT_DOUBLE_EQ(1.7001614900423999e-129, res.delta(rix-1,0));
    EXPECT_DOUBLE_EQ(0.993940157374373, res.delta(0,cix-1));
    EXPECT_DOUBLE_EQ(0.99394015737439212, res.delta(rix-1,cix-1));
    EXPECT_DOUBLE_EQ(0.99996823414579061, res.delta(rix-10,cix-10));

    EXPECT_DOUBLE_EQ(0., res.gamma(0,0));
    EXPECT_DOUBLE_EQ(0., res.gamma(10,10));
    EXPECT_DOUBLE_EQ(0., res.gamma(rix-1,0));
    EXPECT_DOUBLE_EQ(-9.0719544277341756e-08, res.gamma(0,cix-1));
    EXPECT_DOUBLE_EQ(-8.9136842273589639e-08, res.gamma(rix-1,cix-1));
    EXPECT_DOUBLE_EQ(-9.4706563941069902e-08, res.gamma(rix-10,cix-10));

    EXPECT_DOUBLE_EQ(0., res.theta(0,0));
    EXPECT_DOUBLE_EQ(0., res.theta(10,10));
    EXPECT_DOUBLE_EQ(2.0519842169441393e-131, res.theta(rix-1,0));
    EXPECT_DOUBLE_EQ(0., res.theta(0,cix-1));
    EXPECT_DOUBLE_EQ(-0.048582790254290131, res.theta(rix-1,cix-1));
    EXPECT_DOUBLE_EQ(-0.00018032376364196523, res.theta(rix-10,cix-10));

    EXPECT_DOUBLE_EQ(0., res.vega(0,0));
    EXPECT_DOUBLE_EQ(0., res.vega(10,10));
    EXPECT_DOUBLE_EQ(-2.3347557314554352e-130, res.vega(rix-1,0));
    EXPECT_DOUBLE_EQ(0., res.vega(0,cix-1));
    EXPECT_DOUBLE_EQ(0.5535164647594808, res.vega(rix-1,cix-1));
    EXPECT_DOUBLE_EQ(0.0027293413578043961, res.vega(rix-10,cix-10));
}