VPFPMLS2005

Lecture 1

Reasons

4 pillars

Numerical calculations

1 + 1

(2^100)

π

N[π, 20]

Symbolic calculations

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

Factor[e1]
Factor[x + (y^4) + (x^3) (y^6)]
TeXForm[Factor[e1]/(x + 2  y)]

Graphics

Plot[Sin[x + (x^2)], {x, 0, π}]

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

x* y
xx

?Global`*

Programming

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

Factor[h[e1]]

?h

Clear[h]
?h

Remove[h]
?h

Hilbert matrix

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

h3m = h[3] // MatrixForm

Inverse[h[3]]
Inverse[h3m]

?h3m

Inverse[h[3]] . h[3]

{h[3][[3, 
        1]], h[
            3][[3]][[1]]}

h[3][[All, 2]]

Technicalities, basic input

Attributes[Eigensystem]

Attributes[h]

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

Závorky [ ], ( ), { },[[ ]]

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

x = 5;

h15 = h[15];

h[15][[6, 6]]

xx
+yy

xx + 
yy

Clear[x, y, xy];
x = 5; y = 6;

{xy, x*y, x y}

{x, y}

EUR =  30 CZK
7 EUR

l1 = {3.14, 56.12}

e3 = (e2 = {"14", "d"}) // ToExpression

FullForm[e2]
FullForm[e3]

e2[[1]]^2

e2[[1]]^2 // FullForm

e3[[1]]^2

Range["a", "d"]

l1 // FullForm

((2^30)/(3^45))

((2^30)/(3^45)) // N

((2^30.) /(3^45))

Eigenvalues[h[2]] // FullSimplify

Eigenvalues[h[4] // N]

Det[h[10]]

Det[h[10] // N]

Det[h[10]] // N

Timing[Det[h[60]]] // N
Timing[Det[h[60] // N]]

Lecture 2

Numbers and mathematical notation, ...

Real (π,e),rational (23/27),complex (1+2*I),integers (2345)

Mějme integrál (∫_a^b) 
              (x^2) dx nebo ...

ev2 = Eigenvalues[h[20] // N]

ev2 // Chop

1.1 + (10^(-14)) I // Chop

Capabilities of the system

Timing[Eigenvalues[h[3] // N]]

t1 = Timing[Eigenvalues[h[50] // N]];

{t1[[2, 1]], 
      t1[[1]]}

t1[[1]]/Second

tab1 = 
      Table[Timing[
            Eigenvalues[
              h[j] // N]][[1]], {j, 
          40, 65, 5}]

ListPlot[Range[5]^2 + 100]

Options[ListPlot]

ListPlot[Range[5]^2 + 100, PlotJoined -> True]

ListPlot[tab1, PlotJoined -> True]

ListPlot[tab1/Second, PlotJoined -> True]

Attributes[Sin]

Sin[(π/4) Range[7]]

Sin[(π/4) Range[7]] // N

inv1 = Inverse[h[5]]

inv1*h[5] == IdentityMatrix[5]

(inv1*h[5])[[1]]

inv1 . h[5] == IdentityMatrix[5]

Inverse[h[5] // N] . h[5] == IdentityMatrix[5]

Chop[Inverse[h[5] // N] . h[5]] == IdentityMatrix[5]

Chop[Inverse[h[5] // N] . h[5]]

(Chop[Inverse[h[5] // N] . h[5]] // 
          N) == (IdentityMatrix[5] // N)

IdentityMatrix[5] // N

(Chop[Inverse[h[5] // N] . h[5]] // N)

Solving equations

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

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

Solve[(E^x) == x, x]

Solve[(E^x) == Sin[x], x]

Solve[x^6 + 1 == 0, x]

Solve[(x^6) + (x^5) + 1 == 0, x]

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

Curiosity

mat1 = Table[Random[], {200}, {200}];

ListPlot[e2 = Eigenvalues[mat1]]

Take[e2, {1, 5}]

ListPlot[e2 = Abs[Eigenvalues[mat1]], PlotRange -> All]

e2[[1]]

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

Simplify[f1]

i1 = Integrate[(Sqrt[x]) ArcTan[x], x]

d1 = D[i1, x]

d1 // Simplify

Numerics

NSolve[(x^6) + (x^5) + 1 == 0, x]

Plot[(E^x) - Sin[x], {x, -10, 1}]

NSolve[eq1 = (E^x) == Sin[x], x]

Table[FindRoot[eq1, {x, j}], {j, -15, 1}]

FindRoot[eq1, {x, 0}]

(E^x) - 
        Sin[x] /. {x -> -131.94689145077132`}

HoldForm[Sin[0]] === 0

{a, b}
a == b
a === b

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

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

Assuming[
      x ∈ 
        Reals, ∫Sign[
                x]/(x^(1/3)) \[DifferentialD]x]

(∫_0^∞) Exp[a x] \[DifferentialD]x

Integrate[(E^(a x)), {x, 0, ∞}, 
      Assumptions -> Re[a] < 0]

s1 = Series[f2, {x, 0, 10}]

Series[Gamma[n], {n, ∞, 6}]

s1*s1

n1 = Normal[s1]

i4 = Integrate[n1, {x, 0, 2}]
i3 - i4

i5 = NIntegrate[n1, {x, 0, 2}]
i4 - i5

i6 = Integrate[n1, x]

i6 /. x -> 0

i60 = Limit[i6, x -> 0]
i62 = i6 /. x -> 2
i62 - i60

fa1 = FactorInteger[(2^100) - 1]

fa1tr = Transpose[fa1]

Range[3]^Range[3]

pow = 
      fa1tr[[1]]^
        fa1tr[[2]]

fa3 = Apply[Times, pow] // FullSimplify

{(2^100) - 1 == fa3, (2^100) - 1 === fa3}

mat2 = Table[(x^i) + (y^j), {i, 3}, {j, 4}]

MatrixRank[mat2]

nu2 = NullSpace[mat2]

mat2 . Transpose[nu2] // Simplify

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

FullForm[sol1]

Plot[
      x[t] /. sol1[[1]], {t, 0, 
        50}]

ParametricPlot[{x[t], x'[t]} /. sol1, {t, 0, 50}]

DSolve[{x''[t] + x[t]^3 == Sin[x[t]], 
        x[0] == x'[0] == 0}, x[t], t]

Integrate[Log[x] Exp[-(x^3)], x]

i7 = 
      Integrate[
        Log[x] Exp[-(x^3)], {x, 0, Infinity}]
i7 // N

EulerGamma // N

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

Table[{j, (∑_(k = 1)^∞) (1/(k^j))}, {j, 6}] // 
      TableForm

{ex1 = (∏_(n = 1)^5) Gamma[2  (n/5)], ex1 // Simplify, 
      ex1 // FullSimplify}

Simplify[Cos[x]^4]

TrigReduce[Cos[x]^4]

{Prime[4], Prime[(10^9)]}

{PrimeQ[22801763489], PrimeQ[22801763489 - 1]}

Lecture 3

Clear[b];
e1 = Eigenvalues[{{3, 1}, {2, b}}]

Plot[e1, {b, -10, 10}]

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

sol1 = Solve[First[e1] == 0, b]

b1 = sol1

b1 = b /. sol1

(b1 = 
          b /. sol1)[[1]]

Integrate[e1 // First, b] // FullSimplify

?b
?b1

t1 = Table[Prime[i], {i, 40}];
p1 = ListPlot[t1];
Last[t1]

fit1 = Fit[t1, {Log[x], x, (x^2)}, x]
p2 = Plot[fit1, {x, 1, 40}];

Show[p1, p2];

{Sum[fit1, {x, 40}], Apply[Plus, t1]}

{fit1 /. x -> 1, fit1[1]}

Elementary symbolic calculations

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

e13 // FullSimplify

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

e14

e13 == e14
e13 == e14 // Simplify

{Timing[D[e13, {x, 140}]] // First, 
      Timing[D[e14, {x, 140}]] // First}

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

Coefficient[yx, y^2]
CoefficientList[yx, y]

Substitution

{yx /. x -> y, yx /. x -> 2  y, yx /. x -> 3}

{x, y}

{(Sin^2)[x], Sin^2[x]}

Sin[x]^2 // TraditionalForm

i2 = (∫_(-∞)^z) 
              (1/(Sqrt[2  π])) (E^(-((x^2)/2))) \[DifferentialD]x

i2=Φ(x)

Multiple integrals, beware of the integration order!!

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

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

Integration of an ugly function

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

i0 the exact value, i1 its numerical value, i2 the result of a numerical integration

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

Options[NIntegrate]

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

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

Solving equations

p1 = Expand[(((x + 1))^9)]

Solve[p1 == 0, x]

Solve[yx == 0, x]

Clear[a, b, c];
sol2 = 
      Solve[{x + y == 1, a*x + b*y + c == 0}, {x, y}]

sol3 = Flatten[sol2]

x + y /. sol2 // Simplify
x + y /. sol3 // Simplify

Nonlinear equations

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

NSolve[eq3, {x, y}]

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

Lists

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

li1 // TreeForm

li1 // FullForm

t4 = Table[n!, {n, 1, 10}]

Attributes[Factorial]

Factorial[Range[10]]

Range[2, 7]

Range[0, 1, 0.2]

GCD[56, 74]

fun1[i_, j_] := 
      If[GCD[i, j] == 1, 1, 0]
fun2[i_, j_] := 
      Which[
GCD[i, j] == 1, 
        1, 
True, 0]

array1[n_] := 
      Table[fun1[i, j], {i, n}, {j, n}]
array2[n_] := Table[fun2[i, j], {i, n}, {j, n}]

array1[4] == array2[4]

ListDensityPlot[array1[5]];

ListDensityPlot[array1[50]]

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

Range[5]*Range[5]
Range[5]^2
Range[5] . Range[5]

xx = {(Subscript[x, 1]), (Subscript[x, 2])}; yy = {(Subscript[y, 1]), (Subscript[y, 2])};
xx . yy

Transpose[xx]

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

xxx = {xx}
yyy = {yy}

Transpose[xxx] . yyy // MatrixForm

Outer[Times, xx, yy] // MatrixForm

{f1, f2}

Clear[g1, g2, g3];
Outer[D, {g1[x1, x2, x3], g2[x1, x2, x3], g3[x1, x2, x3]}, {x1, x2}] // 
      MatrixForm

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

{Apply[Times, e10], Det[tab10]}

ListPlot[Transpose[{r10 = Re[e10], i10 = Im[e10]}], 
        PlotStyle -> PointSize[0.02]];

ListPlot[Transpose[{r10 = Re[e10], i10 = Im[e10]}], 
        PlotStyle -> PointSize[0.02], 
        PlotRange -> All];

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

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

pos10 = Position[r10, Max[r10]]

{e10 // First, 
      e10[[pos10 // First // 
          First]]}

Take[e10, 3]

Lecture 4

e10 = Eigenvalues[tab10 = Table[Random[], {100}, {100}]];
{Apply[Times, e10], Det[tab10]}
ListPlot[Transpose[{r10 = Re[e10], i10 = Im[e10]}], 
        PlotStyle -> PointSize[0.02]];
ListPlot[Transpose[{r10 = Re[e10], i10 = Im[e10]}], 
        PlotStyle -> PointSize[0.02], 
        PlotRange -> All];
ListPlot[Transpose[{r10 = Re[e10], i10 = Im[e10]}], 
        PlotStyle -> PointSize[0.05], 
        PlotRange -> {{49, 51}, Automatic}];
{Max[r10], Min[r10], Max[i10], Min[i10]}
pos10 = Position[r10, Max[r10]]
{e10 // First, 
      e10[[pos10 // First // 
          First]]}

r10

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

Take[Sort[r10, Greater], 10]

Listablity may improve efficiency:

tim1 = 
        Timing[Table[
            Sin[x], {x, 
              0, π, (π/10000)}]];
tim2 = 
        Timing[Sin[
            Range[0, π, (π/10000)]]];
{tim1[[1]], 
      tim2[[1]]}
mp1 = Map[First, {tim1, tim2}]
mp2 = Apply[First, {tim1, tim2}]

Head[{tim1, tim2}]

mp2 = Apply[First, {{tim1, tim2}}] // Short

Map[Sin, mp1]
Map[#1^2 &, mp1]

Attributes[Second]

Clear[a, b];
{a, b} // FullForm

Map[Sin, {a, b}]

Apply[Plus, {a, b}]
Head[a*b]
a*b // FullForm
Apply[Plus, a*b]
Head[a^b]
Head[2/3]
Head[2. /3]
{Head[2/3] == Rational, Head[2. /3] == Rational, 
      Head[2. /3] === Rational}

x = 5;
Clear[x1, y1];
hu[x_, y_] := x + y
hu[{x1, y1}]
Apply[hu, {x1, y1}]

User-defined functions

Next prime

?*Q

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

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

Attributes[PrimeQ]
PrimeQ[Range[6]]

?Module

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

{nextprime[100000], nextprime[{1, 2}]}
{nextprime2[100000], nextprime2[{1, 2}]}

?nextprime2

?k

?x

Mean, Variance

Mean[{a, b, c}]

VectorQ[{}]

mean1[list_] := 
      Apply[Plus, list]/Length[list]
mean2[list_] := 
      Apply[Plus, list]/Length[list] /; 
        VectorQ[list] ∧ Length[list] > 0

{mean1[{}], mean1[{1}, 3, 4], mean1[Range[5]]}

{mean2[{}], mean2[{1}, 3, 4], mean2[Range[5]]}

\(mean3[list_] := 
      Module[{n}, n = Length[list]; 
        If[VectorQ[list] && n > 0, Apply[Plus, list]/n, 
          Print["Argument is not a nonempty vector"]]\ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

{mean3[{}], mean3[{1}, 3, 4], mean3[Range[5]]}

mean3[nextprime[Range[100000, 110000]]]

Attributes[nextprime] = Listable;

Attributes[nextprime]

nextprime[{20, 30}]

Map[
      Timing, {mean3[nextprime[Range[100000, 140000]]], 
        mean3[nextprime[Range[100000, 140000]] // N]}]

nextprime[10. ]

{PrimeQ[11. ], PrimeQ[11]}

var3[list_] := (((list - mean3[list]))^2) // mean3

var3[Range[5]]

{var3[{a, b}], var3[{a, b}] // Simplify}

Random walk

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

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

randomwalk[10]

ListPlot[randomwalk[10000], PlotJoined -> True];

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

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

Lecture 5

Random walk I

One-dimensional random walk

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

randomwalk[10]

ListPlot[randomwalk[10000], PlotJoined -> True];

Two-dimensional random walk

Table[{Random[], Random[]}, {3}] - 0.5

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

randomwalk2[3]

ListPlot[randomwalk2[10000], PlotJoined -> True]

Random walk II

Delayed evaluation

addone[x_] := x + Random[] -  .5

Immediate evaluation

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

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

Another useful function

Clear[x, a, b];
NestList[Sin, 
      x, 3]
NestList[Times, x, 
      3]
NestList[Times, x, 3] // 
      FullForm
NestList[Times, x, 3] // 
      Head
NestList[Times, x, 
        3][[0]]
Times[x]
Times[x, 
      x]
FullForm[a*b*c]

walkone[n_] := NestList[addone, 0, n]
walktwo[n_] := NestList[addtwo, 0, n]

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

Is a given expression a permutation?

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

permutationQ2[
        list_?VectorQ] := (Sort[
            list] == (Range[Length[list]] + Min[list] - 
              1))
permutationQ2[{5, 8, 6, 7}]
permutationQ2[{5, 8, 6, 7, 3}]

Needs["DiscreteMath`Combinatorica`"]

Options[PermutationQ]

PermutationQ[{5, 8, 6, 7}]
PermutationQ[{3, 1, 2}]

Define your own logarithmic function

log[x_*y_] := log[x] + log[y]
log[x_^n_] := n*log[x]
log[a*b^3*z^g/c^5]

?log

log[1] = 0; log[E] = 1;

?log

log[(E^2) a*b^3*z^g/c^5]

log = 6

?log

Pure functions, Map, and Aply

Usual definition of a function

Clear[g];
g[x_] := (x^2)
g[3]

The same, using the pure function

pf = Function[x, (x^2)]
pf[3]
Function[x, (x^2)][3]

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

Clear[a, b, c, 
        d];

    sol = Solve[(x^4) + (x^3) + (x^2) + 1 == 0, x]

Random walk revisited

walk3[n_] := 
      NestList[(# + Random[] -  .5) &, 0, n]
ListPlot[walk3[10000]];

Other examples using pure functions

Function Select[list, criterion] returns the elements of list satisfying the criterion:

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

Union[{a, b, c}, {f, e, d, c, c}]

Select all elements of the list greater than 10

Select[Range[14], #1 > 10]
Select[Range[14], #1 > 10 &]

From the pairs below select all with the second part greater than 15

mat15 = {Range[5], Range[5]^2} // MatrixForm

FullForm[mat15]

mat15 = Apply[Identity, mat15]

mat15 . (mat15tr = Transpose[mat15])

Select[
      mat15tr, #[[2]] > 15 &]

Sort the rows of a matrix by the second column in the descending order

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

Lecture 6

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

Clear[hh, x, y];
hh[x_, y_] := 
      x + y
hh[{x, 
        y}]
hh[{x, y}, {y, 
        z}]
Apply[hh, {x, y}]

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

log[{(x^2), y*z}]

SetAttributes[log, Listable]
log[{(x^2), y*z}]

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

Attributes[IntegerQ]

IntegerQ[{1, 1/2}]
Map[IntegerQ, {1, 1/2}]

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

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

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

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

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

f1 = 
      FactorInteger[
        c1 = 1132420058460472342*2*16699]
c1

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

Map[First[#]^Last[#] &, f1]

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

f2 = Reverse[f1 // Transpose] // Transpose

Map[#[[1]] &, f1]^
      Map[#[[1]] &, f2]

Packages

Standard Add-On Packages

?LocationReport

Needs["Statistics`Master"]

Needs["Statistics`"]

?LocationReport

l1 = LocationReport[r1 = Range[10]]

m1 = Mean /. l1

{Mean[r1], Variance[r1], VarianceMLE[r1], Quantile[r1, 0.5], 
      InterpolatedQuantile[r1, 0.5], Median[r1]}

DispersionReport[r1]

Seed[13]
r2 = Table[Random[], {(10^6)}];
LocationReport[r2]

Options[MeanCI]

MeanCI[r2]
MeanCI[r2, KnownVariance -> (1/12)]

{Skewness[r2], Kurtosis[r2]}

Kurtosis[du = UniformDistribution[0, 1]] // N

Clear[a, b];
duab = UniformDistribution[a, b]
PDF[duab, x]

Mean[duab]

?ExpectedValue

ExpectedValue[(x^4), duab, x] // Simplify

ExpectedValue[(Sqrt[x]), duab, x]

\({ExpectedValue[\@x, du, x], ExpectedValue[\@# &, du]}\)

dg = GammaDistribution[α, β];

PDF[dg, x]

CharacteristicFunction[dno = NormalDistribution[μ, σ], 
      t]

\(ExpectedValue[#\^10 &, dno]\)

RandomArray[du, {5}]
RandomArray[du, {2, 3, 4}]

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

Visualization of the data

Needs["Graphics`Master"]

Needs["Graphics`"]

?BarChart

BarChart[f4];

Histogram[t4]

r3 = RandomArray[du, {500}];
BinCounts[r3, {0, 1, 0.1}]

Smoothing techniques

l3 = CharacterRange["a", "e"]
FullForm[l3]

MovingAverage[l4 = {a, b, c, d, e}, 3]

FullForm[l4]

(l3[[1]]*
              l3[[2]] // Sin // 
          FullForm) /. a -> 3
(l3[[1]]*
              l3[[2]] // Sin // 
          FullForm) /. "a" -> 3

MovingAverage[r3, 50] // ListPlot

MovingMedian[r3, 20] // ListPlot

User's package

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

l4 = {3, 4};
AppendTo[l4, 5]

$Path

AppendTo[$Path, "C:\Dokumenty\VPFPM"]

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

avg[{a, b, c}]

Financial applications

Calendar

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

Needs["Miscellaneous`Calendar`"]

DayOfWeek[
      today = Date[][[
          Range[3]]]]

{DaysBetween[today, {2005, 6, 31}], 
DaysBetween[
        today, {2005, 7, 1}], 
DaysBetween[{2005, 7, 1}, 
        today]}

DaysPlus[today, 50]

Present value

cf = {cf0, cf1, cf2, cf3};

discountfactors = 1/(1 + r)^(Range[Length[cf]] - 1)

cf . discountfactors

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

Lecture 7

FAQ

s1 = Sum[(v^t), {t, n}]

Solve[pv == Sum[cf (v^t), {t, 0, n}], cf] // Simplify

Continuing

today = Date[][[{1, 2, 
          3}]]

DaysBetween[today, {2005, 4, 21}]

Needs["Miscellaneous`Calendar`"]

DaysBetween[today, {2005, 4, 21}]

Remove["DaysBetween"]

DaysBetween[today, {2005, 4, 21}]

Present value revisited

dates = Map[DaysPlus[today, #] &, 50*Range[5]]

Attributes[DaysBetween]

Map[DaysBetween[today, #] &, dates]

discountfactors = 
      1/(1 + r)^(Map[DaysBetween[today, #] &, dates]/365)

Clear[pv2]; 
    pv2[cashflow_, startdate_, dates_, r_] := 
      cashflow . (1/(1 + r)^(Map[DaysBetween[startdate, #] &, dates]/
                  365))

?pv2

pv2[{cf1, cf2, cf3, cf4, cf5}, today, dates, r]

Clear[pv3];
\(pv3[cashflow_, startdate_, dates_, 
        r_] := 
Module[
{discountfactors, 
          dif}, 
dif = 
          Map[DaysBetween[startdate, #] &, 
            dates]; 
discountfactors = 
          1/\((1 + r)\)^\((dif/365)\); 
cashflow . 
          discountfactors
	]\)

pv3[{cf1, cf2, cf3, cf4, cf5}, today, dates, r]

NestList[Function[x, x/(1 + r)], c, 4]

FoldList[Function[{x, r}, (x/(1 + r))], c, {r1, r2, r3}]

Plot[
        pv2[cf = {1000, -900, -50, -20, -40}, today, dates, 
          r], {r, 0, 0.2}];

IRR

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

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

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

?Cases

r /. NSolve[pv[cf2, r] == 0, r]

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

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

?pv2

NSolve[pv2[cf2, today, dates2, r] == 0, r]

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

NSolve[pv2[cf2, today, dates3, r] == 0, r]

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

FindRoot[pv2[cf2, today, dates3, r] == r, {r, -2}]

Table[
      FindRoot[pv2[cf2, today, dates3, r] == r, {r, 
          0.1*(j - 1)}], {j, 10}]

Import and Export

$ImportFormats

stream = 
      OpenRead["C:\Dokumenty\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]

Drop[{a, b, c}, 1]

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

Expected returns and covariance matrix

Needs["Statistics`MultiDescriptiveStatistics`"]

expectedreturns = Mean[returns]

(covmatrix = CovarianceMatrix[returns]) // MatrixForm

Returns based on prices

Clear[ret];
ret[
        ss_] := ((Drop[ss, 1] - Drop[ss, -1])/Drop)[ss, -1]
Clear[s1, s2, s3, s4, 
        s5];
ss = {s1, s2, s3, s4, 
          s5};
ret[ss]

?Mod

Clear[x];
x[n_] := x[n] = Mod[11*x[n - 1], (2^5)]
x[0] = 3;
Table[x[n], {n, 20}]

x[1]

Needs["Statistics`"]
Needs["Graphics`"]

?Histogram

tab7 = Table[Sum[Random[], {12}] - 6, {100000}];

{Mean[tab7], Variance[tab7]}

Histogram[tab7]

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

tab8 = 
        Table[ListPlot[NestList[# + Random[] - 0.5 &, 10, 50], 
            PlotJoined -> True], {3}];

Show[tab8]

tab9 = tab8 = 
          Table[ListPlot[NestList[# + Random[] - 0.5 &, 10, 50], 
              PlotJoined -> True, PlotRange -> All, 
              DisplayFunction -> Identity], {50}];

Show[tab9, DisplayFunction -> $DisplayFunction]

Show[
        tab8 = Table[
            ListPlot[NestList[# + Random[] - 0.5 &, 10, 50], 
              PlotJoined -> True], {3}]];

Lecture 8

Financial Derivatives (Options)

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]]

?Put

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}]

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];

Actual data from the German stock exchange

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]

{k1, T, r}

Definition of the object OPTION:
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};

optionDB46c[-1, s]

Portfolio of options

Butterfly spread (long CALL k1, long PUT k3, short CALL k2, short PUT k2)

valuepof portfolio[s_] := -c1 - p3 + c2 + p2 + Max[s - k1, 0] + 
        Max[k3 - s, 0] - Max[s - k2, 0] - Max[k2 - s, 0]

There is a space in the definition. Corrected version is

valuepofportfolio[s_] := -c1 - p3 + c2 + p2 + Max[s - k1, 0] + 
        Max[k3 - s, 0] - Max[s - k2, 0] - Max[k2 - s, 0]

Plot[valuepofportfolio[s], {s, 65, 76}];

valueBS[s_, t_] := -c1 - p3 + c2 + p2  + 
        call[s, k1 = 67.5, 0.2, T = 58/360, t, r = 0.05] + 
        put[s, 72.5, 0.2, 58/360, t, 0.05] - 
        call[s, k2 = 70, 0.2, 58/360, t, 0.05] - 
        put[s, 70, 0.2, 58/360, t, 0.05]

valueBS[s, t] // Simplify

Plot[valueBS[s, 30/360], {s, 55, 86}];

Plot[valueBS[s, (58 - 0.01)/360], {s, 55, 86}, 
        PlotRange -> All];

tabBS = 
        Table[Plot[valueBS[s, (t)/360], {s, 55, 86}, 
            PlotRange -> All, DisplayFunction -> Identity], {t, 0, 
            55, 5}];

Show[tabBS, DisplayFunction -> $DisplayFunction]

Illustration

Illustration

(* After marking the cell *)
Illustration

Sensitivities (The greeks)

Clear[s, k, σ, T, t, r];
call[s, k, σ, T, t, r] // 
      Simplify

Delta

(Subscript[Δ, c]) = 
      D[call[s, k, σ, T, t, r], s] // Simplify

Plot[(Subscript[Δ, c]) /. {k -> 50, σ -> 0.2, T -> 20/52, 
          t -> 0, r -> 0.05}, {s, 20, 70}]

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

Yield Curves

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}]

A mistake: even in version 5.1, NonlinearFit is not a built-in function!

Remove[NonlinearFit]
Needs["Statistics`"]

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];

5475/365

2. /3

1.15

2. /3

After copying it becomes

0.6666666666666666`

Remark (Linear Regression)

data = Table[{x, (x^2)}, {x, 6}];
ListPlot[data, PlotJoined -> True];

reg = 
      Regress[data, {1, x}, x, 
        RegressionReport -> {SinglePredictionCITable, 
            ParameterConfidenceRegion}]

Needs["Graphics`"]
Show[Graphics[ParameterConfidenceRegion /. reg], Axes -> True, 
      AxesLabel -> {"Constant", "x squared"}]

Lecture 8

Lecture 8

Lecture 9

FAQ

Remarks (Additional features and applications of Map and Apply and another useful tricks)

Pure function (Ryzi funkce)

Clear[f1, f2];
f1 = 
      Function[x, x^2]

{f1[a], f1[3]}

f2 = #^2 &

{f2[a], f2[3]}

f3 = (#1 + #2^2) &

{f3[a, b], f3[1, 2]}

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

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

Dalsi priklady uziti ryzich funkci
Funkce Select[list, kriterium] vybere ty prvky listu, pro ktere kriterium vraci hodnotu True:

{a, , 23, 12, 0, 3, 5}
Select[{a, , 23, 12, 0, 3, 5}, EvenQ]

Cisla vetsi nez 10:

Select[Range[20], # > 10 &]
Select[Range[20], Function[x, x > 10]]

Cisla  suda a vetsi nez 10:

Select[Range[20], And[# > 10, EvenQ[#]] &]

Select[Range[20], And[# > 10, EvenQ] &]
Select[Range[20], And[# > 10 &, EvenQ]]

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

(mat1 = {Range[5], Range[5]^2} // Transpose) // 
      MatrixForm

Z matice mat1 vybrat radky, ktere maji druhou slozku > 15:
(Pozor, matice mat1 je matice 2x5, takze pozadovane radky jsou {4,5} a {16,25}!)

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

Trideni

SeedRandom[13];
(mat4 = Table[{Random[], Random[]}, {5}]) // MatrixForm

Sort[mat4] // MatrixForm

Sort[{"\:017eu\:017eu", "zuzu"}]
Sort[{"bca", "Acb", "acb", "a2", "a10"}]

Trideni radku matice sestupne podle druheho sloupce:

Sort[
        mat4, #1[[2]] > #2[[2]] &]//MatrixForm

Vzestupne:

Sort[
        mat4, #1[[2]] <= #2[[2]] &]//MatrixForm

Pozor, neco jineho je (nastaveno implicitne):

Sort[mat4]//MatrixForm

{1, 2} ≤ {3, 4}
Thread[{1, 2} ≤ {3, 4}]

(aa = {a1, a2, a3}) == (bb = {b1, b2, 
            b3})
(aa = {a1, a2, a3}) === (bb = {b1, b2, 
            b3})
Thread[aa == bb]

Sort[mat4, 
        OrderedQ[{#1, #2}] &]//MatrixForm

#1 <= #2 &

A mistake

Sort[mat4, #1 <= #2 &]//MatrixForm

Thread[#1 <= #2 &]

Another mistake

Sort[mat4, 
        Thread[#1 <= #2 &]]//MatrixForm

Sort[
        mat4, #1[[1]] <= #2[[1]] &]//MatrixForm

Jiny priklad:

Sort[{-8, -5, 6, 9}, #(1^2) < #(2^2) &]

Dalsi priklady viz Help

Derivace mocniny jako ryzi funkce

(#1^n) &'

Clear[n];
{(#1^n) &'[x], (#1^n) &'''[
        x], ((∂_x) (#1^n) &)[
        x], ((∂_x) (#1^n) &)[x], 
      D[(#1^n &), {#1, 3}] &}

D[(x^n), {x, 4}]

D[(x^n), {x, 3}]

D[((#1^n) &, {#1 &, 4})]

Map and Apply

{Map[f, {a, b, c}], f /@ {a, b, c}}
{Apply[f, {a, b, c}], f @@ {a, b, c}}
Apply[f, g[h[x]]]

{mat4[[0]], g[
          h[x]][[0]]}

{FullForm[{a, b, c}], FullForm[a + b + c]}

Apply[Times, a + b + c]
Map[f, a + b + c]
Map[Times, a + b + c]

list1 = {1. , 2, 3, 3  π};
Sin[list1]
IntegerQ[list1]
Map[IntegerQ, list1]
Map[IntegerQ[#] &, list1]

Duvod spociva v tom, ze funkce IntegerQ neni "listovatelna":

{Attributes[Sin], Attributes[IntegerQ]}

mat1 // Transpose

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

Apply[{#1, 2*#2} &, mat1 // Transpose]

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

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

Reverse[{a, b, c, d, e}]
Reverse[{{a, b}, {c, d}, {e, f}}]
Map[Reverse, {{a, b}, {c, d}, {e, f}}]

(m6 = Map[Reverse, mat4]) // MatrixForm   
Map[Reverse, Sort[m6]] // MatrixForm

Sorting with positions - trideni s puvodnim poradim

sorpos[observ_?VectorQ] := 
      Module[{au, n}, n = Length[observ]; 
        Print[au = Transpose[{observ, Range[n]}], "\nn= ", n]; 
        Sort[au]]

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

y = 3

fun1[x_] := 
      Module[{y}, fun2[y_] := (y^2); fun2[x] + 1]
fun1[5]

sorpos[Table[Random[], {5}]]

Transpose[{1, 2, 3}, {3, 4, 5}]

mat55 = Map[{Random[], #} &, {b, d, a, c}]

Sort[mat55]

If  statement

Clear[a, b];
a == 
      b

If[a == b, c, d]

If[a == b, c, d, e]

If[a === b, c, d]

{a == b, a === b}

Which[7 < 5, a, 
            3 < 2, 
      b, 
            True, aaa = 14; 
      5]
Which[7 < 5, a, 
3 < 2, b, 
True, 
      aaa = 14, 
True, 5]

aaa

Which[7 < 5, a, 3 < 2, b, False, c]

Clear[ff];
ff[x_] := (x^2)
ff[2] = 2
?ff
{ff[5], ff[2]}

Function[x, D[x^5, {x, 3}]]

\(\(Function[x, ∂\_{x, 3}x\^5]\)[x]\)

\(∂\_{#, 3}\ #\^5 &\)

Lecture 10

Patterns, Rules, and Definitions (Vzorky, pravidla a definice)

Patterns 
Blank[]               _    (just one argument) 
double blank       __  (one or more arguments)
triple blank         ___ (arbitrary number of arguments, often used for Options which may be missing)

function[x_,y___]

Clear[hu];
\(hu[x___] := 
      Which[Length[{x}] == 0, 
        Print["No arguments today: ", 
          Date[]], 
			\ \ \ \ \ \ True, 
        If[Length[{x}] == 1, Print[x], 
          Plus[x]]]\)
hu[]
hu[2]
hu[
      r]
hu[r, g, 2005, 12, 
      a]

?MatchQ

MatchQ[f[a, b, h[x]], _]

{MatchQ[h[a, b], _h], MatchQ[h[a, b], _f], 
      MatchQ[h[a, b], _a]}

?Cases

Cases[{1, x, f[a], f[a, b], g[e1, e2]}, _]

Cases[{1, x, f[a], f[a, b], g[e1, e2]}, _f]

Cases[{1, x, f[a], f, f[], f[a, b], g[e1, e2]}, f[_]]

Cases[{1, x, f[a], f, f[], f[a, b], g[e1, e2]}, f[__]]

Cases[{1, x, f[a], f, f[], f[a, b], g[e1, e2]}, f[___]]

Cases[{1, x, f[a], f, f[], f[a, b], g[e1, e2]}, _[___]]

Cases[{f[1], f[0], f[-1], f[a]}, 
      f[x_] /; x ≥ 0]

Cases[{f[1], f[0], f[-1], f[a]}, 
      f[x_ /; x ≥ 0]]

Cases[{f[1], f[0], f[-1], f[a]}, 
      f[x_ /; NonNegative[x]]]

Cases[{f[1], f[0], f[-1], f[a]}, f[x_?NonNegative]]

Cases[{f[1], f[0], f[-1], f[a]}, f[_?NonNegative]]

Substitution

ff[1] /. {1 -> 2}

{f[0], f[1], g[0]} /. {1 -> 2}

{f[0], f[1], g[0]} /. {List -> h}

Srovnej:

Apply[h, {f[0], f[1], g[0]}]

list1 = {f[0], f[1], g[0]};

list1 /. List -> h

list1 /. {1 -> 2}

(list2 = {f[0], f[2], g[0]}) /. {f[x_] -> (x^2)}

Okamzita a odlozena substituce:

{x, x, x} /. x -> Random[]

{x, x, x} /. x :> Random[]

z = 17;

{f[z_] :> (z^3), f[z_] -> (z^3)}

f[5] /. f[z_] :> (z^3)

f[5] /. f[z_] -> (z^3)

Opakovana substituce - Repeated substitution

fac3 = 
      fact[3] /. {fact[1] -> 1, 
          fact[n_] :> n*fact[n - 1]}

fac6 = 
      fac3 /. {fact[1] -> 1, 
          fact[n_] :> n*fact[n - 1]}

fac6 /. {fact[1] -> 1, 
        fact[n_] :> n*fact[n - 1]}

fac3 = 
      fact[3] //. {fact[1] -> 1, 
          fact[n_] :> n*fact[n - 1]}

fact[13] //. {fact[1] -> 1, 
        fact[n_] :> n*fact[n - 1]}

fact[-3] //. {fact[1] -> 1, 
        fact[n_] :> n*fact[n - 1]}

?Thread

Thread[h[a + b + c + d], Plus]

?Select

list22 = Range[6] - 3

Select[list22, Positive]

Cases[list22, n_ /; n > 0]

We can not change "criterion" and "pattern":

{Select[list22, n_ /; n > 0], Cases[list22, Positive]}

Find the number of positive numbers in list22:

Length[Cases[list22, n_ /; n > 0]]

Count[list22, n_ /; n > 0]

Count[list22, _?Positive]

Zahadna funkce - A Mystery Function

Clear[Mystery]
Mystery[list_List] := 
      mystery[list, {}]
mystery[{}, r_] := 
      r
mystery[{e_, list___}, {r___}] := 
      mystery[{list}, {e, r}]

list33 = {a, b, c, d};

Mystery[list33]

Reverse[list33]

Markowitz portfolio

stream = 
      OpenRead["C:\Dokumenty\VPFPM\Returns.dat"]
nn = Read[stream, Number]
tickers = 
      Read[stream, Table[Word, {nn + 1}]]
(datareturns = ReadList[stream, Table[Number, {nn + 1}]]) // 
      TableForm
Close[stream]

returns = 
        Map[Drop[#, 1] &, datareturns];
Needs["Statistics`MultiDescriptiveStatistics`"]
expected = Mean[returns]
cov = CovarianceMatrix[returns]

n = Length[expected]
ones = Table[1, {n}]
xx = Table[(Subscript[x, i]), {i, n}]

lagrange = (1/2) 
            xx . cov . 
              xx + (Subscript[λ, 1])*(1 - ones . xx) + (Subscript[λ, 2])*(μ - 
                expected . xx) // Simplify

Map[D[lagrange, #] &, xx]

?x

?y

xxxx = Table[Unique[], {n}]

x5 = Table[Unique[], {n}]

$1*(xxxx . x5) // Expand

eq1 = 
      Map[D[lagrange, #] &, {xx, (Subscript[λ, 1]), (Subscript[λ, 2])} // 
            Flatten] == 0

Thread[eq1]

Clear[optimalportfolio];
optimalportfolio[expected_, cov_, μ_]:=Module[{n, ones, xx, λ1, λ2, eq1, rp, risk, 
            portf},
n = Length[expected];
ones = Table[1, {n}];
xx = Table[Unique[], {n}];
  
lagrange = (1/2) 
                    xx . cov . 
                      xx + λ1*(1 - 
                        ones . xx) + λ2*(μ - 
                        expected . xx) // Simplify;
  
eq1 = 
              Map[D[lagrange, #] &, {xx, λ1, λ2} // 
                    Flatten] == 0;
eq1 = Thread[eq1];
portf = (xx /. 
                    Solve[eq1, {xx, λ1, λ2} // Flatten]) // 
                Flatten;
Print["Optimal portfolio:  ", 
              portf, "  Expected return = ", 
              portf . expected, "   Risk = ", 
              Sqrt[portf . cov . portf]];
portf]

Priklad: pozadovany vynos 4

μ = 4;
por1 = optimalportfolio[expected, cov, μ];

Apply[Plus, por1]
por1 . expected
por1 . cov . por1 // Sqrt

Markowitz Mean-Variance Portfolio

We assume that N assets have expected returns given by the N×1vector r and the covariance matrix of returns V. We will seek the optimal portfolio in the sense that the Sharpe's ratio (also called Sharpe's measure of portfolio) is maximum:

max(((x^T)r)/(Sqrt[(x^T) V x]))subject to conditionsx≥0and(∑_(n = 1)^N)(Subscript[x, n])=1,

so that short sales are not allowed.

We have 8 assets traded on Prague Stock Exchange with the expected returns

returns = {-9.92488533927622*10^-04, 
          5.72524775783378*10^-03, 
          9.04818959590086*10^-03, -1.44968966764537*10^-03, -1.74045537023702*10^-02, 5.70188394543127*10^-03, 
          3.33476815091662*10^-03, 
          2.15075479228363*10^-03};

and the covariance matrix of returns

covariancematrix = {{1.35640293495777*10^-03, 
            2.56330097904098*10^-04, 3.33037161321956*10^-04, 
            1.80223328810003*10^-04, 2.81325992644196*10^-04, 
            5.93263545509886*10^-04, 7.13610711759175*10^-04, 
            3.20932535820455*10^-04}, {2.56330097904098*10^-04, 
            2.35439634554186*10^-03, 4.42407197078591*10^-04, 
            6.57567906193998*10^-06, 4.01068448708833*10^-04, 
            5.93785356343409*10^-04, 1.30961261626908*10^-04, 
            2.93338399396577*10^-04}, {3.33037161321956*10^-04, 
            4.42407197078591*10^-04, 
            1.57895651834543*10^-03, -1.16804549527263*10^-04, 
            1.85336104268171*10^-04,  .000176059157865, 
            3.41400743910472*10^-04, 
            1.86624552570071*10^-04}, {1.80223328810003*10^-04, 
            6.57567906193998*10^-06, -1.16804549527263*10^-04, 
            6.16786555694974*10^-03, 4.84828341822248*10^-04, 
            8.50333422445824*10^-04, 7.74803938060802*10^-05, 
            4.39412705238888*10^-04}, {2.81325992644196*10^-04, 
            4.01068448708833*10^-04, 1.85336104268171*10^-04, 
            4.84828341822248*10^-04, 9.45147795376653*10^-03, 
            1.78331918975426*10^-03, 5.0987274991634*10^-05, 
            1.49823427374484*10^-03}, {5.93263545509886*10^-04, 
            5.93785356343409*10^-04,  .000176059157865, 
            8.50333422445824*10^-04, 1.78331918975426*10^-03, 
            7.38754498647722*10^-03, 1.03084013803777*10^-03, 
            9.34434830622224*10^-04}, {7.13610711759175*10^-04, 
            1.30961261626908*10^-04, 3.41400743910472*10^-04, 
            7.74803938060802*10^-05, 5.0987274991634*10^-05, 
            1.03084013803777*10^-03, 2.1036055256453*10^-03, 
            2.40818723671803*10^-04}, {3.20932535820455*10^-04, 
            2.93338399396577*10^-04, 1.86624552570071*10^-04, 
            4.39412705238888*10^-04, 1.49823427374484*10^-03, 
            9.34434830622224*10^-04, 2.40818723671803*10^-04, 
            2.77726103166242*10^-03}};

We see that the above matrix is badly conditioned:

Det[covariancematrix]

Next we make use of the procedure NMaximize implemented in the Standard Numerical Package Mathematica, version 4.2. First we create the vector of unknown proportions in the portfolio:

xx = Table[(Subscript[x, i]), {i, 1, Length[returns]}]

The criterion is given by

criterion = (((xx . returns))/Sqrt)[xx . covariancematrix . xx];

This loads the package:

Needs["NumericalMath`NMinimize`"]

The constraints are

constraints1 = 
        Apply[And, 
          Table[(Subscript[x, i]) >= 0, {i, 1, 
              Length[returns]}]];
constraints2 = 
        Sum[(Subscript[x, i]), {i, 1, Length[returns]}] == 
          1;
constraints = constraints1 && (constraints2)

Next we solve the problem using the option for the method:

soldif = 
      NMaximize[{criterion, constraints}, xx, 
        Method -> "DifferentialEvolution"]

We see some problems. Nevertheless, we have got a solution satisfying the constraints:

Apply[Plus, xx /. soldif[[2]]]
Print["Expected return of the optimal portfolio = ", (xx /. 
            soldif[[2]]) . returns]

Looking on the criterion, we see that the solution is independent of a nonzero multiple. So we can suspend the second constraint and make the

constraints3 = 
      Apply[And, Table[(Subscript[x, i]) >= 0, {i, 1, Length[returns]}]]

soldif = NMaximize[{criterion, constraints3}, xx ]

We see that the optimal value is slightly higher than in the above case. The optimal portfolio together with its expected return (substantially higher than in the previous case) is simply obtained:

sol = xx /. soldif[[2]];
sol = sol/Apply[Plus, sol]
sol . returns

0.007859525763768638`

Dalsi poznamky k charakteristikam peneznich toku

npv[cf_, i_] := 
        cf . (1/(1 + i))^(Range[Length[cf]] - 
                1)
dur[cf_, 
        i_] := (cf*(Range[Length[cf]] - 
                  1)) . (1/(1 + i))^(Range[Length[cf]] - 1)/
        npv[cf, i]
con[cf_, i_] := 
      Module[{times = 
              Range[Length[cf]] - 1}, (cf*
                times*(times + 1)) . (1/(1 + i)^times)]/
        npv[cf, i]
altIRR[cf_] := 
      Module[{a}, a = N[Solve[npv[cf, x] == 0, x]]; 
        Cases[x /. a, z_Real /; z ≥ 0]]

Priklad

rancf5 = 
      Table[Random[Real, {-1, 1}], {6}]
Plot[npv[rancf5, i], {i, 0, 0.5}]

Solve[npv[rancf5, i] == 0, i]

altIRR[rancf5]

rancf6 = 
      Table[Random[Real, {-1, 1}], {16}]
Plot[npv[rancf6, i], {i, 0, 0.5}]

Solve[npv[rancf6, i] == 0, i]

altIRR[rancf6]

Priklad

Plot[{npv[rancf5, i], npv[rancf6, i]}, {i, 0, 0.5}];

Plot[
      Map[npv[#, i] &, {rancf5, rancf6}] // Evaluate, {i, 0,  .5}]

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 + r < 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]) = (e^(-r Δt)) 
          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

$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, 100,  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, 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, 50, 0.1, 0, 0.4, 
      0.4167]
put[50, 50, 0.4, 0.4167, 0, 0.1]

Recurssion (Fibonacci)

?Fibonacci

{Fibonacci[1], Fibonacci[2]}

fibA[1] = fibA[2] = 1;
fibA[n_] := fibA[n - 1] + fibA[n - 2]

(tabA = 
          Table[{n, Timing[
                    fibA[n]][[1]]/
                Second}, {n, 20, 28}]) // TableForm

ListPlot[tabA, PlotJoined -> True]

?fibA

fibB[1] = fibB[2] = 1;
fibB[n_] := fibB[n] = fibB[n - 1] + fibB[n - 2]

(tabB = 
          Table[{n, Timing[
                    fibB[n]][[1]]/
                Second}, {n, 100, 105}]) // 
      TableForm
ListPlot[tabB]

?fibB

fibC[n_]:=First[MatrixPower[(({{1, 1}, {1, 0}})),n - 2].{{1}, {1}}]

fibC[100]