Merge a9e32f0 into d867391

ql/experimental/finitedifferences/fdmhestonfwdop.cpp

@@ -48,6 +48,8 @@ namespace QuantLib {
       sigma_(process->sigma()),
       rho_  (process->rho()),
       v0_   (process->v0()),
+      mixingFactor_ (process->mixingFactor()),
+      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));
 
@@ -125,13 +127,13 @@ namespace QuantLib {
             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_->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_,
@@ -142,10 +144,10 @@ namespace QuantLib {
         }
         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

@@ -67,7 +67,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_, mixingFactor_, mixedSigma_;
 
         const ext::shared_ptr<YieldTermStructure> rTS_;
         const ext::shared_ptr<YieldTermStructure> qTS_;

ql/experimental/models/hestonslvfdmmodel.cpp

@@ -316,7 +316,9 @@ 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 mixingFactor = hestonProcess->mixingFactor();
+        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 +368,7 @@ namespace QuantLib {
             = localVolRND.rescaleTimeSteps();
 
         const SquareRootProcessRNDCalculator squareRootRnd(
-            v0, kappa, theta, sigma);
+            v0, kappa, theta, mixedSigma);
 
         const FdmSquareRootFwdOp::TransformationType trafoType
           = params_.trafoType;

ql/experimental/processes/hestonslvprocess.cpp

@@ -44,6 +44,8 @@ namespace QuantLib {
         theta_ = hestonProcess_->theta();
         sigma_ = hestonProcess_->sigma();
         rho_   = hestonProcess_->rho();
+        mixingFactor_ = hestonProcess_->mixingFactor();
+        mixedSigma_ = mixingFactor_ * sigma_;
     }
 
     Disposable<Array> HestonSLVProcess::drift(Time t, const Array& x) const {
@@ -69,7 +71,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 +88,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 +116,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

@@ -58,6 +58,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 +74,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

@@ -62,17 +62,18 @@ namespace QuantLib {
         : Fdm1dMesher(size) {
 
         std::vector<Real> vGrid(size, 0.0), pGrid(size, 0.0);
+        const Real mixedSigma = process->sigma()*process->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 +115,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 +130,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)(

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),
@@ -139,10 +139,10 @@ namespace QuantLib {
         const ext::shared_ptr<FdmQuantoHelper>& quantoHelper,
         const ext::shared_ptr<LocalVolTermStructure>& leverageFct)
     : correlationMap_(SecondOrderMixedDerivativeOp(0, 1, mesher)
-                        .mult(hestonProcess->rho()*hestonProcess->sigma()
+                        .mult(hestonProcess->rho()*hestonProcess->mixedSigma()
                                 *mesher->locations(1))),
       dyMap_(mesher, hestonProcess->riskFreeRate().currentLink(),
-              hestonProcess->sigma(), 
+              hestonProcess->mixedSigma(),
               hestonProcess->kappa(), 
               hestonProcess->theta()),
       dxMap_(mesher,

ql/methods/finitedifferences/operators/fdmhestonop.hpp

@@ -70,7 +70,7 @@ namespace QuantLib {
         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);
 
         void setTime(Time t1, Time t2);
         const TripleBandLinearOp& getMap() const;

ql/models/equity/hestonmodel.cpp

@@ -23,7 +23,7 @@
 namespace QuantLib {
 
     HestonModel::HestonModel(const ext::shared_ptr<HestonProcess> & process)
-    : CalibratedModel(5), process_(process) {
+    : CalibratedModel(6), process_(process) {
         arguments_[0] = ConstantParameter(process->theta(),
                                           PositiveConstraint());
         arguments_[1] = ConstantParameter(process->kappa(),
@@ -34,6 +34,8 @@ namespace QuantLib {
                                           BoundaryConstraint(-1.0, 1.0));
         arguments_[4] = ConstantParameter(process->v0(),
                                           PositiveConstraint());
+        arguments_[5] = ConstantParameter(process->mixingFactor(),
+                                          PositiveConstraint());
         HestonModel::generateArguments();
 
         registerWith(process_->riskFreeRate());
@@ -46,7 +48,8 @@ namespace QuantLib {
                                          process_->dividendYield(),
                                          process_->s0(),
                                          v0(), kappa(), theta(),
-                                         sigma(), rho()));
+                                         sigma(), rho(), HestonProcess::QuadraticExponentialMartingale,
+                                         mixingFactor()));
     }
 
 }

ql/models/equity/hestonmodel.hpp

@@ -54,6 +54,8 @@ namespace QuantLib {
         // spot variance
         Real v0()    const { return arguments_[4](0.0); }
 
+        Real mixingFactor() const { return arguments_[5](0.0); }
+
         // underlying process
         ext::shared_ptr<HestonProcess> process() const { return process_; }
 
@@ -71,8 +73,9 @@ namespace QuantLib {
                 const Real theta = params[0];
                 const Real kappa = params[1];
                 const Real sigma = params[2];
+                const Real mixingFactor = params[5];
 
-                return (sigma >= 0.0 && sigma*sigma < 2.0*kappa*theta);
+                return (mixingFactor*sigma >= 0.0 && mixingFactor*mixingFactor*sigma*sigma < 2.0*kappa*theta);
             }
         };
       public:

ql/pricingengines/vanilla/fdhestonhullwhitevanillaengine.cpp

@@ -192,7 +192,7 @@ namespace QuantLib {
             Real analyticNPV = option.NPV();
 
             ext::shared_ptr<FdHestonVanillaEngine> fdEngine(
-                    new FdHestonVanillaEngine(*model_, tGrid_, xGrid_, 
+                    new FdHestonVanillaEngine(*model_, tGrid_, xGrid_,
                                               vGrid_, dampingSteps_, 
                                               schemeDesc_));
             fdEngine->enableMultipleStrikesCaching(strikes_);

ql/pricingengines/vanilla/fdhestonvanillaengine.hpp

@@ -82,7 +82,7 @@ namespace QuantLib {
         const FdmSchemeDesc schemeDesc_;
         const ext::shared_ptr<LocalVolTermStructure> leverageFct_;
         const ext::shared_ptr<FdmQuantoHelper> quantoHelper_;
-        
+
         std::vector<Real> strikes_;
         mutable std::vector<std::pair<DividendVanillaOption::arguments,
                                       DividendVanillaOption::results> >