capletcoterminalswaptioncalibration.cpp

有很多的函数库

            // now compute contribution from lower right corner
            Real extraConstantPart =0.0;
            // part of caplet approximation formula coming from later rates
            for (Size k = i+1; k < numberOfSteps; ++k)
                for (Size l = i+1; l < numberOfSteps; ++l)
                    extraConstantPart+=invertedZedMatrix[i-1][k] *
                                         CovarianceSwapCovs[i-1][k][l] *
                                         invertedZedMatrix[i-1][l];

            // now compute contribution from top row excluding first two columns and its transpose
            // we divide into two parts as one part is multiplied by a[i-1] and the other by b[i-1]
            // a lot could be precomputed when we need to optimize
            for (Size k = i+1; k < numberOfSteps; ++k)
            {
                if (i > 1)
                {
                    extraConstantPart+=invertedZedMatrix[i-1][i-1] *
                        CovarianceSwapCovs[i-2][i-1][k] *
                        invertedZedMatrix[i-1][k]*a[i-1];

                    extraConstantPart+=invertedZedMatrix[i-1][k] *
                        CovarianceSwapCovs[i-2][k][i-1] *
                        invertedZedMatrix[i-1][i-1]*a[i-1];
                }

                extraConstantPart+=invertedZedMatrix[i-1][i-1] *
                    CovarianceSwapMarginalCovs[i-1][i-1][k] *
                    invertedZedMatrix[i-1][k]*b[i-1];

                extraConstantPart+=invertedZedMatrix[i-1][k] *
                    CovarianceSwapCovs[i-1][k][i-1] *
                    invertedZedMatrix[i-1][i-1]*b[i-1];

            }

            // we also get an extra linear part, this corresponds to row i, and columns j>i+1, and transpose
            double extraLinearPart=0.0;
            for (Size k = i+1; k < numberOfSteps; ++k)
            {
                    extraLinearPart+=invertedZedMatrix[i-1][k] *
                        CovarianceSwapCovs[i-1][k][i] *
                        invertedZedMatrix[i-1][i];

                    extraLinearPart+=invertedZedMatrix[i-1][i] *
                        CovarianceSwapCovs[i-1][i][k] *
                        invertedZedMatrix[i-1][k];
            }



            Real constantPart = w0*w0*totVariance[i-1] +
                                                    extraConstantPart*extraMultiplier-v1t1;
            Real linearPart = -2*w0*w1*(a[i-1]*almostTotCovariance[i] +
                                                b[i-1]*leftCovariance[i]) +extraLinearPart*extraMultiplier ;
            Real quadraticPart = w1*w1*almostTotVariance[i];
            Real disc = linearPart*linearPart-4.0*constantPart*quadraticPart;

            Real root, minimum = -linearPart/(2.0*quadraticPart);
            bool rightUsed = false;
            if (disc <0.0) {
                ++negativeDiscriminants;
                // pick up the minimum vol for the caplet
                root = minimum;
            } else if (lowestRoot) {
                root = (-linearPart-sqrt(disc))/(2.0*quadraticPart);
            } else {
                if (minimum>1.0)
                    root = (-linearPart-sqrt(disc))/(2.0*quadraticPart);
                else {
                    rightUsed = true;
                    root = (-linearPart+sqrt(disc))/(2.0*quadraticPart);
                }
            }

            Real varianceFound = root*root*almostTotVariance[i];
            Real varianceToFind = totVariance[i]-varianceFound;
            Real mult = varianceToFind/swapTimeInhomogeneousVariances[i][i];
            if (mult<=0.0 && rightUsed) {
                root = (-linearPart-sqrt(disc))/(2.0*quadraticPart);
                varianceFound = root*root*almostTotVariance[i];
                varianceToFind = totVariance[i]-varianceFound;
                mult = varianceToFind/swapTimeInhomogeneousVariances[i][i];
            }

            if (mult<0.0) // no solution...
                return false;

            QL_ENSURE(root>=0.0,
                      "negative root -- it should have not happened");

            a[i]=root;
            b[i]=std::sqrt(mult);
        }

        {
            Integer i = numberOfSteps;
            Integer j=0;
            for (; j <= static_cast<Integer>(i)-2; j++)
                swapTimeInhomogeneousVariances[j][i-1]*= a[i-1]*a[i-1];
            swapTimeInhomogeneousVariances[j][i-1]*= b[i-1]*b[i-1];
        }

        // compute the results
        swapCovariancePseudoRoots.resize(numberOfSteps);
        for (Size k=0; k<numberOfSteps; ++k) {
            swapCovariancePseudoRoots[k] = corrPseudo[k];
            for (Size j=0; j<numberOfRates; ++j) {
                Real coeff = std::sqrt(swapTimeInhomogeneousVariances[k][j]);
                 for (Size i=0; i<numberOfFactors; ++i)
                    swapCovariancePseudoRoots[k][j][i]*=coeff;
            }
            QL_ENSURE(swapCovariancePseudoRoots[k].rows()==numberOfRates,
                      "step " << k
                      << " abcd vol wrong number of rows: "
                      << swapCovariancePseudoRoots[k].rows()
                      << " instead of " << numberOfRates);
            QL_ENSURE(swapCovariancePseudoRoots[k].columns()==numberOfFactors,
                      "step " << k
                      << " abcd vol wrong number of columns: "
                      << swapCovariancePseudoRoots[k].columns()
                      << " instead of " << numberOfFactors);
        }

        return true;
    }

    bool CapletCoterminalSwaptionCalibration::calibrate(
                            Size numberOfFactors,
                            const std::vector<Real>& alpha,
                            bool lowestRoot,
                            bool useFullApprox,
                            Size maxIterations,
                            Real tolerance) {

        const std::vector<Time>& rateTimes = evolution_.rateTimes();
        Size numberOfRates = evolution_.numberOfRates();
        negDiscr_ = 0;
        error_ = 987654321; // a positive large number
        calibrated_ = false;
        bool success = true;
        std::vector<Volatility> targetCapletVols(mktCapletVols_);
        for (Size iterCounter=1;
             iterCounter<=maxIterations && success && negDiscr_==0 && error_>tolerance;
             ++iterCounter) {
            success = calibrationFunction(evolution_,
                                          *corr_,
                                          displacedSwapVariances_,
                                          targetCapletVols,
                                          *cs_,
                                          displacement_,
                                          numberOfFactors,
                                          alpha,
                                          lowestRoot,
                                          useFullApprox,
                                          swapCovariancePseudoRoots_,
                                          negDiscr_);
            if (!success)
                return success;
            boost::shared_ptr<MarketModel> smm(new
                PseudoRootFacade(swapCovariancePseudoRoots_,
                                 rateTimes,
                                 cs_->coterminalSwapRates(),
                                 std::vector<Spread>(numberOfRates, displacement_)));
            CotSwapToFwdAdapter flmm(smm);
            // avoid this copy
            Matrix capletTotCovariance = flmm.totalCovariance(numberOfRates-1);
            error_ = 0.0;
            for (Size i=0; i<numberOfRates-1; ++i) {
                Real capletVol = std::sqrt(capletTotCovariance[i][i]/rateTimes[i]);
                Real diff = mktCapletVols_[i]-capletVol;
                error_ += diff*diff;
                targetCapletVols[i] *= mktCapletVols_[i]/capletVol;
            }
            error_ = std::sqrt(error_/(numberOfRates-1));
        }
        calibrated_ = true;
        return success;
    }

    Size CapletCoterminalSwaptionCalibration::negativeDiscriminants() const {
        QL_REQUIRE(calibrated_, "not calibrated yet");
        return negDiscr_;
    }

    Real CapletCoterminalSwaptionCalibration::rmsError() const {
        QL_REQUIRE(calibrated_, "not calibrated yet");
        return error_;
    }

    const std::vector<Matrix>&
    CapletCoterminalSwaptionCalibration::swapPseudoRoots() const {
        QL_REQUIRE(calibrated_, "not calibrated yet");
        return swapCovariancePseudoRoots_;
    }

    const Matrix&
    CapletCoterminalSwaptionCalibration::swapPseudoRoot(Size i) const {
        QL_REQUIRE(calibrated_, "not calibrated yet");
        QL_REQUIRE(i<swapCovariancePseudoRoots_.size(),
                   "invalid index");
        return swapCovariancePseudoRoots_[i];
    }

}