Assembla home | Assembla project page
 

Changeset 220

Show
Ignore:
Timestamp:
06/28/08 13:41:06 (1 year ago)
Author:
snovik
Message:

Fix: ConundrumPricer? (please have a look)

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved
  • trunk/QLNet/QLNet/Cashflows/ConundrumPricer.cs

    r216 r220  
    11/* 
    22 Copyright (C) 2008 Toyin Akin (toyin_akin@hotmail.com) 
     3 Copyright (C) 2008 Siarhei Novik (snovik@gmail.com) 
    34   
    45 This file is part of QLNet Project http://www.qlnet.org 
     
    2324 
    2425namespace 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        } 
    2952        } 
    3053 
    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            } 
    45130                } 
    46131 
    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); 
    171157 
    172158                List<CashFlow> fixedLeg = new List<CashFlow>(swap.fixedLeg()); 
    173                                int n = fixedLeg.Count; 
     159                int n = fixedLeg.Count; 
    174160                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) { 
    311276                        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)) 
    336300                    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)) 
    344308                    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)) 
    383346                    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)) 
    392355                    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; 
    483430 
    484431                shapedSwapPaymentTimes_ = new List<double>(); 
     
    486433                accruals_ = new List<double>(); 
    487434 
    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.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); 
    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        } 
    520465        } 
    521466 
    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) { 
    594526            if (meanReversion_ != null) 
    595527                meanReversion_.unregisterWith(update); 
     
    597529            if (meanReversion_ != null) 
    598530                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; 
    611540 
    612541            if (meanReversion_ != null) 
    613542                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: 
    672598                        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) { 
    770759                a_ = a; 
    771760                b_ = b; 
     
    774763                k_ = k; 
    775764            } 
    776             public double value(double x) 
    777            
     765 
     766            public double value(double x)
    778767                double newVar; 
    779768                double temp = width_; 
    780                 for (int i = 1; i < k_; ++i) 
    781                 { 
     769                for (int i = 1; i < k_; ++i) { 
    782770                    temp *= x; 
    783771                } 
     
    785773                return f_(newVar) * k_ * temp; 
    786774            } 
    787             private double a_; 
    788             private double b_; 
    789             private double width_; 
     775        } 
     776 
     777        public class Spy { 
    790778            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) { 
    798783                f_ = f; 
    799784            } 
    800             public double value(double x) 
    801             { 
     785            public double value(double x) { 
    802786                abscissas.Add(x); 
    803787                double value = f_(x); 
     
    805789                return value; 
    806790            } 
    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(); 
    1020873            double N32 = cumulativeOfNormal.value(((int)optionType) * d32); 
    1021874            double N12 = cumulativeOfNormal.value(((int)optionType) * d12); 
    1022875            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    } 
    1052900}