Merge 07546f4 into d867391

ql/experimental/finitedifferences/fdmhestonfwdop.cpp

@@ -41,13 +41,15 @@ namespace QuantLib {
             const ext::shared_ptr<FdmMesher>& mesher,
             const ext::shared_ptr<HestonProcess>& process,
             FdmSquareRootFwdOp::TransformationType type,
-            const ext::shared_ptr<LocalVolTermStructure>& leverageFct)
+            const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
+            const Real mixingFactor)
     : type_(type),
       kappa_(process->kappa()),
       theta_(process->theta()),
       sigma_(process->sigma()),
       rho_  (process->rho()),
       v0_   (process->v0()),
+      mixedSigma_ (mixingFactor*sigma_),
       rTS_  (process->riskFreeRate().currentLink()),
       qTS_  (process->dividendYield().currentLink()),
       varianceValues_(0.5*mesher->locations(1)),
@@ -59,13 +61,13 @@ namespace QuantLib {
           ))),
       boundary_(ext::make_shared<ModTripleBandLinearOp>(TripleBandLinearOp(SecondDerivativeOp(0, mesher).mult(Array(mesher->locations(0).size(), 0.0))))),
       mapX_  (ext::make_shared<TripleBandLinearOp>(0, mesher)),
-      mapY_  (ext::make_shared<FdmSquareRootFwdOp>(mesher,kappa_,theta_,sigma_, 1, type)),
+      mapY_  (ext::make_shared<FdmSquareRootFwdOp>(mesher,kappa_,theta_,mixedSigma_, 1, type)),
       correlation_(ext::make_shared<NinePointLinearOp>(
           type == FdmSquareRootFwdOp::Log ?
               SecondOrderMixedDerivativeOp(0, 1, mesher)
-              .mult(Array(mesher->layout()->size(), rho_*sigma_))
+              .mult(Array(mesher->layout()->size(), rho_*mixedSigma_))
             : SecondOrderMixedDerivativeOp(0, 1, mesher)
-              .mult(rho_*sigma_*mesher->locations(1))
+              .mult(rho_*mixedSigma_*mesher->locations(1))
            )),
 	   leverageFct_(leverageFct),
 	   mesher_(mesher)
@@ -79,8 +81,8 @@ namespace QuantLib {
         const Real logFacLow = type == FdmSquareRootFwdOp::Log ? exp(mapY_->v(0)) : 1.0;
         const Real logFacUpp = type == FdmSquareRootFwdOp::Log ? exp(mapY_->v(n+1)) : 1.0;
 
-        const Real alpha = -2*rho_/sigma_*lowerBoundaryFactor*logFacLow; 
-        const Real beta  = -2*rho_/sigma_*upperBoundaryFactor*logFacUpp; 
+        const Real alpha = -2*rho_/mixedSigma_*lowerBoundaryFactor*logFacLow;
+        const Real beta  = -2*rho_/mixedSigma_*upperBoundaryFactor*logFacUpp;
 
         ModTripleBandLinearOp fDx(FirstDerivativeOp(0, mesher));
 
@@ -122,30 +124,30 @@ namespace QuantLib {
         if (leverageFct_) {
             L_ = getLeverageFctSlice(t1, t2);
             Array Lsquare = L_*L_;
-            if (type_ == FdmSquareRootFwdOp::Plain) { 
+            if (type_ == FdmSquareRootFwdOp::Plain) {
                 mapX_->axpyb( Array(1, -r + q), *dxMap_,
                     dxxMap_->multR(Lsquare).add(boundary_->multR(L_))
-                    .add(dxMap_->multR(rho_*sigma_*L_))
-                    .add(dxMap_->mult(varianceValues_).multR(Lsquare)), 
+                    .add(dxMap_->multR(rho_*mixedSigma_*L_))
+                    .add(dxMap_->mult(varianceValues_).multR(Lsquare)),
                               Array());
             } else if (type_ == FdmSquareRootFwdOp::Power) {
                 mapX_->axpyb( Array(1, -r + q), *dxMap_,
                     dxxMap_->multR(Lsquare).add(boundary_->multR(L_))
-                    .add(dxMap_->multR(rho_*2.0*kappa_*theta_/(sigma_)*L_))
+                    .add(dxMap_->multR(rho_*2.0*kappa_*theta_/(mixedSigma_)*L_))
                     .add(dxMap_->mult(varianceValues_).multR(Lsquare)), Array());
             } else if (type_ == FdmSquareRootFwdOp::Log) {
                 mapX_->axpyb( Array(1, -r + q), *dxMap_,
                     dxxMap_->multR(Lsquare).add(boundary_->multR(L_))
-                    .add(dxMap_->mult(0.5*Exp(2.0*varianceValues_)).multR(Lsquare)), 
+                    .add(dxMap_->mult(0.5*Exp(2.0*varianceValues_)).multR(Lsquare)),
                               Array());
             }
         }
         else {
             if (type_ == FdmSquareRootFwdOp::Plain) {
-                mapX_->axpyb( - r + q + rho_*sigma_ + varianceValues_, *dxMap_,
+                mapX_->axpyb( - r + q + rho_*mixedSigma_ + varianceValues_, *dxMap_,
                         *dxxMap_, Array());
             } else if (type_ == FdmSquareRootFwdOp::Power) {
-                mapX_->axpyb( - r + q + rho_*2.0*kappa_*theta_/(sigma_) + varianceValues_, 
+                mapX_->axpyb( - r + q + rho_*2.0*kappa_*theta_/(mixedSigma_) + varianceValues_,
                               *dxMap_, *dxxMap_, Array());
             } else if (type_ == FdmSquareRootFwdOp::Log) {
                 mapX_->axpyb( - r + q + 0.5*Exp(2.0*varianceValues_), *dxMap_,

ql/experimental/finitedifferences/fdmhestonfwdop.hpp

@@ -47,7 +47,8 @@ namespace QuantLib {
             FdmSquareRootFwdOp::TransformationType type 
                 = FdmSquareRootFwdOp::Plain,
             const ext::shared_ptr<LocalVolTermStructure> & leverageFct
-                = ext::shared_ptr<LocalVolTermStructure>());
+                = ext::shared_ptr<LocalVolTermStructure>(),
+            const Real mixingFactor = 1.0);
 
         Size size() const;
         void setTime(Time t1, Time t2);
@@ -67,7 +68,7 @@ namespace QuantLib {
       private:
         Disposable<Array> getLeverageFctSlice(Time t1, Time t2) const;
         const FdmSquareRootFwdOp::TransformationType type_;
-        const Real kappa_, theta_, sigma_, rho_, v0_;
+        const Real kappa_, theta_, sigma_, rho_, v0_, mixedSigma_;
 
         const ext::shared_ptr<YieldTermStructure> rTS_;
         const ext::shared_ptr<YieldTermStructure> qTS_;

ql/experimental/models/hestonslvfdmmodel.cpp

@@ -273,12 +273,14 @@ namespace QuantLib {
         const Date& endDate,
         const HestonSLVFokkerPlanckFdmParams& params,
         const bool logging,
-        const std::vector<Date>& mandatoryDates)
+        const std::vector<Date>& mandatoryDates,
+        const Real mixingFactor)
     : localVol_(localVol),
       hestonModel_(hestonModel),
       endDate_(endDate),
       params_(params),
       mandatoryDates_(mandatoryDates),
+      mixingFactor_(mixingFactor),
       logging_(logging) {
 
         registerWith(localVol_);
@@ -316,7 +318,8 @@ namespace QuantLib {
         const Real kappa = hestonProcess->kappa();
         const Real theta = hestonProcess->theta();
         const Real sigma = hestonProcess->sigma();
-        const Real alpha = 2*kappa*theta/(sigma*sigma);
+        const Real mixedSigma = mixingFactor_ * sigma;
+        const Real alpha = 2*kappa*theta/(mixedSigma*mixedSigma);
 
         const Size xGrid = params_.xGrid;
         const Size vGrid = params_.vGrid;
@@ -366,7 +369,7 @@ namespace QuantLib {
             = localVolRND.rescaleTimeSteps();
 
         const SquareRootProcessRNDCalculator squareRootRnd(
-            v0, kappa, theta, sigma);
+            v0, kappa, theta, mixedSigma);
 
         const FdmSquareRootFwdOp::TransformationType trafoType
           = params_.trafoType;
@@ -434,7 +437,7 @@ namespace QuantLib {
             new FixedLocalVolSurface(referenceDate, times, vStrikes, L, dc));
 
         ext::shared_ptr<FdmLinearOpComposite> hestonFwdOp(
-            new FdmHestonFwdOp(mesher, hestonProcess, trafoType, leverageFct));
+            new FdmHestonFwdOp(mesher, hestonProcess, trafoType, leverageFct, mixingFactor_));
 
         Array p = FdmHestonGreensFct(mesher, hestonProcess, trafoType, lv0)
             .get(timeGrid->at(1), params_.greensAlgorithm);

ql/experimental/models/hestonslvfdmmodel.hpp

@@ -78,7 +78,8 @@ class SimpleQuote;
             const Date& endDate,
             const HestonSLVFokkerPlanckFdmParams& params,
             const bool logging = false,
-            const std::vector<Date>& mandatoryDates = std::vector<Date>());
+            const std::vector<Date>& mandatoryDates = std::vector<Date>(),
+            const Real mixingFactor = 1.0);
 
         ext::shared_ptr<HestonProcess> hestonProcess() const;
         ext::shared_ptr<LocalVolTermStructure> localVol() const;
@@ -100,6 +101,7 @@ class SimpleQuote;
         const Date endDate_;
         const HestonSLVFokkerPlanckFdmParams params_;
         const std::vector<Date> mandatoryDates_;
+        const Real mixingFactor_;
 
         mutable ext::shared_ptr<LocalVolTermStructure> leverageFunction_;
 

ql/experimental/models/hestonslvmcmodel.cpp

@@ -44,13 +44,15 @@ namespace QuantLib {
         Size timeStepsPerYear,
         Size nBins,
         Size calibrationPaths,
-        const std::vector<Date>& mandatoryDates)
+        const std::vector<Date>& mandatoryDates,
+        const Real mixingFactor)
     : localVol_(localVol),
       hestonModel_(hestonModel),
       brownianGeneratorFactory_(brownianGeneratorFactory),
       endDate_(endDate),
       nBins_(nBins),
-      calibrationPaths_(calibrationPaths) {
+      calibrationPaths_(calibrationPaths),
+      mixingFactor_(mixingFactor) {
 
         registerWith(localVol_);
         registerWith(hestonModel_);
@@ -120,7 +122,7 @@ namespace QuantLib {
             vStrikes, L, dc);
 
         const ext::shared_ptr<HestonSLVProcess> slvProcess
-            = ext::make_shared<HestonSLVProcess>(hestonProcess, leverageFunction_);
+            = ext::make_shared<HestonSLVProcess>(hestonProcess, leverageFunction_, mixingFactor_);
 
         std::vector<std::pair<Real, Real> > pairs(
                 calibrationPaths_, std::make_pair(spot->value(), v0));

ql/experimental/models/hestonslvmcmodel.hpp

@@ -50,7 +50,8 @@ namespace QuantLib {
             Size timeStepsPerYear = 365,
             Size nBins = 201,
             Size calibrationPaths = (1 << 15),
-            const std::vector<Date>& mandatoryDates = std::vector<Date>());
+            const std::vector<Date>& mandatoryDates = std::vector<Date>(),
+            const Real mixingFactor = 1.0);
 
         ext::shared_ptr<HestonProcess> hestonProcess() const;
         ext::shared_ptr<LocalVolTermStructure> localVol() const;
@@ -65,6 +66,7 @@ namespace QuantLib {
         const ext::shared_ptr<BrownianGeneratorFactory> brownianGeneratorFactory_;
         const Date endDate_;
         const Size nBins_, calibrationPaths_;
+        const Real mixingFactor_;
         ext::shared_ptr<TimeGrid> timeGrid_;
 
         mutable ext::shared_ptr<FixedLocalVolSurface> leverageFunction_;

ql/experimental/processes/hestonslvprocess.cpp

@@ -31,9 +31,11 @@ namespace QuantLib {
 
     HestonSLVProcess::HestonSLVProcess(
         const ext::shared_ptr<HestonProcess>& hestonProcess,
-        const ext::shared_ptr<LocalVolTermStructure>& leverageFct)
+        const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
+        const Real mixingFactor)
     : hestonProcess_(hestonProcess),
-      leverageFct_(leverageFct) {
+      leverageFct_(leverageFct),
+      mixingFactor_(mixingFactor) {
         registerWith(hestonProcess);
         update();
     };
@@ -44,6 +46,7 @@ namespace QuantLib {
         theta_ = hestonProcess_->theta();
         sigma_ = hestonProcess_->sigma();
         rho_   = hestonProcess_->rho();
+        mixedSigma_ = mixingFactor_ * sigma_;
     }
 
     Disposable<Array> HestonSLVProcess::drift(Time t, const Array& x) const {
@@ -69,7 +72,7 @@ namespace QuantLib {
         const Real vol =
             std::min(1e-8, std::sqrt(x[1]*leverageFct_->localVol(t, s, true)));
 
-        const Real sigma2 = sigma_ * vol;
+        const Real sigma2 = mixedSigma_ * vol;
         const Real sqrhov = std::sqrt(1.0 - rho_*rho_);
 
         Matrix tmp(2,2);
@@ -86,8 +89,8 @@ namespace QuantLib {
         const Real ex = std::exp(-kappa_*dt);
 
         const Real m  =  theta_+(x0[1]-theta_)*ex;
-        const Real s2 =  x0[1]*sigma_*sigma_*ex/kappa_*(1-ex)
-                       + theta_*sigma_*sigma_/(2*kappa_)*(1-ex)*(1-ex);
+        const Real s2 =  x0[1]*mixedSigma_*mixedSigma_*ex/kappa_*(1-ex)
+                       + theta_*mixedSigma_*mixedSigma_/(2*kappa_)*(1-ex)*(1-ex);
         const Real psi = s2/(m*m);
 
         if (psi < 1.5) {
@@ -114,7 +117,7 @@ namespace QuantLib {
         const Real v_0 = 0.5*(x0[1]+retVal[1])*l_0*l_0;
 
         retVal[0] = x0[0]*std::exp(mu*dt - 0.5*v_0*dt
-            + rho_/sigma_*l_0 * (
+            + rho_/mixedSigma_*l_0 * (
                   retVal[1] - kappa_*theta_*dt
                   + 0.5*(x0[1]+retVal[1])*kappa_*dt - x0[1])
             + rho1*std::sqrt(v_0*dt)*dw[0]);

ql/experimental/processes/hestonslvprocess.hpp

@@ -34,7 +34,8 @@ namespace QuantLib {
       public:
         HestonSLVProcess(
             const ext::shared_ptr<HestonProcess>& hestonProcess,
-            const ext::shared_ptr<LocalVolTermStructure>& leverageFct);
+            const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
+            const Real mixingFactor = 1.0);
 
         Size size()    const { return Size(2); }
         Size factors() const { return Size(2); }
@@ -58,6 +59,7 @@ namespace QuantLib {
         Real kappa() const { return kappa_; }
         Real theta() const { return theta_; }
         Real sigma() const { return sigma_; }
+        Real mixingFactor() const { return mixingFactor_; }
         ext::shared_ptr<LocalVolTermStructure> leverageFct() const {
             return leverageFct_;
         }
@@ -73,7 +75,7 @@ namespace QuantLib {
         Time time(const Date& d) const { return hestonProcess_->time(d); }
 
       private:
-        Real kappa_, theta_, sigma_, rho_, v0_;
+        Real kappa_, theta_, sigma_, rho_, v0_, mixingFactor_, mixedSigma_;
 
         const ext::shared_ptr<HestonProcess> hestonProcess_;
         const ext::shared_ptr<LocalVolTermStructure> leverageFct_;

ql/methods/finitedifferences/meshers/fdmhestonvariancemesher.cpp

@@ -58,21 +58,23 @@ namespace QuantLib {
     FdmHestonVarianceMesher::FdmHestonVarianceMesher(
         Size size,
         const ext::shared_ptr<HestonProcess> & process,
-        Time maturity, Size tAvgSteps, Real epsilon)
+        Time maturity, Size tAvgSteps, Real epsilon,
+        Real mixingFactor)
         : Fdm1dMesher(size) {
 
         std::vector<Real> vGrid(size, 0.0), pGrid(size, 0.0);
+        const Real mixedSigma = process->sigma()*mixingFactor;
         const Real df  = 4*process->theta()*process->kappa()/
-                            square<Real>()(process->sigma());
+                            square<Real>()(mixedSigma);
         try {
             std::multiset<std::pair<Real, Real> > grid;
 
             for (Size l=1; l<=tAvgSteps; ++l) {
                 const Real t = (maturity*l)/tAvgSteps;
                 const Real ncp = 4*process->kappa()*std::exp(-process->kappa()*t)
-                    /(square<Real>()(process->sigma())
+                    /(square<Real>()(mixedSigma)
                     *(1-std::exp(-process->kappa()*t)))*process->v0();
-                const Real k = square<Real>()(process->sigma())
+                const Real k = square<Real>()(mixedSigma)
                     *(1-std::exp(-process->kappa()*t))/(4*process->kappa());
 
                 const Real qMin = 0.0; // v_min = 0.0;
@@ -114,7 +116,7 @@ namespace QuantLib {
         } 
         catch (const Error&) {
             // use default mesh
-            const Real vol = process->sigma()*
+            const Real vol = mixedSigma*
                 std::sqrt(process->theta()/(2*process->kappa()));
 
             const Real mean = process->theta();
@@ -129,7 +131,7 @@ namespace QuantLib {
         }
 
         Real skewHint = ((process->kappa() != 0.0) 
-                ? std::max(1.0, process->sigma()/process->kappa()) : 1.0);
+                ? std::max(1.0, mixedSigma/process->kappa()) : 1.0);
 
         std::sort(pGrid.begin(), pGrid.end());
         volaEstimate_ = GaussLobattoIntegral(100000, 1e-4)(
@@ -159,11 +161,12 @@ namespace QuantLib {
         Size size,
         const ext::shared_ptr<HestonProcess>& process,
         const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
-        Time maturity, Size tAvgSteps, Real epsilon)
+        Time maturity, Size tAvgSteps, Real epsilon,
+        Real mixingFactor)
      : Fdm1dMesher(size) {
 
         const FdmHestonVarianceMesher mesher(
-            size, process, maturity, tAvgSteps, epsilon);
+            size, process, maturity, tAvgSteps, epsilon, mixingFactor);
 
         for (Size i=0; i < size; ++i) {
             dplus_[i] = mesher.dplus(i);

ql/methods/finitedifferences/meshers/fdmhestonvariancemesher.hpp

@@ -36,7 +36,8 @@ namespace QuantLib {
         FdmHestonVarianceMesher(
             Size size,
             const ext::shared_ptr<HestonProcess> & process,
-            Time maturity, Size tAvgSteps = 10, Real epsilon = 0.0001);
+            Time maturity, Size tAvgSteps = 10, Real epsilon = 0.0001,
+            Real mixingFactor = 1.0);
 
         Real volaEstimate() const { return volaEstimate_; }
 
@@ -53,7 +54,8 @@ namespace QuantLib {
             Size size,
             const ext::shared_ptr<HestonProcess>& process,
             const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
-            Time maturity, Size tAvgSteps = 10, Real epsilon = 0.0001);
+            Time maturity, Size tAvgSteps = 10, Real epsilon = 0.0001,
+            Real mixingFactor = 1.0);
 
         Real volaEstimate() const { return volaEstimate_; }
 

ql/methods/finitedifferences/operators/fdmhestonop.cpp

@@ -115,9 +115,9 @@ namespace QuantLib {
     FdmHestonVariancePart::FdmHestonVariancePart(
         const ext::shared_ptr<FdmMesher>& mesher,
         const ext::shared_ptr<YieldTermStructure>& rTS,
-        Real sigma, Real kappa, Real theta)
+        Real mixedSigma, Real kappa, Real theta)
     : dyMap_(SecondDerivativeOp(1, mesher)
-                .mult(0.5*sigma*sigma*mesher->locations(1))
+                .mult(0.5*mixedSigma*mixedSigma*mesher->locations(1))
              .add(FirstDerivativeOp(1, mesher)
                   .mult(kappa*(theta - mesher->locations(1))))),
       mapT_(1, mesher),
@@ -137,12 +137,14 @@ namespace QuantLib {
         const ext::shared_ptr<FdmMesher>& mesher,
         const ext::shared_ptr<HestonProcess> & hestonProcess,
         const ext::shared_ptr<FdmQuantoHelper>& quantoHelper,
-        const ext::shared_ptr<LocalVolTermStructure>& leverageFct)
+        const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
+        const Real mixingFactor)
     : correlationMap_(SecondOrderMixedDerivativeOp(0, 1, mesher)
                         .mult(hestonProcess->rho()*hestonProcess->sigma()
+                                *mixingFactor
                                 *mesher->locations(1))),
       dyMap_(mesher, hestonProcess->riskFreeRate().currentLink(),
-              hestonProcess->sigma(), 
+              hestonProcess->sigma()*mixingFactor,
               hestonProcess->kappa(), 
               hestonProcess->theta()),
       dxMap_(mesher,