Serving the Quantitative Finance Community

AboodiDaQuant
Topic Author
Posts: 5
Joined: September 11th, 2013, 3:40 pm

### Brigo/Mercurio Theorem 4.2.3

Hi all,I've a question about the theorem 4.2.3 in Brigo/Mercurio. They have a semi-analytic expression for the valuation of european swaptions within the one-factor model. It is a very lengthy expression, however, among others there is a paramter y(x) (x is the integration variable in the swaption price expression), which shall be given as the (unique) solution of: $\sum c_i A(T,t_i)\exp(-B(a,T,t_i)x-B(b,T,t_i)y)=1$ I really have no clue, how to resolve this, since the coefficients A and B are also dependant on the fitting parameters of the underlying model (drift and volatility of the one-factor model). How did they fit this to given ATM-swaption volatilities?! Any hint would be highly appreciated.Thank you very much indeed in advance.Kindest regardsAboodi

daveangel
Posts: 17031
Joined: October 20th, 2003, 4:05 pm

### Brigo/Mercurio Theorem 4.2.3

iteratively ?
knowledge comes, wisdom lingers

AboodiDaQuant
Topic Author
Posts: 5
Joined: September 11th, 2013, 3:40 pm

### Brigo/Mercurio Theorem 4.2.3

Thank you so much for your answer. Could you please a little more specific. I really don't know, since the expressions B and A are also depending on my (not yet known) parameters of my 1-factor model, how to start And I cannot really resolve the above sum analytically. Sorry if I'm really a bit slow to understand in how to tackle it.greetingsaboodi

Gamal
Posts: 2362
Joined: February 26th, 2004, 8:41 am

### Brigo/Mercurio Theorem 4.2.3

QuoteOriginally posted by: AboodiDaQuantAnd I cannot really resolve the above sum analytically. You wanted to have a closed formula? Aha.

AboodiDaQuant
Topic Author
Posts: 5
Joined: September 11th, 2013, 3:40 pm

### Brigo/Mercurio Theorem 4.2.3

No, I did not want to have a closed formula, sorry if I am not clear with my problem. The y to be determined will depend on x and the unknown parameters.... and I don't know how to approach this optimization problem. I would somehow need an expression for y (obviously not analytic) and use it during the fitting procedure. Problem: How can I teach, e.g. matlab, how to handle this problem....

Gamal
Posts: 2362
Joined: February 26th, 2004, 8:41 am

### Brigo/Mercurio Theorem 4.2.3

Use the first optimization procedure in your library.

ChristianR
Posts: 1
Joined: November 13th, 2015, 11:29 am

### Re: Brigo/Mercurio Theorem 4.2.3

Hi folks,
I tried to implement the same calibration procedure and apparently we run into similar difficulties. I will attach my code in python; maybe someone capable of the language can tell me why the code goes bust during the root-search algorithm to determine the y(x). The market data are taken from Brigo/Mercurios book (volatilites and yield curve). I found the actual yield curve here in the forum (alternatively one could check in Bloomberg).

Any help is highly appreciated as this problem bothers me for weeks....

Thanx a lot.

Cheerio
#!/usr/bin/python
2 import scipy.optimize as sp
3 import numpy as np
4 import math
5 import argparse
6 import csv
7 import pprint as pp
8 from utility import Phi, Volatility as V
9

10 def params(n,a,b,sx,sy,rho):
11     n+=1
12     mu_x = -(sx**2/a**2+rho*sx*sy/(a*b))*(1-np.exp(-a*n))+sx**2/(2*a**2)\
13         *(1-np.exp(-2*a*n))+rho*sx*sy/(b*(a+b))*(1-np.exp(-(a+b)*n))
14     mu_y = -(sy**2/b**2+rho*sx*sy/(a*b))*(1-np.exp(-b*n))+sy**2/(2*b**2)\
15         *(1-np.exp(-2*b*n))+rho*sx*sy/(a*(a+b))*(1-np.exp(-(a+b)*n))
16     eta_x = sx*np.sqrt((1-np.exp(-2*a*n))/(2*a))
17     eta_y = sy*np.sqrt((1-np.exp(-2*b*n))/(2*b))
18     eta_xy = rho*sx*sy/((a+b)*eta_x*eta_y)*(1-np.exp(-(a+b)*n))
19     return mu_x, mu_y, eta_x, eta_y, eta_xy
20

21 def h(y,u,n,p,a,b,sx,sy,rho):
22     _,mu_y,_,eta_y,eta_xy = params(n,a,b,sx,sy,rho)
23     h1=(y-mu_y)/(eta_y*np.sqrt(1-eta_xy**2))-np.sqrt(2)*eta_xy*u/np.sqrt(1-eta_xy**2)
24     #print "y=%f,u=%f,n=%f,p=%f,a=%f,b=%f,sx=%f,sy=%f,rho=%f" % (y,u,n,p,a,b,sx,sy,rho)
25     h2=h1+BT(n,p,b)*eta_y*np.sqrt(1-eta_xy**2)
26     return (h1,h2)
27

28 def A(n,p,a,b,sx,sy,rho,B):
29     return B[p]/B[n]*np.exp(1/2*(V(n,p,1,sx,sy,a,b,rho) \
30                                      -V(0,p,1,sx,sy,a,b,rho)
31                                      +V(0,n,1,sx,sy,a,b,rho)))
32

33 def lambdakappa(u,n,p,a,b,sx,sy,rho,K,B):
34     mu_x,_,eta_x,eta_y,eta_xy = params(n,a,b,sx,sy,rho)
35     _,sfwd = SFWD(B)
36     sfwd = sfwd[n][p]
37     if p==n: sfwd+=1
38     lambda1=sfwd*A(n,p,a,b,sx,sy,rho,B)*np.exp(-BT(n,p,a)*(eta_x*np.sqrt(2)*u+mu_x))
39     kappa=BT(n,p,b)*(mu_x-1/2*(1-eta_xy**2)*eta_y**2 \
40             *BT(n,p,b)+eta_xy*eta_y*np.sqrt(2)*u)
41     return lambda1*np.exp(kappa)
42

43 def BT(n,p,a):
44     if not a==0:
45         res = (1-np.exp(a*(n-p)))/a
46     else:
47         res = n-p
48     return res
49

50 def SFWD(B):
51     sumbd = [sum(B[:i+1]) for i in range(len(B))]
52     sfwd = [[(b1-b2)/(sb2-sb1) if i<j else 0 for ((j,b1),sb1) in \
53             zip(enumerate(B),sumbd)] for ((i,b2),sb2) in zip(enumerate(B),sumbd)]
54     return sumbd, sfwd
55

56 def theoreticalSwaption(n,p,a,b,sx,sy,rho,K,B):
57     '''
58     Theoretical price for a swaption with strike K, n the number of years till
59     expiry and p the swap length and discount factors B
60     '''
61     xH, wH = np.polynomial.hermite.hermgauss(10)
62     mu_x,_,eta_x,eta_y,eta_xy = params(n,a,b,sx,sy,rho)
63     _,sfwd = SFWD(B)
64     sfwd = sfwd[n][p]
65     c = [sfwd for i in range(p-n-1)]
66     c.append(sfwd+1)
67     y = []
68     for x in xH:
69         v = [ci*A(n,i+n+1,a,b,sx,sy,rho,B)*np.exp(-BT(n,i+n+1,a)*(np.sqrt(2)*eta_x*x+mu_x)) \
70             for (i,ci) in enumerate(c)]
71         #print " + ".join(["%f*Exp[%f*y]" % (vi,Bi) for (vi,Bi) in zip(v,[-BT(n,i+n+1,b) for (i,_) in enumerate(v)])]) + " == 1"
72         sol = sp.newton(lambda yy: sum((vi*np.exp(-BT(n,i+n+1,b)*yy)) \
73                                     for (i,vi) in enumerate(v))-1.0,0.0, \
74                         fprime=lambda yy: sum((-BT(n,i+n+1,b)*vi*np.exp(-BT(n,i+n+1,b)*yy)) \
75                                                     for (i,vi) in enumerate(v)),maxiter=5000)
76         y.append(sol)
77     #print "Solutions=" + str(y)
78     g1 = [h(yi,xi,n,p,a,b,sx,sy,rho) for (xi,yi) in zip(xH,y)]
79     g1 = map(lambda u: (Phi(-u[0]),Phi(-u[1])),g1)
80     g2 = sum(map(lambda u: lambdakappa(u,n,p,a,b,sx,sy,rho,K,B),xH))
81     integrand = [w*(k1-k2*g2) for (w,(k1,k2)) in zip(wH,g1)]
82     integral = B[n]/np.sqrt(math.pi)*sum(integrand)
83     return integral
84

85 def swaptionFromMarketData(B,Vm):
86     sumbd,sfwd = SFWD(B)
87     d=[[2*Phi(0.5*v*np.sqrt(i+1))+1 for v in row] for (i,row) in enumerate(Vm)]
88     swap = [[s*d1*(sumbd[j]-sumbd[i]) for ((j,s),d1) in zip(enumerate(row1),row2)]
89             for (i,(row1,row2)) in enumerate(zip(sfwd,d))]
90     return swap
91

93     '''
94     Residual sum of squares
95     '''
96     v1 = swaptionFromMarketData(B,Vm)
97     N = len(v1)
98     v2 = [[theoreticalSwaption(n,p,a,b,sx,sy,rho,K,B) if n<p else 0 for \
99             p in range(N)] for n in range(N)]
100     diff = sum(sum(abs(vv1-vv2) if not vv2==0 else 0 for (vv1,vv2) in zip(row1,row2))
101             for (row1,row2) in zip(v1,v2))
102     #print v1
103     #print v2
105     print "a=%f,b=%f,sx=%f,sy=%f,rho=%f" % (a,b,sx,sy,rho)
106     print diff
108

109 def calibrate(K,B,Vm):
110     '''
111     Calibrate parameters
112     '''
113     z = [0.136, -0.02, 0.033, 0.00587, -0.835]
114     sol = sp.minimize(lambda (a,b,sx,sy,rho):
116                                   bounds=((-2.0,2.0),(-2.0,2.0),(0.001,1.0),(0.001,1.0),(-1.0,1.0)))
117     return sol.x
118

120     '''
121     Read market data from a csv file with separator ';' and header line
122     'B;V1;...;Vn', which has exactly as many columns as rows
123     '''
124     csv.register_dialect("mine",delimiter=";")
125     B = []
126     Vm = []
127     with open(file,'rb') as f:
131             B.append(float(row[0]))
132             vrow = []
133             for v in row[1:]:
134                 vrow.append(float(v))
135             Vm.append(vrow)
136     return B, Vm
137

138 def DoIt(file,K=0.01):
139     '''
140     Read 'file' and fit sx,sy,a and b
141     '''
143     res = calibrate(K,B,Vm)
144     result = zip(["a = ","b = ","sx = ","sy = ","rho = "],map(str,res))
145     result = ', '.join(map(lambda (x,y): x+y, result))
146     print result
147

148 if __name__=="__main__":
149     DoIt("marketdata_swaptions.csv")



krs
Posts: 20
Joined: August 20th, 2019, 12:40 pm

### Re: Brigo/Mercurio Theorem 4.2.3

I encountered the same problem, use scipy.optimize.fsolve, it should work...

cheers
K

katastrofa
Posts: 10069
Joined: August 16th, 2007, 5:36 am
Location: Alpha Centauri

### Re: Brigo/Mercurio Theorem 4.2.3

ChristianR, did you try to plot the function to see what you're dealing with? Maybe Newton-Raphson is not a good choice? Try Brent instead or the above.