Financial Modeling with Mathematica and QuantLib

Thomas Weber

Weber & Partner GbR

Introduction

About QuantLib

QuantLib is a free/open-source library for quantitative finance written in C++. It offers tools for practical implementation and for advanced modeling, with features such as market conventions, yield curve models, solvers, PDEs, Monte Carlo, exotic options, VAR, and so on, see http://www.quantlib.org.

About us

Since 1997 we are consulting banks, energy companies and other financial institutions in the area of financial modeling. For that we use quite often Mathematica

as a prototyping environment

for model testing and validation

data analysis and exploration tool

Mathematica is very powerful for that, but some extensions are helpful. That is why we have done it with QuantLib again.

Getting Started

Get the package from http://www.weber-und-partner.com/QuantLib/mma. and install it.

The first step is to load the QuantLib package.

Needs["QuantLib`"]

Getting help >>

Names["QuantLib`*"] // Length

QLCategories // Shallow

?PricingEngines

?BondEngine

What You Get

Boring Stuff, You Expect But You Are Not Interested In

Extended Date functionality (DateEndOfMonth,...,DateNextWeekday, DateNthWeekday, ... ,IMMNextDate,...)

Calendar functionailty (CalendarAddHoliday, ...,CalendarAdvance  ..., CalendarIsBusinessDay,..., CalendarIsHoliday, ...)

Extended DayCounting (30/360,...,30/360 (Italian), Actual/360,Actual/365,..., Actual/Actual (Euro),Actual/Actual (ISDA),Actual/Actual (ISMA)...)

More Boring Stuff, You Might Have Some Interest

Basic Instruments like Bonds, Futures, Deposits, Swaps, Swaptions, Vanilla Options ...

Indices

Curves and Surfaces (Yield Curve construction, Volatility surfaces

More Stuff

Standard models

Standard Pricing Engines

...

The idea is to make everything as low level as possible to avoid black boxes. We will imporve on that.

A Simple Example as a Starter

Keep It Simple: Just use a formula

?BlackFormula

With most functions there are options associated. This simplifies the calling of functions.

Options[BlackFormula]

σ=0.20; (* diffusion parameter *)
S=36;  (* Underlying *)
K=40; (* Strike *)
r =0.06;(* interest rate *)
T = 1;  (* Time to Maturity *)
BlackFormula["Call",K,S,σ(Sqrt[T]),Deflator->(E^(-r T))]

Make It Complicated

Defining the Driving Process

settlementDate =DateList[];
maturityDate =CalendarAdvance[settlementDate,{1,"Year"}]

yieldcurve = FlatForward[r] 
vol =BlackConstantVol[settlementDate,SimpleQuote[σ]]
underlying = SimpleQuote[S]

?BlackProcess

process =BlackProcess[underlying,yieldcurve,vol];

Defining the Instrument

Payoff functions are collected in the category Payoffs.

Category[Payoffs]

In this example a simple striked payoff is assumed.

StrikedTypePayoff[PayoffID,"Vanilla","Call",K];

Available functions within the

Plot[PayoffValue[PayoffID,S],{S,10,60},Frame->True]

ExerciseID =EuropeanExercise[maturityDate]

europeanOption =VanillaOption[PayoffID,ExerciseID];

Doing the Valuation

pricingEngine = AnalyticEuropeanEngine[process]

InstrumentSetPricingEngine[europeanOption,pricingEngine]

InstrumentNPV[europeanOption]

Some More Advanced Examples

Pricing Bermudan Swaption

Describing the market conditions

settlementDate=CalendarAdvance[EvaluationDate[],{2,"Day"}]
yieldcurve =FlatForward[settlementDate,0.05,DayCounter->"Actual/365 (Fixed)"];

swapLengths={1,2,3,4,5};
swaptionVols =(({{0.1490, 0.1340, 0.1228, 0.1189, 0.1148}, {0.1290, 0.1201, 0.1146, 0.1108, 0.1040}, {0.1149, 0.1112, 0.1070, 0.1010, 0.0957}, {0.1047, 0.1021, 0.0980, 0.0951, 0.1270}, {0.1000, 0.0950, 0.0999, 0.1230, 0.1160}}));

ListPointPlot3D[swaptionVols,ColorFunction->"Rainbow",PlotStyle->PointSize[Large],Filling->Bottom]

Constructing the swaptions

indexSixMonths = Euribor["6M",YieldCurve->yieldcurve];
swaptions= Table[SwaptionHelper[{swapLengths[[i]],"Year"},{swapLengths[[-i]],"Year"},swaptionVols[[i,-i]],indexSixMonths,yieldcurve],{i,Length[swapLengths]},
{j,Length[swapLengths],1,-1}] ;

Defining Short Rate Models

One Factor Hull/White Model

modelHW =HullWhite[yieldcurve];
Map[CalibrationHelperSetPricingEngine[#,JamshidianSwaptionEngine[modelHW]]&,swaptions,{2}];
CalibrateModel[modelHW,swaptions // Flatten]

error = Map[CalibrationHelperImpliedVolatility[#,CalibrationHelperModelValue[#]]&,swaptions,{2}]-swaptionVols // Abs;
Total[error// Flatten]

One Factor Black/Karasinski Model

modelBK = BlackKarasinski[yieldcurve];(* One Factor Black/Karasinski *)
times = Map[CalibrationHelperTimes,swaptions,{2}] // Flatten;
timeGrid=TimeGrid[times,30];
Map[CalibrationHelperSetPricingEngine[#,TreeSwaptionEngine[modelBK,timeGrid]]&,swaptions,{2}];
CalibrateModel[modelBK,swaptions // Flatten]
Map[CalibrationHelperImpliedVolatility[#,CalibrationHelperModelValue[#]]&,swaptions,{2}]-swaptionVols // Abs // Flatten // Total

Gaussian two-factor Model

modelG2 = G2[yieldcurve]; (* Gaussian two-factor model *)
Map[CalibrationHelperSetPricingEngine[#,G2SwaptionEngine[modelG2,6,20]]&,swaptions,{2}];
CalibrateModel[modelG2,swaptions // Flatten]
Map[CalibrationHelperImpliedVolatility[#,CalibrationHelperModelValue[#]]&,swaptions,{2}]-swaptionVols // Abs // Flatten // Total

Constructing the bermudan swaption

startDate = CalendarAdvance[settlementDate,{1,"Year"},BusinessDayConvention-> "Modified Following"];
maturityDate =CalendarAdvance[startDate,{5,"Year"},BusinessDayConvention-> "Modified Following"];

fixedSchedule = Schedule[startDate,maturityDate,"1Y",GenRule->"Forward",BusinessDayConvention->  "Unadjusted"];
floatSchedule = Schedule[startDate,maturityDate,"6M",GenRule->"Forward",BusinessDayConvention->  "Modified Following"];

swap = VanillaSwap[1000 ,fixedSchedule,0.05,floatSchedule,indexSixMonths,"Payer",FixedLegDayCounter->"Actual/360"];

InstrumentSetPricingEngine[swap,DiscountingSwapEngine[yieldcurve]];
fixedATMRate =VanillaSwapFairRate[swap]

VanillaSwapFixedLegAnalysis[swap] /.( x_/; IntegerQ[x] && 30000 < x < 60000)  :>  DateString[ ToCalendar[x],"DateShort"]  //  TableForm

accrualDates =VanillaSwapFixedLegAnalysis[swap][[All,4]] // Rest // ToCalendar

bermudanSwaption = Swaption[swap,BermudanExercise[accrualDates]];

InstrumentSetPricingEngine[bermudanSwaption,TreeSwaptionEngine[modelHW,50]] ;
InstrumentNPV[bermudanSwaption]

InstrumentSetPricingEngine[bermudanSwaption,TreeSwaptionEngine[modelBK,50]] ;
InstrumentNPV[bermudanSwaption]

InstrumentSetPricingEngine[bermudanSwaption,TreeSwaptionEngine[modelG2,50]] ;
InstrumentNPV[bermudanSwaption]

Map[CalibrationHelperSetPricingEngine[#,G2SwaptionEngine[modelG2,6,20]]&,swaptions,{2}];

{α,σ,β,η,ρ} =CalibratedModelParams[modelG2]

Plot3D[CalibratedModelSetParams[modelG2,{α,σ,β,η,ρ}];Abs[Map[
CalibrationHelperImpliedVolatility[#,CalibrationHelperModelValue[#]]&,swaptions,{2}]-swaptionVols] // Flatten // Total,{σ,0.01,0.1},{η,0.01,0.1}]

Static Replication

The QuantLib example "Replication.cpp" the Mathematica way. For the theory see Mark Joshi’s “The Concepts and Practice of Mathematical Finance”.

Constructing the barrier option to replicate

today = {2006,5,29};
SetEvaluationDate[today];

maturity = CalendarAdvance[today,{12,"Month"}];
strike = 100.0;
barrier = 70.0;
rebate = 0.0;
underlyingValue= SimpleQuote[100.0];

flatRate = FlatForward[0.04];
flatVol = BlackConstantVol[0.2];
exercise =EuropeanExercise[maturity];
payoff = PlainVanillaPayoff["Put",strike];
bsProcess = BlackScholesProcess[underlyingValue,flatRate,flatVol];

referenceOption=BarrierOption["DownOut",barrier,rebate,payoff,exercise];
InstrumentSetPricingEngine[referenceOption,AnalyticBarrierEngine[bsProcess]];

referencePlot =Plot[{If[S>Barrier[referenceOption],PayoffValue[Payoff[referenceOption],S],Rebate[BarrierOption61]]},{S,50,150},AxesLabel->{"S","V"},Filling->Automatic]

Pricing the barrier option

referenceValue=InstrumentNPV[referenceOption]

Constructing the replication portfolio

Constructing the put for the strike

Final payoff first: as shown in Joshi, a put struck at K...

put1 = EuropeanOption[payoff,exercise];
europeanEngine=AnalyticEuropeanEngine[bsProcess];
InstrumentSetPricingEngine[put1,europeanEngine];

pl1 =Plot[PayoffValue[Payoff[put1],S],{S,50,150},Filling->Automatic]

Constructing the digital option

...minus a digital put struck at B of notional K-B...

digitalPayoff = StrikedTypePayoff["CashOrNothing","Put",barrier];
digitalPut = EuropeanOption[digitalPayoff,exercise];
InstrumentSetPricingEngine[digitalPut,europeanEngine];

Plot[PayoffValue[Payoff[put1],S]-(strike-barrier)PayoffValue[Payoff[digitalPut],S],{S,50,150},Filling->Automatic]

Constructing the put for the barrier

... minus a put option struck at B.

lowerPayoff = PlainVanillaPayoff["Put",barrier];
put2 = EuropeanOption[lowerPayoff,exercise];
InstrumentSetPricingEngine[put2,europeanEngine];

Plot[PayoffValue[Payoff[put1],S]-(strike-barrier)PayoffValue[Payoff[digitalPut],S]-PayoffValue[Payoff[put2],S],{S,50,150},Filling->Automatic]

Putting all insrtruments in the portfolio

pfo0 ={put1,digitalPut,put2}
w0 = {1,-(strike-barrier),-1}

Construction puts at points (B,t)

Hedges for the barrier

Now we need puts at the barrier to kill the value of the portfolio at that point. For the first
portfolio, we’ll use 12 dates at one-month’s distance.

innerMaturities =CalendarAdvance[today,Table[{i,"Weeks"},{i,1,52,4}]]

putn = Map[EuropeanOption[lowerPayoff,EuropeanExercise[#]]& ,innerMaturities];Map[InstrumentSetPricingEngine[#,europeanEngine]&,putn];

pfo = Join[pfo0,putn];
w = Join[w0, Table[0,{Length[putn]}]]

killdates =Most[innerMaturities]

SimpleQuoteSetValue[underlyingValue,barrier];

Do[
SetEvaluationDate[killdates[[-i]] ];
portfolioValue= w.Map[InstrumentNPV,pfo];
(* ...finally,we estimate the notional that kills the portfolio value at that point... *)
(* ...and we subtract from the portfolio a put with such notional. *)
    w[[-i]] =-  portfolioValue/InstrumentNPV[putn[[-i]]],{i,1,Length[killdates]}]

SetEvaluationDate[today];
SimpleQuoteSetValue[underlyingValue,100.];
portfolioValue=w.Map[InstrumentNPV,pfo];
{portfolioValue, referenceValue =InstrumentNPV[referenceOption],portfolioValue- referenceValue}

HedgePerformance[k_,Spot_]:= Module[{pfo,w},
innerMaturities =CalendarAdvance[today,Table[{i,"Weeks"},{i,1,52,k}]];
    putn = Map[EuropeanOption[lowerPayoff,EuropeanExercise[#]]& ,innerMaturities];Map[InstrumentSetPricingEngine[#,europeanEngine]&,putn];
     pfo = Join[pfo0,putn];
w = Join[w0, Table[0,{Length[putn]}]];
    SimpleQuoteSetValue[underlyingValue,barrier];
killdates =Most[innerMaturities];
Do[
SetEvaluationDate[killdates[[-i]] ];
portfolioValue= w.Map[InstrumentNPV,pfo];
 w[[-i]] =-  portfolioValue/InstrumentNPV[putn[[-i]]],{i,1,Length[killdates]}]SetEvaluationDate[today];
SimpleQuoteSetValue[underlyingValue,Spot];
portfolioValue=w.Map[InstrumentNPV,pfo];
{portfolioValue, referenceValue =InstrumentNPV[referenceOption],portfolioValue- referenceValue}
]

HedgePerformance[1,100]

ListPlot[Table[HedgePerformance[k,100][[3]],{k,1,4}],Joined->True]

res =Table[HedgePerformance[k,S][[3]],{k,1,4},{S,80,120,10}]

ListPlot3D[res]

w

What is Next?

Generic Pricing Engines

Generic Instruments