#include <ql/quantlib.hpp>
using namespace QuantLib;
class ReplicationError
{
public:
ReplicationError(Option::Type type,
Time maturity,
double strike,
double s0,
double sigma,
Rate r)
: maturity_(maturity), payoff_(type, strike), s0_(s0),
sigma_(sigma), r_(r) {
EuropeanOption option = EuropeanOption(type, s0_, strike, 0.0, r_,
maturity_, sigma_);
std::cout << "Option value: " << option.value() << std::endl;
vega_ = option.vega();
std::cout << std::endl;
std::cout <<
" | | P&L \t| P&L | Derman&Kamal | P&L"
" \t| P&L" << std::endl;
std::cout <<
"samples | trades | Mean \t| Std Dev | Formula |"
" skewness \t| kurt." << std::endl;
std::cout << "---------------------------------"
"----------------------------------------------" << std::endl;
}
void compute(int nTimeSteps, int nSamples);
private:
Time maturity_;
PlainVanillaPayoff payoff_;
double s0_;
double sigma_;
Rate r_;
double vega_;
};
class ReplicationPathPricer : public PathPricer_old<Path>
{
public:
ReplicationPathPricer(Option::Type type,
double underlying,
double strike,
Rate r,
Time maturity,
double sigma)
: PathPricer_old<Path>(1.0, false), type_(type), underlying_(underlying),
strike_(strike), r_(r), maturity_(maturity), sigma_(sigma) {
QL_REQUIRE(strike_ > 0.0,
"ReplicationPathPricer: strike must be positive");
QL_REQUIRE(underlying_ > 0.0,
"ReplicationPathPricer: underlying must be positive");
QL_REQUIRE(r_ >= 0.0,
"ReplicationPathPricer: risk free rate (r) must"
" be positive or zero");
QL_REQUIRE(maturity_ > 0.0,
"ReplicationPathPricer: maturity must be positive");
QL_REQUIRE(sigma_ >= 0.0,
"ReplicationPathPricer: volatility (sigma)"
" must be positive or zero");
}
double operator()(const Path& path) const;
private:
Option::Type type_;
double underlying_, strike_;
Rate r_;
Time maturity_;
double sigma_;
};
int main(int argc, char* argv[])
{
try {
QL_IO_INIT
Time maturity = 1./12.;
double strike = 100;
double underlying = 100;
double volatility = 0.20;
Rate riskFreeRate = 0.05;
ReplicationError rp(Option::Call, maturity, strike, underlying,
volatility, riskFreeRate);
int scenarios = 50000;
int hedgesNum;
hedgesNum = 21;
rp.compute(hedgesNum, scenarios);
hedgesNum = 84;
rp.compute(hedgesNum, scenarios);
return 0;
} catch (std::exception& e) {
std::cout << e.what() << std::endl;
return 1;
} catch (...) {
std::cout << "unknown error" << std::endl;
return 1;
}
}
double ReplicationPathPricer::operator()(const Path& path) const
{
int n = path.size();
QL_REQUIRE(n>0,
"ReplicationPathPricer: the path cannot be empty");
Time dt = maturity_/n;
double stockDividendYield = 0.0;
Time t = 0;
double stock = underlying_;
double stockLogGrowth = 0.0;
double money_account = 0.0;
EuropeanOption option = EuropeanOption(type_, stock, strike_,
stockDividendYield, r_, maturity_, sigma_);
money_account += option.value();
double delta = option.delta();
double stockAmount = delta;
money_account -= stockAmount*stock;
for(int step = 0; step < n-1; step++){
t += dt;
money_account *= QL_EXP( r_*dt );
stockLogGrowth += path[step];
stock = underlying_*QL_EXP(stockLogGrowth);
EuropeanOption option = EuropeanOption(type_, stock, strike_,
stockDividendYield, r_, maturity_-t, sigma_);
delta = option.delta();
money_account -= (delta - stockAmount)*stock;
stockAmount = delta;
}
money_account *= QL_EXP( r_*dt );
stockLogGrowth += path[n-1];
stock = underlying_*QL_EXP(stockLogGrowth);
double optionPayoff = PlainVanillaPayoff(type_, strike_)(stock);
money_account -= optionPayoff;
money_account += stockAmount*stock;
return money_account;
}
void ReplicationError::compute(int nTimeSteps, int nSamples)
{
QL_REQUIRE(nTimeSteps>0,
"ReplicationError::compute : the number of steps must be > 0");
double drift = r_ - 0.5*sigma_*sigma_;
Handle<GaussianPathGenerator_old> myPathGenerator(
new GaussianPathGenerator_old(drift, sigma_*sigma_,
maturity_, nTimeSteps));
Handle<PathPricer_old<Path> > myPathPricer =
Handle<PathPricer_old<Path> >(new
ReplicationPathPricer(payoff_.optionType(), s0_,
payoff_.strike(), r_, maturity_, sigma_));
Statistics statisticsAccumulator;
OneFactorMonteCarloOption_old MCSimulation(myPathGenerator,
myPathPricer, statisticsAccumulator, false);
MCSimulation.addSamples(nSamples);
double PLMean = MCSimulation.sampleAccumulator().mean();
double PLStDev = MCSimulation.sampleAccumulator().standardDeviation();
double PLSkew = MCSimulation.sampleAccumulator().skewness();
double PLKurt = MCSimulation.sampleAccumulator().kurtosis();
double theorStD = QL_SQRT(M_PI/4/nTimeSteps)*vega_*sigma_;
std::cout << nSamples << "\t| "
<< nTimeSteps << "\t | "
<< DoubleFormatter::toString(PLMean, 3) << " \t| "
<< DoubleFormatter::toString(PLStDev, 2) << " \t | "
<< DoubleFormatter::toString(theorStD, 2) << " \t | "
<< DoubleFormatter::toString(PLSkew, 2) << " \t| "
<< DoubleFormatter::toString(PLKurt, 2) << std::endl;
}