Assembla home | Assembla project page
 

Changeset 220

Show
Ignore:
Timestamp:
06/28/08 13:41:06 (3 months 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   &nb