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.
#!/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")