24#include <ql/errors.hpp>
25#include <ql/experimental/varianceoption/integralhestonvarianceoptionengine.hpp>
26#include <ql/functional.hpp>
55 typedef std::complex<Real> Complex;
61 std::unique_ptr<double[]> xiv(
new double[2048*2048+1]);
66 Real option=0, impart=0;
68 std::unique_ptr<Complex[]> ff(
new Complex[2048*2048]);
70 Complex ui,beta,zita,gamma,csum,vero;
71 Complex contrib, caux, caux1,caux2,caux3;
91 pi= 3.14159265358979324;
93 Real s=2.0*chi*theta/(eps*eps)-1.0;
102 QL_FAIL(
"this parameter must be greater than zero-> " << s);
127 nris=std::sqrt(pi2)/dstep;
128 mm=(int)(pi2/(nris*nris));
135 for (j=0;j<=mm-1;j++)
137 xiv[j+1]=(j-mm/2.0)*nris;
140 for (j=0;j<=mm-1;j++)
149 zita=0.5*std::sqrt(caux2);
151 caux1=std::exp(-2.0*tau*zita);
154 beta=beta+caux1*(zita-0.5*chi);
158 caux2=caux*(zita-0.5*chi);
159 caux=ss*std::log(2.0*(zita/beta));
160 caux3=-v0*ui*xi*(gamma/beta);
164 ff[j+1]=std::exp(caux);
165 if(std::sqrt(std::imag(xi)*std::imag(xi)+std::real(xi)*std::real(xi))>1.e-06)
167 contrib=-eprice/(ui*xi);
173 contrib=contrib+caux/caux2;
177 contrib=eprice*eprice*0.5;
179 ff[j+1]=ff[j+1]*contrib;
182 for (j=0;j<=mm-1;j++)
184 caux=std::pow(-1.0,j);
185 caux2=-2.0*pi*(double)mm*(
double)j*0.5/(double)mm;
187 csum=csum+ff[j+1]*caux*std::exp(caux3);
189 csum=csum*std::sqrt(std::pow(-1.0,mm))*nris/pi2;
190 vero=i0-eprice+theta*tau+(1.0-std::exp(-chi*tau))*(v0-theta)/chi;
192 option=std::exp(-rtax*tau)*std::real(csum);
193 impart=std::imag(csum);
194 QL_ENSURE(impart <= 1e-12,
195 "imaginary part option (must be zero) = " << impart);
203 const ext::function<
Real(
Real)>& payoff) {
206 std::unique_ptr<double[]> xiv(
new double[2048*2048+1]);
207 std::unique_ptr<double[]> ivet(
new double[2048 * 2048 + 1]);
220 std::unique_ptr<Complex[]> ff(
new Complex[2048*2048]);
222 Complex ui,beta,zita,gamma,csum;
223 Complex caux,caux1,caux2,caux3;
240 pi= 3.14159265358979324;
243 Real s=2.0*chi*theta/(eps*eps)-1.0;
252 QL_FAIL(
"this parameter must be greater than zero-> " << s);
277 nris=std::sqrt(pi2)/dstep;
278 mm=(int)(pi2/(nris*nris));
286 for (j=0;j<=mm-1;j++)
288 xiv[j+1]=(j-mm/2.0)*nris;
289 ivet[j+1]=(j-mm/2.0)*pi2/((
double)mm*nris);
292 for (j=0;j<=mm-1;j++)
302 zita=0.5*std::sqrt(caux2);
303 caux1=std::exp(-2.0*tau*zita);
306 beta=beta+caux1*(zita-0.5*chi);
311 caux2=caux*(zita-0.5*chi);
312 caux=ss*std::log(2.0*(zita/beta));
313 caux3=-v0*ui*xi*(gamma/beta);
316 ff[j+1]=std::exp(caux);
321 for (k=0;k<=mm-1;k++)
324 payoffval=payoff(ip);
326 dxi=2.0*pi*(double)k/(
double)mm*ui;
328 for (j=0;j<=mm-1;j++)
331 caux=std::pow(-1.0,j);
332 csum=csum+ff[j+1]*caux*std::exp(z);
334 csum=csum*std::pow(-1.0,k)*nris/pi2;
336 sumr=sumr+payoffval*std::real(csum);
342 option=std::exp(-rtax*tau)*sumr;
349 struct payoff_adapter {
350 ext::shared_ptr<QuantLib::Payoff> payoff;
351 explicit payoff_adapter(ext::shared_ptr<QuantLib::Payoff> payoff)
352 : payoff(
std::move(payoff)) {}
361 ext::shared_ptr<HestonProcess> process)
362 : process_(
std::move(process)) {
368 QL_REQUIRE(
process_->dividendYield().empty(),
369 "this engine does not manage dividend yields");
379 Time tau = riskFreeRate->dayCounter().yearFraction(
383 riskFreeRate->dayCounter(),
386 ext::shared_ptr<PlainVanillaPayoff> plainPayoff =
388 if ((plainPayoff !=
nullptr) && plainPayoff->optionType() ==
Option::Call) {
390 Real strike = plainPayoff->strike();
395 results_.
value = IvopTwoDim(epsilon, chi, theta, rho, v0, tau, r,
VarianceOption::results results_
VarianceOption::arguments arguments_
Shared handle to an observable.
ext::shared_ptr< HestonProcess > process_
IntegralHestonVarianceOptionEngine(ext::shared_ptr< HestonProcess >)
void calculate() const override
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
static Settings & instance()
access to the unique instance
ext::shared_ptr< Payoff > payoff
Real Time
continuous quantity with 1-year units