sabrinterpolation.hpp

有很多的函数库

                    y_[0] = std::sqrt(x[0] - eps1_);
                    y_[1] = std::sqrt(-std::log(x[1]));
                    y_[2] = std::sqrt(x[2] - eps1_);
                    {
                        //arcsin expansion
                        const Real z = x[3]/eps2_;
                        const Real z3 = z*z*z;
                        const Real z5 = z3*z*z;
                        y_[3] = z + z3/6 + 3*z5/40;
                    }
                    return y_;
                }
            };


            class SABRError : public CostFunction {
              public:
                SABRError(SABRInterpolationImpl* sabr)
                : sabr_(sabr) {}

                Real value(const Array& x) const {
                    const Array y = sabr_->transformation_->direct(x);
                    sabr_->alpha_ = y[0];
                    sabr_->beta_  = y[1];
                    sabr_->nu_    = y[2];
                    sabr_->rho_   = y[3];
                    return sabr_->interpolationSquaredError();
                }

                Disposable<Array> values(const Array& x) const{
                    const Array y = sabr_->transformation_->direct(x);
                    sabr_->alpha_ = y[0];
                    sabr_->beta_  = y[1];
                    sabr_->nu_    = y[2];
                    sabr_->rho_   = y[3];
                    return sabr_->interpolationErrors(x);
                }

              private:
                SABRInterpolationImpl* sabr_;
            };

            // optimization method used for fitting
            boost::shared_ptr<EndCriteria> endCriteria_;
            boost::shared_ptr<OptimizationMethod> method_;
            std::vector<Real> weights_;
            const Real& forward_;
            bool vegaWeighted_;
            boost::shared_ptr<ParametersTransformation> transformation_;
            NoConstraint constraint_;

          public:
            SABRInterpolationImpl(
                const I1& xBegin, const I1& xEnd,
                const I2& yBegin,
                Time t,
                const Real& forward,
                Real alpha, Real beta, Real nu, Real rho,
                bool alphaIsFixed,
                bool betaIsFixed,
                bool nuIsFixed,
                bool rhoIsFixed,
                bool vegaWeighted,
                const boost::shared_ptr<EndCriteria>& endCriteria,
                const boost::shared_ptr<OptimizationMethod>& method)
            : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
              SABRCoefficientHolder(t, forward, alpha, beta, nu, rho,
                                    alphaIsFixed, betaIsFixed,
                                    nuIsFixed, rhoIsFixed),
              endCriteria_(endCriteria), method_(method),
              weights_(xEnd-xBegin, 1.0/(xEnd-xBegin)), forward_(forward),
              vegaWeighted_(vegaWeighted)
            {
                // if no method is provided we provide one
                if (!method_)
                    method_ = boost::shared_ptr<OptimizationMethod>(new
                        Simplex(0.01));
                if (!endCriteria_) {
                    endCriteria_ = boost::shared_ptr<EndCriteria>(new
                        EndCriteria(60000, 100, 1e-8, 1e-8, 1e-8));
                }
            }

            void update() {
                // forward_ might have changed
                QL_REQUIRE(forward_>0.0, "forward must be positive: " <<
                           io::rate(forward_) << " not allowed");

                // we should also check that y contains positive values only

                // we must update weights if it is vegaWeighted
                if (vegaWeighted_) {
                    std::vector<Real>::const_iterator x = this->xBegin_;
                    std::vector<Real>::const_iterator y = this->yBegin_;
                    std::vector<Real>::iterator w = weights_.begin();
                    Real weightsSum = 0.0;
                    for ( ; x!=this->xEnd_; ++x, ++y, ++w) {
                        Real stdDev = std::sqrt((*y)*(*y)*t_);
                        *w = blackFormulaStdDevDerivative(*x, forward_, stdDev);
                        weightsSum += *w;
                    }
                    // weight normalization
                    w = weights_.begin();
                    for ( ; w!=weights_.end(); ++w)
                        *w /= weightsSum;
                }

                // there is nothing to optimize
                if (alphaIsFixed_ && betaIsFixed_ && nuIsFixed_ && rhoIsFixed_) {
                    error_ = interpolationError();
                    maxError_ = interpolationMaxError();
                    SABREndCriteria_ = EndCriteria::None;
                    return;
                
                } else {

                    SABRError costFunction(this);
                    transformation_ = boost::shared_ptr<ParametersTransformation>(new
                        SabrParametersTransformation);

                    Array guess(4);
                    guess[0] = alpha_;
                    guess[1] = beta_;
                    guess[2] = nu_;
                    guess[3] = rho_;

                    std::vector<bool> parameterAreFixed(4);
                    parameterAreFixed[0] = alphaIsFixed_;
                    parameterAreFixed[1] = betaIsFixed_;
                    parameterAreFixed[2] = nuIsFixed_;
                    parameterAreFixed[3] = rhoIsFixed_;

                    Array inversedTransformatedGuess(transformation_->inverse(guess));

                    ProjectedCostFunction constrainedSABRError(costFunction,
                                    inversedTransformatedGuess, parameterAreFixed);

                    Array projectedGuess
                        (constrainedSABRError.project(inversedTransformatedGuess));

                    NoConstraint constraint;
                    Problem problem(constrainedSABRError, constraint, projectedGuess);
                    SABREndCriteria_ = method_->minimize(problem, *endCriteria_);
                    Array projectedResult(problem.currentValue());
                    Array transfResult(constrainedSABRError.include(projectedResult));

                    Array result = transformation_->direct(transfResult);
                    alpha_ = result[0];
                    beta_  = result[1];
                    nu_    = result[2];
                    rho_   = result[3];

                }
                error_ = interpolationError();
                maxError_ = interpolationMaxError();

            }

            Real value(Real x) const {
                QL_REQUIRE(x>0.0, "strike must be positive: " <<
                                  io::rate(x) << " not allowed");
                return sabrVolatility(x, forward_, t_,
                                      alpha_, beta_, nu_, rho_);
            }
            Real primitive(Real) const {
                QL_FAIL("SABR primitive not implemented");
            }
            Real derivative(Real) const {
                QL_FAIL("SABR derivative not implemented");
            }
            Real secondDerivative(Real) const {
                QL_FAIL("SABR secondDerivative not implemented");
            }
            // calculate total squared weighted difference (L2 norm)
            Real interpolationSquaredError() const {
                Real error, totalError = 0.0;
                std::vector<Real>::const_iterator x = this->xBegin_;
                std::vector<Real>::const_iterator y = this->yBegin_;
                std::vector<Real>::const_iterator w = weights_.begin();
                for (; x != this->xEnd_; ++x, ++y, ++w) {
                    error = (value(*x) - *y);
                    totalError += error*error * (*w);
                }
                return totalError;
            }
            // calculate weighted differences
            Disposable<Array> interpolationErrors(const Array&) const {
                Array results(this->xEnd_ - this->xBegin_);
                std::vector<Real>::const_iterator x = this->xBegin_;
                Array::iterator r = results.begin();
                std::vector<Real>::const_iterator y = this->yBegin_;
                std::vector<Real>::const_iterator w = weights_.begin();
                for (; x != this->xEnd_; ++x, ++r, ++w, ++y) {
                    *r = (value(*x) - *y)* std::sqrt(*w);
                }
                return results;
            }

            Real interpolationError() const {
                Size n = this->xEnd_-this->xBegin_;
                Real squaredError = interpolationSquaredError();
                return std::sqrt(n*squaredError/(n-1));
            }

            Real interpolationMaxError() const {
                Real error, maxError = QL_MIN_REAL;
                I1 i = this->xBegin_;
                I2 j = this->yBegin_;
                for (; i != this->xEnd_; ++i, ++j) {
                    error = std::fabs(value(*i) - *j);
                    maxError = std::max(maxError, error);
                }
                return maxError;
            }
        };

    }

}

#endif