## amortizing interest rate swap valuation example quantib python

Example of valuation of amortizing interest rate swap in Python with quantlib module

for quantlib excel version see Amortizing interest rate swap valuation excel quantlib addin
First , Install QuantLib package for python (use guide here )
to run the program just copy paste it into iPython editor and press enter

#example amortizing swap valuation
#http://www.pricederivatives.com
from  QuantLib import *
import numpy as np
from math import *

todaysDate=Date(31,12,2013)
startDate=todaysDate
Settings.instance().evaluationDate=todaysDate;
crvToday=FlatForward(todaysDate,0.0121,Actual360())
index = Euribor(Period("6m"),forecastTermStructure)
maturity = Date(31,12,2018);
schedule = Schedule(startDate, maturity,Period("6m"),TARGET(),ModifiedFollowing,ModifiedFollowing,DateGeneration.Forward, False)
nominals=[100.0]*10
couponRates=[0.05]*10
floatingleg=IborLeg(nominals,schedule,index,Actual360())
fixedleg=FixedRateLeg(schedule,Actual360(),nominals,couponRates)
firstPeriodDayCount = DayCounter())
swap1=Swap(floatingleg,fixedleg)
swapEngine = DiscountingSwapEngine(discountTermStructure)
swap1.setPricingEngine(swapEngine)
npv=swap1.NPV()
print npv

Posted in quantlib Tagged with: ,

## Derivatives CVA calculation example Monte-Carlo with python

Here we’ll show an example of code for CVA calculation (credit valuation adjustment) using python and Quantlib with simple Monte-Carlo method with portfolio consisting just of a single interest rate swap.It’s easy to generalize code to include more financial instruments , supported by QuantLib python Swig interface.

# CVA calculation algorithm:

1) Simulate yield curve at future dates

2) Calculate your derivatives portfolio NPV (net present value) at each time point for each scenario

3) Calculate CVA as sum of Expected Exposure multiplied by probability of default at this interval

$$CVA = (1-R) \int{DF(t) EE(t) dQ_t}$$

where R is Recovery (normally set to 40%) EE(t) expected exposure at time t and dQ survival probability density , DF(t) discount factor at time t

## Outline

1) In this simple example we will use Hull White model to generate future yield curves. In practice many banks use some yield curve evolution models based on historical simulations.
In Hull White model short rate $$r_t$$ is distributed normally with known mean and variance.

2) For each point of time we will generate whole yield curve based on short rate. Then we will price our interest rate swap on each of these curves.
after

3) to approximate CVA we will use BASEL III formula for regulatory capital charge approximating default probability [or survival probability ] as exp(-ST/(1-R))
so we get
$$CVA= (1-R) \sum \frac{EE^{*}_{T_{i}}+EE^{*}_{T_{i-1}}}{2} (e^{-\frac{ST_{i-1}}{1-R}}-e^{-\frac{ST_{i}}{1-R}})^+$$
where $$EE^{*}$$ is discounted Expected exposure of portfolio

## Details

For this first example we’ll take 5% flat forward yield curve. In the second test example below we’ll take full curve as of 26/12/2013.

### 1) Hull-White model for future yield curve simulations

the model is given by dynamics:
$$dr_t=(\theta_t-ar_t)dt+\sigma dW_t$$
We will use that in Hull White model short rate $$r_t$$ is distributed normally with mean and variance given by
$$E(r_t|r_s)=r_se^{-a(t-s)}+\gamma_t-\gamma_se^{-a(t-s)}$$
$$Var(r_t|r_s)=\frac{\sigma^2}{2a}(1-e^{-2a(t-s)})$$
where $$\gamma_t=f_t(0)+\frac{\sigma^2}{2a}(1-e^{-at})^2$$
and $$f_t(0)$$ is instantaneous forward rate at time t as seen at time 0.
The calculations will not depend on $$\theta_t$$.
To generate future values of $$r_t$$ first we will generate matrix of standard normal variables using numpy library function numpy.random.standard_normal().After that we will apply transform $$mean+\sqrt{variance} Normal(0,1)$$ to get $$r_t$$
after getting matrix of all $$r_t$$ for all simulations and time points we will construct yield curve for each $$r_t$$ using Hull White discount factors $$P_{->T}(t)=A(t,T)e^{-B(t,T)r_t}$$ for 20 years and store each yield curve in matrix crvMat

### 2) Valuing portfolio

On each of this curve we will value swap using QuantLib function. you can add more financial instruments just by adding its value to global NPV

For our exzample portfolio we’ll take one interest rate swap EUR 10MM notional receiving 5% every 6m , TARGET calendar, with 5 years maturity.
Actual/360 daycounter for both legs.

### 3) CVA calculation

After getting matrix of all NPV at each point for each simulation we will replace negative values with 0.
Then we average each column of the matrix (corresponding to averaging for each time point) to get Expected Exposure.
and finally calculate CVA as sum as in BASEL 3 formula.
here we’ll take 500bps flat CDS spread.

# Python Code for CVA calculation

First , Install QuantLib package for python (use guide here )
to run the program just copy paste it into iPython editor and press enter

#(c) 2013 PriceDerivatives.com
# Distributed without any warranty.For educational purposes only.

from  QuantLib import *
import numpy as np
from math import *
from pylab import *

def A(t,T):
forward = crvToday.forwardRate(t, t,Continuous, NoFrequency).rate()
value = B(t,T)*forward - 0.25*sigma*B(t,T)*sigma*B(t,T)*B(0.0,2.0*t);
return exp(value)*crvToday.discount(T)/crvToday.discount(t);

def B(t,T):
return (1.0-exp(-a*(T-t)))/a;

def gamma(t):
forwardRate =crvToday.forwardRate(t, t, Continuous, NoFrequency).rate()
temp = sigma*(1.0 - exp(-a*t))/a
return (forwardRate + 0.5*temp*temp)

def gamma_v(t):
res=np.zeros(len(t))
for i in xrange(len(t)):
res[i]=gamma(t[i])
return res

Nsim=100
a=0.1
sigma=0.02

todaysDate=Date(30,12,2013);
Settings.instance().evaluationDate=todaysDate;
crvToday=FlatForward(todaysDate,0.02,Actual360())

r0=forwardRate =crvToday.forwardRate(0,0, Continuous, NoFrequency).rate()
months=range(1,12*5,1)
sPeriods=[str(month)+"m" for month in months]
Dates=[todaysDate]+[todaysDate+Period(s) for s in sPeriods]
T=+[Actual360().yearFraction(todaysDate,Dates[i]) for i in xrange(1,len(Dates))]
T=np.array(T)
rmean=r0*np.exp(-a*T)+ gamma_v(T) -gamma(0)*np.exp(-a*T)

np.random.seed(1)
stdnorm = np.random.standard_normal(size=(Nsim,len(T)-1))

rmat=np.zeros(shape=(Nsim,len(T)))
rmat[:,0]=r0

for iSim in xrange(Nsim):
for iT in xrange(1,len(T)):
mean=rmat[iSim,iT-1]*exp(-a*(T[iT]-T[iT-1]))+gamma(T[iT])-gamma(T[iT-1])*exp(-a*(T[iT]-T[iT-1]))
var=0.5*sigma*sigma/a*(1-exp(-2*a*(T[iT]-T[iT-1])))
rmat[iSim,iT]=mean+stdnorm[iSim,iT-1]*sqrt(var)

startDate=Date(30,12,2013);

crvMat= [ [ 0 for i in xrange(len(T)) ] for iSim in range(Nsim) ]
npvMat= [ [ 0 for i in xrange(len(T)) ] for iSim in range(Nsim) ]

for row in crvMat:
row=crvToday

for iT in xrange(1,len(T)):
for iSim in xrange(Nsim):
crvDate=Dates[iT];
crvDates=[crvDate]+[crvDate+Period(k,Years) for k in xrange(1,21)]
crvDiscounts=[1.0]+[A(T[iT],T[iT]+k)*exp(-B(T[iT],T[iT]+k)*rmat[iSim,iT]) for k in xrange(1,21)]
crvMat[iSim][iT]=DiscountCurve(crvDates,crvDiscounts,Actual360(),TARGET())

#indexes definitions
index = Euribor(Period("6m"),forecastTermStructure)

#swap 1 definition
maturity = Date(30,12,2018);
fixedSchedule = Schedule(startDate, maturity,Period("6m"), TARGET(),ModifiedFollowing,ModifiedFollowing,DateGeneration.Forward, False)
floatingSchedule = Schedule(startDate, maturity,Period("6m"),TARGET() ,ModifiedFollowing,ModifiedFollowing,DateGeneration.Forward, False)
swap1 = VanillaSwap(VanillaSwap.Receiver, 1000000,fixedSchedule,0.05 , Actual360(),floatingSchedule, index, 0,Actual360())  #0.01215

for iT in xrange(len(T)):
Settings.instance().evaluationDate=Dates[iT]
allDates= list(floatingSchedule)
fixingdates=[index.fixingDate(floatingSchedule[iDate]) for iDate in xrange(len(allDates)) if index.fixingDate(floatingSchedule[iDate])<=Dates[iT]]
if fixingdates:
for date in fixingdates[:-1]:
except:pass
except:pass
swapEngine = DiscountingSwapEngine(discountTermStructure)
swap1.setPricingEngine(swapEngine)

for iSim in xrange(Nsim):
crv=crvMat[iSim][iT]
npvMat[iSim][iT]=swap1.NPV()

npvMat=np.array(npvMat)
npvMat[npvMat<0]=0
EE=np.mean(npvMat,axis=0)

R=0.4  #Recovery rate 40%

sum=0
#calculate CVA
for i in xrange(len(T)-1):
sum=sum+0.5*(EE[i]*crvToday.discount(T[i])+EE[i+1]*crvToday.discount(T[i+1]))*(exp(-S*T[i]/(1.0-R))-exp(-S*T[i+1]/(1.0-R)))
CVA=(1.0-R)*sum

print "\nCVA=",CVA

plot(T,EE)
print "\nEE\n",npvMat
show()


## Output CVA expected exposure swap

for simple current net exposure method see Current Net Exposure excel example

# Tests with 2% swap

For test results lets take 2% EUR 10MM notional interest rate receiver swap with semiannual payments fixed and floating leg.
To check validity of the python code we can check several quantities:

$$E(r_{t_i})$$ for each future date

$$P_{\rightarrow T}(0)= E( e^{-\int_0^Tr_sds})=E( E(e^{-\int_0^Tr_sds}|F_{t_i}))=E(P_{\rightarrow t_i}(0)E(e^{-\int_{t_i}^Tr_sds}|F_{t_i}))$$ for each future time point $$t_i$$ $$T is maturity$$

also we'll output some generated yield curves.

from  QuantLib import *
import numpy as np
from math import *
from pylab import *

def A(t,T):
evaldate=Settings.instance().evaluationDate
forward = crvToday.forwardRate(t, t,Continuous, NoFrequency).rate()
value = B(t,T)*forward - 0.25*sigma*B(t,T)*sigma*B(t,T)*B(0.0,2.0*t);
return exp(value)*crvToday.discount(T)/crvToday.discount(t);

def B(t,T):
return (1.0-exp(-a*(T-t)))/a;

def gamma(t):
forwardRate =crvToday.forwardRate(t, t, Continuous, NoFrequency).rate()
temp = sigma*(1.0 - exp(-a*t))/a
return (forwardRate + 0.5*temp*temp)

def gamma_v(t): #vectorized version of gamma(t)
res=np.zeros(len(t))
for i in xrange(len(t)):
res[i]=gamma(t[i])
return res

Nsim=1000
#a=0.03
#sigma=0.00743
#parameters calibrated with Quantlib to coterminal swaptions on 26/dec/2013
a=0.376739
sigma=0.0209835

todaysDate=Date(26,12,2013);
Settings.instance().evaluationDate=todaysDate;
crvTodaydates=[Date(26,12,2013),
Date(30,6,2014),
Date(30,7,2014),
Date(29,8,2014),
Date(30,9,2014),
Date(30,10,2014),
Date(28,11,2014),
Date(30,12,2014),
Date(30,1,2015),
Date(27,2,2015),
Date(30,3,2015),
Date(30,4,2015),
Date(29,5,2015),
Date(30,6,2015),
Date(30,12,2015),
Date(30,12,2016),
Date(29,12,2017),
Date(31,12,2018),
Date(30,12,2019),
Date(30,12,2020),
Date(30,12,2021),
Date(30,12,2022),
Date(29,12,2023),
Date(30,12,2024),
Date(30,12,2025),
Date(29,12,2028),
Date(30,12,2033),
Date(30,12,2038),
Date(30,12,2043),
Date(30,12,2048),
Date(30,12,2053),
Date(30,12,2058),
Date(31,12,2063)]
crvTodaydf=[1.0,
0.998022,
0.99771,
0.99739,
0.997017,
0.996671,
0.996337,
0.995921,
0.995522,
0.995157,
0.994706,
0.994248,
0.993805,
0.993285,
0.989614,
0.978541,
0.961973,
0.940868,
0.916831,
0.890805,
0.863413,
0.834987,
0.807111,
0.778332,
0.750525,
0.674707,
0.575192,
0.501258,
0.44131,
0.384733,
0.340425,
0.294694,
0.260792
]

crvToday=DiscountCurve(crvTodaydates,crvTodaydf,Actual360(),TARGET())
#crvToday=FlatForward(todaysDate,0.0121,Actual360())

r0=forwardRate =crvToday.forwardRate(0,0, Continuous, NoFrequency).rate()
months=range(3,12*5+1,3)
sPeriods=[str(month)+"m" for month in months]
print sPeriods
Dates=[todaysDate]+[todaysDate+Period(s) for s in sPeriods]
T=+[Actual360().yearFraction(todaysDate,Dates[i]) for i in xrange(1,len(Dates))]
T=np.array(T)
rmean=r0*np.exp(-a*T)+ gamma_v(T) -gamma(0)*np.exp(-a*T)
rvar=sigma*sigma/2.0/a*(1.0-np.exp(-2.0*a*T))
rstd=np.sqrt(rvar)
np.random.seed(1)
stdnorm = np.random.standard_normal(size=(Nsim,len(T)-1))

rmat=np.zeros(shape=(Nsim,len(T)))
rmat[:,0]=r0
for iSim in xrange(Nsim):
for iT in xrange(1,len(T)):
mean=rmat[iSim,iT-1]*exp(-a*(T[iT]-T[iT-1]))+gamma(T[iT])-gamma(T[iT-1])*exp(-a*(T[iT]-T[iT-1]))
var=0.5*sigma*sigma/a*(1-exp(-2*a*(T[iT]-T[iT-1])))
rnew=mean+stdnorm[iSim,iT-1]*sqrt(var)
#if (rnew<0): rnew=0.001
rmat[iSim,iT]=rnew

bonds=np.zeros(shape=rmat.shape)

#E(E(exp(rt)|ti) test
for iSim in xrange(Nsim):
for iT in xrange(1,len(T)):
bonds[iSim,iT]=bonds[iSim,iT-1]+rmat[iSim,iT]*(T[iT]-T[iT-1])

bonds=-bonds;
bonds=np.exp(bonds)

bondsmean=np.mean(bonds,axis=0)
#plot(T,bondsmean)
#plot(T,[crvToday.discount(T[iT]) for iT in xrange(len(T))])
#show()

startDate=Date(26,12,2013);

crvMat= [ [ 0 for i in xrange(len(T)) ] for iSim in range(Nsim) ]
npvMat= [ [ 0 for i in xrange(len(T)) ] for iSim in range(Nsim) ]

for row in crvMat:
row=crvToday

for iT in xrange(1,len(T)):
for iSim in xrange(Nsim):
crvDate=Dates[iT];
crvDates=[crvDate]+[crvDate+Period(k,Years) for k in xrange(1,21)]
rt=rmat[iSim,iT]
#if (rt<0): rt=0
crvDiscounts=[1.0]+[A(T[iT],T[iT]+k)*exp(-B(T[iT],T[iT]+k)*rt) for k in xrange(1,21)]
crvMat[iSim][iT]=DiscountCurve(crvDates,crvDiscounts,Actual360(),TARGET())

bondT=np.zeros(shape=rmat.shape)
for iSim in xrange(Nsim):
for iT in xrange(len(T)):
bondT[iSim,iT]=bonds[iSim,iT]*crvMat[iSim][iT].discount(19-T[iT])

bondTmean=np.mean(bondT,axis=0)
np.set_printoptions(precision=4,suppress=True)
print 'bondTmean-Terminal bond\n',bondTmean-crvToday.discount(19)

plot(range(10),[crvMat.forwardRate(k, k,Continuous, NoFrequency).rate() for k in range(10)])
for i in xrange(min(Nsim,10)):
plot(range(10),[crvMat[i].forwardRate(k, k,Continuous, NoFrequency).rate() for k in range(10)])
title('generated yield curves')
show()

#indexes definitions
index = Euribor(Period("6m"),forecastTermStructure)

#swap 1 definition
maturity = Date(26,12,2018);
fixedSchedule = Schedule(startDate, maturity,Period("6m"), TARGET(),ModifiedFollowing,ModifiedFollowing,DateGeneration.Forward, False)
floatingSchedule = Schedule(startDate, maturity,Period("6m"),TARGET() ,ModifiedFollowing,ModifiedFollowing,DateGeneration.Forward, False)
swap1 = VanillaSwap(VanillaSwap.Receiver, 10000000,fixedSchedule,0.02, Actual360(),floatingSchedule, index, 0,Actual360())  #0.01215

for iT in xrange(len(T)):
Settings.instance().evaluationDate=Dates[iT]
allDates= list(floatingSchedule)
fixingdates=[index.fixingDate(floatingSchedule[iDate]) for iDate in xrange(len(allDates)) if index.fixingDate(floatingSchedule[iDate])<=Dates[iT]]
if fixingdates:
for date in fixingdates[:-1]:
except:pass
except:pass
swapEngine = DiscountingSwapEngine(discountTermStructure)
swap1.setPricingEngine(swapEngine)

for iSim in xrange(Nsim):
crv=crvMat[iSim][iT]
npvMat[iSim][iT]=swap1.NPV()

npvMat=np.array(npvMat)
npv=npvMat[0,0]
#replace negative values with 0
npvMat[npvMat<0]=0
EE=np.mean(npvMat,axis=0)
print '\nEE:\n',EE
#print '\nrmat:\n',rmat
print '\nrmean:\n',rmean
#print '\nrstd:\n',rstd
#print '\n95% are in \n',zip(rmean-2*rstd,rmean+2*rstd)

S=0.05
R=0.4
sum=0
for i in xrange(len(T)-1):
sum=sum+0.5*crvToday.discount(T[i+1])*(EE[i]+EE[i+1])*(exp(-S*T[i]/(1.0-R))-exp(-S*T[i+1]/(1.0-R)))
CVA=(1.0-R)*sum
print "\nnpv=",npv
print "\nCVA=",CVA

plot(T,EE)
title('Expected Exposure')
xlabel('Time in years')
#plot(T,np.mean(rmat,axis=0))
#plot(T,rmean)
#plot(T,[npvMat[0,0]]*len(T))
show()
print "\nnpv",npvMat


output: Expected exposure

bloomberg CVA = 35266

NB: bloomberg's Exposure is our exposure multiplied by (1-R) CVA swap

## Black Scholes calculator online

price of Call option and Option’s delta and gamma

introduce inputs in yellow cells:
dividend yield (all cash dividends divided by spot)

excel :

Advanced Black-Scholes online calculator takes real market interest rate curve and can calculate American options

## formula Black Scholes

Black Scholes formula is

$$Call=df(F_0N(d_+)-KN(d_-))$$

where $$F_0$$ is forward of the stock S at time 0 for maturity T

$$F_0=S_0e^{-qT}/df_T$$

where:
q – dividend yield
df – discount factor for time T

$$d_+=\frac{ln(F_0/K)}{\sigma\sqrt{T}}+0.5\sigma\sqrt{T}$$
$$d_-=\frac{ln(F_0/K)}{\sigma\sqrt{T}}-0.5\sigma\sqrt{T}$$

Posted in OTC derivatives valuation Tagged with: , ,

## how to value Credit Default Swap default leg and default probabilities

how to value CDS (credit default swap) default leg with following time structure:

0—-t1—–t2—-t3—–t4—–….—T

Suppose that default (at time $$\tau$$) can only occur at discrete times t1,t2,t3,..
and Qi=survival probability until time $$t_i$$

then

$$\tau – time of default {=t_1,=t_2,=t_3 …}$$

$$Q_i=P(\tau>t_i) – survival proba for counterparty$$

$$P(\tau=t_3)=1-P(\tau\neq t_3)=1-( P(\tau<3 | \tau>=4) )$$

$$= 1 – ( 1-Q_2 + Q_3)= Q_2-Q_3$$

default leg of Credit Default Swap pays:

$$discounted payoff = (1-R) ( df_1 E{\bf1}(\tau=t_1) + df_2 E{\bf1}(\tau=t_2)+ …)$$
$$=(1-R) ( df_1 P(\tau=t_1)+df_2 P(\tau=t_2)+ …)$$
$$=(1-R) ( df_1(1-Q_1)+df_2(Q_1-Q_2)+df_3(Q_2-Q_3) … )$$

where $$df_1$$ is discount factor for $$t_1$$ and R – recovery rate (normally assumed ot be 40%)

to use continuous time one normally use intensity:
$$P(\tau>t+dt | \tau>t) = \lambda(t)dt$$
where we ssupose that default is the first jump of Poisson process
$$Q(t)=exp( -\int_{0}^{T}\lambda(s)ds)$$
in practice $$\lambda$$ is assumed to be piecewise-constant chaning value on market CDS maturities

useful approximations:
Credit Default Swap spread $$\approx p(1-R) \approx \lambda (1-R)$$
where $$\lambda$$ is intensity and approx. equal to probability of default in first year.

Posted in math, OTC derivatives valuation Tagged with:

## Part 5

c++ finance course

### Objective

based on this you’ll be able to price Autocallable, Himalaya, Spread and similar basket options

Here let’s see how to price an option on 2 stocks with payoff
$$( w_1S_1+w_2S_2-K)^+$$
where $$w_i$$ are basket weights

let’s consider zero interest rate environment and n non-dividend paying stocks with dynamics: (for interest rates and dividends just add $$(r-q)dt$$ terms to this dynamics

$$dS_1=S_1 \sigma_1 dW_1(t)$$
$$dS_2=S_2 \sigma_2 dW_2(t)$$

$$dW_i dW_j=\rho_{i,j} dt$$

where $$\rho_{i,j}$$ is a stock correlation matrix

given the correlation matrix we’ll construct it’s Cholesky decomposition.
Cholesky decomposition can be done with boost but it’s easier to do it with QuantLib implementation.
to store correlation and cholesky matrices we’ll use QuantLib::Matrix class.
we’ll store vector of independent gaussians in the QuantLib::Array class.
we don’t store in std::vector class because to be able to multiply matrix by vector they need to come from the same library.

to obtain correlated gaussians we’ll generate vector of independent gaussians $$\vec{\epsilon}= (\epsilon_1, \epsilon_2)$$ first.
then vector of correlated gaussians is $$M_{cholesky} \vec{\epsilon}$$

here’s the code for wights =1

#include <boost/random/variate_generator.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <iostream>
#include <ql/quantlib.hpp>

double pp(double x) //positive part   x+
{
if(x>0) return x;
return 0;
}

void main()
{

boost::mt19937 mt;
boost::normal_distribution<> distribution;
boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal(new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution));

int numberSimulations=10000;

std::vector<double> S0(2);
S0=100.0;// S1(0) - spot price of first stock
S0=100.0;//S2(0) - spot price of second stock

std::vector<double> sigma(2);
sigma=0.2; //volatility of S1
sigma=0.2;//volatility of S2

double T=1.0;//maturity in years
double K=100.0;//strike of option
double sum=0;//for monte-carlo averaging
double payoff=0;
double NT=100;//number of time intervals
double dt=T/NT;

QuantLib::Matrix rho(2,2);//correlation matrix

rho=1;rho=0.5;
rho=0.5;rho=1;

QuantLib::Matrix cholesky=QuantLib::CholeskyDecomposition(rho);

for(int iSim=0;iSim<numberSimulations;iSim++) //outer loop for simulations
{

std::vector<double> Sprev(S0); //create a vector equal to S0 with its values
std::vector<double> Snext(2);//create vector with 2 elemens

for(double t=0;t<T;t+=dt) //inner loop for constructing one stock price from time 0 to T
{
QuantLib::Array epsilon(2);//create a vector with 2 elements .we'll put here independent gaussians
epsilon=(*normal)();
epsilon=(*normal)();

QuantLib::Array x(2);
x=cholesky*epsilon; //correlated gaussians

for(int iStock=0;iStock<2;iStock++) //loop over all stocks
{
Snext[iStock]=Sprev[iStock]*std::exp(-0.5*sigma[iStock]*sigma[iStock]*dt+sigma[iStock]*std::sqrt(dt)*x[iStock]);
Sprev[iStock]=Snext[iStock];
}
}

payoff=pp(Snext+Snext-K);

sum+=payoff;
}

int dummy;std::cin>>dummy;
}


#### Exercises

– modify the program to calculate asian basket payoff
– modify the program to N stocks
– modify program to include interest rates and dividend yield

Posted in quant finance coding course

# Part 4

c++ finance course

# Monte-Carlo c++ – Path-dependent payoffs

### Objective

– price asian option with c++
– price lookback option

here we’ll show how to extend previous programs to price path dependent payoffs

#### Asian Option

lets consider following asian payoff:

$$payoff=(\frac{1}{T}\int_0^T S_t dt – K)^+$$

1) instead of calculating the integral we’ll calculate average price of Stock S
2) to compute the average we’ll subdivide interval of time [0..T] in N parts each dt apart

0–dt–2td–3dt–…. T

we’ll have to simulate every path time step by time step
to go from previous time step to next one by black scholes dynamics (interest rates=0 for simplicity)

$$S_0=100$$
$$S_{dt}=S_0 exp( -0.5 \sigma^2 dt + \sigma \sqrt{dt} x )$$
$$S_{next}=S_{prev} exp( -0.5 \sigma^2 dt + \sigma \sqrt{dt} x )$$

here comes the code for asian option:

#include <boost/random/variate_generator.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <iostream>

double pp(double x) //positive part   x+
{
if(x>0) return x;
return 0;
}

void main()
{

boost::mt19937 mt;
boost::normal_distribution<> distribution;
boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal(new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution));

int numberSimulations=1000;

double S0=100.0;
double sigma=0.2;
double T=1.0;
double K=100.0;
double sum=0;
double payoff=0;
double NT=100;
double dt=T/NT;

for(int iSim=0;iSim<numberSimulations;iSim++) //outer loop for simulations
{
double Sprev=S0;
double Snext=0;
double averageSum=0;

for(double t=0;t<=T;t+=dt) //inner loop for constructing one stock price path from time 0 to T
{
double x=(*normal)();
Snext=Sprev*std::exp(-0.5*sigma*sigma*dt+sigma*std::sqrt(dt)*x);
averageSum+=Snext;
Sprev=Snext;
//std::cout<<"\n"<<Snext; output to console of stock dynamics
}

payoff=pp(averageSum/NT-K);
sum+=payoff;
}

std::cout<<"\nCall="<<sum/numberSimulations<<std::endl;

int dummy;std::cin>>dummy;
}



in this example code we use a common construct for path-dependent payoffs :2 monte carlo loops
one for generating one path (inner loop) and another, outer loop for Monte carlo averaging of payoff

#### LookBack Option

in lookback option payoff strike is minimum of the stock price path over the period
so let’s change previous program to calculate lookback option

$$payoff=(S_T-\min_{0 #include <boost/random/variate_generator.hpp> #include <boost/random/mersenne_twister.hpp> #include <boost/random/normal_distribution.hpp> #include <iostream> double pp(double x) //positive part x+ { if(x>0) return x; return 0; } void main() { boost::mt19937 mt; boost::normal_distribution<> distribution; boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal(new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution)); int numberSimulations=1000; double S0=100.0; double sigma=0.2; double T=1.0; double K=100.0; double sum=0; double payoff=0; double NT=100; double dt=T/NT; for(int iSim=0;iSim<numberSimulations;iSim++) //outer loop for simulations { double Sprev=S0; double Snext=0; double currentMin=S0; //set current minimum stock level to initial price for(double t=0;t<=T;t+=dt) //inner loop for constructing one stock price from time 0 to T { double x=(*normal)(); Snext=Sprev*std::exp(-0.5*sigma*sigma*dt+sigma*std::sqrt(dt)*x); if(Snext<currentMin) currentMin=Snext; Sprev=Snext; //std::cout<<"\n"<<Snext<<" MIN="<<currentMin;// output to console of stock dynamics } payoff=pp(Snext-currentMin); sum+=payoff; } std::cout<<"\nLookBack Call="<<sum/numberSimulations<<std::endl; int dummy;std::cin>>dummy; } #### Exercises – introduce non-zero interest rates – calculate lookback option where strike is maximum stock price over the period Posted in quant finance coding course ## finance c++ programming course [part 3] – Monte Carlo c++ # Part 3 c++ finance course # Monte-Carlo c++ ### Objective – generate random numbers – Calculate Call option with Monte – Carlo – intro to std::vector here we’ll see how to use containers in C++ in many books on c++ you’ll see that people use c++ arrays ( for example double x) for storing variables.It has quite many drawbacks for beginners , so here we’ll use std::vector to store real numbers. std::vector is an object – it means you add dot “.” after the object name and you will see it has not only one value but many different “sub-variables” and “sub-functions” ### Gaussian random number generation To implement Monte-Carlo engine in C++ we’ll need standard gaussian numbers. C++ has a built-in random number generators , but they are no good for financial engineering (function random()). One of the most used pseudo – random generators with good characteristics is Mersenne Twister algorithm It’s implemented in Boost. After we have uniform random numbers we can generate Gaussian number with a function.It’s also implemented in Boost. here’s the code to generate 100 standard gaussians: #include <boost/random/variate_generator.hpp> #include <boost/random/mersenne_twister.hpp> #include <boost/random/normal_distribution.hpp> #include <iostream> void main() { boost::mt19937 mt; //create object of mersenne twister generator boost::normal_distribution<> distribution; // create distribution object boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal(new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution)); //create variate generator which will generate gaussians for (int i=0;i<100;i++) std::cout<< (*normal)()<<std::endl; // gaussian is retrieved with function (*normal)() asterisk * is because normal variable is a pointer int dummy;std::cin>>dummy; }  here we use a pointer to variable. pointer is basically address of variable in memory. to access variables by pointers one uses asterisk * read more about pointers here ### Storage of data in std::vector to store these numbers we can use class std::vector read more about this class here here we’ll illustrate its common use: #include <boost/random/variate_generator.hpp> #include <boost/random/mersenne_twister.hpp> #include <boost/random/normal_distribution.hpp> #include <iostream> #include <vector> void main() { boost::mt19937 mt; boost::normal_distribution<> distribution; boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal (new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution)); std::vector<double> v; for (int i=0;i<100;i++) v.push_back( (*normal)() ) ; for (int i=0;i<100;i++) std::cout<<v[i]<<std::endl; int dummy;std::cin>>dummy; }  now lets make a function to calculate the mean of elements of vector #include <boost/random/variate_generator.hpp> #include <boost/random/mersenne_twister.hpp> #include <boost/random/normal_distribution.hpp> #include <iostream> #include <vector> double mean(std::vector<double>& v) //here we use & to indicate that we pass value of variable to function by address and not by creating new variable and copying data into it { int N=v.size(); double sum=0; for(int i=0;i<N;i++) sum+=v[i]; return sum/N; } void main() { boost::mt19937 mt; boost::normal_distribution<> distribution; boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal(new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution)); std::vector<double> v; for (int i=0;i<100;i++) v.push_back( (*normal)() ) ; for (int i=0;i<100;i++) std::cout<<v[i]<<std::endl; std::cout<<"\nmean="<<mean(v); int dummy;std::cin>>dummy; }  here we’ve shown how to to pass vectors as function parameters most useful vector expressions are: std::vector<double> v;//create empty vector std::vector<double> v(10);//create vector of 10 elements first element is v v=4;//set value of first vector element; v.resize(100);//resize vector to have 100 elements; v.size();//returns number of elements  if you want to make simple operations on vectors there are chances that these operations are already implemented in STL or in boost. for example our previous program to compute mean of array could be rewritten using STL algorithms as following: sum =std::accumulate(v.begin(),v.end(),0.0); //need to #include <numeric>  so before writing a function to do simple things on vector google first somthing among the lines “std::vector sum of elements” #### Exercise – change previous program to calculate variance of random numbers as we have random numbers now we are ready to implement call option price ### Call option with monte carlo c++ suppose Black-Scholes dynamics with 0 interest rates, so$$S_T=S_0 exp(-0.5 \sigma^2 T + \sigma \sqrt{T} x )$$where x is standard gaussian Monte-Carlo Call price is just average of payoffs$$ (S_T(x)-K)^+$$#include <boost/random/variate_generator.hpp> #include <boost/random/mersenne_twister.hpp> #include <boost/random/normal_distribution.hpp> #include <iostream> double pp(double x) //positive part x+ { if(x>0) return x; return 0; } void main() { boost::mt19937 mt; boost::normal_distribution<> distribution; boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>* normal(new boost::variate_generator<boost::mt19937&,boost::normal_distribution<>>(mt,distribution)); int numberSimulations=100000; double S0=100.0; double sigma=0.2; double T=1.0; double K=100.0; double sum=0; double payoff=0; for(int iSim=0;iSim<numberSimulations;iSim++) { double x=(*normal)(); payoff=pp(S0*std::exp(-0.5*sigma*sigma*T+sigma*std::sqrt(T)*x)-K); sum+=payoff; } std::cout<<"\nCall="<<sum/numberSimulations<<std::endl; int dummy;std::cin>>dummy; }  #### Exercises – change program to include non-zero constant interest rate – change program to calculate Put option price – verify call put parity Posted in quant finance coding course ## finance c++ programming course [part 2] – numerical integration Black Scholes ## part 2 ## black scholes via numerical integration finance c++ course ## objective – calculate Call option in Black Scholes model using c++ C++ and QuantLib have many types of variables in quant finance most useful types are: double – to store real number int – to store integer number QuantLib::Date – to store dates QuantLib::Array – to store array of double QuantLib::Matrix – to store 2d matrix of doubles lets start by doing some numerical integration delete all the contents of source.cpp file and copy paste following code ## Example – numerical integration this code will calculate integral$$ \int_{0}^{1}e^xdx $$#include <iostream> #include <cmath> //for using math functions //everything after // is a comment double f(double x) //declaration of function f(x)=exp(x) { return std::exp(x); } void main() { double dt=0.0001; //integration step double sum=0; //current value of integral for(double t=0;t<1.0;t+=dt) //for loop t+=dt means t=t+dt { sum+=f(t)*dt; // its equivalent to write sum=sum+f(t)*dt; } std::cout<<"integral "<< sum; //output int dummy;std::cin>>dummy;//screen delay }  read about for loop and other control structures here: http://www.cplusplus.com/doc/tutorial/control/ important notes: use 1.0 instead of just 1 in c++ 1.0 and 1 are different numbers 1.0 is real number and 1 is integer number it means that 1/2=0 while 1.0/2 is 0.5 when doing comparisons of doubles don’t use if(x==0) use if( std::abs(x)for loop see “variables scope” here http://www.cplusplus.com/doc/tutorial/variables/ here we used simplest form of numerical integration for a faster/ more precise integratoin you can use other methods many of them are described in free book “numerical recipes in c” this book is for C but the code structure would be the same ## example – call option in Black Scholes model using integration now lets change previous program to calculate Call option in Black Scholes model with zero interest rates risk-neutral dynamics:$$ dS=S*\sigma dWt$$so$$ S_T=S_0 exp(-0.5 \sigma^2+\sigma \sqrt T x )$$where x is a standard gaussian so the call option price will be$$ Call=\int_{-\infty}^{+\infty} (S_0 e^{-0.5 \sigma^2+\sigma \sqrt T x} – K)^+ \phi(x) dx $$where$$\phi(x)=\frac{1}{\sqrt{2\pi}}exp(-0.5 x^2)$$standard gaussian density S0=100 volatility=20% maturity=T=1 year #include <iostream> #include <cmath> #include <boost/math/constants/constants.hpp> double density(double x) //normal density { const double pi = boost::math::constants::pi<double>(); return 1.0/sqrt(2.0*pi)*std::exp(-0.5*x*x); } double pp(double x) //positive part x+ { if(x>0) return x; return 0; } double f(double K,double S0,double T,double sigma,double x) //function to integrate { return pp(S0*std::exp(-0.5*sigma*sigma*T+sigma*std::sqrt(T)*x)-K)*density(x); } void main() { double S0=100.0; double sigma=0.2; double T=1.0; double K=80.0; double dx=0.0001; double sum=0; for(double x=-10.0;x<10.0;x+=dx) { sum+=f(K,S0,T,sigma,x)*dx; } std::cout<<"call= "<< sum; int dummy;std::cin>>dummy; }  Notes: we have to pass arguments again to the integrated function because their original scope is inside main() , so variables K, S0 etc are not visible inside function f(x) we use boost here to extract Pi constant , so it will be easier to port this program to linux for example (standard c++ has no math constants) ### exercises – debug first program by setting break points and watch variables sum and t change value – change first program to integrate via trapezoidal rule , check the results against known answer$$ e^1-1 
– change second program to include non zero interest rates
– change second program to calculate put price
– verify call put parity

### some tricks

if you’ve got an error always start by correcting first compile error (double click on compile error line)
some c++ error messages are quite cryptic and are results of previous errors

for example remove “;” after return statement
press F7 to compile
now double-click on first error:

Posted in quant finance coding course

## C++ finance course for beginners – part 1 – hello world

REQUIREMENTS:
– microsoft visual studio express 2012 (it’s free,older versions are similar. there are slight differences between microsoft and linux c++)
– some programming experience in another language

### Course objectives:

-price various derivatives using c++ code
-use basic quantlib and boost functions useful in financial engineering

### C++ ,Boost and Quantlib install

#### 2) install boost

– unpack it into drive C:
[so its unpacked into c:\boost_1_53_0 ]

– go to the dir c:\boost_1_53_0 and run bootstrap.bat
– run b2.exe
it will take few hours to build!

– go to dir c:\boost_1_53_0\stage
move whole dir “lib” to c:\boost_1_53_0\ so all built libraries are inside c:\boost_1_53_0\lib

#### 3) install quantlib

http://sourceforge.net/projects/quantlib/files/

unpack it to C:\

– open visual studio express for desktop

File->Open project -> QuantLib_vc10.sln

Choose “update” in the next dialog . it will convert files to visual studio 2012 format

in the solution explorer window select QuantLib

change “include directories” variable to include path to boost
change “library directories” variable to include path to boost\lib

now delete all projects except QuantLib from the solution (click on every project in solution explorer and press “del”)

press F7 to build quantlib [or menu->build->build solution]

## Hello world program

choose Visual C++ -> Win32 Console Application

press OK->NEXT (do NOT press FINISH yet)

check “Empty Project”
uncheck “Security…” box

click Finish

right Click on ConsoleApplication1 anc change include and library directory to include QuantLib and boost directories:

in the window copy paste following code:

#include <ql/quantlib.hpp>
#include <iostream>

using namespace QuantLib;

void main()
{

Date date1(30,January,2013);
Date date2;
date2=date1+1*Years;

std::cout<<date2;

int dummy;std::cin>>dummy;
}


press F7 to compile

double click on error :

and change the following line :

#  error "unknown Microsoft compiler"


by

#  define QL_LIB_TOOLSET "vc100"


so it looks like that:

press F7 to build project and F5 to run

quick explanation of the code:

#include <ql/quantlib.hpp>
#include <iostream>


these lines just insert files named quantlib.hpp and iostream into your source file
quantlib.hpp is needed to have quantlib functions and definitions (Date object)
and iostream is used to write to output console

using namespace QuantLib;


without this line you’d had to write “QuantLib::” before every quantlib object

void main()
{

Date date1(30,January,2013);
Date date2;
date2=date1+1*Years;

std::cout<<date2;

int dummy;std::cin>>dummy;
}


these lines define main function
in c++ when console program starts it automatically start executing the function main()

	Date date1(30,January,2013);
Date date2;


here we define two QuantLib objects(variables) date1 and date2 and assign a value for the first variable

date2=date1+1*Years;
std::cout<<date2;


here we assign value for date2 to be 1 year from date1

and std::cout part prints the result to console

int dummy;std::cin>>dummy;


this is needed to pause the screen

### debugging C++ program

to debug program left click where the red dot appears and press F5

the program will start debug and you can hover cursor over the variable to see it's value

as you can see date1 will show some long number instead of date .It's because internally quantlib uses Excel format for dates .
press F10 to execute instructions one by one of F5 again to continue program until next break point

Posted in quant finance coding course

## FX forward valuation excel

Here we’ll show how to value FX forwards (Foreign eXchange) on EURUSD exchange rate in excel with quantlib excel addin

in this example we’ll need following market data:

• EURUSD spot
• EURUSD forward points (can get from bloomberg or reuters)
• EUR discount curve (for example EUR 6m curve)
• In this spreadsheet we first construct syntetic USD yield curve based on EURUSD forward points, so after we could use usual formulas for FX forward valuation

FX forward valuation excel:

Posted in OTC derivatives valuation