Modelling Financial Derivatives with Mathematica

Update, Errata and V5/4 Compatibility Notes

February 2004
William T. Shaw

Quantitative Analysis Group
Nomura International plc
william.shaw@uk.nomura.com

The Mathematical Institute
24-29 St Giles
Oxford 
OX1 3LB
shaww@maths.ox.ac.uk

Read this first

This notebook has recently been updated to improve compatibility with version 5 of Mathematica. Specifically, this first section comprises a list of changes to get things working, that have been tested under version 5.01 of Mathematica, released in January 2004. The second section consists of a list of supplementary observations, while the third consists of a number of known bug fixes to deal with issues raised by the release of version 4. Logically, you should make the version 4 changes and then the Version 5 changes. However, version 5 compatibility with the early chapters is ensured by making some very simple quick fixes, whereas you will not run into much trouble with the version 4 patches until you start to use the compiler. All of these changes will be incorporated into a revision of the book at some point. (After Complex Analysis with Mathematica (CUP 2004/5) is released!)

Version 5 Quick Fixes

Global Changes

Ignore spelling errors!

In running the code in several chapters you will see messages like (this one is from Chapter 2).

expression  =  Exp[-x^2];
StandardForm[expression]

This sort of message is a helpful/annoying reminder that you have made a mistake/invented a symbol with a name similar to an existing one. In this and most cases thsy can be ignored, but check to see if you really have spelt somefing rong. You get these messages in amy recent version of Mathematica.

The Norm function and Packages

In Mathematica 5 the function Norm is a core kernel function:

?Norm

In the book I used this (rather recklessly!) to label the cumulative normal distribution function. The notation commonly used in finance is N(x), but as Mathematica reserves the N function for enforcing numerical output:

?N

I had to pick something else. Henceforth the function will be called Ncdf. This should be applied globally throughout all chapters and the two packages that came with the book. For example, in Chapter 1, the following should be applied:

Ncdf[x_] := 1/2*(1 + Erf[x/Sqrt[2]]); 
done[s_, σ_, k_, t_, r_, q_] := 
((r - q)*t + Log[s/k])/(σ*Sqrt[t]) + (σ*Sqrt[t])/2; 
dtwo[s_, σ_, k_, t_, r_, q_] := 
done[s, σ, k, t, r, q] - σ*Sqrt[t]; 
  
BlackScholesCall[s_, k_, v_, r_, q_, t_] := 
s*Exp[-q*t]*Ncdf[done[s, v, k, t, r, q]] - k*Exp[-r*t]*Ncdf[dtwo[s, v, k, t, r, q]];

The two packages BlackScholes.m and AnalyticAmerican.m should be similarly edited, along with all other chapters. (Check to see if you have received the updated packages with this notebook).

Precision

There has been some tightening up of how this works. For example, in Chapter 2, the sequence of numerical investigations should now be realised as follows:

N[Sqrt[51]]

Precision[%]

and you need to ask about this:

$MachinePrecision

(3/4)^100

N[(3/4)^100,10]

0.75^100

N[0.75^100,10]

New fitting functions

Note the appearance of the FindFit function. This can be used to do fitting based on other objective functions, for example, Least Absolute Deviation fitting (which is MLE based on double-sided exponential rather than Gaussian noise) (More to come on this!)

Integrate

Symbolic integration appears to be rather slower in version 5, and return results that are carefully couched to the point of major irritation. While I have not got to the bottom of this, part of the story is due to much more careful checking of cases and presentation of results.

Options[Integrate]

?Assumptions

The following examples are interesting in this respect:

Timing[
      Integrate[
        1/(x^2  +  a^2),  {x,  -Infinity,  
          Infinity}]]

Timing[
      Integrate[1/(x^2  +  b^2),  {x,  -Infinity,  Infinity},  
        GenerateConditions  ->  False]]

The output of the following is what you really expect to see:

Timing[
      Integrate[1/(x^2  +  c^2),  {x,  -Infinity,  Infinity},  
        Assumptions  ->  c > 0]]

It's almost as if you really want Mathematica to have a DWIM interface (do what I mean), but it is better to have a pedantic system and then you must tell it about your own (often implicit) assumptions. DWIM behaviour can be sort-of achieved by setting  GenerateConditions -> False. These comments supercede those made in the version 4 update below. I do not now recommend the fix used previously!

Other Chapter-specific issues note covered by global points and existing V4 updates

Chapter 2

See the version 4 update below - the Fourier and Laplace transforms are now in the kernel. Chapter 2 needs a different example!

A bug noticed recently is the non-existent listc used in the MultipleListPlot example. Delete it.

Chapter 3

While there are no M5 issues I should note that I am rewriting Section 3.5 as it is confusing. Watch out for an update.

Chapter 4

Lots of changes due to the differing behaviour of Integrate, and it is best to use the Assumptions option. If you need help with this, e-mail me at shaww@maths.ox.ac.uk for a new version of Chapter 4.

Mathematical Updates

1. Introduction

This notebook contains information known as of July 1st 1999 about 

(i) errors in the first printing of the book
(ii) errors in the second printing of the book
(iii) changes associated with Version 4 of Mathematica. 
(iv) other issues raised.

Other comments of potential interest have been added, some as a result of comments from readers, whose contribution is gratefully acknowledged. 

This note inherits all disclaimers of the original test with regard to both the printed material and its electronic form.

With regard to version 4 of Mathematica, this notebook is a first run through by the author designed to highlight differences and fix any compatibility porblems. There may be further revisions to this note

(a) in the light of feedback from Wolfram Research on issues raised;

(b) becoming more V4 savvy with regard to numerical optimization.

Comments on how to best make use of numerical optimization features are most welcome and will be added to this note after moderation. Revisions to this notebook will be posted in the same locations as this original.

The version 4 checks and changes reported here was based on re-running the code in the book on a Power Macintosh G3 running OS8.x. If you get different effects on other hardware/OS, especially Wintel, please let me know and I will update this update.

2. Errata

In the following sub-sections corrections to the first and second printing known as of end-June 1999. These comments concern mistakes in the text/equations/code related to either general development or Mathematica 3.X implementations. See the section on V4 compatibility for further issues raised by changes made to the V4 software. More important items are denoted by *.

The first/second printing is ISBN 0 521 59233 X. They are distinguished by "Reprinted 1999" on the page containing ISBN and disclaimer information, in the second printing.

A special thanks goes to Roger Wilson and David Kelly, who spotted several of the glitches.

Corrections to First Printing Only

*1. p. 407. In the list of parameter definitions for fine path simulation the explanations of n and navg have been swapped. n is the number of paths to be sampled, and navg is the number of sample points on each path.

2. p. 76 equation (39). There is a spurious extra d t attached to the combination (Subscript[q, X]) (Subscript[Δ, 2]) X within the bracketed collection of terms multiplying d t. This equation should read:

d P=(((∂V)/(∂X))-(Subscript[Δ, 1]) S'-(Subscript[Δ, 2]))d X+(((∂V)/(∂S'))-(Subscript[Δ, 1]) X) d S'+ d t ( ((∂V)/(∂t))+(1/2)(Subscript[σ, X]^2)(X^2) (((∂^2) V)/(∂(X^2))) 
				       + X S'ρ (Subscript[σ, S']) (Subscript[σ, X]) (((∂^2) V)/(∂S' ∂X))+(1/2)(Subscript[σ, S']^2) ((S')^2)(((∂^2) V)/(∂((S')^2))) -(Subscript[q, X]) (Subscript[Δ, 2]) X-(Subscript[q, S' X]) (Subscript[Δ, 1]) S' X)

3. p. 102,  middle of the page. The minus sign should not be in the argument of PowerExpand, and PowerExpand has not in fact produced the right answer, for x > a, which should be

1 - erf(((x - a)/(2 (Sqrt[τ]))))

which agrees with the code line after "Picking out..." and with the n = 0 similarity solution. The correct answer under V3.x can be obtained by using instead:

PowerExpand[-% /. (a - x) ->  (x - a)]

4. p. 123, equations (181,184,187). The partial derivatives w.r.t "t" should be w.r.t "tau". The one in (178), before the time-reversed and vol-scaled time tau is introduced, is correct.

5. Also on p. 123, the gremlins had a feast on the brackets. Here is how p. 123 should really look (without equation numbers):

x = log( S/E)
τ = 
      (1/2) (σ^2) (T - t)

If q(τ) is the τ-dependent function representing the continuously compounded dividend yield DY(t) (this may contain discrete pieces),

\(TraditionalForm\`k\_1[τ] = \(\(2\ R\ \((τ)\)\)\/σ\^2 = \(2\ \(r(\ t)\)\)\/σ\^2\)\)
\(TraditionalForm\`k\_2[τ] = \(\(2\ \((R\ \((τ)\) - q(τ))\)\)\/σ\^2 = \(2\ \((r\ \((t)\) - DY\ \((t)\))\)\)\/σ\^2\)\)

Then we obtain

\(TraditionalForm\`∂V\/∂τ + \ k\_1[τ] 
          V = \ \ ∂\^2 V\/∂x\^2\  + \ \ \((k\_2[τ]\  - 1\ )\) ∂V\/∂x\  + \ \(2\/σ\^2\) 
          K\)

Thus far the change of variables has proceeded as for the European options discussed previously. The final part of the previous analysis does not work for the present problem unless (Subscript[k, 2])[τ]  is constant. To treat the time-dependent case, let

z = x + F(τ)

and set

V(
        x, τ) = u(
          z, τ) (e^(-B(τ)))

Making the required change of variable leads to

\(TraditionalForm\`u\ \((k\_1[τ] - ∂B\/∂τ)\) + ∂u\/∂τ = ∂\^2 u\/∂z\^2 + \(2  Ke\^\(B(τ)\)\)\/σ\^2 + \((\(-\(∂F\/∂τ\)\) + k\_2[τ] - 
              1)\)\ ∂V\/∂x\)

so we can reduce the problem by making the choices

(Subscript[k, 1])[τ] = ∂(B/∂)τ

(Subscript[k, 2])[τ] - 
        1 = ∂(F/∂)τ

whence we get

\(TraditionalForm\`∂u\/∂τ = ∂\^2 u\/∂z\^2 + Q(τ)\)

where the diffusion equation now has a "source term" Q representing the coupon payments in these coordinates:

\(TraditionalForm\`Q(τ)\  = \(2  K\ e\^\(B(τ)\)\)\/σ\^2\)

(End of correction to p. 123)

Corrections to First and Second Printing

6. On p. 83, the commutation of the Black-Scholes operator with the (homogeneity) operator

S ((∂ )/(∂S))

relies on the dividend term being of the form D  =  q S. I found out about Peter Carr's related (and much more detailed) work, deriving the rho-delta link and the vega-gamma identity using probabilistic methods, too late for inclusion in the text, but his working paper on the matter is available on the Web. It is "Deriving Derivatives of Derivative Securities", on www.math.nyu.edu/research/carrp/papers. The operator commutation relations are common to both the probabilistic and PDE approaches.

7. The gremlins got into the editing of the supersymmetric examples in Section 19.4 on pp. 377-378. In the definition of SSdataCApp, the function call should be to FastSSEuropeanCallApp, not FastSSSEuropeanCallApp! Also, in the subsequent plots of the convergence, the string "OTS" should be replaced by "SS". So , for example, the first graphic is obtained by

Show[explot, SSPlotC,   PlotRange  ->  {0.69,  0.70},  
        DisplayFunction  ->  $DisplayFunction];

8.  In Chapter 15, the column headings labelling the numerical and exact solution are the wrong way round, in several examples. On p. 285, 287, 289, the headings should be in the order: S, numerical, exact, error.

9. Regarding Section 6.5 on p141-142, the simple rational approximations to the cumulative Normal distribution were in fact first popularized by Mark Rubenstein and John Cox in their wonderful book, "Options Markets", published by Prentice-Hall in 1985. The function I have called HFuncOne appears on p. 343 of that text within the instructions for a programmable calculator, and HFuncTwo is given on pp. 261-262. I apologize to Profs Robenstein and Cox both for this specific mistake and not referring to their text in many other places, especially for discussions of market issues. 

10. In Chapter 23, p. 447, the 2D numerical integration expressed by BasketCallIntegral produces unstable results - this is currently under investigation. 

11. In Chapter 25, despite saying that I have suppressed the dependence on T, it is left in - the last line on this page and first line of p. 475 should read

TraditionalForm[U[r, t]]

At the top of p. 476, equation (25) is (obviously) obtained by a similar application of Traditional Form to beqn.

12. (second printing only) On p. 123 there is still one outer bracket missing in the lower part of equation (180), which should read

\(TraditionalForm\`k\_2[τ] = \(\(2\ \((R\ \((τ)\) - q(τ))\)\)\/σ\^2 = \(2\ \((r\ \((t)\) - DY\ \((t)\))\)\)\/σ\^2\)\)

These are all known corrections as of 1st Juy 1999.

3. Issues Raised by Recent Papers + More on "Implied Vol"

This section makes brief comments regarding a couple of recent publications that either directly refer to the book, or come up with ideas that I would have included if writing now.

Discretely Observed Barriers

The first paper is the rather wonderful "Connecting discrete and continuous path-dependent options", by M. Broadie, P. Glasserman, and S.G. Kou, (BGK) published in Finance and Stochastics, vol 3, p. 55 (1999).  It generalizes work first reported in Math. Finance vol 7, p. 325 (1997).

I particularly like these papers as they present a means of working with simple closed-form analytics, for a large class of path-dependent options, but where the barrier/max/min along the path is observed only at (regular) discrete time points. This eliminates, with a cunning analytical approximation, a large volume of numerical modelling.

You should really read this paper to get the idea in full, and were I re-writing MFDwM now, I would give several examples both in the development of the analytics, and in the comparisons with Monte Carlo later. Their result is very concisely stated for the case of (zero-rebate) barrier options, where it suffices to make a simple adjustment to the location of the barrier as a function of the observation interval.

Simple Analytical Adjustment for Barriers

The simplest form of the BGK correction applies to barrier options, where we consider zero rebate for simplicity. The analytical model in Chapter 8 then applies with zero rebate, for the European Barrier option. The idea is that a discretely-observed barrier at one level can be replaced by a continuously observed barrier that is slighly further away from the current asset price. (This is obviously a very naive description of the BGK analysis, but is at least intuitive). BGK calculate the correction to the barrier location in closed form! Let H denote the location of the barrier that is observed at a constant time interval Δ t  =  t/m where t is the time to maturity and m is the number of observations. 

For an up-option, the effective barrier to be used in a continuous model is

(Subscript[H, u]) = 
      H (e^(β σ (Sqrt[Δ t])))

For a down-option:

\(TraditionalForm\`H\_d = 
      H\ e\^\(\(-β\)\ σ\ \@\(Δ\ t\)\)\)

where

β=-((ζ((1/2)))/((Sqrt[2 π])))

and ζ(x) denotes the Riemann zeta function evaluated at x. The error in the valuation in this approximate treatment is o(1/(Sqrt[m])). This is all easy to put into Mathematica, which does of course have the zeta-function built in. (If you must use another environment β = 0.5826). Here's the new model for an up and out call - the others are corrected similarly. It is assumed that the relevant functions from Chapter 8 have been entered into the kernel.

BarrUp[h_,  σ_,  t_,  navg_]  :=  
      Module[{rootdt  =  
            Sqrt[t/(navg - 1)],  β  =  -Zeta[1/2]/
              Sqrt[2 Pi]}, 
h*
          Exp[β*σ*rootdt]]

So weekly observation of a barrier at 120 gives an effective barrier at

BarrUp[120,  0.2,  1,  53]

The (zero-rebate) Up and Out Call is then modelled in the BGK approximation by

AnalyticDiscreteUpAndOutCall[p_, k_, h_, sd_, r_, q_, t_,  
        m_] :=
Module[{H  =  
            BarrUp[h,  sd,  t,  m]}, 
{H,  
          UpAndOutCall[0,  p, k, H, sd, r, q, t]}]

TableForm[
      Table[Flatten[{10^k,  
            AnalyticDiscreteUpAndOutCall[100,  100,  120,  0.2,  0.05,  
              0.02,  1,  10^k]}],  {k,  1,  9}]]

The continuum limit is just

UpAndOutCall[0,  100,  100,  120,  0.2,  0.05,  0.02,  
      1]

That is all there is to it, except that you need to be careful with vega and theta, and this simple model is less accurate when you are close to the barrier.

Asian Options

The following is a reply to a paper by Fu, Madan and Wang published in the Journal of Computational Finance early in 1999. It is self-explanatory.  This reply was submitted to the journal - I have to edit out some of the Mathematica details - subject to this my reply is probably going to be published by the journal later in 1999.

Pricing Asian Options with Mathematica:
A Reply to:
Pricing Continuous Asian Options: A Comparison of Monte Carlo and Laplace Transform Inversion Methods, by Fu et al, JCF, 1999

Introduction

In a recent paper, Fu, Madan and Wang (1999) have made a very interesting comparison of several semi-analytical and numerical methods for the pricing of continously averaged Asian Options. One of the methods cited (the "Shaw" method) is based on a presentation in Chapter 10 of the book (Shaw, 1998), where Mathematica is used to simplify and evaluate the Geman-Yor (1993) description of Asian options. In their paper, Fu et al suggest that the "Shaw" approach can give answers that are in error by an order of magnitude, and/or is highly inefficient. Here these claims are rebutted, and we show how-precision answers may be obtained from Mathematica 3.0 when the inversion method is properly applied. This is in a computing framework where the full model takes no more than fourteen lines of code. Furthermore, a straightforward but time-consuming Monte-Carlo analysis suggests that the Mathematica approach may be highly accurate.

2. Recap of the Geman-Yor model

An excellent description of this model is given by Fu et al, so only the briefest of summaries will be given. We are investigating the model of the continuously and arithmetically averaged Asian Call given by Geman and Yor (1993) (henceforth abbreviated as GY), as also summarized by Eydeland and Geman (1995). Suppose that the current time is t, and that the option matures at a time T > t. The averaging is arithmetic, continuous, and began at a time (Subscript[t, 0]) ≤ t. Suppose that the known average value of the underlying over the time interval [(Subscript[t, 0]),  t] is ES. GY define the following changes of variables, as also given by Shaw (1998), and presented here as inputs in Mathematica's StandardForm:

τ[T_,  
          t_,  σ_]  :=  σ^2 (T - t)/4;
ν[r_,  q_,  σ_]  :=  
        2 (r - q)/σ^2  -  1;
α[S_,  ES_,  K_,  σ_,  T_,  t_,  
          to_] :=  σ^2/(4*S) (K*(T - to)  -  (t - to)*
                ES);
μ[ν_,  p_]  :=  Sqrt[ν^2  +  2*p]

The value of the average price option is then given by

\(TraditionalForm\`\(e\^\(\(-r\)\ \((T - t)\)\)\ 4\ S\ \(C(τ, \ ν, α)\)\)\/\(\((T - t\_0)\)\ σ\^2\)\)

The remaining function C[ν,  τ,  α] is not given explicitly, but GY give its Laplace transform,

\(TraditionalForm\`U(
        p, ν, α) = ∫\_0\%∞\( 
            C(τ, \ ν, α)\)\ \(e\^\(\(-p\)\ τ\)\) dτ\)

as an integral:

\(TraditionalForm\`U(
        p, ν, α) = \(∫\_0\%\(1\/\(2\ α\)\)\(x\^\(\(μ - ν\)\/2 - 2\)\ \((1 - 2\ α\ x)\)\^\(\(μ + ν\)\/2 + 1\)\ e\^\(-x\)\) dx\)\/\(p\ \((p - 2\ ν - 2)\)\ \(Γ(\(μ - ν\)\/2 - 1)\)\)\)

where μ is as given above as a function of p and ν. GY develop a series description of the transform and show how it can be inverted. We shall now explore how this can be managed and simplified in Mathematica.

3. Mathematica Implementation of Arithmetic Asian

We have already done the first part of the translation to software - we first entered the definitions of the various basic functions. Now we enter the definition of the integral that is part of the transform, and request immediate evaluation. The use of the [[2]] notation is to extract the answer relevant to parameters of interest:

F[
          p_, μ_,  ν_, α_] = (∫_0^((1/(2 α))))(x^(((μ - ν)/2) - 2)) (((1 - 2 α x))^(((μ + ν)/2) + 1)) Exp[-x] \[DifferentialD]x;
G[p_, μ_,   ν_,  α_]  =  F[
          p, μ,  ν, α][[2]]

There are further cancellations when we insert the other terms that make up the transform - we also present the result in traditional mathematical notation:

U[p_, μ_, ν_, α_] = 
      Simplify[G[p, μ,  ν, α]/(p (p - 2 ν - 2) Gamma[((μ - ν)/2) - 1])]; 
    TraditionalForm[U[p, μ, ν, α]]

We now have the ingredients to build the Mathematica model of the arithmetic average price Asian Call. In the following the Laplace transform inversion is done by direct numerical integration along the truncated Bromwich contour. This contour is a simple vertical line to the right of any finite singularities. In Shaw (1998), the truncation is set initially at a value of 500. However, I issued the following warning: "readers may need to experiment and increase the integration range for peculiar parameter values". That is, the integration range must be taken to be long enough that the results have stabilized. The inaccurate results cited by Fu et al in their Table 4 arise solely from using a truncation that is too short - the default value of 500 can easily be seen to be too small for low volatility cases by plotting the absolute value of the transform. To emphasize the role of the integration truncation parameter, we here make it an explicit argument of the option valuation function, and also parametrize the number of recursive subdivisions of the range needed - we may also need to increase this as the integration range is stretched.

Off[NIntegrate::slwcon]

\(AriAsianPriceCall[S_, \ ES_, \ K_, \ r_, \ q_, \ σ_, \ T_, \ 
        t_, \ to_, \ truncation_, \ recursion_] := 
	
      Module[{ti = τ[T, \ t, \ σ], n = ν[r, \ q, \ σ], 
          a\  = \ α[S, ES, K, σ, T, t, to], \ contour}, 
	
        contour\  = \ 2  n + 3; 
	
        ac\  = \ Re[\ 
            1/\((2\ Pi)\)*
              NIntegrate[	
                U[\((contour + I\ p)\), μ[n, \((contour + I\ p)\)], n, a]*
                  Exp[\((contour + I\ p)\)*ti], {p, \ \(-truncation\), \ 
                  truncation}, \ MaxRecursion\  -> \ recursion]]; 
\ 
        Exp[\(-r\)*\((T\  - \ t)\)]*4*S/\((\((T\  - \ to)\)*σ^2)\)*
          ac]\)

In their paper, Fu et al present, in Table 4, a list of evaluations for a set of parameter values with a common strike set at K = 2. Here are the first three rows with a truncation of 500, and then 750 and 1000.

Table[
      AriAsianPriceCall[1.9, 0, 2, 0.05,  0,  0.5,  1, 0, 0,  trunc,  
        11],  {trunc,  500,  1000,  250}]

Table[
      AriAsianPriceCall[2.0, 0, 2, 0.05,  0,  0.5,  1, 0, 0,  trunc,  
        11],  {trunc,  500,  1000,  250}]

Table[
      AriAsianPriceCall[2.1, 0, 2, 0.05,  0,  0.5,  1, 0, 0,  trunc,  
        11],  {trunc,  500,  1000,  250}]

The values of 0.193, 0.246 and 0.306 quoted by Fu et al as arising from this method therefore represent outputs that are stable as a function of increasing truncation. These same values were given by myself in Shaw (1998). The remaining four rows of Table 4 correspond to (in some cases a computation timing is reported)

Table[
      AriAsianPriceCall[2.0, 0, 2, 0.02,  0,  0.1,  1, 0, 0,  trunc,  
        11],  {trunc,  500,  1000,  250}]

Table[
      AriAsianPriceCall[2.0, 0, 2, 0.18,  0,  0.3,  1, 0, 0,  trunc,  
        11],  {trunc,  500,  1000,  250}]

Table[
      AriAsianPriceCall[2.0, 0, 2, 0.0125,  0,  0.25,  2, 0, 0,  trunc,  
        11],  {trunc,  500,  1000,  250}]

Timing[
      Table[AriAsianPriceCall[2.0, 0, 2, 0.05,  0,  0.5,  2, 0, 0,  
          trunc,  11],  {trunc,  500,  1000,  250}]]

It is clear that for most of these new examples the default integration range is insufficient. Working in reverse order, the last case is clearly stable at 0.350. Row 6 is almost stable. We consider the next few truncations, adjusting the subdivision parameter accordingly as we go.

Table[
      AriAsianPriceCall[2.0, 0, 2, 0.0125,  0,  0.25,  2, 0, 0,  trunc,  
        14],  {trunc,  1250,  2000,  250}]

So row 6 of the Shaw method indeed yields the quoted value of 0.172. Turning our attention to row 5, we consider

Table[
      AriAsianPriceCall[2.0, 0, 2, 0.18,  0,  0.3,  1, 0, 0,  trunc,  
        14],  {trunc,  1250,  2000,  250}]

so that the truncation-stable result is 0.218(4), rather than 0.217 quoted by Fu et al. It is row four that is most problematic. For this case, within the text, Fu et al do comment on the effect of re-running the integration with larger range. Some experimentation leads to rather larger truncations, with stability achieved a truncation of 40000, with the following result:

Timing[
      AriAsianPriceCall[2.0, 0, 2, 0.02,  0,  0.1,  1, 0, 0,  40000,  
        18]]

This timing is for my current "office computer", which is a PowerMacintosh G3 at 433MHz. The times recorded are therefore rather more reasonable than the 16200 CPU seconds reported by Fu et al. Nevertheless, this does raise an efficiency question in relation to other implementations. However, my view is that direct brute-force integration along the original Bromwich contour is completely inappropriate for the low volatility case (σ = 0.1) considered here - indeed it is probably sub-optimal for most cases. This particular example says rather more about the unsuitability of direct integration along the Bromwich contour, as σ → 0, than it does about any alleged inefficiencies in the Mathematica integrator. The right way of dealing with the low volatility case is probably to deform the contour to one with faster convergence. This is under investigation and will be reported elsewhere.

I also have to point out that the matter of "efficiency" has as much to do with the time taken to code up the model as it does with execution time. In a financial environment, the cost of a programmer for one year is now about two orders of magnitude greater than the one-off cost of a state-of-the-art Pentium III computer. Development time is a key cost element, and the implementation here took about 6 hours of my time to put together. Furthermore, we should appreciate that the continuous-time Asian is not a regular desk-traded instrument, rather, it is a limiting case that is useful for testing the integrity of Monte Carlo simulations and PDE approaches to this problem. It's execution time is therefore rather unimportant.

On a related matter, assessing the relative accuracy of various approaches, the Monte Carlo simulations presented in Table 4 of Fu et al are of little use, as the number of samples was far too small for the standard error to be sufficiently low to allow any useful conclusions to be drawn. I re-ran some of the cases myself, using 3650 sampling points per year (i.e., this is my own version of their case "MC10", or ten fixes per day), and several millions of (antithetic) paths. This was also done in Mathematica using a batch processing version of the simulation algorithms given in Shaw (1998). The cases re-run and values obtained were:

Case 1: 3.52 million paths, Value 0.19317, SE < 0.0002;
Case 4: 2.68 million paths. Value 0.055976, SE < 0.00005;
Case 5: 3.24 million paths. Value 0.21823, SE < 0.00015.

which compare very well with the Mathematica continuum values of 0.19317,  0.055986, and 0.21839.  (A copy of the notebook used to generate these results is available on request from the author.) This comparison, and an inspection of the other results in Fu's Table 4, does hint that the Mathematica implementation might be more accurate than any of the other schemes presented by Fu et al, though it would require a further investigation of the MC100 case, and, preferably, an understanding of an analogue of the "Euler-Maclaurin" error formula, for integrating random walks, to nail this down. But there is certainly no evidence to suggest that the Mathematica implementation of the Geman-Yor model is anything other than a high-precision implementation.

References

Eydeland, A. and Geman, H., 1995, Asian options revisited: inverting the Laplace transform, RISK Magazine, March.

Fu, M.C., Madan, D.B. and Wang, T., 1999, Pricing Continuous Asian Options: A Comparison of Monte Carlo and Laplace Transform Inversion Methods, Journal of Computational Finance.

Geman, H. and Yor, M., 1993, Bessel processes, Asian options, and perpetuities, Mathematical Finance, 3, p. 349.

Shaw, W.T., 1998, Modelling Financial Derivatives with Mathematica, Cambridge University Press.

Implied Volatility - yet more nonsense...

I have been continuing to amuse myself in my spare time exploring the bad behaviour of the "implied volatility" construction. In the book I had a good rant in Chapter 1 about the unstabled/multivalued nature of the valuation as a function of σ, and highlighted the particularly bad behaviour of barrier options in this respect. It is a fascinating exercise to pursue this question for simple European binary options, for in this case the calculation (including the multi-valuedness) is analytically tractable using only the simple theory of quadratic equations. (No FindRoot or Newton-Raphson needed). In the book the reader was encouraged to pursue the corresponding analytics for power and log options. The binary analysis goes as follows.

The value of a European Binary Call option is

(Subscript[C, B]) = 
      B (e^(-r t)) Norm ((Subscript[d, 2]))

and that for a Put is

(Subscript[P, B]) = 
      B (e^(-r t)) Norm(- 
            (Subscript[d, 2]))

where (Subscript[d, 2])is the usual thing, which we write in the form

\(TraditionalForm\`d\_2 = α\/\(σ\ \@t\) - \(σ\ \@t\)\/2\)

where

\(TraditionalForm\`α = 
      log(\(S\ e\^\(\(-q\)\ t\)\)\/\(K\ e\^\(\(-r\)\ t\)\))\)

Now suppose we are given market prices (Subscript[C, M]) or (Subscript[P, M]) for such instruments. Provided, in each case

0 < (Subscript[C, M]) < 
      B (e^(-r t))

0 < (Subscript[P, M]) < 
      B (e^(-r t))

the pricing formula can be partially inverted to find (Subscript[d, 2])as a function of market price. We just have, respectively

(Subscript[d, 2]) = 
      (Norm^(-1)) ((Subscript[C, M]) (e^(r t))/
              B)

and, for a Put:

\(TraditionalForm\`d\_2 = \(-Norm\^\(-1\)\)\ \((P\_M\ e\^\(r\ t\)/B)\)\)

The Mathematica implementation of the inverse error function is (Ch. 21, p. 401):

Convert[x_] := (Sqrt[2]) InverseErf[2 x - 1]

Plot[Convert[x],  {x,  0,  1},  
        PlotRange  ->  {-3,  3}];

So (Subscript[d, 2]) is given, and what we want is σ. We have a relation between (Subscript[d, 2]),  σ,  α which is equivalent to a quadratic equation:

t (σ^2) + 2 (Sqrt[t]) (Subscript[d, 2]) σ - 
        2 α = 0

This has (possibly complex) solutions of the form

\(TraditionalForm\`σ = \(\(-d\_2\) ± \@\(d\_2\%2 + 2\ α\)\)\/\@t\)

and the product of the two solutions is -2 α/t. We observe immediately that provided α > 0, there is a unique positive solution of the quadratic for σ, given by taking the upper (+) sign of the square root in the solution. The positivity condition on α is equivalent to the condition

S (e^(-q t)) > 
      K (e^(-r t))

When this condition fails (e.g. out of the money call) there is trouble. Expressed as a function of σ

((∂(Subscript[d, 2]))/(∂σ))=-(α/((σ^2) (Sqrt[t]))) - ((Sqrt[t])/2)

and this will equal zero when

α = ((-1)/2) (σ^2) t

This is not possible when α > 0, but when α < 0 there is a critical value of volatility:

\(TraditionalForm\`σ\_C\^2 = \(\(-\(\(2\ α\)\/t\)\) = \(-\(\(2\ \(log(\(S\ e\^\(\(-q\)\ t\)\)\/\(K\ e\^\(\(-r\)\ t\)\))\)\)\/t\)\)\)\)

where there is a maximum value of (Subscript[d, 2]), and hence a maximum possible Call price, which we call (Subscript[C, crit]), and a minimum possible Put price, (Subscript[P, crit]). We proceed further for the case of a Call - the Put goes similarly.

If in the case α < 0 the market call price exceeds that given by (Subscript[σ, c]), there is no implied volatility consistent with the given market price. If it is less than this critical value, there are two implied volatilities. 

So in total there are three cases for a Call provided 0 < (Subscript[C, M]) < B (e^(-r t)):

(i) α > 0: a unique real and positive implied volatility;
(ii) α<0 and (Subscript[C, M]) >  (Subscript[C, crit]): no implied volatility;
(ii) α<0 and (Subscript[C, M]) <  (Subscript[C, crit]): two implied volatilitities.

This can all be implemented very easily:

InverseN[x_] := Sqrt[2]*InverseErf[2*x - 1];

CallDTwo[CallMkt_, B_, r_, t_] := InverseN[CallMkt*Exp[r*t]/B]; 
PutDTwo[PutMkt_, B_, r_, t_] := -InverseN[PutMkt*Exp[r*t]/B];

Alpha[S_, K_, r_, q_, t_] := Log[S*Exp[(r - q)*t]/K];

CallMax[B_, S_, K_, r_, q_, t_] := Module[{a = Alpha[S, K, r, q, t]}, 
    If[a >= 0, B*Exp[-r*t], B*Exp[-r*t]*Norm[-Sqrt[-2*a]]]];

PutMin[B_, S_, K_, r_, q_, t_] := Module[{a = Alpha[S, K, r, q, t]}, 
    If[a >= 0, 0, B*Exp[-r*t]*Norm[Sqrt[-2*a]]]];

BinaryCallImpVol[CallMkt_, B_, S_, K_, r_, q_, t_] := 
Which[
CallMkt >= B*Exp[-r*t], 
    "Market Price Far Too High", 
CallMkt <= 0, 
    "Market Price Far Too Low", 
True, 
	Module[{dtwo, alpha, root, rawvol, realvol}, 
	dtwo = CallDTwo[CallMkt, B, r, t]; 
	alpha = Alpha[S, K, r, q, t]; 
	root = Sqrt[dtwo^2 + 2*alpha]; 
	Which[
	alpha >= 0, 
	-dtwo/Sqrt[t] + root/Sqrt[t], 
	True, 
	If[CallMkt > CallMax[B, S, K, r, q, t], 
	"Market Price Too High", 
	-dtwo/Sqrt[t] + root/Sqrt[t]*{1, -1}]]]
	];

BinaryPutImpVol[PutMkt_, B_, S_, K_, r_, q_, t_] := 
Which[
PutMkt >= B*Exp[-r*t], 
   "Market Price Far Too High", 
PutMkt <= 0, 
   "Market Price Far Too Low", 
True, 
    Module[{dtwo, alpha, root, rawvol, realvol}, 
    dtwo = PutDTwo[PutMkt, B, r, t]; 
    alpha = Alpha[S, K, r, q, t];
    root = Sqrt[dtwo^2 + 2*alpha]; 
    Which[
    alpha >= 0, 
    -dtwo/Sqrt[t] + root/Sqrt[t], 
    True, 
    If[PutMkt < PutMin[B, S, K, r, q, t], 
    "Market Price Too Low", 
    -dtwo/Sqrt[t] + root/Sqrt[t]*{1, -1}]]]
    ];

Here are some examples with a Call:

BinaryCallImpVol[10,  10,  100,  100,  0.05,  0,  1]

BinaryCallImpVol[5,  10,  100,  100,  0.05,  0,  1]

Going further in the money is fine:

BinaryCallImpVol[5,  10,  105,  100,  0.05,  0,  1]

When we come out of the money a quite plausible market price is has no matching implied vol within the model:

BinaryCallImpVol[4,  10,  90,  100,  0.05,  0,  1]

And which of the following implied vols is the "right" one?

BinaryCallImpVol[3.5,  10,  90,  100,  0.05,  0,  1]

If you do not believe this makes sense, try calculating and plotting the call value as a function of σ with these parameters, using the following code lifted from Chapter 6:

dtwo[s_, v_, k_, t_, r_, q_] := ((r - q)*t + Log[s/k])/(v*Sqrt[t]) - (v*Sqrt[t])/2;

BinaryCall[B_, S_, K_, v_, r_, q_, t_] := B*Exp[-r*t]*Norm[dtwo[S, v, K, t, r, q]];

{BinaryCall[10,  90,  100,  0.282064,  0.05,  0,  1],  
      BinaryCall[10,  90,  100,  0.392539,  0.05,  0,  1]}

Plot[
        BinaryCall[10,  90,  100,  σ,  0.05,  0,  
          1],  {σ,  0.01,  0.5}];

The extreme instability and 2-valued-ness of the inverse mapping is manifest, and has emerged through the simple use of the analytical theory of quadratic equations. As with barrier options, the problem arises when vega changes sign or is just very small. A common mistake is to assume that there is no problem with the vanilla calls and puts. Although vega is technically positive, it is very small for options that are deep in the money, and the region over which it is very small grows as maturity approaches. So inversion of the mapping is meaningless in these regions, and it follows that such portions of the vol surface, calculated off market prices, don't mean much either.

4. Reviews

You might like to look at the March 99 issue of RISK magazine, which contained a very positive review of this book. There are some interesting, and mostly very positive, comments on Amazon.com, which I note is also happy to publish anonymous comments without moderation. Feel free to post your own views here. (Do a book search with author = Shaw, title = Derivatives).

Version 4 Changes

Before going into a chapter by chapter analysis of the impact of V4, it is helpful to note that there are a few global issues which impact several chapters. First you should note that most of the analytical code in Chapters 5-13 works without modification. The main issues relate to a compiler change that affects some code blocks for the FD, tree and Monte Carlo numerical tools in later chapters, and some changes to transform calculus conventions which affect the general theoretical development in the very early chapters. There is also what in my view is a bug in the formatting of integration output which I hope will be fixed. There are other minor issues here and there, of no great signficance - this notebook is only rather long as I have tried to trap all points everywhere they appear.

Overview of Main Points

Global Change - Impact of Compiler Bug Fix

I am grateful to Tom Wickham-Jones (TWJ) for explaining an undocumented change in the way that Compile works that breaks several of the compiled numerical routines (for FD, trees, Monte Carlo). It is easy to fix the code. 

The changes eliminate a type of bug than can arise when a variable is only initialized once inside a conditional construct and then used outside (if it is initialized inside a conditional statement and only used inside there is no problem). The bug can produce gibberish, as evidenced by this simple example from TWJ, which compiles happily under V3.X:

$Version

f  =  
      Compile[{x},  Module[ {t},  If[x > 5,  t = 5.5]; t]]

The following nonsense, or something equivalent, will  then result:

{f[1],  f[6],  f[1]}

The version 4 compiler does not allow this, and returns an error:

$Version

f  =  
      Compile[{x},  Module[ {t},  If[x > 5,  t = 5.5]; t]]

To get this through the compiler, we need to do something like:

g  =  
      Compile[{x},  
        Module[ {t = 1.0},  If[x > 5,  t = 5.5]; t]]

and then all is good:

{g[1],  g[6],  g[1]}

Although none of the examples in the book can generate this type of nonsense, they do involve local Module variables, some of which only get initialized inside a function within the Module, but are reported in the output. An example that is closest in spirit to the example supplied by TWJ is our compiled function for making samples from the Normal distribution quickly. Under V4 we get a cascade of complaints starting with the same initvar message:

\(\(norm = \ 
        Compile[{mu, \ 
            sigma}, \ 
		Module[{va, \ vb, \ rad\  = \ 2.0, \ 
              den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; 
					vb\  = \ 
                  2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			mu + 
              sigma*va*den]];\)\)

The problem is that va is in the output line, but is not initialized (but unlike the TWJ example, will be when the function is called). The problem is fixed by initializing va as follows:

\(\(norm = \ 
        Compile[{mu, \ 
            sigma}, \ 
		Module[{va = 1. , \ vb, \ rad\  = \ 2.0, \ 
              den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; 
					vb\  = \ 
                  2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			mu + 
              sigma*va*den]];\)\)

All of the problems with the compiler in V4 are of this type, and are fixed in exactly the same way. You might find it instructive to do so yourself, but I have supplied fixed code in the Sections below, for every function where it is to be applied. If you have a copy of the electronic version on disk, I recommend you paste these changes into the notebooks for the relevant chapter.

Global Change -Transform Calculus now in kernel; new conventions

The default conventions have now changed, and the manner in which phase and normalizations are controlled are different. Furthermore, the  Fourier and Laplace transform functions are now in core kernel.  Loading a transform package as before does no harm, and you are reminded that they are now in the kernel:

Needs["Calculus`LaplaceTransform`"]

and similarly with Fourier transforms. For example, with the new default conventions, we have

FourierTransform[Exp[-x^2/(4*t)]/Sqrt[4 Pi t],  x,  
      w]

If you want the results to all look like the V3 defaults in the book, you need to use the
FourierParameters -> {1,1} 
option setting, as follows:

FourierTransform[Exp[-x^2/(4*t)]/Sqrt[4 Pi t],  x,  w,  
      FourierParameters  ->  {1,  1}]

Let's go back, just to check:

InverseFourierTransform[Exp[-t w^2],  w,  
      x, 
FourierParameters  ->  {1,  
          1}]

This setting preserves the default conventions used in the V3 packages which were mostly employed in the text (except as noted). A list of parameter settings is given on p. 854 of the V4 book. My own preference is in fact for the "Modern Physics" setting of {0, 1}, which happens to be Mathematica's new default, but in order to simplify V3-V4 cross-compatibility, I suggest the {1,1} setting.

Global Change - Integration Formatting

Well, yours truly definitely thinks this is a bug.

Version 3.0 integration makes use of a certain approach to formatting output. This will be illustrated by use of an integral of the probability density function for the gamma distribution. This is chosen as the integral is manifestly unity. Under Version 3.X, we got

basic  =  
      Integrate[λ^p x^(p - 1) Exp[-λ x]/
            Gamma[p],  {x,  0,  Infinity}]

This is a sensible output, which also ensures that the following simple answer extraction works:

basic[[2]]

Version 4.0 integration gives a rather different default output.

basic  =  
      Integrate[λ^p x^(p - 1) Exp[-λ x]/
            Gamma[p],  {x,  0,  Infinity}]

My view is the format of this last output is just plain stupid, and it also causes the following to fail:

basic[[2]]

However, there are two ways around this, (which also work in V3.X). You can use instead, either

Integrate[λ^p x^(p - 1) Exp[-λ x]/
          Gamma[p],  {x,  0,  Infinity}, 
      GenerateConditions  ->  False]

fix  =  If[_,  expr_,  _]  ->  expr

basic  /.  fix

I have explained below, chapter by chapter, where the use of the result[[2]] construction must be avoided. The raw output of basic under V4 is ridiculuous, IMHO. I suggest that you avoid overusing the GenerateConditions -> False, (a) because you may miss an important parameter restriction,  (b) it is a lot of typing! I now use /. fix as it is as much typing as before, and you can do a default integration, then, e.g.

Integrate[λ^p x^(p - 1) Exp[-λ x]/
            Gamma[p],  {x,  0,  Infinity}];

new  =  %  /.  fix

to carry on to the next stage of your analysis without the If wrapping. If you have finalized a sequence of calculations, this can be abbreviated to

Integrate[λ^p x^(p - 1) Exp[-λ x]/
            Gamma[p],  {x,  0,  Infinity}] /.  fix

and this is the approach suggested for the book. One might argue that my previous approach, which relied on the particular format of an output, was sloppy. Perhaps, but I never imagined software would produce output of the form exemplified above by the V4 output for basic!

Issues with NonLinearFit

This only comes up in Chapter 26 in this book but I use it a lot elsewhere. The following code, which appeared originally in "Applied Mathematica", and which played nicely with version 2.2, now fails! If you make up some data as follows:

exptdata  =  
        Table[{x,  
            1/((x - 10)^2  +  0.1)  +  
              1/((x - 13)^2  +  0.5)  +  0.3 *Random[]},  {x,  
            0,  20,  0.1}];

This obviously has peaks at 10 and 13.

ListPlot[exptdata,  PlotJoined  ->  True,  
        PlotRange  ->  All];

Needs["Statistics`NonlinearFit`"]

You are supposed to be able to fit a function to this, using any available information to supply guesses about some of the parameters. Given a quick look at the data, we'll try to give very rough guesses for the peak locations, at 11 and 14 (we do not want to help too much!):

NonlinearFit[exptdata, 
        a  +  b/((x - c)^2  +  d)  +  
          e/((x - f)^2  +  g),  
        x,  {a, b, {c, 11}, d,  e,  {f,  14},  g},  
        MaxIterations  ->  40];

It can't cope with a starting suggestion or two. It does run, with no suggestions for starting values, but gives answers that are complete rubbish, as in the following:

NonlinearFit[exptdata, 
      a  +  b/((x - c)^2  +  d)  +  
        e/((x - f)^2  +  g),  x,  {a, b, c, d,  e,  f,  g},  
      MaxIterations  ->  40]

This problem now has a history - it worked in V2.2, failed in V3.0.0, was fixed and a MathSource patch supplied. However, I can no longer find a bug fix on MathSource and having reported and obtained a reply the problem, it turns out that only way to get this to work is to supply starting guesses for ALL the parameters - the following is not too bad:

NonlinearFit[exptdata, 
      a  +  b/((x - c)^2  +  d)  +  
        e/((x - f)^2  +  g),  
      x,  {{a, 1}, {b, 1}, {c, 11}, {d, 1},  {e, 1},  {f,  14},  {g, 
          1}},  MaxIterations  ->  40]

The view from WRI at the time of writing is that this function was never intended to cope with only partial suggestions for parameters, and that it worked "by accident" in V2.2. 

Meantime, if you use this function, you must either supply starting values for all parameters, or supply no values. Answers from the latter can sometimes be nonsense.

Other Chapter by Chapter Comments

Chapter 1

No known V4 issues.

Chapter 2

There are no major problems here but a large list of minor changes, as follows.

Standard Form now uses double struck characters for certain symbols. For example, on p. 17:

expression  =  Exp[-x^2];
StandardForm[expression]

On page 21, the operation of N on a previous numerical computation is slightly different. Now we have

0.75^100

N[%,  10]

This may be a bug - it look's rather odd to me, especially as does

N[(3/4)^100,  10]

On p. 26, NSolve returns the results in a different order.

On p. 35, the short form of an expression continues to be unpredictable, depending on the version. I even got two different answers running the same thing in two windows. There must be an option setting I do not know about. On p. 44/45 and 46 there is more of the influence of minor typesetting changes, for example:

DSolve[{y'[x] + x y[x] == 0, y[0] == 1}, y[x], x]

qsoln = Solve[3 x^2 - 4 x + 6 == 0, x]

The double struck character is no longer reserved for Traditional Form. There are many more such trivial changes throughout the book - this will be the last mention of them. 

On p. 49 we get more functions returned thus:

?*Sin*

The transform functions have now been moved into the core kernel, whereas previously they were in a Package. So they show up in the result of this query.

On p. 55/56, you always get different timings depending on the version and your hardware. Note that the new numerical features do NOT alter the importance of remembering stored values, as in the Fibonacci sequence. 

Again, on p. 59-60, the fact that the transform packages are now in the kernel renders part of this discussion obsolete. The functions themselves behave in essentially the same way, but are now normalized diferently.

On p.64 the horizontal plot range now inherits the length of the longest list (c).

Finally, this is not a correction but an observation that some of the algebra tools support assumptions. For example:

Simplify[Sqrt[x^2]]

V4 now supports

Simplify[Sqrt[x^2],  x  ∈  Reals]

Chapter 3

No known V4 issues.

Chapter 4

There are only two real issues in this Chapter. The first concerns the conventions that apply to Fourier transforms, which have changed. This is always a pain, with mathematicians, physicists and engineers all having their own pet views! The second is the implication of the format change for integration output, as noted in the Main Points note above.

First though, on p. 95, two results come out in slightly different form

h[z]  /.  DSolve[
          D[h[z],  {z,  2}]  +  z/2 D[h[z],  z]  + 1/2 h[z]  == 0,  
          h[z],  z][[1]]

The second superficial change is rather irritating - we now have

Integrate[
      A/Sqrt[τ] Exp[-x^2/(4 τ)],  {x,  -Infinity,  
        Infinity}]

instead of the preferable V3 answer in the printed text.

On p. 96, it is no longer necessary to load the Laplace Transform package as it is already in the kernel. The Fourier transforms are also preloaded, and the conventions have changed, we now have, by default:

FourierTransform[Exp[-x^2/(4*t)]/Sqrt[4 Pi t],  x,  
      w]

If you want the results to all look like the V3 defaults in the book, you need to use the
FourierParameters options thus:

FourierTransform[Exp[-x^2/(4*t)]/Sqrt[4 Pi t],  x,  w,  
      FourierParameters  ->  {1,  1}]

Let's go back, just to check:

InverseFourierTransform[Exp[-t w^2],  w,  
      x, 
FourierParameters  ->  {1,  
          1}]

There are other changes, such as the abolition of FourierOverallConstant - these things do not materially affect any of the financial derivatives analysis, but you should watch out in you regularly use these things. In some ways the new default normalization, with the root 2 Pi stuff arranged symmetrically, is nicer, but the fact that things have changed is a real nuisance. In V3 the default forward FT was straight integration:

Integrate[
      Exp[-x^2/(4*t)]/Sqrt[4 Pi t] Exp[
          I w x],  {x,  -Infinity,  Infinity}]

The comment on p 99 about FourierOverallConstant being one no longer applies. The FourierParameters setting given above causes all V4 FTs to look like those in the book.

Regarding p. 102, see the errata also, as this interacts with V3-V4 issue. I do not yet have a clean way of cleaning up this integration under V4.

Matters on p. 103 come out more straightforwardly, but no implications. On p. 104, the change in transform operation has some minor but irritating effects, due to the convention shift. The comparison of direct integration and the use of a transform is now best done with

FourierTransform[-2 α UnitStep[-x] Exp[α x] ,  x,  ω,  FourierParameters  ->  {1,  -1}]

Integrate[-2 α Exp[
          x (α  -  I ω)],  {x,  -Infinity,  
        0}]

There are several places where overall multiplicative constants within an integral have now been brought outside the answer. This causes problems when the integration results in an If statement. Subsequent manipulation with %[[2]] to exrtact the result of the integration then do not work the way it was intended. For example, on p.104, the following calculations no longer work.

Clear[α]

Integrate[-2 α Exp[α p] 1/Sqrt[4 Pi τ] 
        Exp[-(p - x)^2/(4 τ)],  {p,  -Infinity,  
        0}]

As noted in the Main Points section, the following leads to rubbish...

GG[x_,  τ_,  α_]  =  %[[2]]

We use instead,

GG[x_,  τ_,  α_]  = %%  /. fix

with, as before

fix  =  If[_,  expr_,  _]  ->  expr;

All very logical and correct. At the top of p. 106 the preceding sequence of manipulations gives now

Series[%,  {α,  Infinity,  1}]

On p. 102 the integration formatting problem raises its ugly head again. The sequence beginning with

value  =  
      Exp[-r (T - t)] B K^λ Integrate[
            1/(Σ Sqrt[
                      2 Pi]) Exp[λ y  -  (y  -  
                          m)^2/(2 Σ^2)],  {y,  0,  
              Infinity}][[2]]

and then using PowerExpand, which is now useless, needs to be replaced by corresponding fixed integrations of the form

value  =  
      Exp[-r (T - t)] B K^λ Integrate[
            1/(Σ Sqrt[
                      2 Pi]) Exp[λ y  -  (y  -  
                          m)^2/(2 Σ^2)],  {y,  0,  
              Infinity}]  /. fix

There are more changes to do with integration and Fourier conventions starting on p. 115.

The remainder of Chapter 4 requires a great many adjustments to the outputs of Integrate to get around the fact that [[2]] no longer reliably extracts the relevant part of an If-based answer. These adjustments have no financial implications - just replace the [[2]] with /. fix.

Chapter 5

No known V4 issues. Some standard form outputs will now use double-struck characters.

Chapter 6

No significant V4 issues, but just an observation about a point of slight disappointment. The support for continued fractions (in a Package previously, but brought to my attention by its relocation to the kernel) still does not allow one to write the Normal distribution any more elegantly, due to the constraint that the running numerator is unity, whereas the running numerator here is a simple, but non-constant sequence. Consider the following

?FromContinuedFraction

expr  =  a  + 1/(b + 1/(c + 1/d)) ;

TraditionalForm[expr]

Simplify[
      expr  -  FromContinuedFraction[{a,  b,  c,  d}]]

The first approximation, given on p. 148, is therefore equivalent to (after repeated division to get the unit running numerator, and just showing the first few terms)

ContinuedFractionApprox[x_]  :=  
      FromContinuedFraction[{0,  x,  x,  x/2,  2 x/3,  
          3 x/(2*4),  4*2*x/(3*5)}]

so precomputation is necessary just to make a list to feed into FromContinuedFraction, which is disappointing, and not to my mind worth pursuing.

Chapter 7

This all works without modification.

Chapter 8

No significant V4 issues here. On p.186, the plot of machine-precision rumblings of Gamma now looks different, which is OK. On p. 187, as there is no normalization stuff left outside the integral, the integration for UpAndOutPutVega and Gamma use of [[2]] still works, but use of the /. fix construction is now preferable, as explained in the Main Points section.

Chapter 9

This all works without modification.

Chapter 10

First of all, see also the "Recent Papers" section for comments on the Geman-Yor implementation of continuously and arithmetically averaged Asian options. Please note the comment on peculiar parameters. In particular, low vols need a much longer Bromwich contour.

On p. 207, it is recommended that you change

G[p_, μ_,   ν_,  α_]  =  F[
          p, μ,  ν, α][[2]]

to

G[p_, μ_,   ν_,  α_]  =  
      F[p, μ,  ν, α] /.  fix

in line with the new general management of integration output. Next, on p. 209, the penultimate line of code needs to be changed to

Apart[
      Residue[U[p, ν, ν, α],  {p, 0}]  /.  
        Gamma[1/2 (4 + 2  ν)] -> (1 + ν)*
            Gamma[1 + ν]]

as a result of the numerator coming out formatted differently.

Chapter 11

This all works with V4.0.

Chapter 12

No V4 issues in this Chapter.

Chapter 13

No version 4 problems known.

Chapter 14

We need to initialize one more variable - with the existing code, we get

\(\(ExplicitSolver\  = 
	
        Compile[{{initial, _Real, 1}, \ {lower, _Real, 1}, \ {upper, _Real, 
              1}, alpha}, 
		Module[{wold = initial, \ wnew = initial, \ 
              wvold, m, k, \ tsize\  = \ Length[lower], 
              xsize = Length[initial]}, 
			For[m = 2, \ 
              m\  <= \ tsize, \ \(m++\), 
				\((wvold\  = \ wold; \ 
                wold = wnew; 
					For[k = 2, \ 
                  k\  < \ 
                    xsize, \ \(k++\), 
					\((wnew[\([k]\)]\  = \ 
                      alpha\ \((wold[\([k - 1]\)] + 
                              wold[\([k + 1]\)])\) + \((1 - 
                              2  alpha)\)\ wold[\([k]\)])\)]; 
						wnew[\([1]\)]\  = \ lower[\([m]\)]; 
						wnew[\([xsize]\)]\  = \ 
                  upper[\([m]\)])\)]; 
			{wvold, wold, 
              wnew}]
		];\)\)

So initialize wvold as requested in the following, which works fine:

\(\(ExplicitSolver\  = 
	
        Compile[{{initial, _Real, 1}, \ {lower, _Real, 1}, \ {upper, _Real, 
              1}, alpha}, 
		Module[{wold = initial, \ wnew = initial, \ 
              wvold = initial, m, k, \ tsize\  = \ Length[lower], 
              xsize = Length[initial]}, 
			For[m = 2, \ 
              m\  <= \ tsize, \ \(m++\), 
				\((wvold\  = \ wold; \ 
                wold = wnew; 
					For[k = 2, \ 
                  k\  < \ 
                    xsize, \ \(k++\), 
					\((wnew[\([k]\)]\  = \ 
                      alpha\ \((wold[\([k - 1]\)] + 
                              wold[\([k + 1]\)])\) + \((1 - 
                              2  alpha)\)\ wold[\([k]\)])\)]; 
						wnew[\([1]\)]\  = \ lower[\([m]\)]; 
						wnew[\([xsize]\)]\  = \ 
                  upper[\([m]\)])\)]; 
			{wvold, wold, 
              wnew}]
		];\)\)

The timing is now improved, without any further changes - we have in V4, on a 266MHz G3 Power Mac

Timing[soln  =  
          ExplicitSolver[initial, lower, upper,  alpha];]

In V3, with either version, we had:

Timing[soln  =  
          ExplicitSolver[initial, lower, upper,  alpha];]

So that's a handy factor of two, without any attempt to figure out the implications of packed arrays. Everything else in Section 14.5, and indeed, the rest of the Chapter then proceeds as before, only rather faster.

Chapter 15

Once the change to ExplicitSolver inherited from Chapter 14 is made this proceeds as before, only faster. The tridiagonal solver compiles as before.

Chapter 16

No V4 issues here - it all works. The integrator continues to complain about underflow when evaluating the CJM integrals, but the analysis is OK.

Chapter 17

No V4 issues here - it all works, but the naive LP solution appears to be slower (to be investigated).

Chapter 18

No significant V4 issues here. The outputs will look slightly different due to the use of double-struck characters in Standard Form output. On p.353 of Section 18.3 the Series function now supplies extra terms in the expansion of p and q. 

On p. 359 the calculation of provisionalsoln returns just two solution sets, the second of which is still the correct one, rather tidied up now.

Chapter 19

No V4 issues.

Chapter 20

This just needs the by-now familiar initialization adjustments to code that is fed to compile. Specifically, on p. 388, we now use

\(\(ManyTimeStepEuro\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {n, \ _Integer, \ 0}}, 

          Module[{new\  = \ initial, \ k, \ old\  = \ initial, \ 
              vold\  = \ initial}, 

            For[k\  = \ 1, \ 
              k <= n - 1, \ \(k++\), \ 
\((old\  = \ new; \ 
                new\  = \ OneTimeStepEuro[old, p, q])\)]; 
vold\  = \ old; \ 
            old\  = \ new; 

            new\  = \ \ OneTimeStepEuro[old, p, q]; 
{Join[new, old], \ 
              vold}]];\)\)

On pp. 391-392, there are similar changes for the American style, and we give the four explicit modifications for Put and Call, and procedural and functional variants. It is sufficient to simply add the initialization of the local Module variables old and vold.

\(\(ManyTimeStepAmerCall\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {S, \ _Real, \ 0}, \ {u, \ _Real, \ 0}, \ {d, \ _Real, \ 
              0}, {n, \ _Integer, \ 0}, \ {K, \ _Real, \ 0}}, 

          Module[{new\  = \ initial, \ k, \ r, \ old = \ initial, \ 
              vold = \ initial, \ payoff, \ finpay\  = \ {Max[S - K, 0]}}, 

            For[k\  = \ 1, \ 
              k <= n - 1, \ \(k++\), \ 
\((payoff\  = \ 
                  Table[Max[S\ d^r\ u^\((n - k - r)\) - K, 0], \ {r, \ 0, 
                      n - k}]; 
old\  = \ new; \ 
                new\  = \ OneTimeStepAmer[old, payoff, \ p, q])\)]; 

            vold\  = \ old; \ old\  = \ new; 

            new\  = \ \ OneTimeStepAmer[old, finpay, p, q]; 
{Join[new, 
                old], \ vold}]];\)\)

\(\(ManyTimeStepAmerCallFun\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {S, \ _Real, \ 0}, \ {u, \ _Real, \ 0}, \ {d, \ _Real, \ 
              0}, {n, \ _Integer, \ 0}, \ {K, \ _Real, \ 0}}, 

          Module[{new\  = \ initial, \ k, \ r, \ old = \ initial, \ 
              vold = \ initial, \ payoff, \ finpay\  = \ {Max[S - K, 0]}}, 

            For[k\  = \ 1, \ 
              k <= n - 1, \ \(k++\), \ 
\((payoff\  = \ 
                  Table[Max[S\ d^r\ u^\((n - k - r)\) - K, 0], \ {r, \ 0, 
                      n - k}]; 
\ \ \ old\  = \ new; \ 
	
                new\  = \ OneTimeStepAmerFun[old, payoff, \ p, q])\)]; 

            vold\  = \ old; \ old\  = \ new; 

            new\  = \ \ OneTimeStepAmerFun[old, finpay, p, q]; 
{Join[new, 
                old], \ vold}]];\)\)

\(\(ManyTimeStepAmerPut\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {S, \ _Real, \ 0}, \ {u, \ _Real, \ 0}, \ {d, \ _Real, \ 
              0}, {n, \ _Integer, \ 0}, \ {K, \ _Real, \ 0}}, 

          Module[{new\  = \ initial, \ k, \ r, \ old = \ initial, \ 
              vold = \ initial, \ payoff, \ finpay\  = \ {Max[K - S, 0]}}, 

            For[k\  = \ 1, \ 
              k <= n - 1, \ \(k++\), \ 
\((payoff\  = \ 
                  Table[Max[K - S\ d^r\ u^\((n - k - r)\), 0], \ {r, \ 0, 
                      n - k}]; 
old\  = \ new; \ 
                new\  = \ OneTimeStepAmer[old, payoff, \ p, q])\)]; 

            vold\  = \ old; \ old\  = \ new; 

            new\  = \ \ OneTimeStepAmer[old, finpay, p, q]; 
{Join[new, 
                old], \ vold}]];\)\)

\(\(ManyTimeStepAmerPutFun\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {S, \ _Real, \ 0}, \ {u, \ _Real, \ 0}, \ {d, \ _Real, \ 
              0}, {n, \ _Integer, \ 0}, \ {K, \ _Real, \ 0}}, 

          Module[{new\  = \ initial, \ k, \ r, \ old = \ initial, \ 
              vold = \ initial, \ payoff, \ finpay\  = \ {Max[K - S, 0]}}, 

            For[k\  = \ 1, \ 
              k <= n - 1, \ \(k++\), \ 
\((payoff\  = \ 
                  Table[Max[K - S\ d^r\ u^\((n - k - r)\), 0], \ {r, \ 0, 
                      n - k}]; 
old\  = \ new; \ 
                new\  = \ OneTimeStepAmerFun[old, payoff, \ p, q])\)]; 

            vold\  = \ old; \ old\  = \ new; 

            new\  = \ \ OneTimeStepAmerFun[old, finpay, p, q]; 
{Join[new, 
                old], \ vold}]];\)\)

When we come to the trinomial model, we use

\(\(ManyTimeStepEuroTri\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {r, \ _Real, \ 0}, {n, \ _Integer, \ 0}}, 

          Module[{new\  = \ initial, \ k, \ old\  = \ initial}, 

            For[k\  = \ 1, \ 
              k <= n, \ \(k++\), \ 
\((old\  = \ new; \ 
                new\  = \ OneTimeStepEuroTri[old, p, q, r])\)]; 

            Join[new, old]]];\)\)

Chapter 21

The timings of Convert are better under V4, but you should still avoid this method! The function norb on p. 404 still takes 0.25 seconds on the same machine as before.

The change in the compiler requires several changes in this Chapter. As noted in the summary section, the definition on p. 404 gives, now,

\(\(norm = \ 
        Compile[{mu, \ 
            sigma}, \ 
		Module[{va, \ vb, \ rad\  = \ 2.0, \ 
              den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; 
					vb\  = \ 
                  2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			mu + 
              sigma*va*den]];\)\)

and we replace it by:

\(\(norm = \ 
        Compile[{mu, \ 
            sigma}, \ 
		Module[{va = 1. , \ vb, \ rad\  = \ 2.0, \ 
              den = 1. }, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; 
					vb\  = \ 
                  2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			mu + 
              sigma*va*den]];\)\)

A similar fix applies to normpair: before:

\(\(normpair = \ 
        Compile[{mu, \ sigma}, \ 
          Module[{va, \ vb, \ rad\  = \ 2.0, \ den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; \ 
                vb\  = \ 2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			{mu + sigma*va*den, \ 
              mu + sigma*vb*den}]];\)\)

This time we need to initialize both va and vb:

\(\(normpair = \ 
        Compile[{mu, \ sigma}, \ 
          Module[{va = 1. , \ vb = 1. , \ rad\  = \ 2.0, \ 
              den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; \ 
                vb\  = \ 2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			{mu + sigma*va*den, \ 
              mu + sigma*vb*den}]];\)\)

Antithetic Variations

The initialized versions of these to be used with V4 are just:

\(\(normat = \ 
        Compile[{mu, \ 
            sigma}, \ 
		Module[{va = 1.0, \ vb = 1.0, \ rad\  = \ 2.0, \ 
              den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; 
					vb\  = \ 
                  2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			{mu + sigma*va*den, \ 
              mu - sigma*va*den}]];\)\)

\(\(normpairat = \ 
        Compile[{mu, \ 
            sigma}, \ 
		Module[{va = 1.0, \ vb = 1.0, \ rad\  = \ 2.0, \ 
              den}, 
			While[
              rad\  >= 1.00, 
					\((va\  = \ 
                  2.0*Random[]\  - \ 1.0; 
					vb\  = \ 
                  2.0*Random[]\  - 1.0; 
					rad\  = \ 
                  va*va + vb*vb)\)]; 
			den\  = \ 
              Sqrt[\(-2.0\)*Log[rad]/rad]; 
			{mu + sigma*va*den, \ 
              mu - sigma*va*den, \ mu + sigma*vb*den, \ 
              mu - sigma*vb*den}]];\)\)

For reasons that are not yet understood, in Section 21.5, the EuroCallB (and put) versions now run slightly slower than the type A version, under V4. This will be investigated. (As noted in the introduction, this notebook is about making sure code works under V4 - fully optimized material is under consideration).

Chapter 22

The modification to normpairat made in Chapter 21 must also be made here. There are no other V4 issues.

Chapter 23

No V4 issues here - it all works.

Chapter 24

Use the new compiled version of the norm function discussed previously.

\(\(norm = \ 
        Compile[{mu, \ sigma}, \ 

          Module[{va\  = \ 1. , \ vb, \ rad\  = \ 2.0, \ den}, 
            While[rad\  >= 1.00, 
			\((va\  = \ 2.0*Random[]\  - \ 1.0; 
                vb\  = \ 2.0*Random[]\  - 1.0; 
			rad\  = \ 
                  va*va + vb*vb)\)]; 
            den\  = \ Sqrt[\(-2.0\)*Log[rad]/rad]; 
			mu + 
              sigma*va*den]];\)\)

Then everything works.

Chapter 25

The numerical factors of 2 and 1/2 arrange themselves differently under V4 in the construction of equations (24) and (25) by the use of Traditional Form applied to aeqn and beqn. But there is no logical change. At the top of p. 477, the post-application of [[1]] is now un-necessary. Under V4, there are various re-organizations of numerical factors in outputs throughout this Chapter that are of no consequence.

Chapter 26

See the Main Points section for a discussion of NonlinearFit. On p. 486, the following code fails

NonlinearFit[exptdata, 
      a  +  b/((x - c)^2  +  d)  +  
        e/((x - f)^2  +  g),  
      x,  {a, b, {c, 11}, d,  e,  {f, 14},  g},  
      MaxIterations  ->  40]

It does run, with no suggestions for starting values, but gives answers that are rubbish.

NonlinearFit[exptdata, 
      a  +  b/((x - c)^2  +  d)  +  
        e/((x - f)^2  +  g),  x,  {a, b, c, d,  e,  f,  g},  
      MaxIterations  ->  40]

The only way to get anywhere is to guess all the parameters - the following does work:

NonlinearFit[exptdata, 
      a  +  b/((x - c)^2  +  d)  +  
        e/((x - f)^2  +  g),  
      x,  {{a, 1}, {b, 1}, {c, 11}, {d, 1},  {e, 1},  {f,  14},  {g, 
          1}},  MaxIterations  ->  40]

No other V4 issues noted in Chapter 26.

Chapter 27

No V4 issues.

Chapter 28

No V4 issues. The series solution still takes an eternity.