1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2010 Adrian O' Neill
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy of the license along with this program; if not, please email
12<quantlib-dev@lists.sf.net>. The license is also available online at
13<http://quantlib.org/license.shtml>.
14
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include <ql/exercise.hpp>
21#include <ql/experimental/variancegamma/fftengine.hpp>
22#include <ql/math/fastfouriertransform.hpp>
23#include <ql/math/interpolations/linearinterpolation.hpp>
24#include <complex>
25#include <utility>
26
27namespace QuantLib {
28
29 FFTEngine::FFTEngine(ext::shared_ptr<StochasticProcess1D> process, Real logStrikeSpacing)
30 : process_(std::move(process)), lambda_(logStrikeSpacing) {
31 registerWith(h: process_);
32 }
33
34 void FFTEngine::calculate() const
35 {
36 QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
37 "not an European Option");
38
39 ext::shared_ptr<StrikedTypePayoff> payoff =
40 ext::dynamic_pointer_cast<StrikedTypePayoff>(r: arguments_.payoff);
41 QL_REQUIRE(payoff, "non-striked payoff given");
42
43 auto r1 = resultMap_.find(x: arguments_.exercise->lastDate());
44 if (r1 != resultMap_.end())
45 {
46 auto r2 = r1->second.find(x: payoff);
47 if (r2 != r1->second.end())
48 {
49 results_.value = r2->second;
50 return;
51 }
52 }
53
54 // Option not precalculated - do entire FFT for one option. Not very efficient - call precalculate!
55 calculateUncached(payoff, exercise: arguments_.exercise);
56 }
57
58 void FFTEngine::update()
59 {
60 // Process has changed so cached values may no longer be correct
61 resultMap_.clear();
62
63 // Call base class implementation
64 VanillaOption::engine::update();
65 }
66
67 void FFTEngine::calculateUncached(const ext::shared_ptr<StrikedTypePayoff>& payoff,
68 const ext::shared_ptr<Exercise>& exercise) const {
69 ext::shared_ptr<VanillaOption> option(new VanillaOption(payoff, exercise));
70 std::vector<ext::shared_ptr<Instrument> > optionList;
71 optionList.push_back(x: option);
72
73 ext::shared_ptr<FFTEngine> tempEngine(clone().release());
74 tempEngine->precalculate(optionList);
75 option->setPricingEngine(tempEngine);
76 results_.value = option->NPV();
77 }
78
79 void FFTEngine::precalculate(const std::vector<ext::shared_ptr<Instrument> >& optionList) {
80 // Group payoffs by expiry date
81 // as with FFT we can compute a bunch of these at once
82 resultMap_.clear();
83
84 typedef std::vector<ext::shared_ptr<StrikedTypePayoff> > PayoffList;
85 typedef std::map<Date, PayoffList> PayoffMap;
86 PayoffMap payoffMap;
87
88 for (const auto& optIt : optionList) {
89 ext::shared_ptr<VanillaOption> option = ext::dynamic_pointer_cast<VanillaOption>(r: optIt);
90 QL_REQUIRE(option, "instrument must be option");
91 QL_REQUIRE(option->exercise()->type() == Exercise::European,
92 "not an European Option");
93
94 ext::shared_ptr<StrikedTypePayoff> payoff =
95 ext::dynamic_pointer_cast<StrikedTypePayoff>(r: option->payoff());
96 QL_REQUIRE(payoff, "non-striked payoff given");
97
98 payoffMap[option->exercise()->lastDate()].push_back(x: payoff);
99 }
100
101 std::complex<Real> i1(0, 1);
102 Real alpha = 1.25;
103
104 for (PayoffMap::const_iterator payIt = payoffMap.begin(); payIt != payoffMap.end(); ++payIt)
105 {
106 Date expiryDate = payIt->first;
107
108 // Calculate n large enough for maximum strike, and round up to a power of 2
109 Real maxStrike = 0.0;
110 for (const auto& payoff : payIt->second) {
111 if (payoff->strike() > maxStrike)
112 maxStrike = payoff->strike();
113 }
114 Real nR = 2.0 * (std::log(x: maxStrike) + lambda_) / lambda_;
115 Size log2_n = (static_cast<Size>((std::log(x: nR) / std::log(x: 2.0))) + 1);
116 Size n = static_cast<std::size_t>(1) << log2_n;
117
118 // Strike range (equation 19,20)
119 Real b = n * lambda_ / 2.0;
120
121 // Grid spacing (equation 23)
122 Real eta = 2.0 * M_PI / (lambda_ * n);
123
124 // Discount factor
125 Real df = discountFactor(d: expiryDate);
126 Real div = dividendYield(d: expiryDate);
127
128 // Input to fourier transform
129 std::vector<std::complex<Real> > fti;
130 fti.resize(new_size: n);
131
132 // Precalculate any discount factors etc.
133 precalculateExpiry(d: expiryDate);
134
135 for (Size i=0; i<n; i++)
136 {
137 Real v_j = eta * i;
138 Real sw = eta * (3.0 + ((i % 2) == 0 ? -1.0 : 1.0) - ((i == 0) ? 1.0 : 0.0)) / 3.0;
139
140 std::complex<Real> psi = df * complexFourierTransform(u: v_j - (alpha + 1)* i1);
141 psi = psi / (alpha*alpha + alpha - v_j*v_j + i1 * (2 * alpha + 1.0) * v_j);
142
143 fti[i] = std::exp(z: i1 * b * v_j) * sw * psi;
144 }
145
146 // Perform fft
147 std::vector<std::complex<Real> > results(n);
148 FastFourierTransform fft(log2_n);
149 fft.transform(inBegin: fti.begin(), inEnd: fti.end(), out: results.begin());
150
151 // Call prices
152 std::vector<Real> prices, strikes;
153 prices.resize(new_size: n);
154 strikes.resize(new_size: n);
155 for (Size i=0; i<n; i++)
156 {
157 Real k_u = -b + lambda_ * i;
158 prices[i] = (std::exp(x: -alpha * k_u) / M_PI) * results[i].real();
159 strikes[i] = std::exp(x: k_u);
160 }
161
162 for (const auto& payoff : payIt->second) {
163 Real callPrice = LinearInterpolation(strikes.begin(), strikes.end(),
164 prices.begin())(payoff->strike());
165 switch (payoff->optionType())
166 {
167 case Option::Call:
168 resultMap_[expiryDate][payoff] = callPrice;
169 break;
170 case Option::Put:
171 resultMap_[expiryDate][payoff] = callPrice - process_->x0() * div + payoff->strike() * df;
172 break;
173 default:
174 QL_FAIL("Invalid option type");
175 }
176 }
177 }
178 }
179
180}
181
182

source code of quantlib/ql/experimental/variancegamma/fftengine.cpp