APPENDIX B:
The following pseudo-code (C++) demonstrates how to implement the model
to price a callable bond. For the purpose of an easy illustration, we
choose the same settings (the number of nodes and the grid space) across
the lattice and use the Trapezoidal Rule for numerical integration.
// 2*numNodes = 2*mNumNodes = the number of nodes (S); gap = mGap = the
grid space (Phi)
double priceCallableBond (BondTrade* bd, CallableBond* cb, int numNodes,
double gap) {
double pv;
cb->fillLattice();
// The last exercise
CallSchedule& cs = bd->callSch[numCallSch-1];
if (cs.term == bd->cFlow[numCashFlow-1].endDate) // The
last exercise is at maturity
for (int i= -numNodes; i <= numNodes; i++)
cs.reducedValue[i+numNodes] = min (cs.callPrice,
bd->cFlow[numCashFlow-1].reducedPayoff[i+numNodes]);
else { // The last exercise is before maturity
for (int i= -numNodes; i <= numNodes; i++) {
pv = 0;
for (int j = bd->numCF-1;
(bd->cFlow[j].endDate >= cs.term) && (j
>= 0); j–) {
CashFlow& cf = bd->cFlow[j];
(cf.endDate == cs.term) ? pv += cf.reducedPayoff[i+numNodes]
: pv += exp(-bondSpread*(cf.endDate-cs.term)) *
cb->integral(i,
cs.vol, cf.vol, cf.endDate, cs.term, cf.reducedPayoff);
}
cs.reducedValue[i+numNodes] = min
(cs.callPrice/cs.df[i+numNodes], pv);
}
}
if (numCallSch > 1) { // The remaining exercises
for (int i = numCallSch - 2; i>=0; i–) {
CallSchedule& cs = bd->callSch[i];
CallSchedule& preCs = bd->callSch[i+1];
for (int j = -numNodes; j <= numNodes; j++) {
pv = exp(-bondSpread * (preCs.term - cs.term))
* cb->integral (j, cs.vol, preCs.vol, preCs.term, cs.term,
preCs.reducedValue);
for (int k=bd->numCF-1; k >= 0; k–) //
Count intermediate coupons
if ((bd->cFlow[k].endDate < preCs.term) &&
(bd->cFlow[k].endDate >= cs.term))
pv += bd->cFlow[k].reducedPayoff[j+numNodes]
* exp (-bondSpread*(bd->cFlow[k].endDate - cs.term));
cs.reducedValue[j+numNodes] = min
(cs.callPrice/cs.df[j+numNodes], pv);
}
}
}
// The final integral
CallSchedule& preCs = bd->callSch[0];
pv = cb->integral (0, 0, preCs.vol, preCs.term, 0,
preCs.reducedValue) *exp(-bondSpread*(preCs.term));
pv *= bd->cFlow[bd->numCF-1].endDf; //
endDf: discount factor from 0 to the end date
for (int k=bd->numCF-1; k >= 0; k–) //
Count intermediate coupons
if ((bd->cFlow[k].endDate < preCs.term))
pv += bd->cFlow[k].coupon *
bd->cFlow[k].endDf * exp(-bondSpread *
bd->cFlow[k].endDate);
return pv;
}
void CallableBond::fillLattice() {
for (int i = mTrade->numCF-1; i>=0; i–) {
CashFlow& cf = mTrade->cFlow[i];
if (cf.endDate < mTrade->callSch[0].term)
break;
for (int j = -mNumNodes; j <= mNumNodes; j++)
fillNode(i, j, cf.startDate, mDrift);
}
}
void CallableBond::fillNode(int cI, int nI, double vT, DriftAppx flag)
{
int numCF = mTrade->numCF;
double avgF, expon, fwdt, drift = 0;
CashFlow& fl = mTrade->cFlow[cI];
if (cI == numCF-1) { // At maturity
fl.df[nI + mNumNodes] = 1.0;
fl.reducedPayoff[nI + mNumNodes] = fl.notional + fl.coupon;
}
else if (fl.startDate <= 0) // Starting before valuation date)
fl.reducedPayoff[nI + mNumNodes] = fl.coupon * fl.endDf /
mTrade->cFlow[numCF-1].endDf;
else {
fl.df[nI + mNumNodes] = 1.0;
for (int i = numCF - 1; i > cI; i–) {
CashFlow& cf = mTrade->cFlow[i];
expon = (cf.vol * cf.vol * vT / 2) + cf.vol * nI * mGap;
fwdt = cf.fwd0 * exp(-drift + expon);
switch (flag) { // The other cases are similar to either AAFR or CEFR
case AAFR: // Arithemic Average Fwd Rate
avgF = 0.5 * (cf.fwd0 + fwdt);
drift += vT * fl.vol * cf.vol * cf.delta * avgF / (1 + cf.delta * avgF);
break;
case CEFR: // Conditional Expectation of Fwd Rate
drift += fl.vol * cf.vol * integralFwd(cf.fwd0, fwdt, 0, vT, cf.vol,
cf.delta);
break;
default:
break;
}
fl.df[nI + mNumNodes] /= (1 + fwdt * cf.delta); // df: discount
factor maturing at maturity
}
fl.reducedPayoff[nI + mNumNodes] = fl.coupon / fl.df[nI +
mNumNodes];
}
}
// Gauss-Legendre integration for drift
const double xArray[] = {0, 0.1488743389, 0.4333953941,
0.6794095682, 0.8650633666, 0.9739065285};
const double wArray[] = {0, 0.2955242247, 0.2692667193,
0.2190863625, 0.1494513491, 0.0666713443};
double CallableBond::integralFwd(double F0, double Ft, double a, double
b, double vol, double delta) {
double xm = 0.5 * (b + a);
double xr = 0.5 * (b - a);
double ss = 0, dx = 0;
for (int j = 1; j <= 5; j++) {
dx = xr * xArray[j];
ss += wArray[j] * (expectFwd(F0, Ft, (xm + dx), b, vol, delta)
+ expectFwd(F0, Ft, (xm - dx), b, vol, delta));
}
return ss * xr;
}
double CallableBond::expectFwd(double F0, double Ft, double s, double t,
double vol, double delta) {
double mean = F0 * pow ((Ft / F0), (s / t)) * exp(0.5 * vol * vol * s *
(t - s) / t);
return delta * mean / (1 + delta * mean);-
}
// Trapezoidal Rule Integration
double CallableBond::integral (int curPos, double curVol, double preVol,
double preTerm,
double curTerm, double* value){
double diffPos, tmpV, sum = 0;
for (int k = -mNumNodes; k <= mNumNodes; k++) {
diffPos = k*mGap - curPos*mGap + preVol * preTerm - curVol * curTerm;
tmpV = value[k+mNumNodes] * exp (-diffPos * diffPos/(2*(preTerm -
curTerm)));
((k == -mNumNodes) || (k == mNumNodes)) ? sum += 0.5 *
tmpV : sum += tmpV;
}
return sum * mGap / sqrt(2 * PI * (preTerm - curTerm));
}