Serving the Quantitative Finance Community

 
User avatar
AboodiDaQuant
Topic Author
Posts: 0
Joined: September 11th, 2013, 3:40 pm

Brigo/Mercurio Theorem 4.2.3

October 20th, 2014, 11:16 pm

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
 
User avatar
daveangel
Posts: 5
Joined: October 20th, 2003, 4:05 pm

Brigo/Mercurio Theorem 4.2.3

October 21st, 2014, 10:12 am

iteratively ?
knowledge comes, wisdom lingers
 
User avatar
AboodiDaQuant
Topic Author
Posts: 0
Joined: September 11th, 2013, 3:40 pm

Brigo/Mercurio Theorem 4.2.3

October 21st, 2014, 12:54 pm

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
 
User avatar
Gamal
Posts: 1283
Joined: February 26th, 2004, 8:41 am

Brigo/Mercurio Theorem 4.2.3

October 21st, 2014, 1:15 pm

QuoteOriginally posted by: AboodiDaQuantAnd I cannot really resolve the above sum analytically. You wanted to have a closed formula? Aha.
 
User avatar
AboodiDaQuant
Topic Author
Posts: 0
Joined: September 11th, 2013, 3:40 pm

Brigo/Mercurio Theorem 4.2.3

October 21st, 2014, 1:30 pm

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....
 
User avatar
Gamal
Posts: 1283
Joined: February 26th, 2004, 8:41 am

Brigo/Mercurio Theorem 4.2.3

October 21st, 2014, 1:39 pm

Use the first optimization procedure in your library.
 
User avatar
ChristianR
Posts: 1
Joined: November 13th, 2015, 11:29 am

Re: Brigo/Mercurio Theorem 4.2.3

October 11th, 2016, 1:26 pm

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 
 
92 def RSS(K,B,Vm,sx,sy,a,b,rho): 
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 
104     rss = diff 
105     print "a=%f,b=%f,sx=%f,sy=%f,rho=%f" % (a,b,sx,sy,rho) 
106     print diff 
107     return rss 
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): 
115                                   RSS(K,B,Vm,sx,sy,a,b,rho),z,method='SLSQP', 
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 
 
119 def loadData(file): 
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: 
128         reader = csv.reader(f,dialect="mine") 
129         reader.next() 
130         for row in reader: 
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     ''' 
142     B, Vm = loadData(file) 
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: 38
Joined: August 20th, 2019, 12:40 pm

Re: Brigo/Mercurio Theorem 4.2.3

February 5th, 2021, 8:53 am

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

cheers
K
 
User avatar
katastrofa
Posts: 7431
Joined: August 16th, 2007, 5:36 am
Location: Alpha Centauri

Re: Brigo/Mercurio Theorem 4.2.3

February 6th, 2021, 1:12 am

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.