smilesectionutils.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2013, 2015, 2018 Peter Caspers
 
 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/
 
 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.
 
 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/
 
#include <ql/termstructures/volatility/smilesectionutils.hpp>
#include <ql/math/comparison.hpp>
#include <algorithm>
 
namespace QuantLib {
 
    SmileSectionUtils::SmileSectionUtils(const SmileSection &section,
                                         const std::vector<Real> &moneynessGrid,
                                         const Real atm,
                                         const bool deleteArbitragePoints) {
 
        if (!moneynessGrid.empty()) {
            QL_REQUIRE(
                section.volatilityType() == Normal || moneynessGrid[0] >= 0.0,
                "moneyness grid should only contain non negative values ("
                    << moneynessGrid[0] << ")");
            for (Size i = 0; i < moneynessGrid.size() - 1; i++) {
                QL_REQUIRE(moneynessGrid[i] < moneynessGrid[i + 1],
                           "moneyness grid should contain strictly increasing "
                           "values ("
                               << moneynessGrid[i] << ","
                               << moneynessGrid[i + 1] << " at indices " << i
                               << ", " << i + 1 << ")");
            }
        }
 
        if (atm == Null<Real>()) {
            f_ = section.atmLevel();
            QL_REQUIRE(f_ != Null<Real>(),
                       "atm level must be provided by source section or given "
                       "in the constructor");
        } else {
            f_ = atm;
        }
 
        std::vector<Real> tmp;
 
        static const Real defaultMoney[] = { 0.0,  0.01, 0.05, 0.10, 0.25, 0.40,
                                             0.50, 0.60, 0.70, 0.80, 0.90, 1.0,
                                             1.25, 1.5,  1.75, 2.0,  5.0,  7.5,
                                             10.0, 15.0, 20.0 };
        static const Real defaultMoneyNormal[] = {
            -0.20,  -0.15,  -0.10,  -0.075,  -0.05,   -0.04,   -0.03,
            -0.02,  -0.015, -0.01,  -0.0075, -0.0050, -0.0025, 0.0,
            0.0025, 0.0050, 0.0075, 0.01,    0.015,   0.02,    0.03,
            0.04,   0.05,   0.075,  0.10,    0.15,    0.20
        };
 
        if (moneynessGrid.empty()) {
            tmp = section.volatilityType() == Normal
                      ? std::vector<Real>(defaultMoneyNormal,
                                          defaultMoneyNormal + 27)
                      : std::vector<Real>(defaultMoney, defaultMoney + 21);
        }
        else
            tmp = std::vector<Real>(moneynessGrid);
 
        Real shift = section.shift();
 
        if (section.volatilityType() == ShiftedLognormal && tmp[0] > QL_EPSILON) {
            m_.push_back(0.0);
            k_.push_back(-shift);
        }
 
        bool minStrikeAdded = false, maxStrikeAdded = false;
        for (Real& i : tmp) {
            Real k = section.volatilityType() == Normal ? Real(f_ + i) : Real(i * (f_ + shift) - shift);
            if ((section.volatilityType() == ShiftedLognormal && i <= QL_EPSILON) ||
                (k >= section.minStrike() && k <= section.maxStrike())) {
                if (!minStrikeAdded || !close(k, section.minStrike())) {
                    m_.push_back(i);
                    k_.push_back(k);
                }
                if (close(k, section.maxStrike()))
                    maxStrikeAdded = true;
            } else { // if the section provides a limited strike range
                     // we put the respective endpoint in our grid
                     // in order to not loose too much information
                if (k < section.minStrike() && !minStrikeAdded) {
                    m_.push_back(section.volatilityType() == Normal
                                     ? Real(section.minStrike() - f_)
                                     : Real((section.minStrike() + shift) / f_));
                    k_.push_back(section.minStrike());
                    minStrikeAdded = true;
                }
                if (k > section.maxStrike() && !maxStrikeAdded) {
                    m_.push_back(section.volatilityType() == Normal
                                     ? Real(section.maxStrike() - f_)
                                     : Real((section.maxStrike() + shift) / f_));
                    k_.push_back(section.maxStrike());
                    maxStrikeAdded = true;
                }
            }
        }
 
        // only known for shifted lognormal vols, otherwise we include
        // the lower strike in the loop below
        if(section.volatilityType() == ShiftedLognormal)
            c_.push_back(f_ + shift);
 
        for (Size i = (section.volatilityType() == Normal ? 0 : 1);
             i < k_.size(); i++) {
            c_.push_back(section.optionPrice(k_[i], Option::Call, 1.0));
        }
 
        Size centralIndex =
            std::upper_bound(m_.begin(), m_.end(),
                             (section.volatilityType() == Normal ? 0.0 : 1.0) -
                                 QL_EPSILON) -
            m_.begin();
        QL_REQUIRE(centralIndex < k_.size() - 1 && centralIndex > 1,
                   "Atm point in moneyness grid ("
                       << centralIndex << ") too close to boundary.");
 
        // shift central index to the right if necessary
        // (sometimes even the atm point lies in an arbitrageable area)
 
        while (!af(centralIndex, centralIndex, centralIndex + 1) &&
               centralIndex < k_.size() - 1)
            centralIndex++;
 
        QL_REQUIRE(centralIndex < k_.size(),
                   "central index is at right boundary");
 
        leftIndex_ = centralIndex;
        rightIndex_ = centralIndex;
 
        bool done = false;
        while (!done) {
 
            bool isAf = true;
            done = true;
 
            while (isAf && rightIndex_ < k_.size() - 1) {
                rightIndex_++;
                isAf = af(leftIndex_, rightIndex_, rightIndex_) &&
                       af(leftIndex_, rightIndex_ - 1, rightIndex_);
            }
            if (!isAf)
                rightIndex_--;
 
            isAf = true;
            while (isAf && leftIndex_ > 1) {
                leftIndex_--;
                isAf = af(leftIndex_, leftIndex_, rightIndex_) &&
                       af(leftIndex_, leftIndex_ + 1, rightIndex_);
            }
            if (!isAf)
                leftIndex_++;
 
            if (rightIndex_ < leftIndex_)
                rightIndex_ = leftIndex_;
 
            if (deleteArbitragePoints && leftIndex_ > 1) {
                m_.erase(m_.begin() + leftIndex_ - 1);
                k_.erase(k_.begin() + leftIndex_ - 1);
                c_.erase(c_.begin() + leftIndex_ - 1);
                leftIndex_--;
                rightIndex_--;
                done = false;
            }
            if (deleteArbitragePoints && rightIndex_ < k_.size() - 1) {
                m_.erase(m_.begin() + rightIndex_ + 1);
                k_.erase(k_.begin() + rightIndex_ + 1);
                c_.erase(c_.begin() + rightIndex_ + 1);
                rightIndex_--;
                done = false;
            }
        }
 
        QL_REQUIRE(rightIndex_ > leftIndex_,
                   "arbitrage free region must at least contain two "
                   "points (only index is "
                       << leftIndex_ << ")");
 
    }
 
    std::pair<Real, Real> SmileSectionUtils::arbitragefreeRegion() const {
        return {k_[leftIndex_], k_[rightIndex_]};
    }
 
    std::pair<Size, Size> SmileSectionUtils::arbitragefreeIndices() const {
        return {leftIndex_, rightIndex_};
    }
 
    bool SmileSectionUtils::af(const Size i0, const Size i,
                               const Size i1) const {
        if (i == 0)
            return true;
        Size im = i - 1 >= i0 ? i - 1 : 0;
        Real q1 = (c_[i] - c_[im]) / (k_[i] - k_[im]);
        if (q1 < -1.0 || q1 > 0.0)
            return false;
        if (i >= i1)
            return true;
        Real q2 = (c_[i + 1] - c_[i]) / (k_[i + 1] - k_[i]);
        return q1 <= q2 && q2 <= 0.0;
    }
}