An Introduction to Mathematica

A System for Doing Mathematics by Computer
(based on Ch 2 of revised MFDwM)
including data acquisition and analysis; parallel methods.

2.1 Systems, Kernels and Front Ends

The first thing that you need to do is to learn how to communicate with the computer and the Mathematica software system. It is important to realize that the software is in two pieces. The calculations are done by a computational "kernel", but you communicate with the kernel through a front end. The kernel and the front end may be running on the same computer, or on different computers linked by a network. This introduction is targeted at people using a windowing operating system, and the front end is a Notebook interface. This interface creates and uses Mathematica Notebook files - you can exchange Notebooks between computer systems such as MacOS, Windows XP/NT and UNIX. The Notebooks are just text files, and can be passed between systems using only text-translation utilities. 

If you have not used Mathematica before, how to proceed depends on whether you have Mathematica already installed on your computer. If it is not installed at all, you may need some help if you are working in a multi-computer networked environment running under UNIX, for example - under these circumstances, it is almost certainly best to request help from your System Administrator. If you are running on a solitary system, or otherwise doing an installation local to your computer, it is usually sufficient to follow the instructions in the Mathematica installation instructions. These are contained in the booklet "Getting started with Mathematica on Macintosh/Windows (etc.) systems", specific to your own computer. If you have never used Mathematica before, it will be helpful to run through the examples contained in that booklet. The discussion below is designed to complement that booklet, but is essentially self-contained.

Begin by starting Mathematica. You can do this on any system by navigating to the icon for Mathematica, and double-clicking on it. Do not use the "MathKernel" program icon for this. On particular systems there will be other routes to access Mathematica. Under Windows, the program can be found under the "Start" task. On both Windows  and MacOS, it is convenient to make short-cuts/aliases and place them on your desktop - it is really up to you to customize your route to access the system once you have figured out what is most convenient. Mathematica will start up, giving the version number, your licence information, and warnings about software piracy.

Input and Output

When a new Notebook is opened you can type straight into the blank window. Type "3 + 3" into the system (not the quotes, just the expression within the quotes). To enter this expression hold down the shift key and hit the "return" key.  Mathematica labels both your input and its output with numbers, in this case [1], if you are starting a fresh session.  You should see the following, with the input prefixed by In[1], the output prefixed by Out[1]:

3+3

Now type the command  "Expand[(2 x + y)^3]", and enter this by SHIFT-RETURN (ENTER, not to be confused with RETURN, and usually located at the bottom right of keyboards, may also work and is a useful short-cut on some computers).

Expand[(2 x + y)^3]

Throughout this chapter, whenever we present examples like this, expressions in bold type are what you should type, followed by SHIFT-RETURN or ENTER.  To see some of the capabilities of Mathematica try the following sequences of statements.  When you enter the commands in bold, they will acquire  "In[n] := " labels.  Also, if you make a mistake, just re-type the correct expression. We will come back to how to edit expressions shortly; the following are all short enough to make re-typing painless.

The first example uses the N function to ask for π, denoted in Mathematica's Input Form by Pi, to 100-digit precision.

N[Pi, 100]

The next example introduces the plotting capabilities of Mathematica by using Plot.  Given the function and the range over which to draw the graph, the Plot function will produce the appropriate graph.

Plot[x^2 Sin[x], {x, -3, 3}]

In this simple program we define the function f that gives the fourth power of its argument and then apply f twice to y.  Use an ordinary RETURN to get to the second line for typing, and then SHIFT-RETURN to execute the expression. Note the colon and semicolon which we shall examine later.

f[x_] := x^4;  
f[f[y]]

The final example shows the incremental capabilities of Mathematica by using the Table function to produce a list of the cubes of integers up to 10.

cubes = Table[ i^3, {i, 0, 10}]

We can plot given data like the list cubes using the function ListPlot.

ListPlot[cubes]

Don't worry about the details of the positioning of the axes for now. But you should note that if asked to plot a simple one-dimensional list, Mathematica will use the position in the list as the "x" or horizontal coordinate, so that, e.g., the last element is plotted as the point (11, 1000).

Remarks

In our examples we have seen the four basic areas of Mathematica. These are:

Numerical Calculations

For example, we asked for the numerical value of π to 100 digits.

Symbolic or Algebraic Calculations

We asked for the expansion of (((2x+y))^3).

Graphics

We plotted a function, and a sequence of numbers.

Programming

We defined a function, and asked for it to be applied to a specific argument.

We will explore each of these topics in more detail later. Right now it is necessary to do some groundwork on understanding how Mathematica works. From here on, the expressions that we enter may be a little longer, and if we make a mistake we don't want to have to keep re-typing our input. So at this point it will be useful to learn a few more tricks for managing our input.

Using the Cursor

If you move the cursor up and down the Notebook in which you have just done the previous calculations you will see how it changes form dependent on which type of cell it is lying in, or whether it is between cells. A horizontal line under a cell indicates that any text typed will be placed after that cell. To place a horizontal line after any cell you move the cursor so that it indicates that it is lying between cells, i.e. when it is horizontal, and then click. Practice will get you used to being able to place text anywhere in a Notebook.

Opening and Closing Cells

In the Notebook environment all text is held in cells. These can have different styles associated with them, and can be set up to allow different styles of cell to be nested within one another. By double clicking on any cell bracket which is grouping other cells you will be able to hide the text contained in those sub-cells. This can be retrieved by double-clicking on the same cell bracket again. When a cell contains hidden sub-cells you will notice how the lower edge of the bracket changes to a filled triangle. You may find this way of hiding information in the Notebooks useful in organizing a presentation or teaching session. It may take some practice to get used to how these cells are nested.

Editing Text in Cells

Changing the text in a cell is done in a similar way to that used in word-processing applications. The cursor can be placed at any particular point in the text by moving the mouse and clicking on that spot. Selecting text with the mouse allows that text to be removed or typed over. Notice how you can only edit text in certain cells. You will be able to edit your input cells but not in the cells containing the output. This has been set by the way in which the different cells are defined. Here are some handy tips:

(1) use the mouse and scroll bars to get up and down the Notebook;
(2) use an ordinary RETURN to enter the next line of a multi-line expression;
(3) use SHIFT-RETURN (or ENTER) to execute an expression;
(4) use the mouse COPY and PASTE to grab previous formulae.

Quitting and Panicking

Sometimes, say if you make a mistake, Mathematica may go into a long internal session and all you want it to do is stop. If you need to stop Mathematica in the middle of a calculation, go to the menu and choose "Kernel". You then get some options, of which the ones to pick are "Abort Evaluation" or "Interrupt Evaluation...".  The menus also offer keyboard short-cuts specific to your operating system, which you should learn. The kernel can sometimes take time to abort a calculation. It is better to wait for it to finish sorting itself out than to keep thumping the keyboard! If you lose patience waiting, or otherwise wish to exit Mathematica completely, you can click on "Quit" having interrupted the session, or type Quit at an "In" prompt. There are other options in the "Kernel" menu for quitting and re-starting.

Brackets Galore

Our first examples show that we need to understand how Mathematica reads expressions, and especially different types of brackets: [ ], { }, ( ). Sometimes people use square brackets [], round brackets () and curly brackets {} more or less interchangeably. In order to ensure that the computer understands unambiguously what we mean, it is necessary to introduce some conventions.

Arguments of Functions

Arguments of functions are always enclosed in square brackets. This applies just as much to standard mathematical functions like Sin[x] as to built-in Mathematica graphics functions like Plot3D[], which take another function and some ranges as a minimum set of arguments, and user-defined functions like f above.

Grouping

When we wish to group together some expressions for applying a subsequent operation, round brackets are employed. For example (x + y)^3 means the cube of all of x + y.  Note that this distinction gives the expressions f[x+y] and f(x+y) two quite different meanings, and at the same time we have acquired the freedom to use a space between adjacent quantities to denote multiplication. If you wish to you can explicitly type a * to denote multiplication. For example 2 3 is identical in meaning to 2*3.

Note that you don't always need to use round brackets for grouping, since the usual rules of precedence apply. For example x^3 + 4 means (x^3) + 4 and not x^(3 + 4). A good policy is to use grouping brackets in any situation where there is a potential for ambiguity.

List Brackets

Curly brackets { } and double square brackets [[ ]] are also used in Mathematica. These are used with list structures to denote a list, and an indexed member of a list, respectively. We shall look specifically at these structures later on. For now it will suffice to appreciate that a list consists of a number of elements, separated by commas, with a curly bracket at each end.

Input/Output, Standard and Traditional Forms

Older versions of Mathematica prior to version 3.0 used styles where inputs and outputs were both based on pure text representations. Mathematica nwo routinely uses forms called "Standard Form" and "Traditional Form". These are much more like traditional mathematical notation - indeed, this is exactly the idea of "Traditional Form". "Standard Form" is a half-way form that goes as far as possible towards traditional mathematical notation without introducing any ambiguities. Throughout this chapter the default output format type is "Standard". You can ask for explicit conversion as follows:

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

TraditionalForm[expression]

InputForm[expression]

These conversions can be done using the menus, and defaults can also be set. It is impossible to give a hard and fast rule that has been applied throughout this book, for how expressions are represented, as the choice is rather sensitive to the context, but we try to adhere to the "WYSIWMAMA"  (what you see is what Mathematica made) principle - note the partially double-struck notation. Generally, chapters involving more mathematical development will make more use of the traditional form for outputs and equations.

2.2 Basic Graphics

Before we launch into investigating each of the main areas of Mathematica in detail, it will be useful to look at some of the basic graphic commands used in Mathematica. While working in Mathematica it is very often helpful to involve graphics in a calculation. This is never more true than when looking at new ideas and will be used extensively through this chapter and the rest of the book, to illustrate points more clearly.  For this reason, it will be useful if you feel reasonably comfortable with the most basic graphics commands at the outset. Do not be too concerned about the details of their operation at this point. We shall look more carefully at the range of graphics capabilities available in Mathematica later.

Producing a Plot

We have already seen the function Plot which will display the plot for a given function. This needs to have the range of values over which the plot is to be made specified. You supply Plot with two arguments - the first is the function, and the second is a list in the form {variable, start, end}.

Plot[x^3 - 8 x^2 + 19 x - 8, {x, 0, 5}]

Given a set of values expressed as a list data, the function ListPlot will plot data assuming that data is regularly spaced, and will use the position in the list as the horizontal coordinate.  The command which changes the PlotStyle plots the points large so that we can see them clearly.

data = {3,-0.5,2,0,1,0.5};
points = ListPlot[data,PlotStyle -> PointSize[0.03]]

Displaying a Plot Again

All graphics can be displayed again using the function Show. We have given our last plot the name points.  We can now ask to see that particular plot again.

Show[points]

All of the graphics displayed can be altered to suit particular needs.  In our plot for ListPlot we specified a particular size for the points to be plotted. Similarly, we can ask for those points to be joined by a continuous line.

pointsb=ListPlot[data,Joined->True]

If you have executed the previous few commands on a computer you may have received a spelling warning. Here and elsewhere in this text we have suppressed messages from Mathematica warning us about possible spelling mistakes suggested by a new variable name being very much like one previously introduced. The Show command can be used to combine plots:

Show[points,pointsb]

2.3 Numerical Calculations

Exact Arithmetic

Degrees of Precision

With most computer languages you are restricted to certain levels of precision in the values of numbers.  This is not the case when you are using Mathematica. Try the following series of inputs.

30!

{N[30!], N[30!, 23]}

The first result is exact, the second is approximate and has used the degree of precision given by default. In the last expression we have specified a precision of 23 digits.

Exact and Approximate Numbers

When we ask for the square root of a number, for example, Mathematica will return the most exact version.

Sqrt[51]

Note that if you are using a much older version of Mathematica, the output here will differ - you will still see the letters "Sqrt" in the output. Only when we ask explicitly for that to be given in a numerical form will Mathematica return a number which approximates this square root.

N[Sqrt[51]]

In evaluating the last expression Mathematica displays the answer to six significant figures. It will have been calculated, in  its internal numerical form, to a higher precision, which will be dependent on the machine you are using.  We can find the degree of precision by asking for the value of Precision on that particular output line,

Precision[%]

So in this case, the result will have been calculated to 16 SF, but only displayed to six figures.  If that result were to be used in any further calculations, the full 16 figures would be used.

How Mathematica Decides on Exact or Approximate

This following example highlights the different routes Mathematica will take when dealing with numerical expressions.

(3/4)^100

In this expression the answer has been returned in the most exact form possible.  Mathematica will keep to that form unless explicitly asked otherwise, by the use of N[ ], for example.  An alternative way of presenting the above calculation is

0.75^100

This second expression has numbers given in their floating point or approximate form, the expression is then kept in this approximate form and the answer given in that form. Any number with a decimal point in it will be taken to be approximate.  Of course, we can always ask for that approximation to be more precise:

N[%,10]

A Practical Example with Exact and Approximate Numbers

For another illustration of these two paths down which numbers may be sent, we can look at solutions to equations. Giving the equation in exact form:

Solve[{x^2 + 3 x + 11/10 == 0}, x]

N[%]

Alternatively, in an approximate form:

Solve[{x^2 + 3 x + 1.1 == 0}, x]

Again, we chose the way in which numbers would be treated by the way we presented the equation - the presence of a decimal point implied the numerical treatment. This example raises some other issues that it is a good idea to explore right away.

A Note on "Equals"

You will have noticed that in the last section we used a double equals sign to denote that an expression was equal to zero.  The use of the single equals sign will be examined more carefully shortly, where its role in assignment will be illustrated. For the moment it will be sufficient to highlight that when defining a symbolic equation the double equals sign should be used. This does not involve assigning any particular values to any variables.

x == 0

If actually setting a quantity equal to some value for future substitution the single equals sign is used.  This assigns that value to the variable.

x = 0

(Having set this value for x we must clear it for future calculations)

Clear[x];
x

Extracting Results from the Solve Function

It can be frustrating for new users that Mathematica returns results from Solve as some replacement rules indicated by the arrow:

Solve[{x^2 + 3 x + 11/10 == 0}, x]

What we really want are a list of values of x given these rules. In Mathematica the term "given" is invoked by the use of "/.", as shown here:

x/. %

You can do the whole thing at once:

x/. Solve[{x^2 + 3 x + 11/10 == 0}, x]

Some Useful Numerical Functions

The usual arithmetic functions are as you would expect for +, -, *, / etc. Note that multiplication is denoted by a star or a space. A cross or "times" symbol is not used. The square root of a number is calculated by the Sqrt function.

Sqrt[78543.]

Similar functions exist for exponentials, logarithms and many other numerical functions. You should note that sometimes long words are used, and there may be capitals within the name. For example, The function FactorInteger produces a list of the prime factors of a number and the power to which they are raised - useful in cryptography!

FactorInteger[15931157760]

This is telling us that 15931157760=(2^8)× (3^8)× 5× 7× 271.

From Data to Functions

Simple Function Fitting

Suppose that we have a data set and wish to represent it in terms of functions by a least-squares fit. This can be accomplished in a very simple way in Mathematica. Let's invent some data:

fitstuff = {6, 3, 5, 8, 7}

The elements of fitstuff are interpreted as being the values of a function at equally spaced value points 1, 2, 3, 4, 5.  Now we can plot those points.  Calling the plot dataplot will allow us to use it again later.

dataplot = ListPlot[fitstuff, PlotStyle-> PointSize[0.03]]

Now type the following, exactly as written. Do not be too concerned about the syntax used to define the function psi, we will investigate this later - you are most concerned about the function Fit.

psi[x_] = Fit[fitstuff, {1, x, x^2, x^3, x^4}, x]

psiplot = Plot[psi[x], {x, 1, 5}]

Now we can show the original points of the data set as well as the fitting function, and see that the function passes through them. At the same time we can see how to override Mathematica's conventions on the vertical plot range:

Show[psiplot, dataplot, PlotRange -> {0, 10}]

Fitting Functions to an Irregular Array

This approach can be applied to lists of pairs of (x, y) values also, and any set of functions can be used to construct the function, provided they appear with unknown coefficients in a linear combination. For example for the set of data tStuff, we can fit a trigonometric function to the data.

tStuff={{0,2.8},{2,11.},{3,5.3},{5,-4.4},{8,11.5}};

mi[x_] = Fit[tStuff, {1, Sin[x], Cos[x]}, x]

miPlot = Plot[mi[x],{x,0,8}]

tStuffPlot = ListPlot[tStuff,PlotStyle->PointSize[0.03]]

Show[miPlot,tStuffPlot]

You should appreciate that using this ordinary Fit function you can fit any model that is linear in the unknown parameters. The basis functions can be anything you like.

Other Numerical Operations

Numerical Solution of Equations

We have already seen how to solve a quadratic. Mathematica can cope with all polynomials up to degree 4 analytically. For degree 5 and above, usually a numerical approach is required. This involves the NSolve function.

NSolve[5 x^5 - 15 x^3 + 10 x + 1 == 0, x]

NSolve and its symbolic counterpart can be used to solve sets of equations also:

NSolve[{y + 3 x == x^3, 4 y^2 - x == 4}, {x, y}]

To obtain the results as a list of (x,y) pairs, you use /. for given, as before:

{x, y} /. %

For functions other than linear and polynomial systems, the FindRoot function is better. You give it a starting point and it iterates from there to a root:

FindRoot[Cos[x] == 0, {x, 1.5}]

The result is a numerical approximation to π/2:

N[Pi/2]

However, this function should be used with some care.  (The function FindRoot operates using Newton's method which, if you have experience with its use, you will know is most effective when supplied with an approximation very near to a root.)  For example if we supply a first approximation of 3.1, the root finally given is nowhere near 3.1.

FindRoot[Cos[x] == 0, {x, 3.1}]

The tangent to the curve y = cos(x), at x = 3.1, is nearly horizontal, so that the first iteration takes us to a point far away from the starting value. The process then converges! In financial applications FindRoot is particularly useful for solving non-linear equations, such as arise in calculating implied volatility, solving for an internal rate of return, or obtaining a critical price defined by a non-linear equation. We shall see its use many times, and consider in more detail later.

2.4 Lists

There are some recurrent patterns in Mathematica of which you will need to be aware, no matter what the type of calculation you are doing. The use of list structures in Mathematica occurs in many diverse areas. We have already used them without too much fuss - now is a good time to go into a little more detail.

One-Dimensional Lists

Lists are used when dealing with data, or tables of values.  A list is denoted by using curly brackets to enclose the elements of the list, which are separated from each other by commas. We have already used some lists in our examples.  An example of a simple list is

data = {0,1,4,9,16,25,36,49,64,81,100}

We can select individual elements of such a list, or index elements of the list using the number which shows their position in the list. This involves the fourth type of bracketing, that of double square brackets.  For example we can select the fourth element of the list data:

data[[4]]

There are built-in functions for re-arranging the elements of lists including functions such as Reverse, which re-orders the list in the opposite direction, and RotateRight, which takes the last element of the list and brings it to the front.

Reverse[data]

RotateRight[data]

The functions Prepend and Append can be used to add elements to lists.

Prepend[data, 1001]

Append[data, 25000]

Note that neither of these changes the meaning of data:

data

The functions PrependTo and AppendTo have the same effect as Prepend and Append, but also update their argument:

PrependTo[data, 1001]

data

We select particular parts using functions, First, Rest and Last.

First[data]

Rest[data]

Last[data]

Our lists may well be generated within Mathematica by use of the Table command (examined in more detail later), or they may be generated outside Mathematica, imported and then manipulated, possibly to obtain  a graphic representation.

Arrays, Matrices and Lists of Lists

Next up from a list is the concept of a list of lists. We use this idea in several ways. Here is an example of a list of lists:

listb = {{1,2,1,2}, {4,3,4,3}, {7,8,9,6}}

We have already seen how the function FactorInteger supplies its answer in a list of lists form.  Also, lists of lists are how Mathematica represents matrices, and lists of (x,y) pairs.

MatrixForm[listb]

We can ask about the dimensions of this matrix.

Dimensions[listb]

A list of lists is the way we represent any two-dimensional array of numbers. It could be an array for plotting or a matrix, or a rank-2 tensor. Matrix multiplication is specified by a "dot" (.), to distinguish it from element by element multiplication of lists.  The simple square matrix SqMatrix will illustrate this.

SqMatrix={{1,6},{4,2}}

SqMatrix.SqMatrix

SqMatrix SqMatrix

Selecting out specific rows is done in a similar way to that in which we select individual elements from a list. The second row of SqMatrix is therefore denoted by

SqMatrix[[2]]

Similarly we can select an individual element in a matrix by specifying its position in the matrix by row and then column, just as when referring to an entry in the form (Subscript[A, i j]) :

SqMatrix[[1,2]]

Various other standard matrix operations may be carried out using specific Mathematica functions. For example we can transpose a matrix,

MatrixForm[Transpose[listb]]

or find the inverse of a square matrix. To demonstrate this we first add a row to listb to make a square matrix:

invertible = Prepend[listb, {4, 2, 6, 1}];
inverse = Inverse[invertible];
MatrixForm[inverse]

MatrixForm[inverse.invertible]

Flatten and Partition

One of the manipulations which it may be necessary to perform on a list is the conversion from an n-dimensional list to a one-dimensional list and back again.  This is easily accomplished using Flatten and Partition.

listc = Flatten[listb]

Partition[listc, 4]

Note that the second argument to Partition represents the length of the sub-lists generated.

Applying Functions to Lists

There are two different situations under which you may want to apply a function to a list.  Either the elements of the list are to be used as the arguments to a function, or a function is to be applied to each element of the list individually, the result being another list. These are two similar processes, but should not be confused. For functions which act on a sequence of arguments we use the Apply function to supply the argument set.

data

Apply[Plus,data]

When a function is to be applied to each element of the set we can use the function Map, or for certain numerical functions we can use the whole set as the argument.

Sqrt[data]

Map[Sqrt, data]

In Mathematica you can use lists of arbitrary dimension without making further definitions. So e.g tensors of arbitrary rank may be defined. People doing General Relativity can work with the rank 4 Riemann tensor and its derivatives (rank 5!). In finance you should compare this with the utter chaos that pervades the C++ world, where the basic language does not include proper definitions even of matrices. While that language is powerfully extendable to define anything you like, and indeed many authors of books and software now have invented their own structures for arrays, there are now umpteen Classes for matrices lying around!

Incrementation and Nesting

Incrementation

To perform a simple incremental function there are a few Mathematica functions which can be used.

Do[Print[{n, n^2, n^3}], {n, 1, 4}]

In this example the command Print[{n, n^2, n^3}] is executed for n taking values starting at 1 and being incremented by one each time until n equals 4. The function Table, used to produce lists, is of a similar form.

Table[i^2 +3 i,{i,1,15}]

In this case the resulting list or table is made up of each of the results obtained when the quantity (i^2)+3i is calculated for i running from 1 to 15 in steps of 1.

Iterative Application of Functions

Iteration is at the heart of many modern areas of mathematics including dynamical systems, chaos and fractals. Mathematica has some simple tools for dealing with iterated systems, based on Nest and related functions. Nest in particular is a handy way of simulating random walks. We shall give some examples shortly. Before giving some simpler examples, we make sure that we have removed any previous definitions of several handy symbols:

Clear[x, y, f, g, h, z]

The Nest function applies a given function (its first argument) to its second argument a certain number of times, specified by its last argument:

Nest[g, x, 3]

An interesting related function is NestList.

NestList[Sin, x, 3]

There is also a neat pair of generalizations to functions of two variables:

Fold[h, x, {a, b, c}]

FoldList[h, x, {a, b, c}]

An Example Using Nest

For example, what does Sin[Sin[Sin ... (100 times)[x]...]] look like? Ask Mathematica:

Plot[Evaluate[Nest[Sin, x, 100]], {x, -2 Pi, 2 Pi}]

The Evaluate function has been used to make the plot easy for Mathematica to draw. It prevents Mathematica repeating the Nest operation at every point at which the function is sampled. With Evaluate in place it will only do that once. The example shown, and why it looks like a square wave, are discussed in detail by Gray and Glynn (1990).

A Few Useful Observations

Once or twice we have made use of a couple of features that deserve explicit comment.

Referencing Earlier Results

During a sequence of calculations you will probably want to refer to the results of previous calculations.  This is achieved easily using a per cent sign, %, immediately followed by the number of the output corresponding to the result. Recall the definition of data:

data

A special shorthand for the last computed result is %; %% denotes the next-to-last result.

%

%%[[3]]

Suppressing Data Output

Our list data is short, so there is no disadvantage in having Mathematica output the full result.  If you are working with large data sets, however, you will almost certainly not want to see the results at each stage.  Consequently, it may sometimes be desirable to suppress large results of purely symbolic operations.  In the first of the following operations, the data set bigdata is created without being displayed at all. The semicolon is the key ingredient for suppressing output completely.

bigdata=Table[i^2,{i,0,10000}];

The next example shows a way of producing a shortened form of output.  Remember this is only shortening what is displayed. The actual result is still a list of 10001 terms.

Short[bigdata]

You may get different results depending on the exact version of Mathematica you are using. 

n.b. An understanding of the suppression of the output with a semi-colon is important in Mathematica technologies beyond version 5.2. If you get a more recent version and find that graphics are not appearing, make sure you delete any semi-colons after the graphics command. In version 5.2 or earlier, the graphic would appear even when  a semi-colon was present - this is changing!

2.5 Algebra

Simple Algebraic Manipulation

We'll begin by doing some simple algebra. Enter the following statements executing each one as you go.  You will see that Mathematica has its own ideas about how to order the terms in a polynomial – do not let this worry you.

y^4 + 4 y^3 + 6 y^2

Now, add 4y to this last result:

% + 4 y

Add 1 to this last result, and factorize this:

Factor[% + 1]

Undo the factorization by expanding this last result:

Expand[%]

Assign x to be this last quartic expression in y:

x = %

Find the solutions for the quartic expression to be zero, using x as a shorter version of the quartic expression:

Solve[x == 0, y]

Assign the value 2 to y, and without displaying that result ask for the resulting value of x.

y = 2;
x

Now cancel any values associated with y, and ask for the value of x.

Clear[y];x

The three things you should appreciate from this little session are the use of Expand, Factor and Clear.  Also, we used Solve as encountered in the Numerical Calculations section, to give us solutions to an algebraic equation.  It is a good idea to get into the habit of using Clear frequently, to make sure that you do not cause problems simply by trying to define quantities that conflict with earlier definitions.  You should clear the definitions of x and y now, if you are running this session on a computer.

Clear[x,y]

Simplifying Rational Expressions

It is essential that you should have cleared the values for x and y from the last sequence of calculations, if you are trying out this chapter on a computer. If you have not done this Mathematica will go into an internal loop and will complain about exceeding its recursion depth.  If this happens you will want to abort the calculation, which as you will remember from the introductory section is done by selecting "Interrupt Evaluation..."  or "Abort Evaluation"  from the Kernel menu.

First, we will input this simple rational expression involving x and y.

y = 1/(x^3 - 4) - 1/(x^2 - 2)

Now, we can ask for all of the terms to be gathered together over a common denominator, and this expression to be assigned to be u.

u = Together[y]

Using u, the shorthand for the expression, we can select out the denominator and numerator of the expression.

Denominator[u]

Numerator[u]

Expansion can be aimed only at the denominator of the expression.

ExpandDenominator[u]

With rational, rather than polynomial expressions, it may help to separate out the top and bottom for manipulation.  In most cases you will find Expand, Factor and Simplify sufficient.  You should take care in the use of Simplify. It can sometimes get into quite a lengthy process, and even have unpredictable results.  It is much easier to use more functions but keep control over what changes are taking place.

Applying Simplification Rules Selectively

The following example will give some idea of how to apply operations selectively.  We begin by defining the algebraic expression expr.

expr = (x^3 - 216) (x - 3)^2/((x - 2)(x - 4)(x - 6))

Applying Expand to this expression will expand each term in the numerator.

Expand[expr]

ExpandAll, however, will take every part of the expression.

ExpandAll[expr]

You have already used Together, which will gather all terms over a common denominator.

Together[%]

Splitting this expression into partial fractions uses Apart.

Apart[%]

The Factor function will implicitly gather the terms together and factorize the resulting expression:

Factor[%]

Another route to this simplified version of expr would have been using the function Cancel.  This will make sure that any possible cancellations are carried out.

Cancel[expr]

Other Useful Functions

Another useful pair of operations are Collect and Coefficient.  We shall use another polynomial expression to illustrate their use.

Clear[x, y];
stuff = (2 x + 3)^2 (x + y)^2

Expand[stuff]

Collect[%, x]

So Collect literally collects together terms with the same power of x.  We can pick out the coefficients by using Coefficient.  Here we ask for the coefficient of x^2 (all the terms, including those in y) in the expression.

Coefficient[Expand[stuff], x^2]

This type of operation can be very useful in doing perturbation theory, and can also be used with the Series function which we will come to shortly. An alternative syntax is as follows - this must be used to extract the "constant" term:

Coefficient[Expand[stuff], x, 2]

Coefficient[Expand[stuff], x, 0]

Making Substitutions

Substitutions are a common requirement in manipulation. One way of doing this is to make a simple assignment  and then ask for the new value of an expression.

x=3 y;

stuff

This will mean however that from now on x will always be read as 3y and stuff will always equal 16(((y^2)(3+6y))^2). To prevent this we will have to explicitly remove the assignment of x.

Clear[x]

The technique to use to make a substitution which will not make a permanent change is to use a replacement rule. These are the same as are given as solutions to NSolve and FindRoot. To make the same substitution define the rule that x goes to 3y, and apply that using /. to stuff.

stuff /. x -> 3 y

stuff

You can see that this temporary substitution does not alter the stored definition of stuff. This technique can be used for any temporary substitutions whether they are numerical or algebraic.

2.6 Calculus

We now come to a topic that is fundamental to appreciating how useful Mathematica can be in modelling financial derivatives. We are interested in finding the values of certain options, together with the derivatives of those values with respect to variables such as the underlying price, its volatility, interest rates, and so on. When there is some form of analytical prescription available, it is very convenient, not to say less error-prone, to perform the computation of such "Greeks" symbolically. Even when one is working with an essentially numerical scheme, we shall still want to test it against known analytic solutions for special cases, or we may wish to differentiate a function based on interpolated numerical values.

Differentiation

Differentiation and integration work in a straightforward way.  There are various notations for differentiation, of which the D operator is perhaps the most flexible.

Clear[x];
D[x^n, x]

Integrate[%, x]

Partial differentiation can also be carried out by stipulating the variable which you wish to differentiate by.

D[x^4 y^2 + x^2, y]

Mixed and multiple derivatives are also supported.  We can differentiate once with respect to y and once with respect to x.

D[x^3 y^2 + x^2, y, x]

Higher derivatives are specified using the list structure to supply the variable and the number of times by which to differentiate with respect to that variable.

D[x^3 y^2 + x^2, {x, 3}, {y, 2}]

Integration - Symbolic and Numeric

Integrating with respect to y will reverse differentiation up to an arbitrary additive function of x.

Integrate[%%, y]

Similarly multiple integration can be carried out.  As in the examples for differentiation, to integrate with respect to x and y:

Integrate[x^2 y,x,y]

Definite integration needs a list supplying the upper and lower limits of the integration variable.

Integrate[1/(1-x),{x,2,4}]

There is also the facility to calculate integrals numerically.  This is available for any functions, including those which Mathematica cannot do symbolically. As in the case of definite symbolic integration, we must specify the range of integration. This can frequently be used to compute expected values of complicated payoffs, if they do not work symbolically. This next example shows how symbolic integration is not possible, but numerical integration yields a result.

Integrate[Sin[Sin[z]],{z,0,2}]

NIntegrate[Sin[Sin[z]],{z,0,2}]

This is another function under which problems can occur, and you should take care.  For example:

Integrate[Sin[1000000 x],{x,0,2}]

N[%]

NIntegrate[Sin[1000000 x],{x,0,2}]

In trickier cases you might get a warning. Older versions can be fooled by highly oscillatory integrals but it usually works now!

Limits

Of course, differentiation is defined by a limiting procedure. Mathematica can take limits explicitly:

Limit[((x + h)^3 - x^3)/h, h -> 0]

Series

Mathematica also knows how to use Taylor's Theorem through the Series command. The use of Series generates a certain type of expression truncated by an "O" expression. These can be manipulated according to the usual rules of asymptotic series, or converted to a normal form without the "O". Note that the next power in the "O" term is not necessarily computed exactly:

Series[Sin[Sin[x]], {x, 0, 5}]

%^2

Normal[%%]

Simple Ordinary Differential Equations

A large part of this book is concerned with various methods for solving partial differential equations, so we shall explore that topic in some detail later. In our introduction we shall take a brief look at ODEs.

Symbolic Solution of Differential Equations

Mathematica is a useful tool for dealing with certain differential equations.  The function DSolve can deal with a large class of equations.  It is used in a similar way to Solve, except that you need to specify both the independent and dependent variables. Notice how the equation is specified by using the double = sign.

DSolve[y''[x] + y[x] == 0, y[x], x]

If you give appropriate initial conditions, the arbitrary constants are removed.

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

To check that this is indeed a solution of the given equation we use the technique of replacement. The solutions to these equations are given as replacement rules presented in a list. We therefore need to extract a particular element of our solution list and apply that rule to get the solution. Setting the solution to be represented by s we select the first element of the list (we only have one in this example), and then apply that rule to y[x] to see the expression which is a solution.

s=y[x]/.%[[1]]

Now we can check if this is indeed a solution to our original equation.

D[s,x] + x s

Numerical Solutions to Differential Equations

In a similar way to numerical integration there is a function NDSolve which will find numerical solutions to differential equations and extends the scope of the equations to which we can find solutions. This function returns an interpolating function, rather than the usual table of values.  This allows you to ask for the values of the solution anywhere in the solution range, for example.  As usual, the use of "/." allows for the extraction of numerical values.

Clear[x, y, t];
soln = NDSolve[ {x'[t] == y[t],y'[t] == -4.0*x[t] - 0.1*y[t],
x[0] == 1,y[0] == 1},{x, y}, {t, 10}]

Note that the answer is supplied as an InterpolatingFunction - this allows a numerical solution computed in steps to be used just like an ordinary function. We can select out the specific solution for any t.  At t = 2,

{x[2], y[2]} /. soln

We shall use the InterpolatingFunction as a means of handling the output of finite-difference solvers for option equations. We shall explore numerical methods for solving PDEs in some detail later in this book, starting in Chapter 13.

Complex Numbers

Built-In Functions for Complex Manipulation

The use of complex numbers and functions is relatively unusual in derivatives modelling, but recent research has introduced Laplace transform methods as a powerful tool for the analytical pricing of options previously thought to be unpriceable except through an essentially numerical scheme (e.g. arithmetically averaged Asian options). So we need to use complex numbers, and to deal with them in Mathematica. The symbol I (capital) is the square root of -1, in Standard Form.  Mathematica operations work on complex numbers in just the same way they work on real numbers. In addition, various special operations for dealing with complex variables are provided.

Let us find the solution set for the following equation.

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

We put that into an approximate numerical form,

N[qsoln]

and then create the set containing the numerical solutions to the original equation:

solutions = x /. %

As with the exponential function, a special typeface is used for I in Traditional Form notation:

TraditionalForm[qsoln]

The function Re will return the real part of any complex number.  We can apply this function to our set solutions to get the real parts of our set.

Re[solutions]

Similarly, the function Im will give the imaginary part of a complex number. But note that symbols are in general treated as complex - for example:

z = x + I y;

Im[z]

The function Conjugate will give the complex conjugate of a complex number:

Conjugate[solutions]

The function Abs gives the absolute value of a complex number (i.e. (Sqrt[(x^2)+(y^2)])):

Abs[solutions]

The function Arg will give the argument of the complex number (i.e. ArcTan(y/x))

Arg[solutions]

There is a useful ComplexExpand function, which carries out an expansion assuming that all variables are real:

ComplexExpand[Cot[x + I y]]

This function is very useful for getting Mathematica to output what you expect to see when your mind assumes things are real, but the computer can't read your mind!

Complex Contour Integration

For integration in the complex plane the function NIntegrate can take complex numbers to specify the boundary around which the function is to be integrated.

NIntegrate[1/z, {z, 2, I, -2, -I, 2}]

We shall use complex contour integration to invert certain types of Laplace transform for option values.

2.7 More about Functions

We have already seen several Mathematica functions in the course of our investigations so far.  These are only a small proportion of those available. Some of the others are obvious.  One example is trigonometric functions.  Since you already know we can use Sin you would expect Cos, Tan and all the usual inverse and hyperbolic versions. These can all be used in Mathematica in the same way as Sin, as can the inverses of these trigonometric functions, as illustrated below.

Plot[ArcSin[x], {x, -1, 1}]

You will have already spotted that all Mathematica's built-in functions begin with a capital letter. Names which can be thought of as composed of two or more pieces may have capitals in the middle.  This means that when naming your own functions you can use any name starting with lower-case letters and be sure that you are not duplicating any Mathematica built-in function. We'll come back to user-defined functions shortly.

Getting Information about Built-In Functions

As you start to explore the capabilities of Mathematica you will need to get information about the functions available.  There is the facility to get information about built-in functions from within Mathematica.  To ask about a function just type a "?" before the function name.

?Sqrt

The wild card "*" can also be used to obtain a list of functions, with a similar sequence in common.  So we can obtain a list of all built-in functions related to the trigonometric function Sin.

?*Sin*

Many of the built-in functions that you will use have options associated with them.  These set default values for various variables associated with the execution of that function.  To see how the options are set for a function, Plot for example,

Options[Plot]

Changing any of these options can be done when executing the function by explicitly stating the desired option value as the final argument to the function.  For example, we can change the label points on the x-axis for a sine curve plot, and add suitable grid lines.

Plot[Sin[x],{x,0,4 Pi},
		Ticks -> {{0, Pi, 2Pi, 3Pi, 4Pi}, Automatic},
		GridLines -> {Table[k*Pi, {k, 0, 4, 1/2}] ,None} ]

Note that Mathematica allows a slight abuse of notation when it comes to pre-multiplying a symbol by a number. In the above,  2Pi is interpreted as twice π, even in the absence of a star or space. Options may also be used to refine any problems with functions such as NIntegrate:

Options[NIntegrate]

It is especially important that you understand the Assumptions option to Integrate. The following example should help you appreciate that Mathematica will actually analyze a symbolic integration containing parameters to figure out when the integral makes sense. You might well have some assumptions of your own that are actually your internal presumptions - Mathematica cannot read your mind and could well considerably more time figuring out when the integral makes sense than it does working it out.

Integrate[Exp[-a x], {x, 0, Infinity}]

Integrate[Exp[-a x], {x, 0, Infinity}, Assumptions -> a>0]

Integrate[Exp[-a x]Sin[b x], {x, 0, Infinity}]

FullForm[%]

Integrate[Exp[-a x]Sin[b x], {x, 0, Infinity}, Assumptions -> {a>0,b>0}]

Here is the example of the Call payoff from the Heston model:

Integrate[(Exp[x]-K)Exp[I w x], {x, K, Infinity}]

Integrate[(K-Exp[x])Exp[I w x], {x, -Infinity,K}]

Note how you get the correct restrictions on the domain on which the transform is defined!

Built-In Special Functions

Examples of Special Functions

Mathematica also knows a phenomenal number of special functions. If you do the kind of work that employs handbooks of functions and their values and properties, you may be surprised by the breadth of Mathematica's built-in capabilities. As well as a host of number-theoretic and combinatorial functions, there are many "classical" functions. Let us take a look at a few examples.  First of all, we can use Bessel functions:

Plot[{BesselJ[0, x], BesselI[1, x]}, {x, 0, 10},
PlotRange -> {-5, 5}, PlotStyle -> {Thickness[0.015], Thickness[0.01]}]

There is a selection of orthogonal polynomials:

LegendreP[5, x]

SphericalHarmonicY[1, 1, theta, phi]

This is probably a good place to remark on how Greek letters are input - previously we have used π without explaining how we wrote it, and our spherical harmonic function looks rather better written as:

TraditionalForm[SphericalHarmonicY[1, 1, θ, ϕ]]

The trick is to write "ESC\thetaESC" where ESC stands for the escape key. Occurrences of π are generated with "ESC\piESC".

Naming of Special Functions

This should give you a feel for the way in which these functions are named.  Many special functions can be classified under an overall name, like Legendre.  Within this there may be more than one type, like the P and Q Legendre functions.  And finally, these may be indexed by an integer or some other variable. Hence the Bessel function (Subscript[K, n])(x) is denoted by

BesselK[n, x]

Most of these functions are in the core kernel.  However, the statistical distributions are accessed from one of the Packages. We shall look at Packages later. Note that Mathematica knows some particular values of special functions, and also knows series:

Gamma[1/2]

Series[BesselJ[0, x], {x, 0, 6}]

Building Your Own Functions

Defining a Simple Function

When defining a function by using an explicit formula, it is important to specify the variables involved in that function by use of the underscore "_".  The delayed assignment sign ":=" then allows you to define the formula of your function.

h[x_, y_] := x^2 / y^3

This has defined the function h.  We can now ask for that function to be evaluated for specific values of x and y.

h[1, 2]

You can now ask about your function:

?h

If you wish to clear that definition, you use Clear.

Clear[h]

Your definition for h no longer exists.

?h

There are many types of function which can be defined, in the same way as the algebraic one we have just defined. Functions can employ any of the Mathematica built-in functions;

funcSin[x_,y_] := Sin[x] + Cos[y]^2

funcSq[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]

Functions can also be defined which will perform a test on the argument(s) returning the result True or False. testDiv5 will check whether the argument is divisible by 5:

testDiv5[t_] := IntegerQ[t/5]

testDiv5[8]

testDiv5[15]

Defining Functions in Their Pure Form

If you only intend to use the function once, or to use it with operations such as Map and Nest, there is a more economical way to define it. Instead of using a letter or name for the arguments in the function, the syntax for "pure" functions uses the "#" to indicate the arguments. So the expression which defines the action of a function to add 3 to a number is # + 3 instead of x + 3. To actually turn this into a function we use an "&" after the expression of the function. So

f[x_] := x + 3

can be written as

(# + 3) &

without using the definition of f. So if we apply that function to 2,

(# + 3) & [2]

This technique can also be used instead of defining functions of more than one variable. Using the # notation, we now label each of the arguments with a number immediately (i.e. no space) after the #. So instead of

g[x_, y_] := x^2 + y

g[2,3]

we use

(#1^2 + #2) & [2,3]

Notice how the numbers used for labelling directly correspond to the order in which those arguments will be given. The advantage of this technique for defining functions which will only be used once is that the definitions do not need to be stored, they are just used at the appropriate moment.

Using Pure Functions in Simple Programs

Programs can be shortened in length by using pure functions instead of separately defining a general global function which can be used in those functions which specifically ask for function definitions.

NestList[#^2 &, x, 5]

Fold[Sqrt[#1 + #2] &, t, {l, m, n}]

Clear[p,q,r,s]

NestList[Append[Rest[#],1]&, {p,q,r,s}, 3]

This form of writing functions can also be used for defining functions which will perform a test. As an example we can define a function which will return True if a number is divisible by 3, and use this as the criterion in the function Select.

Select[ {1, 2, 3, 4, 5, 6, 7, 8, 9}, IntegerQ[(#/3)] &]

A function which will return True if the second argument is less than the first can be used as the test by which a list can be sorted.

Sort[ {8, 34, 9, 2, 7, 4, 75, 4}, (#2 < #1) &]

This technique can be particularly useful when manipulating lists.

Defining a Function Recursively

There may be functions which you wish to define recursively giving values for the end points. This can be done in the following way.

Clear[f]

f[1] = 1;

f[2] = 1;

f[n_] := f[n-1] + f[n-2]

We can now evaluate f[3].

f[3]

There is a useful modification of this function which makes it much more efficient, that forces the function to remember previously computed values.

Clear[g]

g[1] = 1;

g[2] = 1;

g[n_] := g[n] =g[n-1] + g[n-2]

The following illustration makes it clear that the use of the second function for computing the Fibonacci series is strongly preferred:

Timing[f[35]]

Timing[g[35]]

Such a device will turn out to be useful in a simple implementation of the binomial tree model, where we recurse back through a tree from the payoff. Note that here and throughout this text, the timings reported are based on the computer used for editing. You must make allowances for the clock speed and efficiency of your own system when making comparisons or assessments about what computations are feasible.

The Difference between := and =

Constructing functions is dependent on the use of ":=".  The use of the colon with equals means that the evaluation of the right side of the equation is delayed, and it will only be evaluated when that function is actually called.  As an example let us define some functions which will be dependent on the value of a variable that we call M and to which we assign a value of 2 initially:

M=2

Now if we define funcOne, using :=, the right-hand side will not be evaluated yet.

funcOne[x_] := M x

Defining the function funcTwo using = will mean that the right-hand side is evaluated using the present value of M.

funcTwo[x_] = M x

With the present definitions both funcOne[2] and funcTwo[2] will evaluate to 4.

{funcOne[2],funcTwo[2]}

If we were to change the value of M however, the definition of funcOne would: allow the new value to be used, whereas the definition of funcTwo has already been set by the original definition of M:

M=100

{funcOne[2],funcTwo[2]}

Assignment and Random Walks

Another way to see the difference between = and := is to define two functions from which we can produce random walks and plot those resultant functions. Firstly let us define two functions which will add a random number between -0.5 and +0.5 to their argument:

addrandONE[x_]:=x+RandomReal[{-0.5,0.5}]

addrandTWO[x_]=x+RandomReal[{-0.5,0.5}]

So if we ask for each of these functions evaluated twice on the same value,

{addrandONE[0],addrandONE[0]}

{addrandTWO[0],addrandTWO[0]}

The second function does not act as a random function. Now we can form random walks based on these functions using the function NestList,

walkONE = NestList[addrandONE,0,1000];

walkTWO = NestList[addrandTWO,0,1000];

and now we can display these lists using ListPlot and asking that the points should be joined.

ListPlot[walkONE,Joined->True]

ListPlot[walkTWO,Joined->True]

From this it is obvious what has happened in the function addrandTWO. The function has added a constant number to the argument of addrandTWO, that number being the first random number which was generated when the function was defined. The function addrandTWO was displayed with the Random function replaced by a number when the function definition was entered and from that point on the function definition remained as such. 

Note the ease with which random walks may be constructed - we shall use this facility when carrying out Monte Carlo simulation of paths followed by stock prices.

Functions of Two or More Variables

The basic functions to know about are Plot3D, ContourPlot, DensityPlot and their List analogues. The forms of the command are very similar to those for Plot etc. A few examples will suffice. Let's make a function to use in our examples.

func[x_, y_] := Sin[x^2] Cos[y];
Plot3D[func[x, y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 40]

ContourPlot[func[x, y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 40]

DensityPlot[func[x, y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 60, Mesh -> False]

In the same way as for two-dimensional curves, the functions for three-dimensional surfaces also have options associated with them which allow the alteration of graphics to suit your needs.

Plot3D[func[x, y], {x, -2, 2}, {y, -2, 2},
ViewPoint -> {0, 1, 0.3},
Axes -> None, Boxed -> False, PlotPoints -> 20]

ContourPlot[func[x,y],{x,-2,2},{y,-2,2},Contours->40,ContourStyle->None,ColorFunction->Hue,Frame->False]

Animation

There is also the facility in Mathematica to be able to animate graphics which you have produced. This is done by creating all of the cells which will make up your movie and then asking to see them displayed one after the other. As an example we can see how a three-dimensional surface alters as we step x and y over slightly larger values. Only the first frame is shown in the book. Such techniques can be used to present plots of option values and their Greeks as functions of other variables, using the animation to represent the evolution through time, or perhaps to vary the angle of view, in order to bring out particular features.

Manipulate[Plot3D[Sin[x+t] Cos[y+t], {x, -Pi, Pi}, {y, -Pi, Pi}],
	{t,0,2*Pi-Pi/10,Pi/10}]

Data analysis and distributional estimation

FTSE

series = "^FTSE"

FinancialData[series, "Name"]

ftserets = FinancialData[series,"Return",All];
logftrets = Map[Log[1+#[[2]]]&, ftserets];

First[ftserets]

Last[ftserets]

Skewness[logftrets]

Kurtosis[logftrets]

gaussestimate = EstimatedDistribution[logftrets,NormalDistribution[a,b]]

gausslike = Likelihood[gaussestimate, logftrets]

testimate = EstimatedDistribution[logftrets,StudentTDistribution[a,b, n]]

tlike = Likelihood[testimate, logftrets]

Log[10,tlike/gausslike]

jsu = EstimatedDistribution[logftrets,JohnsonDistribution["SU", γ, δ, μ, σ]]

jsulike = Likelihood[jsu, logftrets]

Log[10,jsulike/gausslike]

See Fergusson and Platen, Applied Mathematical Finance, 2006.

Parallel Programming 
GridMathematica

Introduction

This note gives an introduction to parallel financial computation using Mathematica V8. The standard version will support 4 kernels if your machine has at least 4 cores. In this lecture I will be working off an 8-core box that is also networked.

Although this lecture is targetted at Mathematica, some of the ideas are relevant much more generally and I will flag these as I go along. 

There are ways of configuring parallel kernel configuration and start-up that is more automatic than what I will do. I have stayed rather manual in my own implementation. Once all the needed kernels are up the computational methods will be universal and not dependent on configuration issues. 

My own small grid is also pretty homogenous - all the machines are 8-core Macs running Leopard or Snow Leopard and Mathematica 7.01, clocked at 2.8GHz. You can work much more heterogeneously, with different OS and clock speeds. It is generally a good idea to use the same version of the software, so as to avoid code incompatibility. There are load adjustment tools if you have machines with wildly varying performance.

Local Parallel Start-up

You can trigger this by asking for a parallel computation, but I like to launch my kernels and see what I have to play with before doing calcutions. First of all I will just do a launch of kernels on my own machine (bear in mind I wrote these notes on both a home 8-way machine and revised them on the KCL cluster, so do not worry about variations in the names of machines or Mathematica version/date).

Quit[]

Kernels[]

If these have not started automatically you do this:

LaunchKernels[]

If you go to the parallel kernel status menu you can see something like this:

This reveals that I am running a master kernels and 8 local "slave" kernels. I can do work on either the master or slaves now. There are a pile of basic functions to do with kernels and parallel commands that you can explore.

?*Kernel*

?*Parallel*

For example, I can ask from the notebook some information about my slave kernels.

TableForm[ParallelEvaluate[{$MachineName,$SystemID,$Version,TimeUsed[], MaxMemoryUsed[]}],TableHeadings->{None,{"host", "OS","Mathematica Version", "CPU Time","Memory"}}]

Key Ideas

You can consider various types of calculation for parellelization. Clearly any computation that can be broken down into many identical or similar tasks can benefit in principle. What you do need to watch out for is the balance between doing the task and communication the results back to the master kernel. A demanding task requiring a short result to be sent back is clearly most efficient. You will get much benefit by sending huge lists around the system - this will apply to any type of software. People often spend large amounts of money to increase the inter-node bandwidth. Here I have 1G ethernet between machines, but you can even bog down on-chip communication by being careless with data transfer.

ParallelTable

If you can write what you want to do as Table[something] then there is an easy route to parallelization - ParallelTable!

AbsoluteTiming[Table[Integrate[Sin[x]/Sqrt[x], {x, i Pi, Infinity}], {i, 0,63}];]

AbsoluteTiming[ParallelTable[Integrate[Sin[x]/Sqrt[x], {x, i Pi, Infinity}], {i, 0,63}];]

We have factor of six speedup right away. I have suppressed the actual results in the above - you can delete the semi-colon if you wish to see the output. But you must not do this blindly with huge piles of data. Here is my favourite screw-up! The right way to make lots of random numbers is to use RandomReal with its second argument applied to make lots.

AbsoluteTiming[data = RandomReal[{0,1}, 1000000];]

This is a lot faster than using Table:

AbsoluteTiming[data = Table[RandomReal[{0,1}], {1000000}];]

And the following speaks for itself.

AbsoluteTiming[data = ParallelTable[RandomReal[{0,1}], {1000000}];]

And this does not help that much

AbsoluteTiming[
data = ParallelTable[RandomReal[{0,1}, 25000], {40}];]

The lesson is to simulate what you want and keep the lengthy stuff there, working on it remotely and only sending back what you really need. Let's do a basic expectation summing up locally and then returning that one piece of information.

AbsoluteTiming[
Apply[Plus, ParallelEvaluate[Apply[Plus,RandomReal[{0,1}, 25000]]]]/1000000]

The real benefit comes when you do lots of computation on each slave and send back a simple short result. This is the approach we shall use from now on.

Monte Carlo Parallel foul-up

Another thing NOT to do is to be silly about seeding your Monte Carlo simulations. GridMathematica is set up carefully in such a way that each slave kernel is in a different random seed state.

ParallelTable[RandomReal[], {8}]

Here is how to foul it all up by careless seeding:

ParallelEvaluate[SeedRandom[100]];

ParallelTable[RandomReal[], {8}]

You might well wish to seed each slave in a certain way (e.g. to get reproducible random numbers), but you should do it another way. I tend to feed off the kernel ID, maybe with an integer multiplier.

ParallelEvaluate[SeedRandom[100*$KernelID]];

ParallelTable[RandomReal[], {8}]

Example: Bivariate Exponential Distribution and the DistributeDefinitions command

So what is the BV exponential distribution anyway? Something I have been exploring with my student Asad Munir is the generation of natural multivariate distributions using equilibria of coupled SDEs, each of whose marginals has as an equilibrium a certain target distribution. The maths of this is in our paper that has just appeared:

W.T. Shaw & A. Munir, Dependency without Copulas or Ellipticity, European Journal of Finance, Vol 15, Issues 7&8, October 2009, pp 661-674.

(Subscript[dX, t])=((σ^2)/2) dt (1-λ (Subscript[X, t])) +σ (Subscript[dW, t])

(Subscript[dY, t])=((σ^2)/2) dt (1-λ (Subscript[Y, t])) +σ (Subscript[dZ, t])

and the two Brownian motions are correlated. We get a natural bivariate exponential distribution in equilibrium.

BiExpPathGen[{x0_, y0_}, sig_, rho_, λ_, t_, ts_] := 
Module[{srt = sig*Sqrt[t/ts],nlista = RandomReal[NormalDistribution[0,1], {ts}],
nlistb = RandomReal[NormalDistribution[0,1], {ts}], corrlist, pairlist},
corrlist = rho*nlista + Sqrt[1-rho^2]*nlistb;
pairlist = Transpose[{nlista, corrlist}];
Fold[Abs[#1+1/2 srt^2(1-λ #1)+ srt*Sqrt[#1] #2]&, {x0, y0}, pairlist]]

So that is fine, but the slave kernels do not know the definition - only the master kernel does. We fix this as follows:

DistributeDefinitions[BiExpPathGen]

AbsoluteTiming[bidata = Table[BiExpPathGen[{1,1}, 0.4, 0.5, 0.5, 200, 3200], {k, 1, 4000}];]

AbsoluteTiming[bidata = ParallelTable[BiExpPathGen[{1,1}, 0.4, 0.5, 0.5, 200, 3200], {k, 1, 4000}];]

%%/%

So here is the efficiency of parallelization for this case:

100*%[[1]]/8

Here is a nice bivariate exponential distribution as well!

ListPlot[bidata, PlotRange -> {{-1, 10}, {-1, 10}}, AspectRatio -> 1, PlotStyle ->AbsolutePointSize[2]]

While we are on this topic, note that by varying the SDEs you can go heterogeneous on the choice of marginal distirbutions with this technique.

Bivariate Exotic: "exponential-normal" distribution!

(Subscript[dX, t])=((σ^2)/2) dt (1-λ (Subscript[X, t])) +σ (Subscript[dW, t])

(Subscript[dY, t])=-((σ^2)/2) dt (Subscript[Y, t]) +σ (Subscript[dZ, t])

The equilibrium is exponential in X and Gaussian in Y, and they are correlated:

ExpNorPathGen[{x0_, y0_}, sig_, rho_, λ_, t_, ts_] := 
Module[{srt = sig*Sqrt[t/ts],nlista = RandomReal[NormalDistribution[0,1], {ts}],
nlistb = RandomReal[NormalDistribution[0,1], {ts}], corrlist, pairlist},
corrlist = rho*nlista + Sqrt[1-rho^2]*nlistb;
pairlist = Transpose[{nlista, corrlist}];
Fold[{Abs[#1[[1]]+1/2 srt^2(1-λ #1[[1]])+ srt*Sqrt[#1[[1]]] #2[[1]]],
#1[[2]]-1/2 srt^2#1[[2]]+ srt*#2[[2]]}&, {x0, y0}, pairlist]]

ExpNorPathGen[{1,0}, 0.4, 0.5, 0.5, 1, 10]

DistributeDefinitions[ExpNorPathGen]

bidata = ParallelTable[ExpNorPathGen[{1,1}, 0.4, 0.5, 0.5, 200, 3200], {k, 1, 4000}];

ListPlot[Map[Reverse,bidata], PlotRange -> { {-4, 4},{0,10}}, AspectRatio ->5/4, PlotStyle ->AbsolutePointSize[2]]

This is a naturally-generated bivariate exponential-normal distribution!

Correlation[bidata]

Needs["MultivariateStatistics`"]

?KendallRankCorrelation

SpearmanRankCorrelation[Transpose[bidata][[1]], Transpose[bidata][[2]]]//N

KendallRankCorrelation[Transpose[bidata][[1]], Transpose[bidata][[2]]]//N

We know how to do this with mixtures of normal, exp, gamma, Student and so on. See "Dependency without copulas".

Other Parallel*** Commands

Once you have seen ParallelTable in operation the operation of some of the other commands becomes pretty much self-evident:

?ParallelSum

More about ParallelMap - an example doing many implied vols

This one I find particularly useful. If you can run your computer model essentially in the form

Map[model, lotsofdata]

then provided you have told the slaves what "model" is then you an try instead

ParallelMap[model, lotsofdata]

Here is a classic example working out many implied volatilities. We first define the BS model in the usual way:

Ncdf[x_] := (1+Erf[x/Sqrt[2]])/2

The analytical valuation of the vanilla Put is then given by the following definitions

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

BlackScholesPut[s_, k_, σ_, r_, q_, t_] := 
k*Exp[-r*t]*Ncdf[-dtwo[s, σ, k, t, r, q]] - s*Exp[-q*t]*Ncdf[-done[s, σ, k, t, r, q]]

Now I am going to make up a lot of random option data (normally this would come from markets)

SeedRandom[200]

randomoptionparameters = Table[{RandomReal[{0.001, 30}],RandomReal[{30, 200}], RandomReal[{80, 120}], 0.045, 0, RandomReal[{1, 10}]}, {2000}];

Some of these might, being random, be pathological with regard to an IV compuation,  so I will suppress some complaints here, just for the demo - one would normally trap and report these.

Off[FindRoot::lstol];Off[FindRoot::jsing];
ParallelEvaluate[Off[FindRoot::lstol];Off[FindRoot::jsing]];

test = randomoptionparameters[[1]]

So here is the usual IV computation using FindRoot (which to remind you is a Newton-Raphson solver if you feed it one starting value):

FindRoot[test[[1]]==BlackScholesPut[test[[2]], test[[3]], sigma, 0.045, 0, test[[6]]], {sigma, 0.1}]

ta = AbsoluteTiming[valsa = Map[sigma/.FindRoot[#[[1]]==BlackScholesPut[#[[2]], #[[3]], sigma, 0.045, 0, #[[6]]], {sigma, 0.1}]&, randomoptionparameters];]

Parallelization

Important - you just export the definitions to the grid environment and change Map to ParallelMap:

DistributeDefinitions[Ncdf,done,dtwo,BlackScholesPut];

tb = AbsoluteTiming[valsb = ParallelMap[sigma/.FindRoot[#[[1]]==BlackScholesPut[#[[2]], #[[3]], sigma, 0.045, 0, #[[6]]], {sigma, 0.1}]&, randomoptionparameters];]

Some people call this sort of thing "embarassingly parallel". It's certainly easy! The performance improvement is a factor of:

ta[[1]]/tb[[1]]

Close Down

CloseKernels[]

Quit[]