MarketModels.cpp

This example shows the use of interest-rate market models.

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
#include <ql/qldefines.hpp>
#ifdef BOOST_MSVC
#  include <ql/auto_link.hpp>
#endif
#include <ql/models/marketmodels/marketmodel.hpp>
#include <ql/models/marketmodels/accountingengine.hpp>
#include <ql/models/marketmodels/pathwiseaccountingengine.hpp>
#include <ql/models/marketmodels/products/multiproductcomposite.hpp>
#include <ql/models/marketmodels/products/multistep/multistepswap.hpp>
#include <ql/models/marketmodels/products/multistep/callspecifiedmultiproduct.hpp>
#include <ql/models/marketmodels/products/multistep/exerciseadapter.hpp>
#include <ql/models/marketmodels/products/multistep/multistepnothing.hpp>
#include <ql/models/marketmodels/products/multistep/multistepinversefloater.hpp>
#include <ql/models/marketmodels/products/pathwise/pathwiseproductswap.hpp>
#include <ql/models/marketmodels/products/pathwise/pathwiseproductinversefloater.hpp>
#include <ql/models/marketmodels/products/pathwise/pathwiseproductcallspecified.hpp>
#include <ql/models/marketmodels/models/flatvol.hpp>
#include <ql/models/marketmodels/callability/swapratetrigger.hpp>
#include <ql/models/marketmodels/callability/swapbasissystem.hpp>
#include <ql/models/marketmodels/callability/swapforwardbasissystem.hpp>
#include <ql/models/marketmodels/callability/nothingexercisevalue.hpp>
#include <ql/models/marketmodels/callability/collectnodedata.hpp>
#include <ql/models/marketmodels/callability/lsstrategy.hpp>
#include <ql/models/marketmodels/callability/upperboundengine.hpp>
#include <ql/models/marketmodels/correlations/expcorrelations.hpp>
#include <ql/models/marketmodels/browniangenerators/mtbrowniangenerator.hpp>
#include <ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp>
#include <ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp>
#include <ql/models/marketmodels/evolvers/lognormalfwdrateeuler.hpp>
#include <ql/models/marketmodels/pathwisegreeks/bumpinstrumentjacobian.hpp>
#include <ql/models/marketmodels/utilities.hpp>
#include <ql/methods/montecarlo/genericlsregression.hpp>
#include <ql/legacy/libormarketmodels/lmlinexpcorrmodel.hpp>
#include <ql/legacy/libormarketmodels/lmextlinexpvolmodel.hpp>
#include <ql/time/schedule.hpp>
#include <ql/time/calendars/nullcalendar.hpp>
#include <ql/time/daycounters/simpledaycounter.hpp>
#include <ql/pricingengines/blackformula.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
#include <ql/utilities/dataformatters.hpp>
#include <ql/math/integrals/segmentintegral.hpp>
#include <ql/math/statistics/convergencestatistics.hpp>
#include <ql/termstructures/volatility/abcd.hpp>
#include <ql/termstructures/volatility/abcdcalibration.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/optimization/simplex.hpp>
#include <ql/quotes/simplequote.hpp>
#include <sstream>
#include <iostream>
#include <ctime>
 
using namespace QuantLib;
 
#if defined(QL_ENABLE_SESSIONS)
namespace QuantLib {
 
    ThreadKey sessionId() { return 0; }
 
}
#endif
 
 
std::vector<std::vector<Matrix> >
theVegaBumps(bool factorwiseBumping, const ext::shared_ptr<MarketModel>& marketModel, bool doCaps) {
    Real multiplierCutOff = 50.0;
    Real projectionTolerance = 1E-4;
    Size numberRates= marketModel->numberOfRates();
 
    std::vector<VolatilityBumpInstrumentJacobian::Cap> caps;
 
    if (doCaps)
    {
 
        Rate capStrike = marketModel->initialRates()[0];
 
        for (Size i=0; i< numberRates-1; i=i+1)
        {
            VolatilityBumpInstrumentJacobian::Cap nextCap;
            nextCap.startIndex_ = i;
            nextCap.endIndex_ = i+1;
            nextCap.strike_ = capStrike;
            caps.push_back(nextCap);
        }
 
 
    }
 
 
 
    std::vector<VolatilityBumpInstrumentJacobian::Swaption> swaptions(numberRates);
 
    for (Size i=0; i < numberRates; ++i)
    {
        swaptions[i].startIndex_ = i;
        swaptions[i].endIndex_ = numberRates;
 
    }
 
    VegaBumpCollection possibleBumps(marketModel,
        factorwiseBumping);
 
    OrthogonalizedBumpFinder  bumpFinder(possibleBumps,
        swaptions,
        caps,
        multiplierCutOff, // if vector length grows by more than this discard
        projectionTolerance);      // if vector projection before scaling less than this discard
 
    std::vector<std::vector<Matrix> > theBumps;
 
    bumpFinder.GetVegaBumps(theBumps);
 
    return theBumps;
 
}
 
 
 
int Bermudan()
{
 
    Size numberRates =20;
    Real accrual = 0.5;
    Real firstTime = 0.5;
 
 
    std::vector<Real> rateTimes(numberRates+1);
    for (Size i=0; i < rateTimes.size(); ++i)
        rateTimes[i] = firstTime + i*accrual;
 
    std::vector<Real> paymentTimes(numberRates);
    std::vector<Real> accruals(numberRates,accrual);
    for (Size i=0; i < paymentTimes.size(); ++i)
        paymentTimes[i] = firstTime + (i+1)*accrual;
 
 
 
 
    Real fixedRate = 0.05;
    std::vector<Real> strikes(numberRates,fixedRate);
    Real receive = -1.0;
 
    // 0. a payer swap
    MultiStepSwap payerSwap(rateTimes, accruals, accruals, paymentTimes,
        fixedRate, true);
 
    // 1. the equivalent receiver swap
    MultiStepSwap receiverSwap(rateTimes, accruals, accruals, paymentTimes,
        fixedRate, false);
 
    //exercise schedule, we can exercise on any rate time except the last one
    std::vector<Rate> exerciseTimes(rateTimes);
    exerciseTimes.pop_back();
 
    // naive exercise strategy, exercise above a trigger level
    std::vector<Rate> swapTriggers(exerciseTimes.size(), fixedRate);
    SwapRateTrigger naifStrategy(rateTimes, swapTriggers, exerciseTimes);
 
    // Longstaff-Schwartz exercise strategy
    std::vector<std::vector<NodeData> > collectedData;
    std::vector<std::vector<Real> > basisCoefficients;
 
    // control that does nothing, need it because some control is expected
    NothingExerciseValue control(rateTimes);
 
//    SwapForwardBasisSystem basisSystem(rateTimes,exerciseTimes);
    SwapBasisSystem basisSystem(rateTimes,exerciseTimes);
 
 
 
    // rebate that does nothing, need it because some rebate is expected
    // when you break a swap nothing happens.
    NothingExerciseValue nullRebate(rateTimes);
 
    CallSpecifiedMultiProduct dummyProduct =
        CallSpecifiedMultiProduct(receiverSwap, naifStrategy,
        ExerciseAdapter(nullRebate));
 
    const EvolutionDescription& evolution = dummyProduct.evolution();
 
 
    // parameters for models
 
 
    Size seed = 12332; // for Sobol generator
    Size trainingPaths = 65536;
    Size paths = 16384;
    Size vegaPaths = 16384*64;
 
    std::cout << "training paths, " << trainingPaths << "\n";
    std::cout << "paths, " << paths << "\n";
    std::cout << "vega Paths, " << vegaPaths << "\n";
#ifdef _DEBUG
   trainingPaths = 512;
  paths = 1024;
  vegaPaths = 1024;
#endif
 
 
    // set up a calibration, this would typically be done by using a calibrator
 
 
 
    Real rateLevel =0.05;
 
 
    Real initialNumeraireValue = 0.95;
 
    Real volLevel = 0.11;
    Real beta = 0.2;
    Real gamma = 1.0;
    Size numberOfFactors = std::min<Size>(5,numberRates);
 
    Spread displacementLevel =0.02;
 
    // set up vectors
    std::vector<Rate> initialRates(numberRates,rateLevel);
    std::vector<Volatility> volatilities(numberRates, volLevel);
    std::vector<Spread> displacements(numberRates, displacementLevel);
 
    ExponentialForwardCorrelation correlations(
        rateTimes,volLevel, beta,gamma);
 
 
 
 
    FlatVol  calibration(
        volatilities,
        ext::shared_ptr<PiecewiseConstantCorrelation>(new  ExponentialForwardCorrelation(correlations)),
        evolution,
        numberOfFactors,
        initialRates,
        displacements);
 
    ext::shared_ptr<MarketModel> marketModel(new FlatVol(calibration));
 
    // we use a factory since there is data that will only be known later
    SobolBrownianGeneratorFactory generatorFactory(
        SobolBrownianGenerator::Diagonal, seed);
 
    std::vector<Size> numeraires( moneyMarketMeasure(evolution));
 
    // the evolver will actually evolve the rates
    LogNormalFwdRatePc  evolver(marketModel,
        generatorFactory,
        numeraires   // numeraires for each step
        );
 
    ext::shared_ptr<MarketModelEvolver> evolverPtr(new LogNormalFwdRatePc(evolver));
 
    int t1= clock();
 
    // gather data before computing exercise strategy
    collectNodeData(evolver,
        receiverSwap,
        basisSystem,
        nullRebate,
        control,
        trainingPaths,
        collectedData);
 
    int t2 = clock();
 
 
    // calculate the exercise strategy's coefficients
    genericLongstaffSchwartzRegression(collectedData,
        basisCoefficients);
 
 
    // turn the coefficients into an exercise strategy
    LongstaffSchwartzExerciseStrategy exerciseStrategy(
        basisSystem, basisCoefficients,
        evolution, numeraires,
        nullRebate, control);
 
    //  bermudan swaption to enter into the payer swap
    CallSpecifiedMultiProduct bermudanProduct =
        CallSpecifiedMultiProduct(
        MultiStepNothing(evolution),
        exerciseStrategy, payerSwap);
 
    //  callable receiver swap
    CallSpecifiedMultiProduct callableProduct =
        CallSpecifiedMultiProduct(
        receiverSwap, exerciseStrategy,
        ExerciseAdapter(nullRebate));
 
    // lower bound: evolve all 4 products togheter
    MultiProductComposite allProducts;
    allProducts.add(payerSwap);
    allProducts.add(receiverSwap);
    allProducts.add(bermudanProduct);
    allProducts.add(callableProduct);
    allProducts.finalize();
 
    AccountingEngine accounter(evolverPtr,
        Clone<MarketModelMultiProduct>(allProducts),
        initialNumeraireValue);
 
    SequenceStatisticsInc stats;
 
    accounter.multiplePathValues (stats,paths);
 
    int t3 = clock();
 
    std::vector<Real> means(stats.mean());
 
    for (Size i=0; i < means.size(); ++i)
        std::cout << means[i] << "\n";
 
    std::cout << " time to build strategy, " << (t2-t1)/static_cast<Real>(CLOCKS_PER_SEC)<< ", seconds.\n";
    std::cout << " time to price, " << (t3-t2)/static_cast<Real>(CLOCKS_PER_SEC)<< ", seconds.\n";
 
    // vegas
 
    // do it twice once with factorwise bumping, once without
    Size pathsToDoVegas = vegaPaths;
 
    for (Size i=0; i < 4; ++i)
    {
 
        bool allowFactorwiseBumping = i % 2 > 0 ;
 
        bool doCaps = i / 2 > 0 ;
 
 
 
 
 
        LogNormalFwdRateEuler evolverEuler(marketModel,
            generatorFactory,
            numeraires
            ) ;
 
        MarketModelPathwiseSwap receiverPathwiseSwap(  rateTimes,
            accruals,
            strikes,
            receive);
        Clone<MarketModelPathwiseMultiProduct> receiverPathwiseSwapPtr(receiverPathwiseSwap.clone());
 
        //  callable receiver swap
        CallSpecifiedPathwiseMultiProduct callableProductPathwise(receiverPathwiseSwapPtr,
            exerciseStrategy);
 
        Clone<MarketModelPathwiseMultiProduct> callableProductPathwisePtr(callableProductPathwise.clone());
 
 
        std::vector<std::vector<Matrix> > theBumps(theVegaBumps(allowFactorwiseBumping,
            marketModel,
            doCaps));
 
        PathwiseVegasOuterAccountingEngine
            accountingEngineVegas(ext::make_shared<LogNormalFwdRateEuler>(evolverEuler),
            callableProductPathwisePtr,
            marketModel,
            theBumps,
            initialNumeraireValue);
 
        std::vector<Real> values,errors;
 
        accountingEngineVegas.multiplePathValues(values,errors,pathsToDoVegas);
 
 
        std::cout << "vega output \n";
        std::cout << " factorwise bumping " << allowFactorwiseBumping << "\n";
        std::cout << " doCaps " << doCaps << "\n";
 
 
 
        Size r=0;
 
        std::cout << " price estimate, " << values[r++] << "\n";
 
        for (Size i=0; i < numberRates; ++i, ++r)
            std::cout << " Delta, " << i << ", " << values[r] << ", " << errors[r] << "\n";
 
        Real totalVega = 0.0;
 
        for (; r < values.size(); ++r)
        {
            std::cout << " vega, " << r - 1 -  numberRates<< ", " << values[r] << " ," << errors[r] << "\n";
            totalVega +=  values[r];
        }
 
        std::cout << " total Vega, " << totalVega << "\n";
    }
 
    // upper bound
 
    MTBrownianGeneratorFactory uFactory(seed+142);
 
    ext::shared_ptr<MarketModelEvolver> upperEvolver(new LogNormalFwdRatePc( ext::shared_ptr<MarketModel>(new FlatVol(calibration)),
            uFactory,
            numeraires   // numeraires for each step
            ));
 
    std::vector<ext::shared_ptr<MarketModelEvolver> > innerEvolvers;
 
    std::valarray<bool> isExerciseTime =   isInSubset(evolution.evolutionTimes(),    exerciseStrategy.exerciseTimes());
 
    for (Size s=0; s < isExerciseTime.size(); ++s)
    {
        if (isExerciseTime[s])
        {
            MTBrownianGeneratorFactory iFactory(seed+s);
            ext::shared_ptr<MarketModelEvolver> e =ext::shared_ptr<MarketModelEvolver> (static_cast<MarketModelEvolver*>(new   LogNormalFwdRatePc(ext::shared_ptr<MarketModel>(new FlatVol(calibration)),
                    uFactory,
                    numeraires ,  // numeraires for each step
                    s)));
 
            innerEvolvers.push_back(e);
        }
    }
 
 
 
    UpperBoundEngine uEngine(upperEvolver,  // does outer paths
                             innerEvolvers, // for sub-simulations that do continuation values
                             receiverSwap,
                             nullRebate,
                             receiverSwap,
                             nullRebate,
                             exerciseStrategy,
                             initialNumeraireValue);
 
    Statistics uStats;
    Size innerPaths = 255;
    Size outerPaths =256;
 
    int t4 = clock();
 
    uEngine.multiplePathValues(uStats,outerPaths,innerPaths);
    Real upperBound = uStats.mean();
    Real upperSE = uStats.errorEstimate();
 
    int t5=clock();
 
    std::cout << " Upper - lower is, " << upperBound << ", with standard error " << upperSE << "\n";
    std::cout << " time to compute upper bound is,  " << (t5-t4)/static_cast<Real>(CLOCKS_PER_SEC) << ", seconds.\n";
 
    return 0;
}
 
int InverseFloater(Real rateLevel)
{
 
    Size numberRates =20;
    Real accrual = 0.5;
    Real firstTime = 0.5;
 
    Real strike =0.15;
    Real fixedMultiplier = 2.0;
    Real floatingSpread =0.0;
    bool payer = true;
 
 
    std::vector<Real> rateTimes(numberRates+1);
    for (Size i=0; i < rateTimes.size(); ++i)
        rateTimes[i] = firstTime + i*accrual;
 
    std::vector<Real> paymentTimes(numberRates);
    std::vector<Real> accruals(numberRates,accrual);
    std::vector<Real> fixedStrikes(numberRates,strike);
    std::vector<Real> floatingSpreads(numberRates,floatingSpread);
    std::vector<Real> fixedMultipliers(numberRates,fixedMultiplier);
 
    for (Size i=0; i < paymentTimes.size(); ++i)
        paymentTimes[i] = firstTime + (i+1)*accrual;
 
  MultiStepInverseFloater inverseFloater(
                                                        rateTimes,
                                                        accruals,
                                                         accruals,
                                                        fixedStrikes,
                                                        fixedMultipliers,
                                                        floatingSpreads,
                                                         paymentTimes,
                                                         payer);
 
 
 
 
    //exercise schedule, we can exercise on any rate time except the last one
    std::vector<Rate> exerciseTimes(rateTimes);
    exerciseTimes.pop_back();
 
    // naive exercise strategy, exercise above a trigger level
    Real trigger =0.05;
    std::vector<Rate> swapTriggers(exerciseTimes.size(), trigger);
    SwapRateTrigger naifStrategy(rateTimes, swapTriggers, exerciseTimes);
 
    // Longstaff-Schwartz exercise strategy
    std::vector<std::vector<NodeData> > collectedData;
    std::vector<std::vector<Real> > basisCoefficients;
 
    // control that does nothing, need it because some control is expected
    NothingExerciseValue control(rateTimes);
 
   SwapForwardBasisSystem basisSystem(rateTimes,exerciseTimes);
//    SwapBasisSystem basisSystem(rateTimes,exerciseTimes);
 
 
 
    // rebate that does nothing, need it because some rebate is expected
    // when you break a swap nothing happens.
    NothingExerciseValue nullRebate(rateTimes);
 
    CallSpecifiedMultiProduct dummyProduct =
        CallSpecifiedMultiProduct(inverseFloater, naifStrategy,
        ExerciseAdapter(nullRebate));
 
    const EvolutionDescription& evolution = dummyProduct.evolution();
 
 
    // parameters for models
 
 
    Size seed = 12332; // for Sobol generator
    Size trainingPaths = 65536;
    Size paths = 65536;
    Size vegaPaths =16384;
 
#ifdef _DEBUG
   trainingPaths = 8192;
  paths = 8192;
  vegaPaths = 1024;
#endif
 
 
    std::cout <<  " inverse floater \n";
    std::cout << " fixed strikes :  "  << strike << "\n";
    std::cout << " number rates :  " << numberRates << "\n";
 
    std::cout << "training paths, " << trainingPaths << "\n";
    std::cout << "paths, " << paths << "\n";
    std::cout << "vega Paths, " << vegaPaths << "\n";
 
 
    // set up a calibration, this would typically be done by using a calibrator
 
 
 
    //Real rateLevel =0.08;
 
    std::cout << " rate level " <<  rateLevel << "\n";
 
    Real initialNumeraireValue = 0.95;
 
    Real volLevel = 0.11;
    Real beta = 0.2;
    Real gamma = 1.0;
    Size numberOfFactors = std::min<Size>(5,numberRates);
 
    Spread displacementLevel =0.02;
 
    // set up vectors
    std::vector<Rate> initialRates(numberRates,rateLevel);
    std::vector<Volatility> volatilities(numberRates, volLevel);
    std::vector<Spread> displacements(numberRates, displacementLevel);
 
    ExponentialForwardCorrelation correlations(
        rateTimes,volLevel, beta,gamma);
 
 
 
 
    FlatVol  calibration(
        volatilities,
        ext::shared_ptr<PiecewiseConstantCorrelation>(new  ExponentialForwardCorrelation(correlations)),
        evolution,
        numberOfFactors,
        initialRates,
        displacements);
 
    ext::shared_ptr<MarketModel> marketModel(new FlatVol(calibration));
 
    // we use a factory since there is data that will only be known later
    SobolBrownianGeneratorFactory generatorFactory(
        SobolBrownianGenerator::Diagonal, seed);
 
    std::vector<Size> numeraires( moneyMarketMeasure(evolution));
 
    // the evolver will actually evolve the rates
    LogNormalFwdRatePc  evolver(marketModel,
        generatorFactory,
        numeraires   // numeraires for each step
        );
 
    ext::shared_ptr<MarketModelEvolver> evolverPtr(new LogNormalFwdRatePc(evolver));
 
    int t1= clock();
 
    // gather data before computing exercise strategy
    collectNodeData(evolver,
        inverseFloater,
        basisSystem,
        nullRebate,
        control,
        trainingPaths,
        collectedData);
 
    int t2 = clock();
 
 
    // calculate the exercise strategy's coefficients
    genericLongstaffSchwartzRegression(collectedData,
        basisCoefficients);
 
 
    // turn the coefficients into an exercise strategy
    LongstaffSchwartzExerciseStrategy exerciseStrategy(
        basisSystem, basisCoefficients,
        evolution, numeraires,
        nullRebate, control);
 
 
    //  callable receiver swap
    CallSpecifiedMultiProduct callableProduct =
        CallSpecifiedMultiProduct(
        inverseFloater, exerciseStrategy,
        ExerciseAdapter(nullRebate));
 
     MultiProductComposite allProducts;
    allProducts.add(inverseFloater);
    allProducts.add(callableProduct);
    allProducts.finalize();
 
 
    AccountingEngine accounter(evolverPtr,
        Clone<MarketModelMultiProduct>(allProducts),
        initialNumeraireValue);
 
    SequenceStatisticsInc stats;
 
    accounter.multiplePathValues (stats,paths);
 
    int t3 = clock();
 
    std::vector<Real> means(stats.mean());
 
    for (Size i=0; i < means.size(); ++i)
        std::cout << means[i] << "\n";
 
    std::cout << " time to build strategy, " << (t2-t1)/static_cast<Real>(CLOCKS_PER_SEC)<< ", seconds.\n";
    std::cout << " time to price, " << (t3-t2)/static_cast<Real>(CLOCKS_PER_SEC)<< ", seconds.\n";
 
    // vegas
 
    // do it twice once with factorwise bumping, once without
    Size pathsToDoVegas = vegaPaths;
 
    for (Size i=0; i < 4; ++i)
    {
 
        bool allowFactorwiseBumping = i % 2 > 0 ;
 
        bool doCaps = i / 2 > 0 ;
 
 
        LogNormalFwdRateEuler evolverEuler(marketModel,
            generatorFactory,
            numeraires
            ) ;
 
        MarketModelPathwiseInverseFloater pathwiseInverseFloater(
                                                         rateTimes,
                                                         accruals,
                                                         accruals,
                                                         fixedStrikes,
                                                         fixedMultipliers,
                                                         floatingSpreads,
                                                         paymentTimes,
                                                         payer);
 
        Clone<MarketModelPathwiseMultiProduct> pathwiseInverseFloaterPtr(pathwiseInverseFloater.clone());
 
        //  callable inverse floater
        CallSpecifiedPathwiseMultiProduct callableProductPathwise(pathwiseInverseFloaterPtr,
                                                                                                                                               exerciseStrategy);
 
        Clone<MarketModelPathwiseMultiProduct> callableProductPathwisePtr(callableProductPathwise.clone());
 
 
        std::vector<std::vector<Matrix> > theBumps(theVegaBumps(allowFactorwiseBumping,
            marketModel,
            doCaps));
 
        PathwiseVegasOuterAccountingEngine
            accountingEngineVegas(ext::make_shared<LogNormalFwdRateEuler>(evolverEuler),
   //         pathwiseInverseFloaterPtr,
            callableProductPathwisePtr,
            marketModel,
            theBumps,
            initialNumeraireValue);
 
        std::vector<Real> values,errors;
 
        accountingEngineVegas.multiplePathValues(values,errors,pathsToDoVegas);
 
 
        std::cout << "vega output \n";
        std::cout << " factorwise bumping " << allowFactorwiseBumping << "\n";
        std::cout << " doCaps " << doCaps << "\n";
 
 
 
        Size r=0;
 
        std::cout << " price estimate, " << values[r++] << "\n";
 
        for (Size i=0; i < numberRates; ++i, ++r)
            std::cout << " Delta, " << i << ", " << values[r] << ", " << errors[r] << "\n";
 
        Real totalVega = 0.0;
 
        for (; r < values.size(); ++r)
        {
            std::cout << " vega, " << r - 1 -  numberRates<< ", " << values[r] << " ," << errors[r] << "\n";
            totalVega +=  values[r];
        }
 
        std::cout << " total Vega, " << totalVega << "\n";
    }
 
    // upper bound
 
    MTBrownianGeneratorFactory uFactory(seed+142);
 
 
    ext::shared_ptr<MarketModelEvolver> upperEvolver(new LogNormalFwdRatePc( ext::shared_ptr<MarketModel>(new FlatVol(calibration)),
            uFactory,
            numeraires   // numeraires for each step
            ));
 
    std::vector<ext::shared_ptr<MarketModelEvolver> > innerEvolvers;
 
    std::valarray<bool> isExerciseTime =   isInSubset(evolution.evolutionTimes(),    exerciseStrategy.exerciseTimes());
 
    for (Size s=0; s < isExerciseTime.size(); ++s)
    {
        if (isExerciseTime[s])
        {
            MTBrownianGeneratorFactory iFactory(seed+s);
            ext::shared_ptr<MarketModelEvolver> e =ext::shared_ptr<MarketModelEvolver> (static_cast<MarketModelEvolver*>(new   LogNormalFwdRatePc(ext::shared_ptr<MarketModel>(new FlatVol(calibration)),
                    uFactory,
                    numeraires ,  // numeraires for each step
                    s)));
 
            innerEvolvers.push_back(e);
        }
    }
 
 
 
    UpperBoundEngine uEngine(upperEvolver,  // does outer paths
                             innerEvolvers, // for sub-simulations that do continuation values
                             inverseFloater,
                             nullRebate,
                             inverseFloater,
                             nullRebate,
                             exerciseStrategy,
                             initialNumeraireValue);
 
    Statistics uStats;
    Size innerPaths = 255;
    Size outerPaths =256;
 
    int t4 = clock();
 
    uEngine.multiplePathValues(uStats,outerPaths,innerPaths);
    Real upperBound = uStats.mean();
    Real upperSE = uStats.errorEstimate();
 
    int t5=clock();
 
    std::cout << " Upper - lower is, " << upperBound << ", with standard error " << upperSE << "\n";
    std::cout << " time to compute upper bound is,  " << (t5-t4)/static_cast<Real>(CLOCKS_PER_SEC) << ", seconds.\n";
 
 
    return 0;
 
}
 
int main()
{
    try {
        for (Size i=5; i < 10; ++i)
            InverseFloater(i/100.0);
 
        return 0;
    } catch (std::exception& e) {
        std::cerr << e.what() << std::endl;
        return 1;
    } catch (...) {
        std::cerr << "unknown error" << std::endl;
        return 1;
    }
}