Login  Register

Numerical Greeks

Posted by DPaulino on Sep 18, 2015; 7:59pm
URL: http://quantlib.414.s1.nabble.com/Numerical-Greeks-tp16919.html

Hello,

I'm new to quantlib and have been trying to write a code to calculate a
numerical approximation to the greeks of an one-asset option.

For Delta, Gamma and Vega, what I have this far works great, but both rho
and theta are completely off the mark, even for an European vanilla.

I´ve read both this topic and this topic and as far as I can tell, it should work.

Can anyone give any clue as to what is not working?

The relevant part of the code is:


                Calendar calNACIONAL = Brazil();
                Calendar calBOVESPA = Brazil(Brazil::Exchange);
                Date todaysDate(18, Sep, 2015);
                Date maturity(16, Dec, 2015);
                DayCounter dayCounter = Business252();
                Natural settlementDays = 0;
                Settings::instance().evaluationDate() = todaysDate;

                Option::Type type
                type = Option::Call;

                Real underlying = 4869.00;
                Real strike = 5000.00;
                Spread dividendYield = 0.0;
                Rate riskFreeRate = .1427;
                Volatility volatility = .24679;
                Real dF=0,01
                Real dvol=0,0001
                Real dr=0,0001
 
               boost::shared_ptr<Exercise> europeanExercise(
                        new EuropeanExercise(maturity));

                // Create the quotes

                boost::shared_ptr<SimpleQuote> underlyingQ(new SimpleQuote(underlying));
                boost::shared_ptr<SimpleQuote> volatilityQ(new SimpleQuote(volatility));
                boost::shared_ptr<SimpleQuote> riskFreeRateQ(new SimpleQuote(riskFreeRate));
                boost::shared_ptr<SimpleQuote> dividendYieldQ(new SimpleQuote(dividendYield));
               

                // bootstrap the yield/dividend/vol curves
               
               
                boost::shared_ptr<YieldTermStructure>  flatTermStructure(
                        new FlatForward(settlementDays, calNACIONAL,
                        Handle<Quote>(riskFreeRateQ), dayCounter));
                Handle<YieldTermStructure> flatDividendTS(
                        boost::shared_ptr<YieldTermStructure>(
                        new FlatForward(settlementDays, calNACIONAL,
                        Handle<Quote>(dividendYieldQ), dayCounter)));
                 
                boost::shared_ptr<BlackVolTermStructure> flatVolTS(
                        new BlackConstantVol(settlementDays, calBOVESPA,
                           Handle<Quote>(volatilityQ),dayCounter));
                boost::shared_ptr<StrikedTypePayoff> payoff(
                        new PlainVanillaPayoff(type, strike));
                boost::shared_ptr<BlackProcess> bsmProcess(
                        new BlackProcess(Handle<Quote>(underlyingQ),
                                                         Handle<YieldTermStructure>(flatTermStructure),
                                                         Handle<BlackVolTermStructure>(flatVolTS)));

                // option
                VanillaOption europeanOption(payoff, europeanExercise);

                Real P;
                Real Pu;
                Real Pd;
               
                // Numerical Delta
                underlyingQ->setValue(underlying + dF);
                Pu = europeanOption.NPV();
                underlyingQ->setValue(underlying - dF);
                Pd = europeanOption.NPV();
                std::cout << "Delta = " << underlyingQ->value()*(Pu-Pd)/(2*dF)
                << std::endl; underlyingQ->setValue(underlying);

                // Numerical Gamma
                P = europeanOption.NPV();
                underlyingQ->setValue(underlying + dF);
                Pu = europeanOption.NPV();
                underlyingQ->setValue(underlying - dF);
                Pd = europeanOption.NPV();
                std::cout << "Gamma = " << underlyingQ->value()*underlyingQ->value()*(Pu +
                      Pd - 2 * P) / (100*dF*dF) << std::endl;underlyingQ->setValue(underlying);

                // Numerical Vega
                volatilityQ->setValue(volatility + dvol);
                Pu = europeanOption.NPV();
                volatilityQ->setValue(volatility - dvol);
                Pd = europeanOption.NPV();
                std::cout << "Vega = " << (Pu - Pd) / (200 * dvol) << std::endl;
                volatilityQ->setValue(volatility);
               
                //Numerical Rho
                riskFreeRateQ->setValue(riskFreeRate + dr);
                Pu = europeanOption.NPV();
                riskFreeRateQ->setValue(riskFreeRate - dr);
                Pd = europeanOption.NPV();
                std::cout << "Rho = " << (Pu - Pd) / (200 * dr) << std::endl;
                riskFreeRateQ->setValue(riskFreeRate);
               

                //Numerical Theta
                Settings::instance().evaluationDate() = todaysDate + 1;
                std::cout << "Theta = " << (europeanOption.NPV() - P) << std::endl;
                Settings::instance().evaluationDate() = todaysDate;

                std::cout << std::endl;