Computational Aspects of Financial Mathematics

Open the closed cells by double-clicking them. Evaluate the input cells by <Shift>+<Enter> or just <Enter> on the numeric pad.

Summer School in Financial Mathematics, June 2002, Warsaw, Poland

Jan Hurt
Division of Financial & Insurance Mathematics
Faculty of Mathematics & Physics
CHARLES UNIVERSITY PRAGUE

Home page: http://adela.karlin.mff.cuni.cz/~hurt/

System Mathematica

System for doing mathematics by computer. 
Reference:

www.wolfram.com

Why Mathematica in Financial Mathematics?

According to William Shaw (Modelling Financial Derivatives with Mathematica, Nomura International plc, Cambridge University Press, Cambridge 1998),

The system has the ability to

deal with thousands of special functions and perform symbolic calculations with them

manage advanced numerical computation

allow complex structures to be programmed

visualize functions in 2 or 3 dimensions

4 pillars

Numerical calculation

1 + 1

2^1000

N[Pi, 1000]

Symbolic calculation

e1 = Expand[(x + y)^10]
Factor[e1] // TeXForm

TeXForm[(a + 1)/(x^2 + 2*x - 3)]

(a + 1)/(x^2 + 2*x - 3)

Graphics

Plot[Sin[x + x^2], {x, 0, Pi}];

Plot3D[Sin[x + y^2], {x, 0, Pi}, {y, 0, Pi}];

Programming

Important commands to avoid collisions of possible previous definitions:

h[n_] := n^2
{h[2], h[3]}

?h

Clear[h]
?h

Remove[h]
?h

The Hilbert matrix of order n:

h[n_] := Table[1/(i + j - 1), {i, n}, {j, n}]

?h

(h3 = h[3]) // MatrixForm

?h3

h[4][[3, 1]]
h[4][[3]][[1]]

Technicalities, basic input

Front end and kernel
Menu

Saving, printing, aborting and interrupting evaluation  (Menu -> Kernel -> Abort Evaluation)

Case sensitivity: Built-in functions start with the capitals. Sin and sin are different symbols. It is a good practise not to use your own symbols that begin with capital letters, to prevent confusion with built-in Mathematica functions. Symbol "sin" is not known at the moment so the corresponding expression remains unevaluated. Observe also that Sin[x] may be written as x // Sin:

{Sin[3.14], 3.14 // Sin, sin[3.14]}

Parentheses and brackets:
 [ ] - arguments of functions
 ( ) - priority in algebraic expressions
 { }- lists
  [[ ]] - part of an expression

Operations + - * / ^    <space> = * . Beware of the difference "xy" and "x y"!

Clear[x, y, xy]
x = 5; y = 6;
xy == xy
{xy, x*y, x y}

EUR = 4 PLN;

5*EUR

"_" (underscore):  x_ a dummy variable

<;> behind the expression causes that the calculation is done but not displayed

hh = h[100];
hh[[50, 30]]

N[E, 5]

Notebooks in 7-bit ASCII code

The difference between exact and approximate (rational) numbers

2 == 2

2 == 3

{IntegerQ[2], IntegerQ[2. ]}

But

2 == 2.0

<Enter> = simple linefeed to enter

Entering the calculation: <Shift>+<Enter> or <Enter> on the numeric pad.

Complete expression: <Ctrl>+<K> (Very useful particularly for functions with long names.)

Using 2D mathematical notation either from a palette (similar to MS Word equation editor) or by a key sequence. Return to the line by <Ctrl>+<Space> (sometimes repeatedly).

(∑_(i = 1)^n) (i^4)

π

π // N

N[Eigenvectors[h[3]]] // Chop
Eigenvectors[N[h[3]]]

Capabilities of the system

For a more detailed materials go to  http://library.wolfram.com/tutorials/

This is only a short adapted extract from tour.nb:

Mathematica as a Calculator

6^20

6^200

Eigenvalues[(({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}))]

Eigenvalues[h[15] // N]

An important feature of Mathematica is its ability to handle formulas as well as numbers.

Clear[x, a]
Solve[(Sqrt[x]) + a == 2  x, x]

∫(Sqrt[x]) (Sqrt[a + x]) \[DifferentialD]x

You can use Mathematica to make 2D and 3D graphics.

Plot[Sin[x] + Sin[1.6  x], {x, 0, 40}];

Plot3D[Sin[x y], {x, 0, 4}, {y, 0, 4}, 
        PlotPoints -> 30];

Power Computing with Mathematica

Even though you can use it as easily as a calculator, Mathematica gives you access to immense computational power.

m = Table[Random[], {100}, {100}];

ListPlot[Abs[Eigenvalues[m]], PlotRange -> All];

Mathematica can handle numbers of any size. On most computers Mathematica takes under a second to compute the exact factorial of 1000.

1000!

N[π, 500]

Clear[x, y]
f111 = Factor[(x^99) + (y^99)];

Simplify[f111]

Accessing Algorithms in Mathematica

Whenever you use Mathematica, you are accessing the world's largest collection of computational algorithms.

FindRoot[Cos[x] == x + Log[x], {x, 1}]

Mathematica cannot find the explicit formula for this integral so that it returns the expression unevaluated

Integrate[Log[x + Sin[x]], {x, 0, 2}]

The numerical integration works, however:

NIntegrate[Log[x + Sin[x]], {x, 0, 2}]

We may also find a series expansion with the remainder term of order O((x^10)) and then evaluate the integral:

(s1 = 
          Series[Log[x + Sin[x]], {x, 0, 10}] // 
            Normal);
i1 = Integrate[s1, x]

This is an example of a transformation rule:

Clear[x]
(1/x) /. x -> 11

We calculate the limit of i1:

Limit[i1, x -> 0]

Here is a comparison of the results:

inum1 = 
      NIntegrate[Log[x + Sin[x]], {x, 0, 2}]
inum2 = (i1 /. x -> 2) - Limit[i1, x -> 0] // 
        N
inum1 - inum2

Solution to the algebraic equations:

Solve[x^5 - 6  x^3 + 8  x + 1 == 0, x]
NSolve[x^5 - 6  x^3 + 8  x + 1 == 0, x]

FactorInteger[2^105 - 1]

m = Table[2^i  +  x^j, {i, 3}, {j, 4}]

MatrixForm[m]

NullSpace[m]

matrix1 = {{1, 2, 3}, {4, 5, 6}}

Mathematica can solve differential equations both symbolically and numerically.

\(NDSolve[{\(\(x'\)'\)[t] + x[t]^3 == Sin[t], 
        x[0] == \(x'\)[0] == 0}, 
									x, {t, 0, 
        50}]\)

ParametricPlot[
        Evaluate[{x[t], x'[t]} /. %], {t, 0, 50}];

Mathematical Knowledge in Mathematica

Mathematica incorporates the knowledge from the world's mathematical handbooks — and uses its own revolutionary algorithms to go much further.

Mathematica knows about all the hundreds of special functions in pure and applied mathematics.

LegendreQ[3, x]

Mathematica is now able to do vastly more integrals than were ever before possible for either humans or computers.

∫(Sqrt[x]) ArcTan[x] \[DifferentialD]x

(∫_0^∞) Log[
            x] Exp[-(x^3)] \[DifferentialD]x

(∫_0^∞) Sin[
            (x^2)] Exp[-x] \[DifferentialD]x

(∑_(k = 1)^n) (1/(k^6))

(∑_(k = 1)^∞) (1/(k^6))

(∏_(n = 1)^5) Gamma[((2 n)/5)] // Simplify

FullSimplify[(∏_(n = 1)^5) Gamma[((2 n)/5)]]

TrigReduce[Cos[x]^4]

Prime[(10^9)]

Building Up Computations

Being able to work with formulas lets you easily integrate all the parts of a computation.

Eigenvalues[{{3,  1},  {2,  6}}]

v  =  Eigenvalues[{{3,  1},  {2,  b}}]

Mathematica's functions are carefully designed so that output from one can easily be used as input to others.

Plot[Evaluate[v],  {b,  -10,  10}];

Solve[First[v]  ==  0,  b]

int = (∫_0^c) First[v] \[DifferentialD]b

Series[int,  {c,  0,  5}]

FindRoot[int == 1 + c, {c, 1}]

Being able to work with formulas is also important in summarizing data.

tab1 = Table[Prime[i],  {i,  40}]

ListPlot[tab1];

?Fit

fit11 = Fit[tab1, {Log[x], x, (x^2)}, x]

Sum[fit11,  {x,  1,  40}]

Sum[Prime[i], {i, 1, 40}]

MaxMemoryUsed[]

Etc., see Help Browser, Tour of Mathematica

Basic operations

Recommended textbook on Mathematica for the beginners and for intermediate: Shaw, W. T. and Tigg, J.: Applied Mathematica (Getting Started, Getting It Done). Addison-Wesley. Reading (MA), 1994.

Elementary symbolic computation

Manipulating polynomials:

poly = 9*(x + 2)*(x + y) + (x + y)^2

e1 = Expand[poly]

e2 = Expand[e1^3]

f1 = Factor[e2]

y2 = (x^2 - 
              4)*(x - 
                  3)^2/((x - 2)*(x - 4)*(x - 
                    6));
e12 = Expand[y2]
e13 = ExpandAll[y2]

Traditional form of the above input has been obtained by marking the cell and then Menu Cell -> Convert To -> Traditional Form:

y2 = ((((x^2) - 4) (((x - 3))^2))/((x - 2) (x - 4) (x - 6)));

Apart is a command which is particularly useful for handling probability generating functions or characteristic functions.

{e14 = Apart[e13], Apart[y2]}

tim1 = Timing[D[Sin[x], {x, 4}]]
tim1[[1]]
tim1[[2]]
tim1[[1]]/Second

{Timing[D[e13, {x, 40}]][[1]], Timing[
          D[e14, {x, 40}]][[1]]}

Clear[x, y]
yx = (1 + x)^2*(x - y)^3;
Collect[yx, x]
yxexp = Expand[yx]
Collect[yxexp, x]
Collect[yxexp, y]

Coefficient[(x^4) + 4  (x^3) y + 4  x (y^3) + (y^4), (x^3)]

Here is an equivalent way.

Coefficient[(x^4) + 4  (x^3) y + 4  x (y^3) + (y^4), x, 3]

Substitution: "/." is the replacement operator and "->" means a transformation rule:

yx /. x -> y
yx /. x -> 2  y
yx /. x -> 3

Integration

First, beware of the syntax: The expression

Sin^2[x]

should be written as

Sin[x]^2

which becomes

Sin[x]^2 // TraditionalForm

in the traditional form.

Input form:

i1 = Integrate[x^2*Sin[x]^2, x]

Traditional form:

i1 = ∫(x^2) (sin^2)(
              x) \[DifferentialD]x

Checking the integration process:

d1 = D[i1, x]

This is not the original function. But we get it after simplification:

d1 // Simplify

Another examples of a series expansions:

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

s2 = Series[Sin[x], {x, a, 3}]

s3 = Series[d1, {x, 0, 5}]

Series[(f[x + h] - f[x - h])/(2  h), {h, 0, 6}]

Further examples of integration
Integrals which may be expressed in terms of special functions:

Integrate[Sin[x^2], x]

Integrate[Exp[-x^2], {x, -∞, z}]

Multiple integrals, beware of the order of integration:

Input form:

Integrate[x*y, {x, a, b}, {y, 0, x}]

Traditional form:

(∫_a^b)(∫_0^x) 
            x y \[DifferentialD]y \[DifferentialD]x

Integration of an ugly function:

Plot[Sin[200  x], {x, 0, 10}];

i0 is the exact value of the integral, i1 its numerical value, and i2 is a result of a numerical integration:

i0 = Integrate[Sin[200  x], {x, 0, 10}]
i1 = 
      N[Integrate[Sin[200  x], {x, 0, 10}], 20]
i2 = NIntegrate[Sin[200  x], {x, 0, 10}]
i1 - i2

The result seems not to be completely satisfactory. We may try to change some of the options:

Options[NIntegrate]

i3 = 
      NIntegrate[Sin[200  x], {x, 0, 10}, 
        MaxRecursion -> 16]
i1 - i3

i4 = 
      NIntegrate[Sin[200  x], {x, 0, 10}, 
        Method -> Oscillatory]

End 20020627

Solving equations

Here an equation eq is defined

eq = x^3 - 7  x^2 + 3*a*x == 0

and consequently solved:

s1 = Solve[eq, x]

Another example:

p1 = Expand[(x + 1)^7]
Solve[p1 == 0, x]

Back substitution.  The following construction is of vital importance. Instead of transformation rules we get the list of the roots:

roots = x /. s1

We may extract the second root:

roots[[2]]

A system of two linear equations. Observe that the result is in the form of a double list.

s2 = Solve[{x + y == 0, a*x + b*y == c}, {x, y}]

Here is the method how to overcome the problem. Use the "Flatten" command:

Flatten[s2]

(x + y) /. s2
(x + y) /. (s2 // Flatten)

Numerical solution to algebraic equations. Notice that the definition of the equation may be done within the command.

Solve[eq2 = x^5 + 2  x + 1 == 0, x]

x^5 + 2  x + 1 // Factor

NSolve[eq2, x]

Solving a system of nonlinear equations.

Solve[
      eq3 = {Sin[x] == x - y, Cos[y] == x + y}, {x, 
        y}]
NSolve[{Sin[x] == x - y, Cos[y] == x + y}, {x, 
        y}]

The FindRoot function with the initial values (1,0) does work, however:

FindRoot[eq3, {x, 10}, {y, 10}]

Lists

Lists are general objects taking the form {{...},{...,{...,{}}}}. The special cases are vectors {(Subscript[a, 1]),  ... , (Subscript[a, n])}and matrices {{(Subscript[a, 11]),  ... , (Subscript[a, 1  n])},  ... , {(Subscript[a, m1]),  ... , (Subscript[a, mn])}}.
This is a structured list:

li1 = {1, 2, 3, 4, {5, 6, {7, 8}}}

Its tree form looks like

li1 // TreeForm

Here is another (vector) list:

t1 = Table[n!, {n, 1, 10, 0.5}]

Function Log (as many other built-in functions) is a listable function:

t2 = Log[t1]

ListPlot[t2, PlotJoined -> True];

The built-in function Range returns also a list:

?Range
Range[5]
Range[3, 7]
Range[1, 2, 0.2]

Example of a matrix, Note that GCD[n1, n2, ... ] gives the greatest common divisor of the integers n1, n2, ....

(array = Table[If[GCD[i, j] == 1, 1, 0], {j, 3}, {i, 3}]) // 
      MatrixForm

array

ListDensityPlot[IdentityMatrix[10]];

ListDensityPlot[Table[Random[], {10}, {10}]];

array = 
        Table[If[GCD[i, j] == 1, 1, 0], {j, 50}, {i, 50}];

ListDensityPlot[array];

Let us turn to the Hilbert matrix.

mat = Table[1/(i + j - 1), {i, 7}, {j, 7}];

This is the inverse:

invmat = Inverse[mat]

Check it. Observe that "."  (dot) means the matrix multiplication (in case of vectors always the scalar product), while " * " (asterisk) is the Hadamard product (componentwise multiplication). Also note that up to this point we work with the exact numbers.

mat . invmat

Here is a numerical solution to the problem:

matnum = N[mat];

j1 = Inverse[matnum] . matnum

?Chop
j1 // Chop

e5 = Eigenvalues[mat // N]

d5 = Apply[Times, e5]
d6 = Det[mat]
d5 - d6

Remark on matrix multiplication. If x and y are (nx1) column vectors, then (x^T) y is a scalar product while (xy^T) is an (nxn) matrix. Since the vector is not considered as a matrix in Mathematica, the second case must be handled separately.

Range[5]*Range[5]
Range[5] . Range[5]

MatrixForm[{x, y}]
MatrixForm[{x, y} // Transpose]
MatrixForm[{{x}, {y}}]
MatrixForm[Transpose[{{x}, {y}}]]

{{x}, {y}} . {x, y}
{{x}, {y}} . {{x, y}}

t3 = Transpose[t4 = {{1, 2, 3}}]
t3 . t4
t4 . t3
(t4 . t3)[[1, 1]]

Another method:

Outer[Times, {x, y}, {x, y}] // MatrixForm

Another use of Outer  (Jacobian matrix):

Outer[
      D, {f1[x1, x2, x3], f2[x1, x2, x3], f3[x1, x2, x3]}, {x1, 
        x2}]

Flatten, Partition

t1 = {1, 2, {3, , 4, {5, 6}}};
f10 = Flatten[t1]
f11 = Flatten[t1, 1]
TreeForm[t1]
p1 = Partition[f10, 2]
p2 = Partition[f10, 4]

Random matrix

e10 = Eigenvalues[Table[Random[], {100}, {100}]];

Apply[Times, e10]

Observe a peculiar property of such a random matrix of order n. The eigenvalue with the maximum real part has the absolute value about n/2.

ListPlot[Transpose[{r10 = Re[e10], i10 = Im[e10]}], 
        PlotRange -> {{49, 51}, Automatic}, 
        PlotStyle -> PointSize[0.04]];

{Max[r10], Min[r10], Max[i10], Min[i10]}

We sort the real parts and take 10 greatest of them:

Take[Sort[r10] // Reverse, 10]

We generate 20 random matrices of order 120.

tab10 = 
        Table[Eigenvalues[Table[Random[], {120}, {120}]], {20}];

Next we take maximum real parts of these 20 matrices:

Map[Max, Re[tab10]]

Remark

Listability may be used to improve efficiency. Compare:

tim1 = 
        Timing[Table[
            Sin[x], {x, 
              0, π, π/10000. }]];
tim2 = 
        Timing[Sin[
            Range[0, π, π/10000. ]]];
{tim1[[1]], tim2[[1]]}

User-defined functions

Next prime

Built-in functions ending with "Q" are logical functions which return "True" or "False". The exclamation sign "!" means logical "Not", i.e., negation.

{VectorQ[{Pi, 1/3, 1. /3}], 
      MatrixQ[{Pi, 1/3, 
          1. /3}], 
MatrixQ[{{Pi, 1/3, 
            1. /3}}], 
Not[
        MatrixQ[{{Pi, 1/3, 1. /3}}]]}

Table[! PrimeQ[j], {j, 1, 6}]

{LowerCaseQ["a"], LowerCaseQ["B"]}

In the next example we will find the next prime number which follows after a given natural number. Remain unevaluated if the argument is not a number. The construction uses "Module" is used to form an expression (a program) with more than one command. It consists of two parts, the list of local variables with possible initial values assigned and the body. If there is more than one command in the body, they are separated by the semicolon. Last expression of the body is returned unless there is a "Print" command inside the body.

?Module

nextprime[n_] := 
      Module[{k = n}, While[! PrimeQ[k], k++]; 
        k /; NumberQ[n]]

{nextprime[1000000], nextprime[{2, 1}]}

Mean, variance

obs = Range[6]
Length[obs]

mean[data_] := Apply[Plus, data]/Length[data]

mean[obs]

mean[{a, b, c, d}]

mean[{}]

A slight improvement. Evaluated only if data is a vector of a positive length.

Clear[mean]
mean[
        data_] := (Apply[Plus, data]/Length[data] // 
            N) /; (VectorQ[data] && Length[data] > 0)

mean[{a, 2, 3}]
mean[{a, 2, {3}}]
mean[{}]

Variance is simply the mean of the squared deviations of the data from their mean. Listability is utilized.

var[data_] := mean[(data - mean[data])^2]

var[{a, b, c}]
var[obs]

Random walk 1

First look on another important built-in function.

FoldList[Plus, 0, {a, b, c}]
FoldList[Times, 1, {a, b, c}]

One-dimensional random walk.

randomwalk[n_] := 
      FoldList[Plus, 0, Table[Random[] - 0.5, {n}]]

SeedRandom[13]
ListPlot[randomwalk[100], 
        PlotJoined -> True];
ListPlot[randomwalk[100], PlotJoined -> True];

Two-dimensional random walk.

randomwalk2[n_] := 
      FoldList[Plus, {0, 0}, 
        Table[{Random[], Random[]} - 0.5, {n}]]

randomwalk2[5]

ListPlot[rw2 = randomwalk2[1000], 
        PlotJoined -> True];

{rw2[[50, 1]], rw2[[50]][[1]]}

Random walk 2

First define the function with delayed evaluation.

addone[x_] := x + Random[] - 0.5

Second, define the function with immediate evaluation. Random number is generated just once and subsequently used every time the function is called.

addtwo[x_] = x + Random[] -  .5

Look at the difference.

{addone[0], addtwo[0]}
{addone[0], addtwo[0]}

Next we use another useful function.

NestList[Sin, x, 3]

walkone = NestList[addone, 0, 1000];
walktwo = NestList[addtwo, 0, 1000];

ListPlot[walkone, 
        PlotJoined -> True];
ListPlot[walktwo, PlotJoined -> True];

Is a given expression a permutation?

First we sort the list. Let its length be n. If the list is a permutation, then after sorting it should look like {1,2,...,n}. So we compare it with Range[n]. Note that we do not need the explicit value n.

permutationQ[
        list_?VectorQ] := (Sort[list] == 
          Range[Length[list]])
permutationQ[{2, 4, 5, 1, 3}]
permutationQ[{2, 4, 5, 1, 3, 2}]
permutationQ[{2, 4, 5, 1, {3}}]

Define your own logarithmic function

Definition

log[x_*y_] := log[x] + log[y]
log[x_^n_] := n*log[x]

Test

log[a*b^3*z^g/c^5]

Pure functions, Map, and Apply

Pure functions

Usual definition of a function

Clear[g]
g[x_] := x^2
g[3]

The same using a pure function:

pf = Function[w, w^2]
pf[3]

A more economical notation. Observe that both symbols  # and & must appear in the definition!

(#^2 &)
(#^2 &)[3]

Pure functions of more than one argument:

(#1 + #2^2 &)[a, b]

Random walk using a pure function:

SeedRandom[13]
ListPlot[NestList[# + Random[] - 0.5 &, 0, 1000], 
        PlotJoined -> True];

More examples of the use of pure functions

Function Select[list, criterion] picks out the elements of list for which criterion  is True. 
We select even numbers of the given list:

Select[{a, 23, 12, 0, 3.5}, EvenQ]

For selecting numbers greater than 10 in another list we must use the pure function:

Select[{23, 6, 2, 62, 673}, #1 > 10 &]

From the pairs below we select those for which the second part is greater than 15:

mat3 = {Range[5], Range[5]^2} // Transpose
Select[mat3, #[[2]] > 15 &] // Transpose

We sort the rows of the following matrix by the second column in a descending order:

(mat4 = {Table[Random[], {5}], Table[Random[], {5}]} // Transpose) // 
      MatrixForm
Sort[mat4, #1[[2]] > #2[[2]] &] // MatrixForm

Map and Apply

Very important functions are Map and Apply. Look at the difference between Map and Apply:

Map[f, {a, b, c}]
Apply[f, {a, b, c}]

Plus[a, b, c]
Plus[{a, b, c}]
Apply[Plus, {a, b, c}]

Sometimes the result may be obtained without mapping, particularly in case of listable functions:

Map[Sin, {a, b, c}]
Sin[{a, b, c}]

This works without mapping since

Attributes[Sin]

show that Sin is a listable function. But

IntegerQ[{1.1, 2, π}]
Attributes[IntegerQ]

because IntegerQ is not a listable function. We must map this function instead:

Map[IntegerQ, {1.1, 2, π}]

A simple example:

Map[#^3 &, Range[5]]

Exercise

Think what will happen in the following examples:

data = {Range[5], Range[5]^2}

Apply[{#1, 2*#2} &, data]

data2 = {data[[1]], 2*data[[2]]}

Transpose[data]
Map[{#[[1]], 2*#[[2]]} &, data]
Map[{#[[1]], 2*#[[2]]} &, Transpose[data]] // 
      Transpose

Example

By the built-in function FactorInteger we get a list of prime factors together with their exponents in the form of a list {{(Subscript[factor, 1]),(Subscript[exponent, 1])},...,{...,...}}.

f1 = FactorInteger[11059400205764080766512]

We will check the result. Look, how to create (3^2) from the list {3,2}:

First[{3, 2}]^Last[{3, 2}]

Now we just map the above function over the list f1 and multiply the items of the resulting list:

Map[(First[#]^Last[#]) &, f1]
Apply[Times, Map[(First[#]^Last[#]) &, f1]]

Packages

Add-on packages

See Help Browser for the complete list of Add-on packages. The packages are available in the source code of Mathematica developed by the professionals so that they can serve as a good material for further study. Here we give some examples of the use of the statistical package.

Mean[Range[5]]

It was a mistake. The proper package has not been loaded and the function Mean is not known at the moment.

Needs["Statistics`Master`"]

Mean[Range[5]]

The symbol Mean must be removed. To be on the save side, load the package again.

Remove[Mean]
Needs["Statistics`Master`"]
Mean[Range[5]]

Note that loading the Master suffices in many cases. There are some exceptions, like Statistics`MultiDescriptiveStatistics`. In such cases, the more specific loading is necessary. In this case, by the command Needs["Statistics`MultiDescriptiveStatistics`"]. It is recommended to use "Needs" instead of  "Get" (equivalent to "<<").

SeedRandom[13]
data2 = Table[Random[], {100}];
data3 = Table[Random[], {100}];

{Mean[data2], Variance[data2], Median[data2]}

LocationReport[data2]
DispersionReport[data2]

Table[Quantile[data2, j/10], {j, 9}]

TrimmedMean[data2, {0.1, 0.1}]

{Skewness[data2], Kurtosis[data2]}

Correlation[data2, data3]

Variance (with the factor 1/(n-1)) and VarianceMLE (with the factor 1/n) of the data:

{Variance[data2], VarianceMLE[data2]}

Confidence intervals for the mean, default confidence level is 95 per cent.

MeanCI[data2]
MeanCI[data2, ConfidenceLevel -> 0.99]

The confidence interval for the difference in means, options available.

MeanDifferenceCI[data2, data3]
MeanDifferenceCI[data2, data3, 
      KnownVariance -> {1. /12, 1. /12}]

We can handle the distributions involved also in a symbolic way. Continuous distributions (PDF standing for the power density function, CDF standing for the cumulative distribution function).

It is a good idea to check the meaning of the parameters:

distr1 = 
        WeibullDistribution[α, β];
PDF[distr1, x]
CDF[distr1, x]

Quantile[distr1, y]

distr2 = 
        NormalDistribution[μ, σ];
{Mean[distr2], Variance[distr2], Kurtosis[distr2]}

distr2 is defined as N (μ, (σ^2)). In Mathematica, the second parameter is the standard deviation, not the variance!

PDF[distr2, x]

Symbolic expression of the quantile function of  N (μ, (σ^2)):

Quantile[distr2, y]

Our own definition of the cumulative distribution function of a N(0,1):

Φ[x_] := 
      0.5*(1 + Erf[x/(Sqrt[2])])
{Φ[-10], Φ[10]}

Of course, we can use predefined cumulative distribution function of a normal distribution:

Plot[Φ[x] - 
          CDF[NormalDistribution[0, 1], x], {x, -19, 19}, 
        PlotRange -> All];

{Skewness[distr2], Kurtosis[distr2]}

CharacteristicFunction[distr2, z]

Further possibilities

distr3 = NormalDistribution[0, 1];
{ExpectedValue[x^2, distr2, x], 
      ExpectedValue[#^2 &, distr2]}

Table[Random[distr3], {5}]

RandomArray[distr3, 5]

RandomArray[distr3, {5, 5}] // MatrixForm

distr4 = PoissonDistribution[5];
SeedRandom[13]
t4 = RandomArray[distr4, 100]
f4 = Frequencies[t4]
QuantileForm[t4];

Visualization of the data

Needs["Graphics`Master`"]

?BarChart

?Histogram

BarChart[f4];

Histogram[t4, HistogramCategories -> Max[t4], 
        ApproximateIntervals -> False];

Options[Histogram]

t5 = RandomArray[distr3, 100];
BinCounts[t5, {-3, 3, 0.25}]

Hypothesis testing

MeanTest[t5, 0, TwoSided -> True]

MeanTest // Options

Smoothing techniques

MovingAverage[{a, g, b, c, d, e}, 3]
MovingMedian[{a, g, b, c, d, e}, 3]

MovingAverage[t5, 20] // ListPlot;
MovingMedian[t5, 20] // ListPlot;

User's packages

If you have a collection of your own functions which you suppose to be frequently used in your further work, it is a good idea to create your own package. See the Skeleton.nb in the Mathematica directory how to create your own package. See also:
 Roman E. Maeder. Programming in Mathematica, 3rd ed. Addison-Wesley, 1996.

Here is only a sample package with two functions:

BeginPackage["Hurt`fmdm`"]
\(\(\(fmdm::usage = "fmdm.m contains two descriptive statisitcs."\)\(
\)
    \)\)
avg::usage = "avg[data] returns the arithmetic mean of the \
data"

    
\(\(\(sd::usage = "sd returns the standard deviation of the data"\)\(
\)
    \)\)
Begin["`Private`"]
Unprotect[avg, sd]
avg[
        data_] := (Apply[Plus, data]/Length[data] // 
            N) /; (VectorQ[data] && 
            Length[data] > 0)
sd[data_] := 
      avg[(data - avg[data])^2] // Sqrt
Protect[avg, sd]
End[]
EndPackage[]

Select this cell and save it by Menu -> Edit -> Save Selection As -> Package Format ... -> c:\Hurt\fmdm.m, e.g. Look at the file in some text editor. If it contains "(*" at the beginning and "*)", remove them, and save the file again, with the extension *.m. Afterwards (in our case), proceed by

AppendTo[$Path, "C:\\Hurt"]

This is not necessary if you activate (evaluate) your package during the same session of Mathematica. Otherwise, you can open the package, copy  and paste it into a current Mathematica session anytime. And activate it, of course. Then

Needs["Hurt`fmdm`"]

?avg

Financial applications

Measures of cash flows

Calendar

Calendar conventions play a crucial role in time value of money. There is a variety of calendar conventions for  financial instruments. Within Mathematica framework, they are included in Applications Packages like Financial Essentials or Derivatives Expert. They must be bought separately, however. In Mathematica Add-On Packages the only available calendar is Actual/Actual, i.e., actual number of days in a month and actual number of days in a year. For the illustration of the financial calculations we will make use of this one.

The actual date in the form {Year, Month, Day, Hour, Minute, Second}:

Date[]

We may extract just the date:

Date[]
Date[][[Range[3]]]

Needs["Miscellaneous`Calendar`"]

DayOfWeek[today = {2002, 6, 28}]

DaysBetween[today, {2010, 3, 28}]
DaysBetween[{2010, 3, 28}, today]

DaysPlus[today, 50]

Present Value

For a cash flow CF = ((Subscript[CF, 1]),  ... , (Subscript[CF, T]))and a  valuation interest rate r the formula for the present value reads:

PV(CF,r)=(∑_(t = 0)^T) ((Subscript[CF, t])/(((1 + r))^t))

In fact, this is the scalar product of the cash flow and the vector of discount factors. Implementation heavily uses listability:

cf = {cf0, cf1, cf2, cf3};
len = Length[cf]
disc = 1/(1 + r)^(Range[len] - 1)
cf . disc

pv1[cf_, r_] := 
      cf . (1/(1 + r)^(Range[Length[cf]] - 
                  1))
pv1[cf, r]

Modification for the payments going on given different dates and we evaluate the present value with respect to the startday:

dif = 
      Map[DaysBetween[#, today] &, {{2004, 3, 28}, {2003, 1, 13}, {2003, 1, 
            1}}]
discounts = 1/((1 + r)^(dif/360))

pv2[cf_, dates_, r_, startday_] := 
      Module[{diskonty, dif}, 
dif = 
          Map[DaysBetween[#, startday] &, 
            dates]; 
diskonty = 
          1/((1 + r)^(dif/360)); 
cf . 
          diskonty]

dates = {{2004, 3, 28}, {2003, 1, 13}, {2003, 1, 
            1}};
cf = {1000, -900, -50};
pv2[cf, dates, 0.1, today]

Plot[pv2[cf, dates, r, today], {r, 0, 0.1}];

IRR (Internal Rate of Return)

cf2 = {-1000, Table[100, {10}], 1100} // Flatten;

NSolve[pv1[cf2, r] == 0, r]

Cases[r /. NSolve[pv1[cf2, r] == 0, r], 
        z_Real /; z ≥ 0][[1]]

dates2 = 
      Map[DaysPlus[today, #] &, 
        Table[365*(j - 1), {j, Length[cf2]}]]

Plot[pv2[cf2, dates2, r, today], {r, 0.09, 0.11}];

FindRoot[pv2[cf2, dates2, r, today] == 0, {r, 0.01}]

Import and export

The following data  come from Finance Pack (Applications Library) Wolfram Research, 1994:

stream = 
      OpenRead["c:/My Documents/Vpfpm/Data/stocks1.ret"]
nn = Read[stream, Number]
ticker = Read[stream, Table[Word, {nn + 1}]]
(datareturns = ReadList[stream, Table[Number, {nn + 1}]]) // 
      TableForm
Close[stream]

returns = Map[Drop[#, 1] &, datareturns];

Needs["Statistics`MultiDescriptiveStatistics`"]

expectedreturns = Mean[returns]

(covreturns = CovarianceMatrix[returns]) // MatrixForm

Remark. How to get returns (Subscript[r, t]) from prices (Subscript[S, t]):

(Subscript[r, t]) = (((Subscript[S, t + 1]) - (Subscript[S, t]))/(Subscript[S, t]))

ss = {11, 12, 10, 12, 13, 14};
Drop[ss, 1] - Drop[ss, -1]
(Drop[ss, 1] - Drop[ss, -1])/Drop[ss, -1]

Markowitz portfolio

Simulation

ListPlot[NestList[(# + Random[] - 0.5) &, 10, 50], 
        PlotJoined -> True];

Needs["Statistics`Master`"]
tab1 = 
        Table[ListPlot[NestList[(# + Random[] - 0.5) &, 10, 50], 
            PlotJoined -> True, PlotRange -> All, 
            DisplayFunction -> 
              Identity], {100}];
tab1[[1]];

Show[tab1, DisplayFunction -> $DisplayFunction];

s0 = 100; μ = 0.05/365; σ = 
      0.003;
distr = 
        NormalDistribution[μ, σ];
tab2 = 
        Table[ListPlot[NestList[#*(1 + Random[distr]) &, s0, 360], 
            PlotJoined -> True, PlotRange -> All, 
            DisplayFunction -> 
              Identity], {30}];
Show[tab2, DisplayFunction -> $DisplayFunction];

Financial derivatives

Black-Scholes Formula

Notation: s ... current price of underlying asset at time t, k ... strike price, r ..  riskless interest rate, T ... expiry date, σ ... volatility of the underlying asset
call (put) ... the value of a CALL (PUT) option

Clear[Φ]
Φ[x_] := 
      (1/2) (1 + Erf[(x/(Sqrt[2]))])
ϕ[x_] := Exp[-x^2/2]/Sqrt[2*π]

d1 = (Log[s/k] + (r + σ^2/2)*(T - t))/(σ*
              Sqrt[T - t]);
d2 = (Log[s/k] + (r - σ^2/2)*(T - t))/(σ*
              Sqrt[T - t]);
call[s_, k_, σ_, T_, t_, r_] := 
      Module[{d1, d2}, 
        d1 = (Log[s/k] + (r + σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
d2 = (Log[
                  s/k] + (r - σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
s*Φ[d1] - 
          k*Exp[-r*(T - t)]*Φ[d2]]
put[s_, k_, σ_, T_, t_, r_] := 
      Module[{d1, d2}, 
        d1 = (Log[s/k] + (r + σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
d2 = (Log[
                  s/k] + (r - σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
k*
            Exp[-r*(T - t)]*Φ[-d2] - 
          s*Φ[-d1]]

Implied volatility: price ... the market price of an option at time t

impliedcall[s_, k_, T_, t_, r_, price_] := 
      sd /. FindRoot[
          call[s, k, sd, T, t, r] == price, {sd, 
            0.2}]
impliedput[s_, k_, T_, t_, r_, price_] := 
      sd /. FindRoot[
          put[s, k, sd, T, t, r] == price, {sd, 0.2}]

Real data: Underlying Deutsche Bank April 23, 1999 - Price EUR 52.44

optionDB46c[pos_, S_] = {pos, 1, 5.26, S, 
          46};

    optionDB48c[pos_, S_] = {pos, 1, 3.75, S, 48};

    optionDB50c[pos_, S_] = {pos, 1, 2.51, S, 50};

    optionDB46p[pos_, S_] = {pos, -1, 0.74, S, 
          46};

    optionDB48p[pos_, S_] = {pos, -1, 1.39, S, 
          48};

    optionDB50p[pos_, S_] = {pos, -1, 2.33, S, 50};

call[52.44, 50, 0.2, 20/52, 0, 0.05]
put[49, 50, 0.2, 20/52, 0, 0.05]

Plot[call[s, 46, 0.3, 28/360, 0, 0.05], {s, 40, 53}, 
        PlotRange -> All];
Plot[call[s, 50, 0.3, 28/360, 0, 0.05], {s, 40, 53}, 
        PlotRange -> All];

Plot[{call[  s, 50, 0.2, 22/52, 0, 
            0.05], put[s, 50, 0.2, 20/52, 0, 0.05]}, {s, 0, 70}, 
        PlotRange -> All];

Real data: Wolkswagen

{impliedcall[70.72, k1 = 67.5, T = 58/360, 
        t = 0, r = 0.05, c1 = 6.31], 
impliedcall[70.72, 
        k2 = 70, 58/360, 0, 0.05, c2 = 4.92], 
impliedcall[
        70.72, k3 = 72.5, 58/360, 0, 0.05, 
        c3 = 3.37], 
impliedput[70.72, 67.5, 58/360, 0, 
        0.05, p1 = 2.92], 
impliedput[70.72, 70, 58/360, 0, 
        0.05, p2 = 4.08], 
impliedput[70.72, 72.5, 58/360, 
        0, 0.05, p3 = 5.48]}

Portfolio of options

Butterfly spread

valuepor[s_] := -c1 - p3 + c2 + p2 + Max[s - k1, 0] + 
        Max[k3 - s, 0] - Max[s - k2, 0] - 
        Max[k2 - s, 0]
Plot[valuepor[s], {s, 65, 76}];

The Greeks (sensitivities of the option price on the parameters

Clear[s, k, σ, T, t, r]
delta = D[call[s, k, σ, T, t, r], s] // PowerExpand // 
        Simplify

!!! Warning !!! Time consuming calculation ahead:

theta = D[call[s, k, σ, T, t, r], t] // PowerExpand // 
        FullSimplify

rho = D[put[s, k, σ, T, t, r], r] // PowerExpand // 
        Simplify

The following calculation takes  cca 2 minutes.

vegac = D[call[s, k, σ, T, t, r], σ] // 
          PowerExpand // FullSimplify

Binomial trees

In this Section we closely follow Chapter 18, Shaw, W.: Modelling Financial Derivatives with Mathematica.

Introduction

Let us consider a time interval (t, t + Δt). Under the Black-Scholes model the conditional distribution of (Subscript[S, t + Δt]) given (Subscript[S, t]) is log-normal \(TraditionalForm\`LN(
        ln\ S\_t + \((μ - 
                1\/σ\^2)\) Δt, \(σ\^2\) Δt)\). In a risk neutral world μ = r, where r is the riskless rate and we have

E(
          (Subscript[S, t + Δt]) | (Subscript[S, t])) = (Subscript[S, t]) 
            (e^(r Δt)) = :(Subscript[E, c])     and     var(
          (Subscript[S, t + Δt]) | 
            (Subscript[S, t])) = ((Subscript[S, t])^2) (e^(2  r Δt))(
              (e^((σ^2) Δt)) - 
                1) = :(Subscript[V, c]).

Under the discrete binomial random walk we suppose that the stock price (Subscript[S, t]) moves to one of two new values, (Subscript[S, t + Δt]) = (Subscript[S, t]) u and (Subscript[S, t + Δt]) = (Subscript[S, t]) d where d < 1 < u with probabilities p and q := 1 - p, respectively. In this case we have

E(
          (Subscript[S, t + Δt]) | (Subscript[S, t])) = (Subscript[S, t])(
            p u  +  q d) = :(Subscript[E, b])    and     var(
          (Subscript[S, t + Δt]) | (Subscript[S, t])) = ((Subscript[S, t])^2) 
            p ((q(u - d))^2) = (Subscript[V, b]).

While in the continuous model we have two parameters r and (σ^2), ceteris paribus, in the discrete model we have three parameters, i.e., u,  d,  p. In practice we choose the parameters of the binomial model to match the parameters of the continuous one in terms of the moments:

(Subscript[E, b]) = (Subscript[E, c])     and     (Subscript[V, b]) = 
          (Subscript[V, c]).

The reason is simple; usually we can estimate r and (σ^2) directly which is hardly possible for the binomial model. Since there is one more degree of freedom, we have to postulate a reasonable relationship between the parameters involved. The corresponding system of equations now reads:

Clear[p, u, d, r, σ]
eqs = {p*u + (1 - p)*d == Exp[r*Δt], 
          p*(1 - p)*(((u - d))^2) == 
            Exp[2*r*Δt]*(Exp[(σ^2) Δt] - 1)};

Jarrow-Rudd scheme

Under this scheme it is supposed that both up and down movements have the same probability.

Clear[p, u, d, r, σ]
p = (1/2); solJR = Solve[eqs, {u, d}]

Only the largest solution makes sense:

{uJR, dJR} = {u, d} /. solJR[[2]]

This is an approximation:

approxJR = {Series[uJR, {Δt, 0, 1}] // Normal // 
          PowerExpand, 
Series[
              dJR, {Δt, 0, 1}] // Normal // 
          PowerExpand}

which is very close to frequently used

Exp[approxJR]/E

Cox-Ross-Rubinstein scheme

Under this scheme it is supposed that u = 1/d:

Clear[p, u, d, r, σ]
u = (1/d); 
    solCRR = Solve[eqs, {p, d}] // PowerExpand // 
        FullSimplify

After a hand-made simplification (not too difficult) we get

Clear[CRRupBin];
CRRupBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g + Sqrt[(f+g)^2-4])];
Clear[CRRdownBin];
CRRdownBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g - Sqrt[(f+g)^2-4])];
Discount[n_, r_, T_] := N[Exp[-r T/n]];

Binomial tree for American options

Let us consider an n steps' evolution, time to expiry T and the time step Δt = T/n. Denote (i, j) the jth node at time iΔt, i = 1,  ... , n,  j = 0,  ... , i. Let (Subscript[V, i, j]) be the value of the option (either CALL or PUT) at node (i, j). The price of the stock is (Subscript[S, 0]) (u^j) (d^(i - j)) at that node. Since the value for an American PUT is max(K - (Subscript[S, T]), 0) at expiry, the value is

(Subscript[V, n, j]) = 
        max(K - (Subscript[S, 0]) (u^j) (d^(n - j)), 0),      j = 0, 
      1,  ... , n.

With probability p the system moves from node (i, j) at time iΔt to (i + 1, j + 1) at time (i + 1) Δt and with probability 1 - p to (i + 1, j). Without possibility of an early exercise we get

(Subscript[V, i, j]) = (e^(-r Δt))(
          p (Subscript[V, i + 1, j + 1]) + (1 - p) (Subscript[V, i + 1, j])),    i = 
        0,  ... , n - 1,    j = 0,  ... , i.

If an early exercise is taken into account, the value must be compared with its intrinsic value. Thus we finally get

(Subscript[V, i, j]) = 
        max {K - (Subscript[S, 0]) (u^j) 
                (d^(i - j)), (e^(-r Δt))(
              p (Subscript[V, i + 1, j + 1]) + (1 - p) 
                  (Subscript[V, i + 1, j]))},    i = 0,  ... , 
      n - 1,    j = 0,  ... , i.

Starting from level n we go backward through the tree to the desired value (Subscript[V, 0, 0]).

Implementation

Here also a dividend yield q is considered.

$RecursionLimit=5000;

afunc[n_,  r_,  q_,  T_]  :=  
        N[Exp[(r - q) T/n]];
Pfunc[up_, down_, a_, disc_] := 
        N[disc*(a - down)/(up - down)];
Qfunc[up_, down_, a_, disc_] := 
        N[disc - Pfunc[up, down, a, disc]];

Clear[CRRAmericanOption]
CRRAmericanOption[S_, n_, σ_, T_, r_, q_, exercise_Function] :=
Module[{u = CRRupBin[n,r,q,σ,T], d = CRRdownBin[n,r,q,σ,T],
a = afunc[n,r,q,T],
disc = Discount[n,r,T],
P,Q,TreeOp,value,delta,gamma,theta},
P = Pfunc[u,d,a,disc];
Q = Qfunc[u,d,a,disc];
TreeOp[node_,level_] := TreeOp[node,level] = 
If[level==n,exercise[S d^node u^(level-node)],
Max[{P,Q}.{TreeOp[node,level+1],TreeOp[node+1,level+1]},
exercise[S d^node u^(level-node)]]];
value=TreeOp[0,0];
delta = (TreeOp[0,1] - TreeOp[1,1])/((u-d)*S);
gamma = ((TreeOp[0,2] - TreeOp[1,2])/(u^2-1) - 
(TreeOp[1,2]-TreeOp[2,2])/(1-d^2))*2/(S^2 (u^2-d^2));
theta = (TreeOp[1,2] - TreeOp[0,0])/(2*T/n);
Clear[TreeOp];
{value, delta, gamma, theta}];

CRRAmericanPut[S_, K_, n_, r_, q_, σ_, T_] := 
        CRRAmericanOption[S, n, σ, T, r, q, 
          Max[K - #, 0] &];
CRRAmericanCall[S_, K_, n_, r_, q_, σ_, T_] := 
        CRRAmericanOption[S, n, σ, T, r, q, 
          Max[# - K, 0] &];

Recall the meaning of the parameters:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanCall[9,  10, 10,  0.1,  0,  0.2,  
      1]
call[9, 10, 0.2, 1, 0, 0.1]

Look at the convergence:

Table[{m, CRRAmericanCall[9,  10, m,  0.1,  0,  0.2,  
              1][[1]] // Evaluate}, {m, 3, 8}]

The following plot (number of steps equal to 80) shows that the convergence is not too convincing and that the recommended number of steps 30 may not be sufficient.

ListPlot[
        Table[{m, CRRAmericanCall[9,  10, m,  0.1,  0,  0.2,  
                  1][[1]] // Evaluate}, {m, 3, 80}], 
        PlotJoined -> True];

The following Example comes from HULL, J. C.: (2000) Options, Futures, and Other Derivative Securities, p.391, where the approximate values u = (e^(σ (Sqrt[ Δ])t)), d = 1/u are used. His result (after 5 steps) is 4.49. A comparison with usual BS formula for an European CALL is also given:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanPut[50, 50, 5, 0.1, 0, 0.4, 
      0.4167]
put[50, 50, 0.4, 0.4167, 0, 0.1]

Yield Curve

ir3 = {{7, 4.74`}, {14, 4.72833333333333`}, {30, 4.705`}, {60, 
            4.595`}, {90, 4.53`}, {180, 4.42`}, {270, 4.415`}, {360, 
            4.41`}, {730, 4.51`}, {1095, 4.68`}, {1460, 4.85`}, {1825, 
            4.99`}, {2190, 5.09`}, {2555, 5.19`}, {2920, 5.27`}, {3285, 
            5.33`}, {3650, 5.38`}, {4015, 5.42`}, {4380, 5.46`}, {4745, 
            5.48666666666667`}, {5110, 5.51333333333333`}, {5475, 
            5.54`}};

p2ir3 = 
        ListPlot[ir3  ];

Clear[b0, b1, b2, c]
fit3 = 
      NonlinearFit[ir3, 
        b0 + (b1 + b2)*(1 - Exp[-t/c])/(t/c) - 
          b2*Exp[-t/c], t, {b0, b1, b2, c}]

p1ir3 = Plot[fit3, {t, 0, 6230}];

Show[p1ir3, p2ir3];