Skip to content

Commit cc439d3

Browse files
pcaspersjenkins
pcaspers
authored and
jenkins
committed
QPR-12028 add the option to use smoothed indicators / max / min in forward eval
1 parent b249bc0 commit cc439d3

4 files changed

Lines changed: 86 additions & 24 deletions

File tree

‎QuantExt/qle/math/randomvariable.cpp‎

Lines changed: 50 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -26,9 +26,9 @@
2626
#include <boost/math/distributions/normal.hpp>
2727

2828
#include <boost/accumulators/accumulators.hpp>
29+
#include <boost/accumulators/statistics/mean.hpp>
2930
#include <boost/accumulators/statistics/stats.hpp>
3031
#include <boost/accumulators/statistics/variance.hpp>
31-
#include <boost/accumulators/statistics/mean.hpp>
3232

3333
#include <iostream>
3434
#include <map>
@@ -75,6 +75,16 @@ inline void stopCalcStats(const std::size_t n) {}
7575

7676
#endif
7777

78+
double getDelta(const RandomVariable& x, const Real eps) {
79+
Real sum = 0.0;
80+
for (Size i = 0; i < x.size(); ++i) {
81+
sum += x[i] * x[i];
82+
}
83+
Real delta = std::sqrt(sum / static_cast<Real>(x.size())) * eps / 2.0;
84+
// std::cout << "delta = " << delta << std::endl;
85+
return delta;
86+
}
87+
7888
} // namespace
7989

8090
Filter::~Filter() { clear(); }
@@ -891,14 +901,30 @@ RandomVariable indicatorEq(RandomVariable x, const RandomVariable& y, const Real
891901
return x;
892902
}
893903

894-
RandomVariable indicatorGt(RandomVariable x, const RandomVariable& y, const Real trueVal, const Real falseVal) {
904+
RandomVariable indicatorGt(RandomVariable x, const RandomVariable& y, const Real trueVal, const Real falseVal,
905+
const Real eps) {
895906
if (!x.initialised() || !y.initialised())
896907
return RandomVariable();
897908
QL_REQUIRE(x.size() == y.size(), "RandomVariable: indicatorEq(x,y): x size ("
898909
<< x.size() << ") must be equal to y size (" << y.size() << ")");
899910
x.checkTimeConsistencyAndUpdate(y.time());
900911
if (!y.deterministic_)
901912
x.expand();
913+
if (eps != 0.0) {
914+
Real delta = getDelta(x - y, eps);
915+
if (!QuantLib::close_enough(delta, 0.0)) {
916+
// logistic function
917+
x.expand();
918+
Real delta = getDelta(x - y, eps);
919+
resumeCalcStats();
920+
for (Size i = 0; i < x.n_; ++i) {
921+
x.data_[i] = falseVal + (trueVal - falseVal) * 1.0 / (1.0 + std::exp(-x.data_[i] / delta));
922+
}
923+
stopCalcStats(x.n_);
924+
return x;
925+
}
926+
}
927+
// eps == 0.0 or delta == 0.0
902928
if (x.deterministic()) {
903929
x.constantData_ =
904930
(x.constantData_ > y.constantData_ && !QuantLib::close_enough(x.constantData_, y.constantData_)) ? trueVal
@@ -913,14 +939,30 @@ RandomVariable indicatorGt(RandomVariable x, const RandomVariable& y, const Real
913939
return x;
914940
}
915941

916-
RandomVariable indicatorGeq(RandomVariable x, const RandomVariable& y, const Real trueVal, const Real falseVal) {
942+
RandomVariable indicatorGeq(RandomVariable x, const RandomVariable& y, const Real trueVal, const Real falseVal,
943+
const Real eps) {
917944
if (!x.initialised() || !y.initialised())
918945
return RandomVariable();
919946
QL_REQUIRE(x.size() == y.size(), "RandomVariable: indicatorEq(x,y): x size ("
920947
<< x.size() << ") must be equal to y size (" << y.size() << ")");
921948
x.checkTimeConsistencyAndUpdate(y.time());
922949
if (!y.deterministic_)
923950
x.expand();
951+
if (eps != 0.0) {
952+
Real delta = getDelta(x - y, eps);
953+
if (!QuantLib::close_enough(delta, 0.0)) {
954+
// logistic function
955+
x.expand();
956+
Real delta = getDelta(x - y, eps);
957+
resumeCalcStats();
958+
for (Size i = 0; i < x.n_; ++i) {
959+
x.data_[i] = falseVal + (trueVal - falseVal) * 1.0 / (1.0 + std::exp(-x.data_[i] / delta));
960+
}
961+
stopCalcStats(x.n_);
962+
return x;
963+
}
964+
}
965+
// eps == 0.0 or delta == 0.0
924966
if (x.deterministic()) {
925967
x.constantData_ =
926968
(x.constantData_ > y.constantData_ || QuantLib::close_enough(x.constantData_, y.constantData_)) ? trueVal
@@ -1206,33 +1248,27 @@ RandomVariable indicatorDerivative(const RandomVariable& x, const double eps) {
12061248
// Fries, 2017: Automatic Backward Differentiation for American Monte-Carlo Algorithms -
12071249
// ADD for Conditional Expectations and Indicator Functions
12081250

1209-
if (QuantLib::close_enough(eps, 0.0) || x.deterministic())
1210-
return tmp;
1211-
12121251
resumeCalcStats();
12131252

1214-
Real sum = 0.0;
1215-
for (Size i = 0; i < x.size(); ++i) {
1216-
sum += x[i] * x[i];
1217-
}
1218-
1219-
Real delta = std::sqrt(sum / static_cast<Real>(x.size())) * eps / 2.0;
1253+
Real delta = getDelta(x, eps);
12201254

12211255
if (QuantLib::close_enough(delta, 0.0))
12221256
return tmp;
12231257

12241258
// compute derivative
12251259

12261260
for (Size i = 0; i < tmp.size(); ++i) {
1227-
Real ax = std::abs(x[i]);
12281261

12291262
// linear approximation of step
1263+
// Real ax = std::abs(x[i]);
12301264
// if (ax < delta) {
12311265
// tmp.set(i, 1.0 / (2.0 * delta));
12321266
// }
12331267

12341268
// logistic function
1235-
tmp.set(i, std::exp(-1.0 / delta * ax) / (delta * std::pow(1.0 + std::exp(-1.0 / delta * ax), 2.0)));
1269+
// f(x) = 1 / (1 + exp(-x / delta))
1270+
// f'(x) = exp(-x/delta) / (delta * (1 + exp(-x / delta))^2), this is an even function
1271+
tmp.set(i, std::exp(-1.0 / delta * x[i]) / (delta * std::pow(1.0 + std::exp(-1.0 / delta * x[i]), 2.0)));
12361272
}
12371273

12381274
stopCalcStats(x.size() * 8);

‎QuantExt/qle/math/randomvariable.hpp‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -204,8 +204,10 @@ struct RandomVariable {
204204
friend RandomVariable applyInverseFilter(RandomVariable, const Filter&);
205205
friend RandomVariable conditionalResult(const Filter&, RandomVariable, const RandomVariable&);
206206
friend RandomVariable indicatorEq(RandomVariable, const RandomVariable&, const Real trueVal, const Real falseVal);
207-
friend RandomVariable indicatorGt(RandomVariable, const RandomVariable&, const Real trueVal, const Real falseVal);
208-
friend RandomVariable indicatorGeq(RandomVariable, const RandomVariable&, const Real trueVal, const Real falseVal);
207+
friend RandomVariable indicatorGt(RandomVariable, const RandomVariable&, const Real trueVal, const Real falseVal,
208+
const Real eps);
209+
friend RandomVariable indicatorGeq(RandomVariable, const RandomVariable&, const Real trueVal, const Real falseVal,
210+
const Real eps);
209211

210212
void expand();
211213
// pointer to raw data, this is null for deterministic variables
@@ -257,8 +259,10 @@ RandomVariable cos(RandomVariable);
257259
RandomVariable normalCdf(RandomVariable);
258260
RandomVariable normalPdf(RandomVariable);
259261
RandomVariable indicatorEq(RandomVariable, const RandomVariable&, const Real trueVal = 1.0, const Real falseVal = 0.0);
260-
RandomVariable indicatorGt(RandomVariable, const RandomVariable&, const Real trueVal = 1.0, const Real falseVal = 0.0);
261-
RandomVariable indicatorGeq(RandomVariable, const RandomVariable&, const Real trueVal = 1.0, const Real falseVal = 0.0);
262+
RandomVariable indicatorGt(RandomVariable, const RandomVariable&, const Real trueVal = 1.0, const Real falseVal = 0.0,
263+
const Real eps = 0.0);
264+
RandomVariable indicatorGeq(RandomVariable, const RandomVariable&, const Real trueVal = 1.0, const Real falseVal = 0.0,
265+
const Real eps = 0.0);
262266

263267
Filter close_enough(const RandomVariable&, const RandomVariable&);
264268
bool close_enough_all(const RandomVariable&, const RandomVariable&);

‎QuantExt/qle/math/randomvariable_ops.cpp‎

Lines changed: 24 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,8 @@
2525
namespace QuantExt {
2626

2727
std::vector<RandomVariableOp> getRandomVariableOps(const Size size, const Size regressionOrder,
28-
QuantLib::LsmBasisSystem::PolynomialType polynomType) {
28+
QuantLib::LsmBasisSystem::PolynomialType polynomType,
29+
const double eps) {
2930
std::vector<RandomVariableOp> ops;
3031

3132
// None = 0
@@ -65,16 +66,34 @@ std::vector<RandomVariableOp> getRandomVariableOps(const Size size, const Size r
6566
ops.push_back([](const std::vector<const RandomVariable*>& args) { return indicatorEq(*args[0], *args[1]); });
6667

6768
// IndicatorGt = 8
68-
ops.push_back([](const std::vector<const RandomVariable*>& args) { return indicatorGt(*args[0], *args[1]); });
69+
ops.push_back([eps](const std::vector<const RandomVariable*>& args) {
70+
return indicatorGt(*args[0], *args[1], 1.0, 0.0, eps);
71+
});
6972

7073
// IndicatorGeq = 9
71-
ops.push_back([](const std::vector<const RandomVariable*>& args) { return indicatorGeq(*args[0], *args[1]); });
74+
ops.push_back([eps](const std::vector<const RandomVariable*>& args) {
75+
return indicatorGeq(*args[0], *args[1], 1.0, 0.0, eps);
76+
});
7277

7378
// Min = 10
74-
ops.push_back([](const std::vector<const RandomVariable*>& args) { return QuantExt::min(*args[0], *args[1]); });
79+
if (eps == 0.0) {
80+
ops.push_back(
81+
[](const std::vector<const RandomVariable*>& args) { return QuantExt::min(*args[0], *args[1]); });
82+
} else {
83+
ops.push_back([eps](const std::vector<const RandomVariable*>& args) {
84+
return indicatorGt(*args[0], *args[1], 1.0, 0.0, eps) * (*args[1] - *args[0]) + *args[0];
85+
});
86+
}
7587

7688
// Max = 11
77-
ops.push_back([](const std::vector<const RandomVariable*>& args) { return QuantExt::max(*args[0], *args[1]); });
89+
if (eps == 0.0) {
90+
ops.push_back(
91+
[](const std::vector<const RandomVariable*>& args) { return QuantExt::max(*args[0], *args[1]); });
92+
} else {
93+
ops.push_back([eps](const std::vector<const RandomVariable*>& args) {
94+
return indicatorGt(*args[0], *args[1], 1.0, 0.0, eps) * (*args[0] - *args[1]) + *args[1];
95+
});
96+
}
7897

7998
// Abs = 12
8099
ops.push_back([](const std::vector<const RandomVariable*>& args) { return QuantExt::abs(*args[0]); });

‎QuantExt/qle/math/randomvariable_ops.hpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,15 +35,18 @@ namespace QuantExt {
3535

3636
using RandomVariableOp = std::function<RandomVariable(const std::vector<const RandomVariable*>&)>;
3737

38+
// eps determines the smoothing, 0 means no smoothing (default)
3839
std::vector<RandomVariableOp>
3940
getRandomVariableOps(const Size size, const Size regressionOrder = 2,
40-
const QuantLib::LsmBasisSystem::PolynomialType polynomType = QuantLib::LsmBasisSystem::Monomial);
41+
const QuantLib::LsmBasisSystem::PolynomialType polynomType = QuantLib::LsmBasisSystem::Monomial,
42+
const double eps = 0.0);
4143

4244
// random variable gradients
4345

4446
using RandomVariableGrad =
4547
std::function<std::vector<RandomVariable>(const std::vector<const RandomVariable*>&, const RandomVariable*)>;
4648

49+
// eps determines the smoothing, 0 means no smoothing
4750
std::vector<RandomVariableGrad> getRandomVariableGradients(
4851
const Size size, const Size regressionOrder = 2,
4952
const QuantLib::LsmBasisSystem::PolynomialType polynomType = QuantLib::LsmBasisSystem::Monomial,

0 commit comments

Comments
 (0)