Navigation

← Previous Changeset
Next Changeset →

Changeset 220

Timestamp:

06/28/08 13:41:06 (1 year ago)

Author:

snovik

Message:

Fix: ConundrumPricer? (please have a look)

Files:

trunk/QLNet/QLNet/Cashflows/ConundrumPricer.cs (modified) (7 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

trunk/QLNet/QLNet/Cashflows/ConundrumPricer.cs

r216	r220
1	1	/*
2	2	Copyright (C) 2008 Toyin Akin (toyin_akin@hotmail.com)
	3	Copyright (C) 2008 Siarhei Novik (snovik@gmail.com)
3	4
4	5	This file is part of QLNet Project http://www.qlnet.org
…	…
23	24
24	25	namespace QLNet {
25
26		public abstract class VanillaOptionPricer
27		{
28		public abstract double value(double strike, Option.Type optionType, double deflator);
	26	public abstract class VanillaOptionPricer {
	27	public abstract double value(double strike, Option.Type optionType, double deflator);
	28	}
	29
	30	//===========================================================================//
	31	// BlackVanillaOptionPricer //
	32	//===========================================================================//
	33	public class BlackVanillaOptionPricer : VanillaOptionPricer {
	34	private double forwardValue_;
	35	private Date expiryDate_;
	36	private Period swapTenor_;
	37	private SwaptionVolatilityStructure volatilityStructure_;
	38	private SmileSection smile_;
	39
	40	public BlackVanillaOptionPricer(double forwardValue, Date expiryDate, Period swapTenor, SwaptionVolatilityStructure volatilityStructure) {
	41	forwardValue_ = forwardValue;
	42	expiryDate_ = expiryDate;
	43	swapTenor_ = swapTenor;
	44	volatilityStructure_ = volatilityStructure;
	45	smile_ = volatilityStructure_.smileSection(expiryDate_, swapTenor_);
	46	}
	47
	48	public override double value(double strike, Option.Type optionType, double deflator) {
	49	double variance = smile_.variance(strike);
	50	return deflator * Utils.blackFormula(optionType, strike, forwardValue_, Math.Sqrt(variance));
	51	}
29	52	}
30	53
31		public class BlackVanillaOptionPricer : VanillaOptionPricer
32		{
33
34		//===========================================================================//
35		// BlackVanillaOptionPricer //
36		//===========================================================================//
37
38		public BlackVanillaOptionPricer(double forwardValue, Date expiryDate, Period swapTenor, SwaptionVolatilityStructure volatilityStructure)
39		{
40		forwardValue_ = forwardValue;
41		expiryDate_ = expiryDate;
42		swapTenor_ = swapTenor;
43		volatilityStructure_ = volatilityStructure;
44		smile_ = volatilityStructure_.smileSection(expiryDate_, swapTenor_);
	54	public abstract class GFunction {
	55	public abstract double value(double x);
	56	public abstract double firstDerivative(double x);
	57	public abstract double secondDerivative(double x);
	58	}
	59
	60	public class GFunctionFactory {
	61	public enum YieldCurveModel : int {
	62	Standard,
	63	ExactYield,
	64	ParallelShifts,
	65	NonParallelShifts
	66	}
	67	public static GFunction newGFunctionStandard(int q, double delta, int swapLength) {
	68	return new GFunctionStandard(q, delta, swapLength) as GFunction;
	69	}
	70	public static GFunction newGFunctionExactYield(CmsCoupon coupon) {
	71	return new GFunctionExactYield(coupon) as GFunction;
	72	}
	73	public static GFunction newGFunctionWithShifts(CmsCoupon coupon, Handle<Quote> meanReversion) {
	74	return new GFunctionWithShifts(coupon, meanReversion) as GFunction;
	75	}
	76
	77	//===========================================================================//
	78	// GFunctionStandard //
	79	//===========================================================================//
	80	private class GFunctionStandard : GFunction {
	81	// number of period per year
	82	protected int q_;
	83	// fraction of a period between the swap start date and the pay date
	84	protected double delta_;
	85	// length of swap
	86	protected int swapLength_;
	87
	88	public GFunctionStandard(int q, double delta, int swapLength) {
	89	q_ = q;
	90	delta_ = delta;
	91	swapLength_ = swapLength;
	92	}
	93
	94	public override double value(double x) {
	95	double n = swapLength_ * q_;
	96	return x / Math.Pow((1.0 + x / q_), delta_) * 1.0 / (1.0 - 1.0 / Math.Pow((1.0 + x / q_), n));
	97	}
	98
	99	public override double firstDerivative(double x) {
	100	double n = swapLength_ * q_;
	101	double a = 1.0 + x / q_;
	102	double AA = a - delta_ / q_ * x;
	103	double B = Math.Pow(a, (n - delta_ - 1.0)) / (Math.Pow(a, n) - 1.0);
	104
	105	double secNum = n * x * Math.Pow(a, (n - 1.0));
	106	double secDen = q_ * Math.Pow(a, delta_) * (Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0);
	107	double sec = secNum / secDen;
	108
	109	return AA * B - sec;
	110	}
	111
	112	public override double secondDerivative(double x) {
	113	double n = swapLength_ * q_;
	114	double a = 1.0 + x / q_;
	115	double AA = a - delta_ / q_ * x;
	116	double A1 = (1.0 - delta_) / q_;
	117	double B = Math.Pow(a, (n - delta_ - 1.0)) / (Math.Pow(a, n) - 1.0);
	118	double Num = (1.0 + delta_ - n) * Math.Pow(a, (n - delta_ - 2.0)) - (1.0 + delta_) * Math.Pow(a, (2.0 * n - delta_ - 2.0));
	119	double Den = (Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0);
	120	double B1 = 1.0 / q_ * Num / Den;
	121
	122	double C = x / Math.Pow(a, delta_);
	123	double C1 = (Math.Pow(a, delta_) - delta_ / q_ * x * Math.Pow(a, (delta_ - 1.0))) / Math.Pow(a, 2 * delta_);
	124
	125	double D = Math.Pow(a, (n - 1.0)) / ((Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0));
	126	double D1 = ((n - 1.0) * Math.Pow(a, (n - 2.0)) * (Math.Pow(a, n) - 1.0) - 2 * n * Math.Pow(a, (2 * (n - 1.0)))) / (q_ * (Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0));
	127
	128	return A1 * B + AA * B1 - n / q_ * (C1 * D + C * D1);
	129	}
45	130	}
46	131
47		public override double value(double strike, Option.Type optionType, double deflator)
48		{
49		double variance = smile_.variance(strike);
50		return deflator * Utils.blackFormula(optionType, strike, forwardValue_, Math.Sqrt(variance));
51		}
52		private double forwardValue_;
53		private Date expiryDate_;
54		private Period swapTenor_;
55		private SwaptionVolatilityStructure volatilityStructure_;
56		private SmileSection smile_;
57		}
58
59		public abstract class GFunction
60		{
61		public abstract double value(double x);
62		public abstract double firstDerivative(double x);
63		public abstract double secondDerivative(double x);
64		}
65
66		public class GFunctionFactory
67		{
68		public enum YieldCurveModel: int
69		{
70		Standard,
71		ExactYield,
72		ParallelShifts,
73		NonParallelShifts
74		}
75		public static GFunction newGFunctionStandard(int q, double delta, int swapLength)
76		{
77		return new GFunctionStandard(q, delta, swapLength) as GFunction;
78		}
79		public static GFunction newGFunctionExactYield(CmsCoupon coupon)
80		{
81		return new GFunctionExactYield(coupon) as GFunction;
82		}
83		public static GFunction newGFunctionWithShifts(CmsCoupon coupon, Handle<Quote> meanReversion)
84		{
85		return new GFunctionWithShifts(coupon, meanReversion) as GFunction;
86		}
87
88		private class GFunctionStandard : GFunction
89		{
90		public GFunctionStandard(int q, double delta, int swapLength)
91		{
92		q_ = q;
93		delta_ = delta;
94		swapLength_ = swapLength;
95		}
96
97		//===========================================================================//
98		// GFunctionStandard //
99		//===========================================================================//
100
101		public override double value(double x)
102		{
103		double n = swapLength_ * q_;
104		return x / Math.Pow((1.0 + x/q_), delta_) * 1.0 / (1.0 - 1.0 / Math.Pow((1.0 + x/q_), n));
105		}
106		public override double firstDerivative(double x)
107		{
108		double n = swapLength_ * q_;
109		double a = 1.0 + x / q_;
110		double AA = a - delta_/q_ * x;
111		double B = Math.Pow(a,(n - delta_ - 1.0))/(Math.Pow(a,n) - 1.0);
112
113		double secNum = n * x * Math.Pow(a,(n-1.0));
114		double secDen = q_ * Math.Pow(a, delta_) * (Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0);
115		double sec = secNum / secDen;
116
117		return AA * B - sec;
118		}
119		public override double secondDerivative(double x)
120		{
121		double n = swapLength_ * q_;
122		double a = 1.0 + x/q_;
123		double AA = a - delta_/q_ * x;
124		double A1 = (1.0 - delta_)/q_;
125		double B = Math.Pow(a,(n - delta_ - 1.0))/(Math.Pow(a,n) - 1.0);
126		double Num = (1.0 + delta_ - n) * Math.Pow(a, (n-delta_-2.0)) - (1.0 + delta_) * Math.Pow(a, (2.0 *n-delta_-2.0));
127		double Den = (Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0);
128		double B1 = 1.0 / q_ * Num / Den;
129
130		double C = x / Math.Pow(a, delta_);
131		double C1 = (Math.Pow(a, delta_) - delta_ /q_ * x * Math.Pow(a, (delta_ - 1.0))) / Math.Pow(a, 2 * delta_);
132
133		double D = Math.Pow(a, (n-1.0))/ ((Math.Pow(a, n) - 1.0) * (Math.Pow(a, n) - 1.0));
134		double D1 = ((n - 1.0) * Math.Pow(a, (n-2.0)) * (Math.Pow(a, n) - 1.0) - 2 * n * Math.Pow(a, (2 * (n-1.0)))) / (q_ * (Math.Pow(a, n) - 1.0)(Math.Pow(a, n) - 1.0)(Math.Pow(a, n) - 1.0));
135
136		return A1 * B + AA * B1 - n/q_ * (C1 * D + C * D1);
137		}
138		// number of period per year
139		protected int q_;
140		// fraction of a period between the swap start date and
141		// the pay date
142		protected double delta_;
143		// length of swap
144		protected int swapLength_;
145		}
146
147		private class GFunctionExactYield : GFunction
148		{
149
150		//===========================================================================//
151		// GFunctionExactYield //
152		//===========================================================================//
153
154		public GFunctionExactYield(CmsCoupon coupon)
155		{
156
157		SwapIndex swapIndex = coupon.swapIndex();
158		VanillaSwap swap = swapIndex.underlyingSwap(coupon.fixingDate());
159
160		Schedule schedule = swap.fixedSchedule();
161		Handle<YieldTermStructure> rateCurve = swapIndex.termStructure();
162
163		DayCounter dc = swapIndex.dayCounter();
164
165		double swapStartTime = dc.yearFraction(rateCurve.link.referenceDate(), schedule.startDate());
166		double swapFirstPaymentTime = dc.yearFraction(rateCurve.link.referenceDate(), schedule.date(1));
167
168		double paymentTime = dc.yearFraction(rateCurve.link.referenceDate(), coupon.date());
169
170		delta_ = (paymentTime-swapStartTime) / (swapFirstPaymentTime-swapStartTime);
	132	//===========================================================================//
	133	// GFunctionExactYield //
	134	//===========================================================================//
	135	private class GFunctionExactYield : GFunction {
	136	// fraction of a period between the swap start date and the pay date
	137	protected double delta_;
	138	// accruals fraction
	139	protected List<double> accruals_;
	140
	141	public GFunctionExactYield(CmsCoupon coupon) {
	142
	143	SwapIndex swapIndex = coupon.swapIndex();
	144	VanillaSwap swap = swapIndex.underlyingSwap(coupon.fixingDate());
	145
	146	Schedule schedule = swap.fixedSchedule();
	147	Handle<YieldTermStructure> rateCurve = swapIndex.termStructure();
	148
	149	DayCounter dc = swapIndex.dayCounter();
	150
	151	double swapStartTime = dc.yearFraction(rateCurve.link.referenceDate(), schedule.startDate());
	152	double swapFirstPaymentTime = dc.yearFraction(rateCurve.link.referenceDate(), schedule.date(1));
	153
	154	double paymentTime = dc.yearFraction(rateCurve.link.referenceDate(), coupon.date());
	155
	156	delta_ = (paymentTime - swapStartTime) / (swapFirstPaymentTime - swapStartTime);
171	157
172	158	List<CashFlow> fixedLeg = new List<CashFlow>(swap.fixedLeg());
173		int n = fixedLeg.Count;
	159	int n = fixedLeg.Count;
174	160	accruals_ = new List<double>();
175		for (int i =0; i<n; ++i)
176		{
177		Coupon coupon1 = fixedLeg[i] as Coupon;
178		accruals_.Add(coupon1.accrualPeriod());
179		}
180		}
181
182		public override double value(double x)
183		{
184		double product = 1.0;
185		for(int i =0; i<accruals_.Count; i++)
186		{
187		product = 1.0/(1.0+ accruals_[i]x);
188		}
189		return x Math.Pow(1.0+ accruals_[0]x,-delta_)*(1.0/(1.0-product));
190		}
191		public override double firstDerivative(double x)
192		{
193		double c = -1.0;
194		double derC = 0.0;
195		List<double> b = new List<double>();
196		for (int i =0; i<accruals_.Count; i++)
197		{
198		double temp = 1.0/(1.0+ accruals_[i]*x);
199		b.Add(temp);
200		c *= temp;
201		derC += accruals_[i]*temp;
202		}
203		c += 1.0;
204		c = 1.0/c;
205		derC = (c-c c);
206
207		return -delta_ accruals_[0]Math.Pow(b[0],delta_+1.0)x c+ Math.Pow(b[0],delta_)c+ Math.Pow(b[0],delta_)x *derC;
208		//double dx = 1.0e-8;
209		//return (operator()(x+dx)-operator()(x-dx))/(2.0*dx);
210		}
211		public override double secondDerivative(double x)
212		{
213		double c = -1.0;
214		double sum = 0.0;
215		double sumOfSquare = 0.0;
216		List<double> b = new List<double>();
217		for(int i =0; i<accruals_.Count; i++)
218		{
219		double temp = 1.0/(1.0+ accruals_[i]*x);
220		b.Add(temp);
221		c *= temp;
222		sum += accruals_[i]*temp;
223		sumOfSquare += Math.Pow(accruals_[i]*temp, 2.0);
224		}
225		c += 1.0;
226		c = 1.0/c;
227		double derC =sum(c-c c);
228
229		return (-delta_ accruals_[0]Math.Pow(b[0],delta_+1.0)c+ Math.Pow(b[0],delta_)derC)* (-delta_ accruals_[0]b[0]x + 1.0 + x(1.0-c)sum)+ Math.Pow(b[0],delta_)c(delta_ Math.Pow(accruals_[0]b[0],2.0)x - delta_* accruals_[0]b[0] - x derC sum + (1.0-c)sum - x(1.0-c)sumOfSquare);
230		//double dx = 1.0e-8;
231		//return (firstDerivative(x+dx)-firstDerivative(x-dx))/(2.0*dx);
232		}
233		// fraction of a period between the swap start date and
234		// the pay date
235		protected double delta_;
236		// accruals fraction
237		protected List<double> accruals_;
238		}
239
240		private class GFunctionWithShifts : GFunction
241		{
242
243		private double swapStartTime_;
244
245		private double shapedPaymentTime_;
246		private List<double> shapedSwapPaymentTimes_;
247
248		private List<double> accruals_;
249		private List<double> swapPaymentDiscounts_;
250		private double discountAtStart_;
251		private double discountRatio_;
252
253		private double swapRateValue_;
254		private Handle<Quote> meanReversion_;
255
256		private double calibratedShift_;
257		private double tmpRs_;
258		private double accuracy_;
259
260		//* function describing the non-parallel shape of the curve shift*/
261		private double shapeOfShift(double s)
262		{
263		double x = s-swapStartTime_;
264		double meanReversion = meanReversion_.link.value();
265		if(meanReversion>0)
266		{
267		return (1.0-Math.Exp(-meanReversion *x))/meanReversion;
268		}
269		else
270		{
271		return x;
272		}
273		}
274		//* calibration of shift*/
275		private double calibrationOfShift(double Rs)
276		{
277
278		if(Rs!=tmpRs_)
279		{
280		double initialGuess;
281		double N =0;
282		double D =0;
283		for(int i =0; i<accruals_.Count; i++)
284		{
285		N+=accruals_[i]*swapPaymentDiscounts_[i];
286		D+=accruals_[i]swapPaymentDiscounts_[i]shapedSwapPaymentTimes_[i];
287		}
288		N *= Rs;
289		D *= Rs;
290		N += accruals_.Last() * swapPaymentDiscounts_.Last() - objectiveFunction_.gFunctionWithShifts().discountAtStart_;
291		D += accruals_.Last() * swapPaymentDiscounts_.Last()* shapedSwapPaymentTimes_.Last();
292		initialGuess = N/D;
293
294		objectiveFunction_.setSwapRateValue(Rs);
295		Newton solver = new Newton();
296		solver.setMaxEvaluations(1000);
297
298		// these boundaries migth not be big enough if the volatility
299		// of big swap rate values is too high . In this case the G function
300		// is not even integrable, so better to fix the vol than increasing
301		// these values
302		double lower = -20;
303		double upper = 20.0;
304
305		try
306		{
307		calibratedShift_ = solver.solve(objectiveFunction_, accuracy_, Math.Max(Math.Min(initialGuess, upper.99), lower.99), lower, upper);
308		}
309		catch (Exception e)
310		{
	161	for (int i = 0; i < n; ++i) {
	162	Coupon coupon1 = fixedLeg[i] as Coupon;
	163	accruals_.Add(coupon1.accrualPeriod());
	164	}
	165	}
	166
	167	public override double value(double x) {
	168	double product = 1.0;
	169	for (int i = 0; i < accruals_.Count; i++) {
	170	product = 1.0 / (1.0 + accruals_[i] x);
	171	}
	172	return x * Math.Pow(1.0 + accruals_[0] * x, -delta_) * (1.0 / (1.0 - product));
	173	}
	174
	175	public override double firstDerivative(double x) {
	176	double c = -1.0;
	177	double derC = 0.0;
	178	List<double> b = new List<double>();
	179	for (int i = 0; i < accruals_.Count; i++) {
	180	double temp = 1.0 / (1.0 + accruals_[i] * x);
	181	b.Add(temp);
	182	c *= temp;
	183	derC += accruals_[i] * temp;
	184	}
	185	c += 1.0;
	186	c = 1.0 / c;
	187	derC = (c - c c);
	188
	189	return -delta_ * accruals_[0] * Math.Pow(b[0], delta_ + 1.0) * x * c + Math.Pow(b[0], delta_) * c + Math.Pow(b[0], delta_) * x * derC;
	190	//double dx = 1.0e-8;
	191	//return (operator()(x+dx)-operator()(x-dx))/(2.0*dx);
	192	}
	193
	194	public override double secondDerivative(double x) {
	195	double c = -1.0;
	196	double sum = 0.0;
	197	double sumOfSquare = 0.0;
	198	List<double> b = new List<double>();
	199	for (int i = 0; i < accruals_.Count; i++) {
	200	double temp = 1.0 / (1.0 + accruals_[i] * x);
	201	b.Add(temp);
	202	c *= temp;
	203	sum += accruals_[i] * temp;
	204	sumOfSquare += Math.Pow(accruals_[i] * temp, 2.0);
	205	}
	206	c += 1.0;
	207	c = 1.0 / c;
	208	double derC = sum * (c - c * c);
	209
	210	return (-delta_ * accruals_[0] * Math.Pow(b[0], delta_ + 1.0) * c + Math.Pow(b[0], delta_) * derC) * (-delta_ * accruals_[0] * b[0] * x + 1.0 + x * (1.0 - c) * sum) + Math.Pow(b[0], delta_) * c * (delta_ * Math.Pow(accruals_[0] * b[0], 2.0) * x - delta_ * accruals_[0] * b[0] - x * derC * sum + (1.0 - c) * sum - x * (1.0 - c) * sumOfSquare);
	211	//double dx = 1.0e-8;
	212	//return (firstDerivative(x+dx)-firstDerivative(x-dx))/(2.0*dx);
	213	}
	214	}
	215
	216	private class GFunctionWithShifts : GFunction {
	217	private double swapStartTime_;
	218
	219	private double shapedPaymentTime_;
	220	private List<double> shapedSwapPaymentTimes_;
	221
	222	private List<double> accruals_;
	223	private List<double> swapPaymentDiscounts_;
	224	private double discountAtStart_;
	225	private double discountRatio_;
	226
	227	private double swapRateValue_;
	228	private Handle<Quote> meanReversion_;
	229
	230	private double calibratedShift_;
	231	private double tmpRs_;
	232	private double accuracy_;
	233
	234	private ObjectiveFunction objectiveFunction_;
	235
	236	//* function describing the non-parallel shape of the curve shift*/
	237	private double shapeOfShift(double s) {
	238	double x = s - swapStartTime_;
	239	double meanReversion = meanReversion_.link.value();
	240	if (meanReversion > 0) {
	241	return (1.0 - Math.Exp(-meanReversion * x)) / meanReversion;
	242	} else {
	243	return x;
	244	}
	245	}
	246	//* calibration of shift*/
	247	private double calibrationOfShift(double Rs) {
	248	if (Rs != tmpRs_) {
	249	double initialGuess;
	250	double N = 0;
	251	double D = 0;
	252	for (int i = 0; i < accruals_.Count; i++) {
	253	N += accruals_[i] * swapPaymentDiscounts_[i];
	254	D += accruals_[i] * swapPaymentDiscounts_[i] * shapedSwapPaymentTimes_[i];
	255	}
	256	N *= Rs;
	257	D *= Rs;
	258	N += accruals_.Last() * swapPaymentDiscounts_.Last() - objectiveFunction_.gFunctionWithShifts().discountAtStart_;
	259	D += accruals_.Last() * swapPaymentDiscounts_.Last() * shapedSwapPaymentTimes_.Last();
	260	initialGuess = N / D;
	261
	262	objectiveFunction_.setSwapRateValue(Rs);
	263	Newton solver = new Newton();
	264	solver.setMaxEvaluations(1000);
	265
	266	// these boundaries migth not be big enough if the volatility
	267	// of big swap rate values is too high . In this case the G function
	268	// is not even integrable, so better to fix the vol than increasing
	269	// these values
	270	double lower = -20;
	271	double upper = 20.0;
	272
	273	try {
	274	calibratedShift_ = solver.solve(objectiveFunction_, accuracy_, Math.Max(Math.Min(initialGuess, upper * .99), lower * .99), lower, upper);
	275	} catch (Exception e) {
311	276	throw new ApplicationException("meanReversion: " + meanReversion_.link.value() + ", swapRateValue: " + swapRateValue_ + ", swapStartTime: " + swapStartTime_ + ", shapedPaymentTime: " + shapedPaymentTime_ + "\n error message: " + e.Message);
312		}
313		tmpRs_ =Rs;
314		}
315		return calibratedShift_;
316		}
317		private double functionZ(double x)
318		{
319		return Math.Exp(-shapedPaymentTime_ x) / (1.0-discountRatio_ Math.Exp(-shapedSwapPaymentTimes_.Last()*x));
320		}
321		private double derRs_derX(double x)
322		{
323		double sqrtDenominator = 0;
324		double derSqrtDenominator = 0;
325		for(int i =0; i<accruals_.Count; i++)
326		{
327		sqrtDenominator += accruals_[i]swapPaymentDiscounts_[i] Math.Exp(-shapedSwapPaymentTimes_[i]*x);
328		derSqrtDenominator -= shapedSwapPaymentTimes_[i]* accruals_[i]swapPaymentDiscounts_[i] Math.Exp(-shapedSwapPaymentTimes_[i]*x);
329		}
330		double denominator = sqrtDenominator* sqrtDenominator;
331
332		double numerator = 0;
333		numerator += shapedSwapPaymentTimes_.Last()* swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)sqrtDenominator;
334		numerator -= (discountAtStart_ - swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)) derSqrtDenominator;
335		if (!(denominator!=0))
	277	}
	278	tmpRs_ = Rs;
	279	}
	280	return calibratedShift_;
	281	}
	282
	283	private double functionZ(double x) {
	284	return Math.Exp(-shapedPaymentTime_ * x) / (1.0 - discountRatio_ * Math.Exp(-shapedSwapPaymentTimes_.Last() * x));
	285	}
	286
	287	private double derRs_derX(double x) {
	288	double sqrtDenominator = 0;
	289	double derSqrtDenominator = 0;
	290	for (int i = 0; i < accruals_.Count; i++) {
	291	sqrtDenominator += accruals_[i] * swapPaymentDiscounts_[i] * Math.Exp(-shapedSwapPaymentTimes_[i] * x);
	292	derSqrtDenominator -= shapedSwapPaymentTimes_[i] * accruals_[i] * swapPaymentDiscounts_[i] * Math.Exp(-shapedSwapPaymentTimes_[i] * x);
	293	}
	294	double denominator = sqrtDenominator * sqrtDenominator;
	295
	296	double numerator = 0;
	297	numerator += shapedSwapPaymentTimes_.Last() * swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x) * sqrtDenominator;
	298	numerator -= (discountAtStart_ - swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x)) * derSqrtDenominator;
	299	if (!(denominator != 0))
336	300	throw new ApplicationException("GFunctionWithShifts::derRs_derX: denominator == 0");
337		~~return numerator/~~denominator;
338		}
339		private double derZ_derX(double x)
340		{
341		double sqrtDenominator = (1.0-discountRatio_ Math.Exp(-shapedSwapPaymentTimes_.Last()x));
342		~~double denominator = sqrtDenominator~~* sqrtDenominator;
343		~~if (!(denominator!=~~0))
	301	return numerator / denominator;
	302	}
	303
	304	private double derZ_derX(double x) {
	305	double sqrtDenominator = (1.0 - discountRatio_ * Math.Exp(-shapedSwapPaymentTimes_.Last() * x));
	306	double denominator = sqrtDenominator * sqrtDenominator;
	307	if (!(denominator != 0))
344	308	throw new ApplicationException("GFunctionWithShifts::derZ_derX: denominator == 0");
345
346		double numerator = 0;
347		numerator -= shapedPaymentTime_* Math.Exp(-shapedPaymentTime_ x) sqrtDenominator;
348		numerator -= shapedSwapPaymentTimes_.Last()* Math.Exp(-shapedPaymentTime_ x) (1.0-sqrtDenominator);
349
350		return numerator/denominator;
351		}
352		private double der2Rs_derX2(double x)
353		{
354		double denOfRfunztion = 0.0;
355		double derDenOfRfunztion = 0.0;
356		double der2DenOfRfunztion = 0.0;
357		for(int i =0; i<accruals_.Count; i++)
358		{
359		denOfRfunztion += accruals_[i]swapPaymentDiscounts_[i] Math.Exp(-shapedSwapPaymentTimes_[i]*x);
360		derDenOfRfunztion -= shapedSwapPaymentTimes_[i]* accruals_[i]swapPaymentDiscounts_[i] Math.Exp(-shapedSwapPaymentTimes_[i]*x);
361		der2DenOfRfunztion+= shapedSwapPaymentTimes_[i]shapedSwapPaymentTimes_[i] accruals_[i]* swapPaymentDiscounts_[i]Math.Exp(-shapedSwapPaymentTimes_[i]x);
362		}
363
364		double denominator = Math.Pow(denOfRfunztion, 4);
365
366		double numOfDerR = 0;
367		numOfDerR += shapedSwapPaymentTimes_.Last()* swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)denOfRfunztion;
368		numOfDerR -= (discountAtStart_ - swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)) derDenOfRfunztion;
369
370		double denOfDerR = Math.Pow(denOfRfunztion,2);
371
372		double derNumOfDerR = 0.0;
373		derNumOfDerR -= shapedSwapPaymentTimes_.Last()shapedSwapPaymentTimes_.Last() swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)denOfRfunztion;
374		derNumOfDerR += shapedSwapPaymentTimes_.Last()* swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)derDenOfRfunztion;
375
376		derNumOfDerR -= (shapedSwapPaymentTimes_.Last()swapPaymentDiscounts_.Last() Math.Exp(-shapedSwapPaymentTimes_.Last()x)) derDenOfRfunztion;
377		derNumOfDerR -= (discountAtStart_ - swapPaymentDiscounts_.Last()* Math.Exp(-shapedSwapPaymentTimes_.Last()x)) der2DenOfRfunztion;
378
379		double derDenOfDerR = 2 denOfRfunztion derDenOfRfunztion;
380
381		double numerator = derNumOfDerR denOfDerR -numOfDerR derDenOfDerR;
382		if (!(denominator!=0))
	309
	310	double numerator = 0;
	311	numerator -= shapedPaymentTime_ * Math.Exp(-shapedPaymentTime_ * x) * sqrtDenominator;
	312	numerator -= shapedSwapPaymentTimes_.Last() * Math.Exp(-shapedPaymentTime_ * x) * (1.0 - sqrtDenominator);
	313
	314	return numerator / denominator;
	315	}
	316
	317	private double der2Rs_derX2(double x) {
	318	double denOfRfunztion = 0.0;
	319	double derDenOfRfunztion = 0.0;
	320	double der2DenOfRfunztion = 0.0;
	321	for (int i = 0; i < accruals_.Count; i++) {
	322	denOfRfunztion += accruals_[i] * swapPaymentDiscounts_[i] * Math.Exp(-shapedSwapPaymentTimes_[i] * x);
	323	derDenOfRfunztion -= shapedSwapPaymentTimes_[i] * accruals_[i] * swapPaymentDiscounts_[i] * Math.Exp(-shapedSwapPaymentTimes_[i] * x);
	324	der2DenOfRfunztion += shapedSwapPaymentTimes_[i] * shapedSwapPaymentTimes_[i] * accruals_[i] * swapPaymentDiscounts_[i] * Math.Exp(-shapedSwapPaymentTimes_[i] * x);
	325	}
	326
	327	double denominator = Math.Pow(denOfRfunztion, 4);
	328
	329	double numOfDerR = 0;
	330	numOfDerR += shapedSwapPaymentTimes_.Last() * swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x) * denOfRfunztion;
	331	numOfDerR -= (discountAtStart_ - swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x)) * derDenOfRfunztion;
	332
	333	double denOfDerR = Math.Pow(denOfRfunztion, 2);
	334
	335	double derNumOfDerR = 0.0;
	336	derNumOfDerR -= shapedSwapPaymentTimes_.Last() * shapedSwapPaymentTimes_.Last() * swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x) * denOfRfunztion;
	337	derNumOfDerR += shapedSwapPaymentTimes_.Last() * swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x) * derDenOfRfunztion;
	338
	339	derNumOfDerR -= (shapedSwapPaymentTimes_.Last() * swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x)) * derDenOfRfunztion;
	340	derNumOfDerR -= (discountAtStart_ - swapPaymentDiscounts_.Last() * Math.Exp(-shapedSwapPaymentTimes_.Last() * x)) * der2DenOfRfunztion;
	341
	342	double derDenOfDerR = 2 * denOfRfunztion * derDenOfRfunztion;
	343
	344	double numerator = derNumOfDerR * denOfDerR - numOfDerR * derDenOfDerR;
	345	if (!(denominator != 0))
383	346	throw new ApplicationException("GFunctionWithShifts::der2Rs_derX2: denominator == 0");
384		~~return numerator/~~denominator;
385		}
386		private double der2Z_derX2(double x)
387		{
388		double denOfZfunction = (1.0-discountRatio_ Math.Exp(-shapedSwapPaymentTimes_.Last()x));
389		double derDenOfZfunction = shapedSwapPaymentTimes_.Last()discountRatio_ Math.Exp(-shapedSwapPaymentTimes_.Last()*x);
390		double denominator = Math.Pow(denOfZfunction, 4);
391		~~if (!(denominator!=~~0))
	347	return numerator / denominator;
	348	}
	349
	350	private double der2Z_derX2(double x) {
	351	double denOfZfunction = (1.0 - discountRatio_ * Math.Exp(-shapedSwapPaymentTimes_.Last() * x));
	352	double derDenOfZfunction = shapedSwapPaymentTimes_.Last() * discountRatio_ * Math.Exp(-shapedSwapPaymentTimes_.Last() * x);
	353	double denominator = Math.Pow(denOfZfunction, 4);
	354	if (!(denominator != 0))
392	355	throw new ApplicationException("GFunctionWithShifts::der2Z_derX2: denominator == 0");
393
394		double numOfDerZ = 0;
395		numOfDerZ -= shapedPaymentTime_* Math.Exp(-shapedPaymentTime_ x) denOfZfunction;
396		numOfDerZ -= shapedSwapPaymentTimes_.Last()* Math.Exp(-shapedPaymentTime_ x) (1.0-denOfZfunction);
397
398		double denOfDerZ = Math.Pow(denOfZfunction,2);
399		double derNumOfDerZ = (-shapedPaymentTime_* Math.Exp(-shapedPaymentTime_ x) (-shapedPaymentTime_+(shapedPaymentTime_ discountRatio_- shapedSwapPaymentTimes_.Last()discountRatio_)* Math.Exp(-shapedSwapPaymentTimes_.Last()x)) -shapedSwapPaymentTimes_.Last()Math.Exp(-shapedPaymentTime_ x) (shapedPaymentTime_ discountRatio_- shapedSwapPaymentTimes_.Last()discountRatio_)* Math.Exp(-shapedSwapPaymentTimes_.Last()*x));
400
401		double derDenOfDerZ = 2 denOfZfunction derDenOfZfunction;
402		double numerator = derNumOfDerZ denOfDerZ -numOfDerZ derDenOfDerZ;
403
404		return numerator/denominator;
405		}
406
407		private class ObjectiveFunction : ISolver1d
408		{
409		private GFunctionWithShifts o_;
410		private double Rs_;
411		private double derivative_;
412
413		public ObjectiveFunction(GFunctionWithShifts o, double Rs)
414		{
415		o_ = o;
416		Rs_ = Rs;
417		}
418		public override double value(double x)
419		{
420		double result = 0;
421		derivative_ = 0;
422		for(int i =0; i<o_.accruals_.Count; i++)
423		{
424		double temp = o_.accruals_[i]o_.swapPaymentDiscounts_[i] Math.Exp(-o_.shapedSwapPaymentTimes_[i]*x);
425		result += temp;
426		derivative_ -= o_.shapedSwapPaymentTimes_[i] * temp;
427		}
428		result *= Rs_;
429		derivative_ *= Rs_;
430		double temp1 = o_.swapPaymentDiscounts_.Last() * Math.Exp(-o_.shapedSwapPaymentTimes_.Last()*x);
431
432		result += temp1-o_.discountAtStart_;
433		derivative_ -= o_.shapedSwapPaymentTimes_.Last()*temp1;
434		return result;
435		}
436		public override double derivative(double UnnamedParameter1)
437		{
438		return derivative_;
439		}
440		public void setSwapRateValue(double x)
441		{
442		Rs_ = x;
443		}
444		public GFunctionWithShifts gFunctionWithShifts()
445		{
446		return o_;
447		}
448		}
449
450		private ObjectiveFunction objectiveFunction_;
451
452		//===========================================================================//
453		// GFunctionWithShifts //
454		//===========================================================================//
455
456		public GFunctionWithShifts(CmsCoupon coupon, Handle<Quote> meanReversion)
457		{
458		meanReversion_ = meanReversion;
459		calibratedShift_ = 0.03;
460		tmpRs_ = 10000000.0;
461		accuracy_ = 1.0e-14;
462
463		SwapIndex swapIndex = coupon.swapIndex();
464		VanillaSwap swap = swapIndex.underlyingSwap(coupon.fixingDate());
465
466		swapRateValue_ = swap.fairRate();
467
468		objectiveFunction_ = new ObjectiveFunction( this, swapRateValue_);
469
470		Schedule schedule = swap.fixedSchedule();
471		Handle<YieldTermStructure> rateCurve = swapIndex.termStructure();
472		DayCounter dc = swapIndex.dayCounter();
473
474		swapStartTime_ = dc.yearFraction(rateCurve.link.referenceDate(), schedule.startDate());
475		discountAtStart_ = rateCurve.link.discount(schedule.startDate());
476
477		double paymentTime = dc.yearFraction(rateCurve.link.referenceDate(), coupon.date());
478
479		shapedPaymentTime_ = shapeOfShift(paymentTime);
480
481		List<CashFlow> fixedLeg = new List<CashFlow>(swap.fixedLeg());
482		int n = fixedLeg.Count;
	356
	357	double numOfDerZ = 0;
	358	numOfDerZ -= shapedPaymentTime_ * Math.Exp(-shapedPaymentTime_ * x) * denOfZfunction;
	359	numOfDerZ -= shapedSwapPaymentTimes_.Last() * Math.Exp(-shapedPaymentTime_ * x) * (1.0 - denOfZfunction);
	360
	361	double denOfDerZ = Math.Pow(denOfZfunction, 2);
	362	double derNumOfDerZ = (-shapedPaymentTime_ * Math.Exp(-shapedPaymentTime_ * x) * (-shapedPaymentTime_ + (shapedPaymentTime_ * discountRatio_ - shapedSwapPaymentTimes_.Last() * discountRatio_) * Math.Exp(-shapedSwapPaymentTimes_.Last() * x)) - shapedSwapPaymentTimes_.Last() * Math.Exp(-shapedPaymentTime_ * x) * (shapedPaymentTime_ * discountRatio_ - shapedSwapPaymentTimes_.Last() * discountRatio_) * Math.Exp(-shapedSwapPaymentTimes_.Last() * x));
	363
	364	double derDenOfDerZ = 2 * denOfZfunction * derDenOfZfunction;
	365	double numerator = derNumOfDerZ * denOfDerZ - numOfDerZ * derDenOfDerZ;
	366
	367	return numerator / denominator;
	368	}
	369
	370	private class ObjectiveFunction : ISolver1d {
	371	private GFunctionWithShifts o_;
	372	private double Rs_;
	373	private double derivative_;
	374
	375	public ObjectiveFunction(GFunctionWithShifts o, double Rs) {
	376	o_ = o;
	377	Rs_ = Rs;
	378	}
	379	public override double value(double x) {
	380	double result = 0;
	381	derivative_ = 0;
	382	for (int i = 0; i < o_.accruals_.Count; i++) {
	383	double temp = o_.accruals_[i] * o_.swapPaymentDiscounts_[i] * Math.Exp(-o_.shapedSwapPaymentTimes_[i] * x);
	384	result += temp;
	385	derivative_ -= o_.shapedSwapPaymentTimes_[i] * temp;
	386	}
	387	result *= Rs_;
	388	derivative_ *= Rs_;
	389	double temp1 = o_.swapPaymentDiscounts_.Last() * Math.Exp(-o_.shapedSwapPaymentTimes_.Last() * x);
	390
	391	result += temp1 - o_.discountAtStart_;
	392	derivative_ -= o_.shapedSwapPaymentTimes_.Last() * temp1;
	393	return result;
	394	}
	395
	396	public override double derivative(double UnnamedParameter1) { return derivative_; }
	397	public void setSwapRateValue(double x) { Rs_ = x; }
	398	public GFunctionWithShifts gFunctionWithShifts() { return o_; }
	399	}
	400
	401	//===========================================================================//
	402	// GFunctionWithShifts //
	403	//===========================================================================//
	404	public GFunctionWithShifts(CmsCoupon coupon, Handle<Quote> meanReversion) {
	405	meanReversion_ = meanReversion;
	406	calibratedShift_ = 0.03;
	407	tmpRs_ = 10000000.0;
	408	accuracy_ = 1.0e-14;
	409
	410	SwapIndex swapIndex = coupon.swapIndex();
	411	VanillaSwap swap = swapIndex.underlyingSwap(coupon.fixingDate());
	412
	413	swapRateValue_ = swap.fairRate();
	414
	415	objectiveFunction_ = new ObjectiveFunction(this, swapRateValue_);
	416
	417	Schedule schedule = swap.fixedSchedule();
	418	Handle<YieldTermStructure> rateCurve = swapIndex.termStructure();
	419	DayCounter dc = swapIndex.dayCounter();
	420
	421	swapStartTime_ = dc.yearFraction(rateCurve.link.referenceDate(), schedule.startDate());
	422	discountAtStart_ = rateCurve.link.discount(schedule.startDate());
	423
	424	double paymentTime = dc.yearFraction(rateCurve.link.referenceDate(), coupon.date());
	425
	426	shapedPaymentTime_ = shapeOfShift(paymentTime);
	427
	428	List<CashFlow> fixedLeg = new List<CashFlow>(swap.fixedLeg());
	429	int n = fixedLeg.Count;
483	430
484	431	shapedSwapPaymentTimes_ = new List<double>();
…	…
486	433	accruals_ = new List<double>();
487	434
488		for(int i =0; i<n; ++i)
489		{
490		Coupon coupon1 = fixedLeg[i] as Coupon;
491		accruals_.Add(coupon1.accrualPeriod());
492		Date paymentDate = new Date(coupon1.date().serialNumber());
493		double swapPaymentTime = dc.yearFraction(rateCurve.link.referenceDate(), paymentDate);
494		shapedSwapPaymentTimes_.Add(shapeOfShift(swapPaymentTime));
495		swapPaymentDiscounts_.Add(rateCurve.link.discount(paymentDate));
496		}
497		discountRatio_ = swapPaymentDiscounts_.Last()/discountAtStart_;
498		}
499		public override double value(double Rs)
500		{
501		double calibratedShift = calibrationOfShift(Rs);
502		return Rs* functionZ(calibratedShift);
503		}
504		public override double firstDerivative(double Rs)
505		{
506		//double dRs = 1.0e-8;
507		//return (operator()(Rs+dRs)-operator()(Rs-dRs))/(2.0*dRs);
508		double calibratedShift = calibrationOfShift(Rs);
509		return functionZ(calibratedShift) + Rs * derZ_derX(calibratedShift)/derRs_derX(calibratedShift);
510		}
511		public override double secondDerivative(double Rs)
512		{
513		//double dRs = 1.0e-8;
514		//return (firstDerivative(Rs+dRs)-firstDerivative(Rs-dRs))/(2.0*dRs);
515		double calibratedShift = calibrationOfShift(Rs);
516		return 2.0derZ_derX(calibratedShift)/derRs_derX(calibratedShift) + Rs der2Z_derX2(calibratedShift)/Math.Pow(derRs_derX(calibratedShift),2.0)- Rs * derZ_derX(calibratedShift)*der2Rs_derX2(calibratedShift)/ Math.Pow(derRs_derX(calibratedShift),3.0);
517		}
518		}
519
	435	for (int i = 0; i < n; ++i) {
	436	Coupon coupon1 = fixedLeg[i] as Coupon;
	437	accruals_.Add(coupon1.accrualPeriod());
	438	Date paymentDate = new Date(coupon1.date().serialNumber());
	439	double swapPaymentTime = dc.yearFraction(rateCurve.link.referenceDate(), paymentDate);
	440	shapedSwapPaymentTimes_.Add(shapeOfShift(swapPaymentTime));
	441	swapPaymentDiscounts_.Add(rateCurve.link.discount(paymentDate));
	442	}
	443	discountRatio_ = swapPaymentDiscounts_.Last() / discountAtStart_;
	444	}
	445
	446	public override double value(double Rs) {
	447	double calibratedShift = calibrationOfShift(Rs);
	448	return Rs * functionZ(calibratedShift);
	449	}
	450
	451	public override double firstDerivative(double Rs) {
	452	//double dRs = 1.0e-8;
	453	//return (operator()(Rs+dRs)-operator()(Rs-dRs))/(2.0*dRs);
	454	double calibratedShift = calibrationOfShift(Rs);
	455	return functionZ(calibratedShift) + Rs * derZ_derX(calibratedShift) / derRs_derX(calibratedShift);
	456	}
	457
	458	public override double secondDerivative(double Rs) {
	459	//double dRs = 1.0e-8;
	460	//return (firstDerivative(Rs+dRs)-firstDerivative(Rs-dRs))/(2.0*dRs);
	461	double calibratedShift = calibrationOfShift(Rs);
	462	return 2.0 * derZ_derX(calibratedShift) / derRs_derX(calibratedShift) + Rs * der2Z_derX2(calibratedShift) / Math.Pow(derRs_derX(calibratedShift), 2.0) - Rs * derZ_derX(calibratedShift) * der2Rs_derX2(calibratedShift) / Math.Pow(derRs_derX(calibratedShift), 3.0);
	463	}
	464	}
520	465	}
521	466
522		//! CMS-coupon pricer
523		// ! Base class for the pricing of a CMS coupon via static replication
524		// as in Hagan's "Conundrums..." article
525		//
526		public abstract class HaganPricer: CmsCouponPricer
527		{
528		//
529		public override double swapletRate()
530		{
531		return swapletPrice()/(coupon_.accrualPeriod()*discount_);
532		}
533		public override double capletPrice(double effectiveCap)
534		{
535		// caplet is equivalent to call option on fixing
536		Date today = Settings.evaluationDate();
537		if (fixingDate_ <= today)
538		{
539		// the fixing is determined
540		double Rs = Math.Max(coupon_.swapIndex().fixing(fixingDate_)-effectiveCap, 0.0);
541		double price = (gearing_ Rs)(coupon_.accrualPeriod()*discount_);
542		return price;
543		}
544		else
545		{
546		double cutoffNearZero = 1e-10;
547		double capletPrice = 0;
548		if (effectiveCap < cutoffForCaplet_)
549		{
550		double effectiveStrikeForMax = Math.Max(effectiveCap,cutoffNearZero);
551		capletPrice = optionletPrice(Option.Type.Call, effectiveStrikeForMax);
552		}
553		return gearing_ * capletPrice;
554		}
555		}
556		public override double capletRate(double effectiveCap)
557		{
558		return capletPrice(effectiveCap)/(coupon_.accrualPeriod()*discount_);
559		}
560		public override double floorletPrice(double effectiveFloor)
561		{
562		// floorlet is equivalent to put option on fixing
563		Date today = Settings.evaluationDate();
564		if (fixingDate_ <= today)
565		{
566		// the fixing is determined
567		double Rs = Math.Max(effectiveFloor-coupon_.swapIndex().fixing(fixingDate_),0.0);
568		double price = (gearing_ Rs)(coupon_.accrualPeriod()*discount_);
569		return price;
570		}
571		else
572		{
573		double cutoffNearZero = 1e-10;
574		double floorletPrice = 0;
575		if (effectiveFloor > cutoffForFloorlet_)
576		{
577		double effectiveStrikeForMin = Math.Max(effectiveFloor,cutoffNearZero);
578		floorletPrice =optionletPrice(Option.Type.Put, effectiveStrikeForMin);
579		}
580		return gearing_ * floorletPrice;
581		}
582		}
583		public override double floorletRate(double effectiveFloor)
584		{
585		return floorletPrice(effectiveFloor)/(coupon_.accrualPeriod()*discount_);
586		}
587		//
588		public double meanReversion()
589		{
590		return meanReversion_.link.value();
591		}
592		public void setMeanReversion(Handle<Quote> meanReversion)
593		{
	467
	468	//===========================================================================//
	469	// HaganPricer //
	470	//===========================================================================//
	471	//! Base class for the pricing of a CMS coupon via static replication as in Hagan's "Conundrums..." article
	472	public abstract class HaganPricer : CmsCouponPricer {
	473	public override double swapletRate() {
	474	return swapletPrice() / (coupon_.accrualPeriod() * discount_);
	475	}
	476
	477	public override double capletPrice(double effectiveCap) {
	478	// caplet is equivalent to call option on fixing
	479	Date today = Settings.evaluationDate();
	480	if (fixingDate_ <= today) {
	481	// the fixing is determined
	482	double Rs = Math.Max(coupon_.swapIndex().fixing(fixingDate_) - effectiveCap, 0.0);
	483	double price = (gearing_ * Rs) * (coupon_.accrualPeriod() * discount_);
	484	return price;
	485	} else {
	486	double cutoffNearZero = 1e-10;
	487	double capletPrice = 0;
	488	if (effectiveCap < cutoffForCaplet_) {
	489	double effectiveStrikeForMax = Math.Max(effectiveCap, cutoffNearZero);
	490	capletPrice = optionletPrice(Option.Type.Call, effectiveStrikeForMax);
	491	}
	492	return gearing_ * capletPrice;
	493	}
	494	}
	495
	496	public override double capletRate(double effectiveCap) {
	497	return capletPrice(effectiveCap) / (coupon_.accrualPeriod() * discount_);
	498	}
	499
	500	public override double floorletPrice(double effectiveFloor) {
	501	// floorlet is equivalent to put option on fixing
	502	Date today = Settings.evaluationDate();
	503	if (fixingDate_ <= today) {
	504	// the fixing is determined
	505	double Rs = Math.Max(effectiveFloor - coupon_.swapIndex().fixing(fixingDate_), 0.0);
	506	double price = (gearing_ * Rs) * (coupon_.accrualPeriod() * discount_);
	507	return price;
	508	} else {
	509	double cutoffNearZero = 1e-10;
	510	double floorletPrice = 0;
	511	if (effectiveFloor > cutoffForFloorlet_) {
	512	double effectiveStrikeForMin = Math.Max(effectiveFloor, cutoffNearZero);
	513	floorletPrice = optionletPrice(Option.Type.Put, effectiveStrikeForMin);
	514	}
	515	return gearing_ * floorletPrice;
	516	}
	517	}
	518	public override double floorletRate(double effectiveFloor) {
	519	return floorletPrice(effectiveFloor) / (coupon_.accrualPeriod() * discount_);
	520	}
	521	//
	522	public double meanReversion() {
	523	return meanReversion_.link.value();
	524	}
	525	public void setMeanReversion(Handle<Quote> meanReversion) {
594	526	if (meanReversion_ != null)
595	527	meanReversion_.unregisterWith(update);
…	…
597	529	if (meanReversion_ != null)
598	530	meanReversion_.registerWith(update);
599		update();
600		}
601
602		//===========================================================================//
603		// HaganPricer //
604		//===========================================================================//
605		protected HaganPricer(Handle<SwaptionVolatilityStructure> swaptionVol, GFunctionFactory.YieldCurveModel modelOfYieldCurve, Handle<Quote> meanReversion) : base(swaptionVol)
606		{
607		modelOfYieldCurve_ = modelOfYieldCurve;
608		cutoffForCaplet_ = 2;
609		cutoffForFloorlet_ = 0;
610		meanReversion_ = meanReversion;
	531	update();
	532	}
	533
	534	protected HaganPricer(Handle<SwaptionVolatilityStructure> swaptionVol, GFunctionFactory.YieldCurveModel modelOfYieldCurve, Handle<Quote> meanReversion)
	535	: base(swaptionVol) {
	536	modelOfYieldCurve_ = modelOfYieldCurve;
	537	cutoffForCaplet_ = 2;
	538	cutoffForFloorlet_ = 0;
	539	meanReversion_ = meanReversion;
611	540
612	541	if (meanReversion_ != null)
613	542	meanReversion_.registerWith(update);
614		}
615		public override void initialize(FloatingRateCoupon coupon)
616		{
617		coupon_ = coupon as CmsCoupon;
618		gearing_ = coupon_.gearing();
619		spread_ = coupon_.spread();
620
621		fixingDate_ = coupon_.fixingDate();
622		paymentDate_ = coupon_.date();
623		SwapIndex swapIndex = coupon_.swapIndex();
624		rateCurve_ = swapIndex.termStructure().link;
625
626		Date today = Settings.evaluationDate();
627
628		if(paymentDate_ > today)
629		discount_ = rateCurve_.discount(paymentDate_);
630		else
631		discount_ = 1.0;
632
633		spreadLegValue_ = spread_ * coupon_.accrualPeriod()* discount_;
634
635		if (fixingDate_ > today)
636		{
637		swapTenor_ = swapIndex.tenor();
638		VanillaSwap swap = swapIndex.underlyingSwap(fixingDate_);
639
640		swapRateValue_ = swap.fairRate();
641
642		double bp = 1.0e-4;
643		annuity_ = (swap.floatingLegBPS()/bp);
644
645		int q = (int)swapIndex.fixedLegTenor().frequency();
646		Schedule schedule = swap.fixedSchedule();
647		DayCounter dc = swapIndex.dayCounter();
648		//DayCounter dc = coupon.dayCounter();
649		double startTime = dc.yearFraction(rateCurve_.referenceDate(), swap.startDate());
650		double swapFirstPaymentTime = dc.yearFraction(rateCurve_.referenceDate(), schedule.date(1));
651		double paymentTime = dc.yearFraction(rateCurve_.referenceDate(), paymentDate_);
652		double delta = (paymentTime-startTime) / (swapFirstPaymentTime-startTime);
653
654		switch (modelOfYieldCurve_)
655		{
656		case GFunctionFactory.YieldCurveModel.Standard:
657		gFunction_ = GFunctionFactory.newGFunctionStandard(q, delta, swapTenor_.length());
658		break;
659		case GFunctionFactory.YieldCurveModel.ExactYield:
660		gFunction_ = GFunctionFactory.newGFunctionExactYield(coupon_);
661		break;
662		case GFunctionFactory.YieldCurveModel.ParallelShifts:
663		{
664		Handle<Quote> nullMeanReversionQuote = new Handle<Quote>(new SimpleQuote(0.0));
665		gFunction_ = GFunctionFactory.newGFunctionWithShifts(coupon_, nullMeanReversionQuote);
666		}
667		break;
668		case GFunctionFactory.YieldCurveModel.NonParallelShifts:
669		gFunction_ = GFunctionFactory.newGFunctionWithShifts(coupon_, meanReversion_);
670		break;
671		default:
	543	}
	544
	545	public override void initialize(FloatingRateCoupon coupon) {
	546	coupon_ = coupon as CmsCoupon;
	547	gearing_ = coupon_.gearing();
	548	spread_ = coupon_.spread();
	549
	550	fixingDate_ = coupon_.fixingDate();
	551	paymentDate_ = coupon_.date();
	552	SwapIndex swapIndex = coupon_.swapIndex();
	553	rateCurve_ = swapIndex.termStructure().link;
	554
	555	Date today = Settings.evaluationDate();
	556
	557	if (paymentDate_ > today)
	558	discount_ = rateCurve_.discount(paymentDate_);
	559	else
	560	discount_ = 1.0;
	561
	562	spreadLegValue_ = spread_ * coupon_.accrualPeriod() * discount_;
	563
	564	if (fixingDate_ > today) {
	565	swapTenor_ = swapIndex.tenor();
	566	VanillaSwap swap = swapIndex.underlyingSwap(fixingDate_);
	567
	568	swapRateValue_ = swap.fairRate();
	569
	570	double bp = 1.0e-4;
	571	annuity_ = (swap.floatingLegBPS() / bp);
	572
	573	int q = (int)swapIndex.fixedLegTenor().frequency();
	574	Schedule schedule = swap.fixedSchedule();
	575	DayCounter dc = swapIndex.dayCounter();
	576	//DayCounter dc = coupon.dayCounter();
	577	double startTime = dc.yearFraction(rateCurve_.referenceDate(), swap.startDate());
	578	double swapFirstPaymentTime = dc.yearFraction(rateCurve_.referenceDate(), schedule.date(1));
	579	double paymentTime = dc.yearFraction(rateCurve_.referenceDate(), paymentDate_);
	580	double delta = (paymentTime - startTime) / (swapFirstPaymentTime - startTime);
	581
	582	switch (modelOfYieldCurve_) {
	583	case GFunctionFactory.YieldCurveModel.Standard:
	584	gFunction_ = GFunctionFactory.newGFunctionStandard(q, delta, swapTenor_.length());
	585	break;
	586	case GFunctionFactory.YieldCurveModel.ExactYield:
	587	gFunction_ = GFunctionFactory.newGFunctionExactYield(coupon_);
	588	break;
	589	case GFunctionFactory.YieldCurveModel.ParallelShifts: {
	590	Handle<Quote> nullMeanReversionQuote = new Handle<Quote>(new SimpleQuote(0.0));
	591	gFunction_ = GFunctionFactory.newGFunctionWithShifts(coupon_, nullMeanReversionQuote);
	592	}
	593	break;
	594	case GFunctionFactory.YieldCurveModel.NonParallelShifts:
	595	gFunction_ = GFunctionFactory.newGFunctionWithShifts(coupon_, meanReversion_);
	596	break;
	597	default:
672	598	throw new ApplicationException("unknown/illegal gFunction type");
673		}
674		vanillaOptionPricer_ = new BlackVanillaOptionPricer(swapRateValue_, fixingDate_, swapTenor_, swaptionVolatility().link);
675		}
676		}
677
678		protected YieldTermStructure rateCurve_;
679		protected GFunctionFactory.YieldCurveModel modelOfYieldCurve_;
680		protected GFunction gFunction_;
681		protected CmsCoupon coupon_;
682		protected Date paymentDate_;
683		protected Date fixingDate_;
684		protected double swapRateValue_;
685		protected double discount_;
686		protected double annuity_;
687		protected double gearing_;
688		protected double spread_;
689		protected double spreadLegValue_;
690		protected double cutoffForCaplet_;
691		protected double cutoffForFloorlet_;
692		protected Handle<Quote> meanReversion_;
693		protected Period swapTenor_;
694		protected VanillaOptionPricer vanillaOptionPricer_;
695		}
696
697
698		//! CMS-coupon pricer
699		// ! Prices a cms coupon via static replication as in Hagan's
700		// "Conundrums..." article via numerical integration based on
701		// prices of vanilla swaptions
702		//
703		public class NumericHaganPricer : HaganPricer
704		{
705		public NumericHaganPricer(Handle<SwaptionVolatilityStructure> swaptionVol, GFunctionFactory.YieldCurveModel modelOfYieldCurve, Handle<Quote> meanReversion, double lowerLimit, double upperLimit, double precision) : base(swaptionVol, modelOfYieldCurve, meanReversion)
706		{
707		upperLimit_ = upperLimit;
708		lowerLimit_ = lowerLimit;
709		requiredStdDeviations_ = 8;
710		precision_ = precision;
711		refiningIntegrationTolerance_ = 0.0001;
712
713		}
714
715		protected override double optionletPrice(Option.Type optionType, double strike)
716		{
717
718		ConundrumIntegrand integrand = new ConundrumIntegrand(vanillaOptionPricer_, rateCurve_, gFunction_, fixingDate_, paymentDate_, annuity_, swapRateValue_, strike, optionType);
719		stdDeviationsForUpperLimit_= requiredStdDeviations_;
720		double a;
721		double b;
722		double integralValue;
723		if (optionType ==Option.Type.Call)
724		{
725		upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
726		// while(upperLimit_ <= strike){
727		// stdDeviationsForUpperLimit_ += 1.;
728		// upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
729		// }
730		integralValue = integrate(strike, upperLimit_, integrand);
731		//refineIntegration(integralValue, *integrand);
732		}
733		else
734		{
735		a = Math.Min(strike, lowerLimit_);
736		b = strike;
737		integralValue = integrate(a, b, integrand);
738		}
739
740		double dFdK = integrand.firstDerivativeOfF(strike);
741		double swaptionPrice = vanillaOptionPricer_.value(strike, optionType, annuity_);
742
743		// v. HAGAN, Conundrums..., formule 2.17a, 2.18a
744		return coupon_.accrualPeriod() * (discount_/annuity_) * ((1 + dFdK) * swaptionPrice + ((int)optionType) *integralValue);
745		}
746
747		public double upperLimit()
748		{
749		return upperLimit_;
750		}
751		public double stdDeviations()
752		{
753		return stdDeviationsForUpperLimit_;
754		}
755
756		//private:
757		public abstract class Function
758		{
759		public abstract double value(double x);
760		}
761
762		//===========================================================================//
763		// NumericHaganPricer //
764		//===========================================================================//
765
766		public class VariableChange
767		{
768		public VariableChange(ref Func<double, double> f, double a, double b, int k)
769		{
	599	}
	600	vanillaOptionPricer_ = new BlackVanillaOptionPricer(swapRateValue_, fixingDate_, swapTenor_, swaptionVolatility().link);
	601	}
	602	}
	603
	604	protected YieldTermStructure rateCurve_;
	605	protected GFunctionFactory.YieldCurveModel modelOfYieldCurve_;
	606	protected GFunction gFunction_;
	607	protected CmsCoupon coupon_;
	608	protected Date paymentDate_;
	609	protected Date fixingDate_;
	610	protected double swapRateValue_;
	611	protected double discount_;
	612	protected double annuity_;
	613	protected double gearing_;
	614	protected double spread_;
	615	protected double spreadLegValue_;
	616	protected double cutoffForCaplet_;
	617	protected double cutoffForFloorlet_;
	618	protected Handle<Quote> meanReversion_;
	619	protected Period swapTenor_;
	620	protected VanillaOptionPricer vanillaOptionPricer_;
	621	}
	622
	623
	624	//===========================================================================//
	625	// NumericHaganPricer //
	626	//===========================================================================//
	627	// ! Prices a cms coupon via static replication as in Hagan's
	628	// "Conundrums..." article via numerical integration based on
	629	// prices of vanilla swaptions
	630	public class NumericHaganPricer : HaganPricer {
	631	public double upperLimit_;
	632	public double stdDeviationsForUpperLimit_;
	633	public double lowerLimit_;
	634	public double requiredStdDeviations_;
	635	public double precision_;
	636	public double refiningIntegrationTolerance_;
	637
	638	public NumericHaganPricer(Handle<SwaptionVolatilityStructure> swaptionVol, GFunctionFactory.YieldCurveModel modelOfYieldCurve, Handle<Quote> meanReversion, double lowerLimit, double upperLimit, double precision)
	639	: base(swaptionVol, modelOfYieldCurve, meanReversion) {
	640	upperLimit_ = upperLimit;
	641	lowerLimit_ = lowerLimit;
	642	requiredStdDeviations_ = 8;
	643	precision_ = precision;
	644	refiningIntegrationTolerance_ = 0.0001;
	645	}
	646
	647	protected override double optionletPrice(Option.Type optionType, double strike) {
	648	ConundrumIntegrand integrand = new ConundrumIntegrand(vanillaOptionPricer_, rateCurve_, gFunction_, fixingDate_, paymentDate_, annuity_, swapRateValue_, strike, optionType);
	649	stdDeviationsForUpperLimit_ = requiredStdDeviations_;
	650	double a;
	651	double b;
	652	double integralValue;
	653	if (optionType == Option.Type.Call) {
	654	upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
	655	// while(upperLimit_ <= strike){
	656	// stdDeviationsForUpperLimit_ += 1.;
	657	// upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
	658	// }
	659	integralValue = integrate(strike, upperLimit_, integrand);
	660	//refineIntegration(integralValue, *integrand);
	661	} else {
	662	a = Math.Min(strike, lowerLimit_);
	663	b = strike;
	664	integralValue = integrate(a, b, integrand);
	665	}
	666
	667	double dFdK = integrand.firstDerivativeOfF(strike);
	668	double swaptionPrice = vanillaOptionPricer_.value(strike, optionType, annuity_);
	669
	670	// v. HAGAN, Conundrums..., formule 2.17a, 2.18a
	671	return coupon_.accrualPeriod() * (discount_ / annuity_) * ((1 + dFdK) * swaptionPrice + ((int)optionType) * integralValue);
	672	}
	673
	674	public double upperLimit() { return upperLimit_; }
	675	public double stdDeviations() { return stdDeviationsForUpperLimit_; }
	676
	677	public double integrate(double a, double b, ConundrumIntegrand integrand) {
	678	double result = .0;
	679	//double abserr =.0;
	680	//double alpha = 1.0;
	681
	682	//double epsabs = precision_;
	683	//double epsrel = 1.0; // we are interested only in absolute precision
	684	//size_t neval =0;
	685
	686	// we use the non adaptive algorithm only for semi infinite interval
	687	if (a > 0) {
	688	// we estimate the actual boudary by testing integrand values
	689	double upperBoundary = 2 * a;
	690	while (integrand.value(upperBoundary) > precision_)
	691	upperBoundary *= 2.0;
	692	// sometimes b < a because of a wrong estimation of b based on stdev
	693	if (b > a)
	694	upperBoundary = Math.Min(upperBoundary, b);
	695
	696	GaussKronrodNonAdaptive gaussKronrodNonAdaptive = new GaussKronrodNonAdaptive(precision_, 1000000, 1.0);
	697	// if the integration intervall is wide enough we use the
	698	// following change variable x -> a + (b-a)*(t/(a-b))^3
	699	if (upperBoundary > 2 * a) {
	700	VariableChange variableChange = new VariableChange(integrand.value, a, upperBoundary, 3);
	701	result = gaussKronrodNonAdaptive.value(variableChange.value, .0, 1.0);
	702	} else {
	703	result = gaussKronrodNonAdaptive.value(integrand.value, a, upperBoundary);
	704	}
	705
	706	// if the expected precision has not been reached we use the old algorithm
	707	if (!gaussKronrodNonAdaptive.integrationSuccess()) {
	708	GaussKronrodAdaptive integrator = new GaussKronrodAdaptive(precision_, 1000000);
	709	result = integrator.value(integrand.value, a, b);
	710	}
	711	} // if a < b we use the old algorithm
	712	else {
	713	GaussKronrodAdaptive integrator = new GaussKronrodAdaptive(precision_, 1000000);
	714	result = integrator.value(integrand.value, a, b);
	715	}
	716	return result;
	717	}
	718
	719	public override double swapletPrice() {
	720	Date today = Settings.evaluationDate();
	721	if (fixingDate_ <= today) {
	722	// the fixing is determined
	723	double Rs = coupon_.swapIndex().fixing(fixingDate_);
	724	double price = (gearing_ * Rs + spread_) * (coupon_.accrualPeriod() * discount_);
	725	return price;
	726	} else {
	727	double atmCapletPrice = optionletPrice(Option.Type.Call, swapRateValue_);
	728	double atmFloorletPrice = optionletPrice(Option.Type.Put, swapRateValue_);
	729	return gearing_ * (coupon_.accrualPeriod() * discount_ * swapRateValue_ + atmCapletPrice - atmFloorletPrice) + spreadLegValue_;
	730	}
	731	}
	732
	733	public double resetUpperLimit(double stdDeviationsForUpperLimit) {
	734	//return 1.0;
	735	double variance = swaptionVolatility().link.blackVariance(fixingDate_, swapTenor_, swapRateValue_);
	736	return swapRateValue_ * Math.Exp(stdDeviationsForUpperLimit * Math.Sqrt(variance));
	737	}
	738
	739	public double refineIntegration(double integralValue, ConundrumIntegrand integrand) {
	740	double percDiff = 1000.0;
	741	while (Math.Abs(percDiff) < refiningIntegrationTolerance_) {
	742	stdDeviationsForUpperLimit_ += 1.0;
	743	double lowerLimit = upperLimit_;
	744	upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
	745	double diff = integrate(lowerLimit, upperLimit_, integrand);
	746	percDiff = diff / integralValue;
	747	integralValue += diff;
	748	}
	749	return integralValue;
	750	}
	751
	752	#region Nested classes
	753	public class VariableChange {
	754	private double a_, b_, width_;
	755	private Func<double, double> f_;
	756	private int k_;
	757
	758	public VariableChange(Func<double, double> f, double a, double b, int k) {
770	759	a_ = a;
771	760	b_ = b;
…	…
774	763	k_ = k;
775	764	}
776		public double value(double x)
777		{
	765
	766	public double value(double x) {
778	767	double newVar;
779	768	double temp = width_;
780		for (int i = 1; i < k_; ++i)
781		{
	769	for (int i = 1; i < k_; ++i) {
782	770	temp *= x;
783	771	}
…	…
785	773	return f_(newVar) * k_ * temp;
786	774	}
787		~~private double a_;~~
788		private double b_;
789		~~private double width_;~~
	775	}
	776
	777	public class Spy {
790	778	Func<double, double> f_;
791		private int k_;
792		}
793
794		public class Spy
795		{
796		public Spy(Func<double, double> f)
797		{
	779	private List<double> abscissas = new List<double>();
	780	private List<double> functionValues = new List<double>();
	781
	782	public Spy(Func<double, double> f) {
798	783	f_ = f;
799	784	}
800		public double value(double x)
801		{
	785	public double value(double x) {
802	786	abscissas.Add(x);
803	787	double value = f_(x);
…	…
805	789	return value;
806	790	}
807		Func<double, double> f_;
808		private List<double> abscissas = new List<double>();
809		private List<double> functionValues = new List<double>();
810		}
811
812		public class ConundrumIntegrand : Function
813		{
814
815		//===========================================================================//
816		// ConundrumIntegrand //
817		//===========================================================================//
818
819		public ConundrumIntegrand(VanillaOptionPricer o, YieldTermStructure curve, GFunction gFunction, Date fixingDate, Date paymentDate, double annuity, double forwardValue, double strike, Option.Type optionType)
820		{
821		vanillaOptionPricer_ = o;
822		forwardValue_ = forwardValue;
823		annuity_ = annuity;
824		fixingDate_ = fixingDate;
825		paymentDate_ = paymentDate;
826		strike_ = strike;
827		optionType_ = optionType;
828		gFunction_ = gFunction;
829		}
830		public override double value(double x)
831		{
832		double option = vanillaOptionPricer_.value(x, optionType_, annuity_);
833		return option * secondDerivativeOfF(x);
834		}
835		protected double functionF (double x)
836		{
837		double Gx = gFunction_.value(x);
838		double GR = gFunction_.value(forwardValue_);
839		return (x - strike_) * (Gx/GR - 1.0);
840		}
841		public double firstDerivativeOfF (double x)
842		{
843		double Gx = gFunction_.value(x);
844		double GR = gFunction_.value(forwardValue_);
845		double G1 = gFunction_.firstDerivative(x);
846		return (Gx/GR - 1.0) + G1/GR * (x - strike_);
847		}
848		public double secondDerivativeOfF(double x)
849		{
850		double GR = gFunction_.value(forwardValue_);
851		double G1 = gFunction_.firstDerivative(x);
852		double G2 = gFunction_.secondDerivative(x);
853		return 2.0 * G1/GR + (x - strike_) * G2/GR;
854		}
855
856		protected double strike()
857		{
858		return strike_;
859		}
860		protected double annuity()
861		{
862		return annuity_;
863		}
864		protected Date fixingDate()
865		{
866		return fixingDate_;
867		}
868		protected void setStrike(double strike)
869		{
870		strike_ = strike;
871		}
872
873		protected VanillaOptionPricer vanillaOptionPricer_;
874		protected double forwardValue_;
875		protected double annuity_;
876		protected Date fixingDate_;
877		protected Date paymentDate_;
878		protected double strike_;
879		protected Option.Type optionType_;
880		protected GFunction gFunction_;
881		}
882
883		public double integrate(double a, double b, ConundrumIntegrand integrand)
884		{
885		double result =.0;
886		//double abserr =.0;
887		//double alpha = 1.0;
888
889
890		//double epsabs = precision_;
891		//double epsrel = 1.0; // we are interested only in absolute precision
892		//size_t neval =0;
893
894		// we use the non adaptive algorithm only for semi infinite interval
895		if (a>0)
896		{
897
898		// we estimate the actual boudary by testing integrand values
899		double upperBoundary = 2 *a;
900		while(integrand.value(upperBoundary)>precision_)
901		upperBoundary *=2.0;
902		// sometimes b < a because of a wrong estimation of b based on stdev
903		if (b > a)
904		upperBoundary = Math.Min(upperBoundary, b);
905
906		Func<double, double> f;
907		GaussKronrodNonAdaptive gaussKronrodNonAdaptive = new GaussKronrodNonAdaptive(precision_, 1000000, 1.0);
908		// if the integration intervall is wide enough we use the
909		// following change variable x -> a + (b-a)*(t/(a-b))^3
910		if (upperBoundary > 2 *a)
911		{
912		int k = 3;
913		Func<double, double> temp = integrand.value;
914		VariableChange variableChange = new VariableChange(ref temp, a, upperBoundary, k);
915
916		throw new NotImplementedException();
917
918		f = integrand.value;
919		// TODO f = boost.bind(variableChange.value, variableChange, _1);
920
921		result = gaussKronrodNonAdaptive.value(f, .0, 1.0);
922		}
923		else
924		{
925		f = integrand.value;
926		result = gaussKronrodNonAdaptive.value(f, a, upperBoundary);
927		}
928
929		// if the expected precision has not been reached we use the old algorithm
930		if (!gaussKronrodNonAdaptive.integrationSuccess())
931		{
932		GaussKronrodAdaptive integrator = new GaussKronrodAdaptive(precision_, 1000000);
933		result = integrator.value(integrand.value, a, b);
934		}
935
936		} // if a < b we use the old algorithm
937		else
938		{
939		GaussKronrodAdaptive integrator = new GaussKronrodAdaptive(precision_, 1000000);
940		result = integrator.value(integrand.value, a, b);
941		}
942
943		return result;
944		}
945		public override double swapletPrice()
946		{
947
948		Date today = Settings.evaluationDate();
949		if (fixingDate_ <= today)
950		{
951		// the fixing is determined
952		double Rs = coupon_.swapIndex().fixing(fixingDate_);
953		double price = (gearing_ Rs + spread_)(coupon_.accrualPeriod()*discount_);
954		return price;
955		}
956		else
957		{
958		double atmCapletPrice = optionletPrice(Option.Type.Call, swapRateValue_);
959		double atmFloorletPrice = optionletPrice(Option.Type.Put, swapRateValue_);
960		return gearing_ (coupon_.accrualPeriod() discount_ * swapRateValue_ + atmCapletPrice - atmFloorletPrice) + spreadLegValue_;
961		}
962		}
963		public double resetUpperLimit(double stdDeviationsForUpperLimit)
964		{
965		//return 1.0;
966		double variance = swaptionVolatility().link.blackVariance(fixingDate_,swapTenor_,swapRateValue_);
967		return swapRateValue_ * Math.Exp(stdDeviationsForUpperLimit *Math.Sqrt(variance));
968		}
969		public double refineIntegration(double integralValue, ConundrumIntegrand integrand)
970		{
971		double percDiff = 1000.0;
972		while(Math.Abs(percDiff) < refiningIntegrationTolerance_)
973		{
974		stdDeviationsForUpperLimit_ += 1.0;
975		double lowerLimit = upperLimit_;
976		upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
977		double diff = integrate(lowerLimit, upperLimit_, integrand);
978		percDiff = diff/integralValue;
979		integralValue += diff;
980		}
981		return integralValue;
982		}
983
984		public double upperLimit_;
985		public double stdDeviationsForUpperLimit_;
986		public double lowerLimit_;
987		public double requiredStdDeviations_;
988		public double precision_;
989		public double refiningIntegrationTolerance_;
990		}
991
992		//! CMS-coupon pricer
993		public class AnalyticHaganPricer : HaganPricer
994		{
995
996		//===========================================================================//
997		// AnalyticHaganPricer //
998		//===========================================================================//
999
1000		public AnalyticHaganPricer(Handle<SwaptionVolatilityStructure> swaptionVol, GFunctionFactory.YieldCurveModel modelOfYieldCurve, Handle<Quote> meanReversion) : base(swaptionVol, modelOfYieldCurve, meanReversion)
1001		{
1002		}
1003
1004		//Hagan, 3.5b, 3.5c
1005		protected override double optionletPrice(Option.Type optionType, double strike)
1006		{
1007		double variance = swaptionVolatility().link.blackVariance(fixingDate_, swapTenor_, swapRateValue_);
1008		double firstDerivativeOfGAtForwardValue = gFunction_.firstDerivative(swapRateValue_);
1009		double price = 0;
1010
1011		double CK = vanillaOptionPricer_.value(strike, optionType, annuity_);
1012		price += (discount_/annuity_)*CK;
1013		double sqrtSigma2T = Math.Sqrt(variance);
1014		double lnRoverK = Math.Log(swapRateValue_/strike);
1015		double d32 = (lnRoverK+1.5 *variance)/sqrtSigma2T;
1016		double d12 = (lnRoverK+.5 *variance)/sqrtSigma2T;
1017		double dminus12 = (lnRoverK-.5 *variance)/sqrtSigma2T;
1018
1019		CumulativeNormalDistribution cumulativeOfNormal = new CumulativeNormalDistribution();
	791	}
	792
	793	//===========================================================================//
	794	// ConundrumIntegrand //
	795	//===========================================================================//
	796	public class ConundrumIntegrand : IValue {
	797	public ConundrumIntegrand(VanillaOptionPricer o, YieldTermStructure curve, GFunction gFunction, Date fixingDate, Date paymentDate, double annuity, double forwardValue, double strike, Option.Type optionType) {
	798	vanillaOptionPricer_ = o;
	799	forwardValue_ = forwardValue;
	800	annuity_ = annuity;
	801	fixingDate_ = fixingDate;
	802	paymentDate_ = paymentDate;
	803	strike_ = strike;
	804	optionType_ = optionType;
	805	gFunction_ = gFunction;
	806	}
	807
	808	public double value(double x) {
	809	double option = vanillaOptionPricer_.value(x, optionType_, annuity_);
	810	return option * secondDerivativeOfF(x);
	811	}
	812
	813	protected double functionF(double x) {
	814	double Gx = gFunction_.value(x);
	815	double GR = gFunction_.value(forwardValue_);
	816	return (x - strike_) * (Gx / GR - 1.0);
	817	}
	818
	819	public double firstDerivativeOfF(double x) {
	820	double Gx = gFunction_.value(x);
	821	double GR = gFunction_.value(forwardValue_);
	822	double G1 = gFunction_.firstDerivative(x);
	823	return (Gx / GR - 1.0) + G1 / GR * (x - strike_);
	824	}
	825
	826	public double secondDerivativeOfF(double x) {
	827	double GR = gFunction_.value(forwardValue_);
	828	double G1 = gFunction_.firstDerivative(x);
	829	double G2 = gFunction_.secondDerivative(x);
	830	return 2.0 * G1 / GR + (x - strike_) * G2 / GR;
	831	}
	832
	833	protected double strike() { return strike_; }
	834	protected double annuity() { return annuity_; }
	835	protected Date fixingDate() { return fixingDate_; }
	836	protected void setStrike(double strike) { strike_ = strike; }
	837
	838	protected VanillaOptionPricer vanillaOptionPricer_;
	839	protected double forwardValue_;
	840	protected double annuity_;
	841	protected Date fixingDate_;
	842	protected Date paymentDate_;
	843	protected double strike_;
	844	protected Option.Type optionType_;
	845	protected GFunction gFunction_;
	846	}
	847	#endregion
	848	}
	849
	850	//===========================================================================//
	851	// AnalyticHaganPricer //
	852	//===========================================================================//
	853	public class AnalyticHaganPricer : HaganPricer {
	854	public AnalyticHaganPricer(Handle<SwaptionVolatilityStructure> swaptionVol, GFunctionFactory.YieldCurveModel modelOfYieldCurve, Handle<Quote> meanReversion)
	855	: base(swaptionVol, modelOfYieldCurve, meanReversion) {
	856	}
	857
	858	//Hagan, 3.5b, 3.5c
	859	protected override double optionletPrice(Option.Type optionType, double strike) {
	860	double variance = swaptionVolatility().link.blackVariance(fixingDate_, swapTenor_, swapRateValue_);
	861	double firstDerivativeOfGAtForwardValue = gFunction_.firstDerivative(swapRateValue_);
	862	double price = 0;
	863
	864	double CK = vanillaOptionPricer_.value(strike, optionType, annuity_);
	865	price += (discount_ / annuity_) * CK;
	866	double sqrtSigma2T = Math.Sqrt(variance);
	867	double lnRoverK = Math.Log(swapRateValue_ / strike);
	868	double d32 = (lnRoverK + 1.5 * variance) / sqrtSigma2T;
	869	double d12 = (lnRoverK + .5 * variance) / sqrtSigma2T;
	870	double dminus12 = (lnRoverK - .5 * variance) / sqrtSigma2T;
	871
	872	CumulativeNormalDistribution cumulativeOfNormal = new CumulativeNormalDistribution();
1020	873	double N32 = cumulativeOfNormal.value(((int)optionType) * d32);
1021	874	double N12 = cumulativeOfNormal.value(((int)optionType) * d12);
1022	875	double Nminus12 = cumulativeOfNormal.value(((int)optionType) * dminus12);
1023
1024		price += ((int)optionType) * firstDerivativeOfGAtForwardValue * annuity_ * swapRateValue_ * (swapRateValue_ * Math.Exp(variance) * N32- (swapRateValue_+strike) * N12 + strike * Nminus12);
1025		price *= coupon_.accrualPeriod();
1026		return price;
1027		}
1028
1029		//Hagan 3.4c
1030		public override double swapletPrice()
1031		{
1032
1033		Date today = Settings.evaluationDate();
1034		if (fixingDate_ <= today)
1035		{
1036		// the fixing is determined
1037		double Rs = coupon_.swapIndex().fixing(fixingDate_);
1038		double price = (gearing_ Rs + spread_)(coupon_.accrualPeriod()*discount_);
1039		return price;
1040		}
1041		else
1042		{
1043		double variance = swaptionVolatility().link.blackVariance(fixingDate_, swapTenor_, swapRateValue_);
1044		double firstDerivativeOfGAtForwardValue = gFunction_.firstDerivative(swapRateValue_);
1045		double price = 0;
1046		price += discount_ *swapRateValue_;
1047		price += firstDerivativeOfGAtForwardValue annuity_ swapRateValue_* swapRateValue_*(Math.Exp(variance)-1.0);
1048		return gearing_ * price * coupon_.accrualPeriod() + spreadLegValue_;
1049		}
1050		}
1051		}
	876
	877	price += ((int)optionType) * firstDerivativeOfGAtForwardValue * annuity_ * swapRateValue_ * (swapRateValue_ * Math.Exp(variance) * N32 - (swapRateValue_ + strike) * N12 + strike * Nminus12);
	878	price *= coupon_.accrualPeriod();
	879	return price;
	880	}
	881
	882	//Hagan 3.4c
	883	public override double swapletPrice() {
	884	Date today = Settings.evaluationDate();
	885	if (fixingDate_ <= today) {
	886	// the fixing is determined
	887	double Rs = coupon_.swapIndex().fixing(fixingDate_);
	888	double price = (gearing_ * Rs + spread_) * (coupon_.accrualPeriod() * discount_);
	889	return price;
	890	} else {
	891	double variance = swaptionVolatility().link.blackVariance(fixingDate_, swapTenor_, swapRateValue_);
	892	double firstDerivativeOfGAtForwardValue = gFunction_.firstDerivative(swapRateValue_);
	893	double price = 0;
	894	price += discount_ * swapRateValue_;
	895	price += firstDerivativeOfGAtForwardValue * annuity_ * swapRateValue_ * swapRateValue_ * (Math.Exp(variance) - 1.0);
	896	return gearing_ * price * coupon_.accrualPeriod() + spreadLegValue_;
	897	}
	898	}
	899	}
1052	900	}

QLNet

Navigation

Changeset 220

Legend:

trunk/QLNet/QLNet/Cashflows/ConundrumPricer.cs

Download in other formats: