#include <ql/quantlib.hpp>
using namespace QuantLib;
class ReplicationError
{
public:
ReplicationError(Option::Type type,
Time maturity,
Real strike,
Real s0,
Volatility sigma,
Rate r)
: maturity_(maturity), payoff_(type, strike), s0_(s0),
sigma_(sigma), r_(r) {
DiscountFactor rDiscount = QL_EXP(-r_*maturity_);
DiscountFactor qDiscount = 1.0;
Real forward = s0_*qDiscount/rDiscount;
Real variance = sigma_*sigma_*maturity_;
boost::shared_ptr<StrikedTypePayoff> payoff(
new PlainVanillaPayoff(payoff_));
BlackFormula black(forward,rDiscount,variance,payoff);
std::cout << "Option value: " << black.value() << std::endl;
vega_ = black.vega(maturity_);
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(Size nTimeSteps, Size nSamples);
private:
Time maturity_;
PlainVanillaPayoff payoff_;
Real s0_;
Volatility sigma_;
Rate r_;
Real vega_;
};
class ReplicationPathPricer : public PathPricer<Path>
{
public:
ReplicationPathPricer(Option::Type type,
Real underlying,
Real strike,
Rate r,
Time maturity,
Volatility sigma)
: type_(type), underlying_(underlying),
strike_(strike), r_(r), maturity_(maturity), sigma_(sigma) {
QL_REQUIRE(strike_ > 0.0, "strike must be positive");
QL_REQUIRE(underlying_ > 0.0, "underlying must be positive");
QL_REQUIRE(r_ >= 0.0,
"risk free rate (r) must be positive or zero");
QL_REQUIRE(maturity_ > 0.0, "maturity must be positive");
QL_REQUIRE(sigma_ >= 0.0,
"volatility (sigma) must be positive or zero");
}
Real operator()(const Path& path) const;
private:
Option::Type type_;
Real underlying_, strike_;
Rate r_;
Time maturity_;
Volatility sigma_;
};
int main(int, char* [])
{
try {
QL_IO_INIT
Time maturity = 1.0/12.0;
Real strike = 100;
Real underlying = 100;
Volatility volatility = 0.20;
Rate riskFreeRate = 0.05;
ReplicationError rp(Option::Call, maturity, strike, underlying,
volatility, riskFreeRate);
Size scenarios = 50000;
Size 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;
}
}
Real ReplicationPathPricer::operator()(const Path& path) const
{
Size n = path.size();
QL_REQUIRE(n>0, "the path cannot be empty");
Time dt = maturity_/n;
Rate stockDividendYield = 0.0;
Time t = 0;
Real stock = underlying_;
Real stockLogGrowth = 0.0;
Real money_account = 0.0;
DiscountFactor rDiscount = QL_EXP(-r_*maturity_);
DiscountFactor qDiscount = QL_EXP(-stockDividendYield*maturity_);
Real forward = stock*qDiscount/rDiscount;
Real variance = sigma_*sigma_*maturity_;
boost::shared_ptr<StrikedTypePayoff> payoff(
new PlainVanillaPayoff(type_,strike_));
BlackFormula black(forward,rDiscount,variance,payoff);
money_account += black.value();
Real delta = black.delta(stock);
Real stockAmount = delta;
money_account -= stockAmount*stock;
for (Size step = 0; step < n-1; step++){
t += dt;
money_account *= QL_EXP( r_*dt );
stockLogGrowth += path[step];
stock = underlying_*QL_EXP(stockLogGrowth);
rDiscount = QL_EXP(-r_*(maturity_-t));
qDiscount = QL_EXP(-stockDividendYield*(maturity_-t));
forward = stock*qDiscount/rDiscount;
variance = sigma_*sigma_*(maturity_-t);
BlackFormula black(forward,rDiscount,variance,payoff);
delta = black.delta(stock);
money_account -= (delta - stockAmount)*stock;
stockAmount = delta;
}
money_account *= QL_EXP( r_*dt );
stockLogGrowth += path[n-1];
stock = underlying_*QL_EXP(stockLogGrowth);
Real optionPayoff = PlainVanillaPayoff(type_, strike_)(stock);
money_account -= optionPayoff;
money_account += stockAmount*stock;
return money_account;
}
void ReplicationError::compute(Size nTimeSteps, Size nSamples)
{
QL_REQUIRE(nTimeSteps>0, "the number of steps must be > 0");
Date today = Date::todaysDate();
DayCounter dayCount = Actual365Fixed();
Handle<Quote> stateVariable(
boost::shared_ptr<Quote>(new SimpleQuote(s0_)));
Handle<YieldTermStructure> riskFreeRate(
boost::shared_ptr<YieldTermStructure>(
new FlatForward(today, r_, dayCount)));
Handle<YieldTermStructure> dividendYield(
boost::shared_ptr<YieldTermStructure>(
new FlatForward(today, 0.0, dayCount)));
Handle<BlackVolTermStructure> volatility(
boost::shared_ptr<BlackVolTermStructure>(
new BlackConstantVol(today, sigma_, dayCount)));
boost::shared_ptr<StochasticProcess> diffusion(
new BlackScholesProcess(stateVariable, dividendYield,
riskFreeRate, volatility));
PseudoRandom::rsg_type rsg =
PseudoRandom::make_sequence_generator(nTimeSteps, 0);
bool brownianBridge = false;
typedef SingleAsset<PseudoRandom>::path_generator_type generator_type;
boost::shared_ptr<generator_type> myPathGenerator(new
generator_type(diffusion, maturity_, nTimeSteps,
rsg, brownianBridge));
boost::shared_ptr<PathPricer<Path> > myPathPricer(new
ReplicationPathPricer(payoff_.optionType(), s0_,
payoff_.strike(), r_, maturity_, sigma_));
Statistics statisticsAccumulator;
OneFactorMonteCarloOption MCSimulation(myPathGenerator,
myPathPricer,
statisticsAccumulator,
false);
MCSimulation.addSamples(nSamples);
Real PLMean = MCSimulation.sampleAccumulator().mean();
Real PLStDev = MCSimulation.sampleAccumulator().standardDeviation();
Real PLSkew = MCSimulation.sampleAccumulator().skewness();
Real PLKurt = MCSimulation.sampleAccumulator().kurtosis();
Real theorStD = QL_SQRT(M_PI/4/nTimeSteps)*vega_*sigma_;
std::cout << nSamples << "\t| "
<< nTimeSteps << "\t | "
<< DecimalFormatter::toString(PLMean, 3) << " \t| "
<< DecimalFormatter::toString(PLStDev, 2) << " \t | "
<< DecimalFormatter::toString(theorStD, 2) << " \t | "
<< DecimalFormatter::toString(PLSkew, 2) << " \t| "
<< DecimalFormatter::toString(PLKurt, 2) << std::endl;
}