text: /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
  
  /*
   Copyright (C) 2014 Jose Aparicio
  
   This file is part of QuantLib, a free-software/open-source library
   for financial quantitative analysts and developers - http://quantlib.org/
  
   QuantLib is free software: you can redistribute it and/or modify it
   under the terms of the QuantLib license.  You should have received a
   copy of the license along with this program; if not, please email
   <quantlib-dev@lists.sf.net>. The license is also available online at
   <http://quantlib.org/license.shtml>.
  
   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE.  See the license for more details.
  */
  
  #include <ql/qldefines.hpp>
  #ifdef BOOST_MSVC
  #  include <ql/auto_link.hpp>
  #endif
  #include <ql/experimental/credit/randomdefaultlatentmodel.hpp>
  #include <ql/termstructures/credit/flathazardrate.hpp>
  #include <ql/currencies/europe.hpp>
  #include <ql/time/calendars/target.hpp>
  #include <ql/time/daycounters/actual365fixed.hpp>
  
  #include <boost/timer.hpp>
  #include <boost/function.hpp>
  
  #include <iostream>
  #include <iomanip>
  
  using namespace std;
  using namespace QuantLib;
  
  #if defined(QL_ENABLE_SESSIONS)
  namespace QuantLib {
  
      Integer sessionId() { return 0; }
  
  }
  #endif
  
  
  /* This sample code shows basic usage of a Latent variable model.
     The data and correlation problem presented is the same as in:
       'Modelling Dependent Defaults: Asset Correlations Are Not Enough!'
       Frey R., A. J. McNeil and M. A. Nyfeler RiskLab publications March 2001
  */
  int main(int, char* []) {
  
      try {
  
          boost::timer timer;
          std::cout << std::endl;
  
          Calendar calendar = TARGET();
          Date todaysDate(19, March, 2014);
          // must be a business day
          todaysDate = calendar.adjust(todaysDate);
  
          Settings::instance().evaluationDate() = todaysDate;
  
  
          /* --------------------------------------------------------------
                          SET UP BASKET PORTFOLIO
          -------------------------------------------------------------- */
          // build curves and issuers into a basket of three names
          std::vector<Real> hazardRates(3, -std::log(1.-0.01));
          std::vector<std::string> names;
          for(Size i=0; i<hazardRates.size(); i++)
              names.push_back(std::string("Acme") + 
                  boost::lexical_cast<std::string>(i));
          std::vector<Handle<DefaultProbabilityTermStructure> > defTS;
          for(Size i=0; i<hazardRates.size(); i++)
              defTS.push_back(Handle<DefaultProbabilityTermStructure>(
                  ext::make_shared<FlatHazardRate>(0, TARGET(), hazardRates[i], 
                      Actual365Fixed())));
          std::vector<Issuer> issuers;
          for(Size i=0; i<hazardRates.size(); i++) {
              std::vector<QuantLib::Issuer::key_curve_pair> curves(1, 
                  std::make_pair(NorthAmericaCorpDefaultKey(
                      EURCurrency(), QuantLib::SeniorSec,
                      Period(), 1. // amount threshold
                      ), defTS[i]));
              issuers.push_back(Issuer(curves));
          }
  
          ext::shared_ptr<Pool> thePool = ext::make_shared<Pool>();
          for(Size i=0; i<hazardRates.size(); i++)
              thePool->add(names[i], issuers[i], NorthAmericaCorpDefaultKey(
                      EURCurrency(), QuantLib::SeniorSec, Period(), 1.));
  
          std::vector<DefaultProbKey> defaultKeys(hazardRates.size(), 
              NorthAmericaCorpDefaultKey(EURCurrency(), SeniorSec, Period(), 1.));
          // Recoveries are irrelevant in this example but must be given as the 
          //   lib stands.
          std::vector<ext::shared_ptr<RecoveryRateModel> > rrModels(
              hazardRates.size(), ext::make_shared<ConstantRecoveryModel>(
              ConstantRecoveryModel(0.5, SeniorSec)));
          ext::shared_ptr<Basket> theBskt = ext::make_shared<Basket>(
              todaysDate, names, std::vector<Real>(hazardRates.size(), 100.), 
              thePool);
          /* --------------------------------------------------------------
                          SET UP JOINT DEFAULT EVENT LATENT MODELS
          -------------------------------------------------------------- */
          // Latent model factors, corresponds to the first entry in Table1 of the
          //   publication mentioned. It is a single factor model
          std::vector<std::vector<Real> > fctrsWeights(hazardRates.size(), 
              std::vector<Real>(1, std::sqrt(0.1)));
          // --- Default Latent models -------------------------------------
          #ifndef QL_PATCH_SOLARIS
          // Gaussian integrable joint default model:
          ext::shared_ptr<GaussianDefProbLM> lmG(new 
              GaussianDefProbLM(fctrsWeights, 
              LatentModelIntegrationType::GaussianQuadrature,
  			GaussianCopulaPolicy::initTraits() // otherwise gcc screams
  			));
          #endif
          // Define StudentT copula
          // this is as far as we can be from the Gaussian, 2 T_3 factors:
          std::vector<Integer> ordersT(2, 3);
          TCopulaPolicy::initTraits iniT;
          iniT.tOrders = ordersT;
          // StudentT integrable joint default model:
          ext::shared_ptr<TDefProbLM> lmT(new TDefProbLM(fctrsWeights, 
              // LatentModelIntegrationType::GaussianQuadrature,
              LatentModelIntegrationType::Trapezoid,
              iniT));
  
          // --- Default Loss models ----------------------------------------
          // Gaussian random joint default model:
          Size numSimulations = 100000;
          // Size numCoresUsed = 4;
          #ifndef QL_PATCH_SOLARIS
          // Sobol, many cores
          ext::shared_ptr<DefaultLossModel> rdlmG(
              ext::make_shared<RandomDefaultLM<GaussianCopulaPolicy> >(lmG, 
                  std::vector<Real>(), numSimulations, 1.e-6, 2863311530UL));
          #endif
          // StudentT random joint default model:
          ext::shared_ptr<DefaultLossModel> rdlmT(
              ext::make_shared<RandomDefaultLM<TCopulaPolicy> >(lmT, 
              std::vector<Real>(), numSimulations, 1.e-6, 2863311530UL));
  
          /* --------------------------------------------------------------
                          DUMP SOME RESULTS
          -------------------------------------------------------------- */
          /* Default correlations in a T copula should be below those of the 
          gaussian for the same factors.
          The calculations on the MC show dispersion on both copulas (thats
          ok) and too large values with very large dispersions on the T case.
          Computations are ok, within the dispersion, for the gaussian; compare
          with the direct integration in both cases.
          However the T does converge to the gaussian value for large value of
          the parameters.
          */
          Date calcDate(TARGET().advance(Settings::instance().evaluationDate(), 
              Period(120, Months)));
          std::vector<Probability> probEventsTLatent, probEventsTRandLoss;
          #ifndef QL_PATCH_SOLARIS
          std::vector<Probability> probEventsGLatent, probEventsGRandLoss;
          #endif
          //
          lmT->resetBasket(theBskt);
          for(Size numEvts=0; numEvts <=theBskt->size(); numEvts++) {
              probEventsTLatent.push_back(lmT->probAtLeastNEvents(numEvts, 
                  calcDate));
           }
          //
          theBskt->setLossModel(rdlmT);
          for(Size numEvts=0; numEvts <=theBskt->size(); numEvts++) {
              probEventsTRandLoss.push_back(theBskt->probAtLeastNEvents(numEvts, 
                  calcDate));
           }
          //
          #ifndef QL_PATCH_SOLARIS
          lmG->resetBasket(theBskt);
          for(Size numEvts=0; numEvts <=theBskt->size(); numEvts++) {
              probEventsGLatent.push_back(lmG->probAtLeastNEvents(numEvts, 
                  calcDate));
           }
          //
          theBskt->setLossModel(rdlmG);
          for(Size numEvts=0; numEvts <=theBskt->size(); numEvts++) {
              probEventsGRandLoss.push_back(theBskt->probAtLeastNEvents(numEvts, 
                  calcDate));
           }
          #endif
  
          Date correlDate = TARGET().advance(
              Settings::instance().evaluationDate(), Period(12, Months));
          std::vector<std::vector<Real> > correlsGlm, correlsTlm, correlsGrand, 
              correlsTrand;
          //
          lmT->resetBasket(theBskt);
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              std::vector<Real> tmp;
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  tmp.push_back(lmT->defaultCorrelation(correlDate, 
                      iName1, iName2));
              correlsTlm.push_back(tmp);
          }
          //
          theBskt->setLossModel(rdlmT);
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              std::vector<Real> tmp;
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  tmp.push_back(theBskt->defaultCorrelation(correlDate, 
                      iName1, iName2));
              correlsTrand.push_back(tmp);
          }
          #ifndef QL_PATCH_SOLARIS
          //
          lmG->resetBasket(theBskt);
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              std::vector<Real> tmp;
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  tmp.push_back(lmG->defaultCorrelation(correlDate, 
                      iName1, iName2));
              correlsGlm.push_back(tmp);
          }
          //
          theBskt->setLossModel(rdlmG);
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              std::vector<Real> tmp;
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  tmp.push_back(theBskt->defaultCorrelation(correlDate, 
                      iName1, iName2));
              correlsGrand.push_back(tmp);
          }
          #endif
  
  
          std::cout << 
              " Gaussian versus T prob of extreme event (random and integrable)-" 
              << std::endl;
          for(Size numEvts=0; numEvts <=theBskt->size(); numEvts++) {
              std::cout << "-Prob of " << numEvts << " events... " <<
                  #ifndef QL_PATCH_SOLARIS
                  probEventsGLatent[numEvts] << " ** " <<
                  #else
                  "n/a" << " ** " <<
                  #endif
                  probEventsTLatent[numEvts] << " ** " << 
                  #ifndef QL_PATCH_SOLARIS
                  probEventsGRandLoss[numEvts]<< " ** " <<
                  #else
                  "n/a" << " ** " <<
                  #endif
                  probEventsTRandLoss[numEvts] 
              << std::endl;
          }
  
          cout << endl;
          cout << "-- Default correlations G,T,GRand,TRand--" << endl;
          cout << "-----------------------------------------" << endl;
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  cout << 
                  #ifndef QL_PATCH_SOLARIS
                      correlsGlm[iName1][iName2] << " , ";
                  #else
                      "n/a" << " , ";
                  #endif
              cout << endl;
          }
          cout << endl;
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  cout << 
                      correlsTlm[iName1][iName2] << " , ";
              ;
                  cout << endl;
          }
          cout << endl;
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  cout << 
                  #ifndef QL_PATCH_SOLARIS
                      correlsGrand[iName1][iName2] << " , ";
                  #else
                      "n/a" << " , ";
                  #endif
              cout << endl;
          }
          cout << endl;
          for(Size iName1=0; iName1 <theBskt->size(); iName1++) {
              for(Size iName2=0; iName2 <theBskt->size(); iName2++)
                  cout << 
                      correlsTrand[iName1][iName2] << " , ";
              ;
                  cout << endl;
          }
  
  
  
          Real seconds  = timer.elapsed();
          Integer hours = Integer(seconds/3600);
          seconds -= hours * 3600;
          Integer minutes = Integer(seconds/60);
          seconds -= minutes * 60;
          cout << "Run completed in ";
          if (hours > 0)
              cout << hours << " h ";
          if (hours > 0 || minutes > 0)
              cout << minutes << " m ";
          cout << fixed << setprecision(0)
               << seconds << " s" << endl;
  
          return 0;
      } catch (exception& e) {
          cerr << e.what() << endl;
          return 1;
      } catch (...) {
          cerr << "unknown error" << endl;
          return 1;
      }
  }
  
  