VPFPMLS2006

Lecture 1

Motivace

4 pilire

Numericke vypocty

Shift+Enter nebo Enter na numericke klavesnici

1 + 1

(2^1000)

N[π, 1000]

Mozno sdruzit vice vypoctu do jedne bunky:

1 + 1
2*3
2*
3
2
+3

Symbolicke pocitovani

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

e2 = Factor[e1]

TeXForm[e1/(Exp[x] + x)]

Grafika

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

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

Programovani

hu[x_] := x^2

{hu[2], hu[Log[2]], hu[Log[2.0]], hu[2, 3]}

hu[x_, y_] := x^y

{hu[Log[2.0]], hu[2, 3]}

?hu

?Global`*

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

(hil15 = hil1[5]) // MatrixForm

{Det[hil15], N[Det[hil15]]}

N[Det[hil1[45]]]

Inverse[hil1[50]][[30, 30]]

$ImportFormats // Length

$ExportFormats // Length

?$*

$Path

Technicke zalezitosti, vstup a vystup

Front End a Kernel
Menu

Vstup 2D

Vstup 2D pomoci klaves. zkratek nebo z palet (pohyb bud Tabelator nebo ⟶)

(∫_e^∞) 
              (1/((t^(2 - p)) (((Log[t]))^q))) \[DifferentialD]t // 
      Timing

Log[E]

Pozor na syntaxi:

{(Sin^2)[(π/4)], Sin[(π/4)]^2, Sin[(π/4)]}

Clear[cos];
cos[x_] := 
      Sin[x]

Typy cisel

Integer 456, Rational =Integer/Integer (both are precise), Real 3.14, Complex (1+I,1+i)

Zavorky [] (argumenty funkci), () priorita algebraickych operaci, {} listy, seznamy, matice, vektory, iteratory, [[ ]] casti vyrazu, nejcasteji listu

Date[]

list1 = {Sin[Pi], Sin[3.14], Sin@3.14, 3.14 // Sin, sin[3.14], 
        Sin[1001 (π/4)], Sin[1001.  (π/4)]}

list1[[2]]

Priklad s matici

hil13 = hil1[3]

hil13[[3, 3]]

Druhy radek

{hil13[[2, All]], Take[hil13, {2}]}

Treti sloupec

hil13[[All, 3]]

{Take[hil13, 2], Drop[hil13, -1]}

Take[{1, 2, 3, 4, 5}, -2]

Range[5]

Operace +,-,* (ekvivalentni je mezera), /, ^

hil30 = hil1[50];

    tim1 = Timing[Det[hil30] // N]

    Timing[Det[hil30 // N]]

Vypreparovani casu

cas = tim1[[1]]

Sin[cas]

cas2 = tim1[[1]]/Second
Sin[cas2]

$MaxPrecision

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

    ?x
xy == x y
xy === x y
{xy, x*y, x y}

EUR = 28.46 CZK

16 EUR

Lecture 2

Technicke zalezitosti, vstupy - pokracovani

Clear[hil3];
hil3[n_] := 
      Array[(1/(#1 + #2 - 1)) &, {n, n}]

h3 = hil3[3]

Esc[[Esc dava [[

h3[[2]]

h3[[{1, 2}]]

2.`40

2.  == 2

2.  - 2.`40

v1 = Sin[3]
v1num1 = N[Sin[3]]
v1num2 = N[Sin[3], 1000];
N[v1num1 - v1num2, 1000]

Sin[3.`100]

{BaseForm[3. , 2], BaseForm[3.`100, 2]}

{BaseForm[ .2, 2], BaseForm[ .2`100, 2]}

BaseForm[ .2 -  .2`100, 2]

0.1411200080598672`

0.1411200080598672`

Moznosti systemu

nahodnamatice[n_] := Table[Random[], {n}, {n}]

Eigenvalues[hil3[2]]

tab1 = 
      Table[Timing[Eigenvalues[nahodnamatice[j]];], {j, 50, 300, 
          10}]

tab1[[All, 1]]

dava totez co

(tab1 // 
          Transpose)[[1]]

Range[50, 300, 10]

ListPlot[
      Transpose[{Range[50, 300, 10], 
          tab1[[All, 1]]/Second}], 
      PlotJoined -> True]

ListPlot[
      Transpose[{Range[50, 300, 10], 
          tab1[[All, 1]]/Second}], 
      PlotJoined -> True]

Vlastni cislo nahodne matice nejvetsi v vabsolutni hodnote. Vuyzijeme "listovatelnost" funkce Abs:

Attributes[Abs]
Abs[{1 + I, -2, -1, -5, 4 I}]

tab2 = 
      Table[Max[Abs[Eigenvalues[nahodnamatice[j]]]], {j, 50, 300, 
          10}]

plot21 = 
      ListPlot[mat2 = ({Range[50, 300, 10], tab2} // 
              Transpose)]

fun2 = Fit[mat2, {1, x}, x]

plot22 = Plot[fun2, {x, 50, 300}]
Show[plot21, plot22]

tim4 = Timing[mat4 = nahodnamatice[1000];]

Det[mat4] // Timing

inv4 = Inverse[mat4];

Presnost

Abs[inv4 . mat4 - IdentityMatrix[1000]] // Max

Needs["DiscreteMath`Combinatorica`"]

Needs["DiscreteMath`"]

Transpozice matice a vektoru

Transpose[{{a, b}, {c, d}}] // MatrixForm
Transpose[{a, b, c, d}]
{{a, b, c, d}} // Dimensions
Transpose[{{a, b, c, d}}] // MatrixForm

Mathematica vhodne interpretuje vektory jako radkove nebo sloupcove v maticovych operacich podle kontextu:

(({{a, b}, {c, d}})).{e, f}

{e, f}.(({{a, b}, {c, d}}))

{e, f}.(({{a, b}, {c, d}})).{e, f}

{{e}, {f}}.(({{a, b}, {c, d}}))

(r5 = Range[5]) . r5

Transpose[{r5}] . {r5} // MatrixForm

r5*r5

Jiny zpusob maticoveho nasobeni dvou vektoru

Clear[a, b, c, d, e, f];
vec1 = {a, b, c}; 
    vec2 = {d, e, f};
vec1*
      vec2
vec1 . 
      vec2
Outer[Times,  vec1, vec2] // 
      MatrixForm
Outer[Plus,  vec1, 
      vec2]
FullForm[
      a + b + c]
Apply[Times, 
      a + b + c]
Apply[List, 
      a + b + c]
FullForm[(a + b)*
        c]
{Apply[Plus, (a + b)*c], 
      Apply[Times, (a + b)*c]}
Plus[
      Plus[Plus[a, b], c], d]

Upozorneni: MatrixForm nema formu matice:

aa = Table[(Subscript[a, ij]), {i, 2}, {j, 2}]
FullForm[aa]
aa . aa // MatrixForm
aa2 = 
      MatrixForm[Table[(Subscript[a, ij]), {i, 2}, {j, 2}]]
FullForm[aa2]
FullForm[Apply[Identity, aa2]]
aa2 . aa2

(ar3 = Array[a, {3, 3}]) // TableForm
ar3

ar3 . ar3

Z palety (pridavani sloupcu Ctrl ,pridavani radek Ctrl Enter

aa4=(({{a, b}, {c, d}}))
aa4 . aa4 // MatrixForm

Sparse Arrays, PseudoInverse

?Pseudo*

m=(({{1, 2}, {1, 2}}));
MatrixForm[
      p = PseudoInverse[m]]
Simplify[
      Apply[Plus, 
          Apply[Plus, (((m . p - IdentityMatrix[2]))^2)]] ≤ 
        Apply[Plus, 
          Apply[Plus, (((m . {{a, b}, {c, d}} - IdentityMatrix[2]))^2)]], {a, b, c, d} ∈ Reals]

{a ≤ b, a ≤ 1}

{a ≤ b, a ≤ 1, π ≤ 3}

Apply[Plus, 
        Apply[Plus, (((m . p - IdentityMatrix[2]))^2)]] ≤ 
      Apply[Plus, 
        Apply[Plus, (((m . {{a, b}, {c, d}} - IdentityMatrix[2]))^2)]]

Frekvence

tab060308 = 
      RandomArray[PoissonDistribution[1], {100}]
freq = Frequencies[tab060308]

Needs["Graphics`"]

BarChart[freq];

MaxMemoryUsed[]

freq2 = Map[Reverse, freq]

ListPlot[freq2, PlotStyle -> PointSize[0.02]]

Generovani pozitivne defintni matice P = C . (C^T):

posdef[n_] := 
      Module[{c, ctransp}, 
        c = Table[Random[], {n}, {n}]; 
ctransp = 
          Transpose[c]; 
c . 
          ctransp]
posdef[3] // MatrixForm
Eigenvalues[posdef[400]] // Min

Zabudovane algoritmy

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

sol1[[1, 1]]

x /. sol1[[1, 1]]

?x

sol3 = 
        Solve[{a (x^2) + y == x, b x +  c (y^2) == y}, {x, 
            y}];
sol3 // Length

Zkouska

({a (x^2) + y == x, b x +  c (y^2) == y} /. 
          sol3[[2]]) // 
      Simplify

eq4 = 
      Apply[Plus, Table[(j + 1)*(x^j), {j, 0, 6}]] == 0

Apply[Plus, Table[(j + 1)*(x^j), {j, 0, 6}]] // Factor

Solve[eq4, x]

sol4 = NSolve[eq4, x]

Map[(eq4 /. #) &, sol4]

FullForm[eq4]

test4 = 
      Map[(Apply[List, 
                  eq4][[1]] /. #) &, 
        sol4]

Chop[test4]

NSolve[eq5 = Cos[x] == x + Log[x], x]

fr5 = FindRoot[eq5, {x, 1}]

eq5 /. fr5

Plot[Cos[x] - x - Log[x], {x, 0.01, 1000}]

eq5 /. fr5i

Cos[x] - x - Log[x] /. fr5i

frtab = Table[FindRoot[eq5, {x, (10^j)}], {j, 0, 6}]

Map[(Cos[x] - x - Log[x] /. #) &, frtab]

Assuming[x ∈ Reals, FindRoot[eq5, {x, 10000}]]

Solve[
      Floor[(x^2)] == 196 && Mod[x, 10] == 4, {x}]

NMinimize[{1, 
        Floor[(x^2)] == 196 && Mod[x, 10] == 4}, {x}]

Lecture 3

Range[5]
Attributes[Range]
Range[5]^2
Select[Range[10], # > 5 &]
{3 =!= 3, 3 =!= 5}
Select[Range[10], # =!= 5 &]

Vymazani vsech promennych

a = 3; b = 4; c = 5;
{a, b, c}
Clear["`*"]
{a, b, c}

Vymazani vsech promennych z kontextu Global:

a = 3;
f[n_] := (n^2)
?f
a
Clear["Global`*"]
a
?f
Remove["Global`*"]
?a
?f

a = 3

Knihovny programů (Statistika)

Uvod

Needs["Statistics`"]
distr1 = NormalDistribution[μ, σ]

PDF[  distr1, 
      x]
CDF[  distr1, 
      x]
CDF[ExponentialDistribution[λ], x]
Table[
      ExpectedValue[(x^k), distr1, x], {k, 4}]
\(Table[ExpectedValue[#\^k &, distr1], {k, 4}]\)

Plot[CDF[ExponentialDistribution[1], x], {x, -1, 1}]

?Pareto*

PDF[distr2 = ParetoDistribution[k, α], 
      x]
Domain[distr2]

Priklad 1

tab1 = Table[Prime[i], {i, 40}];
pl11 = ListPlot[tab1];

Fit (linearni v parametrech)nam da pribliznou hodnotu

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

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

pl12 = Plot[fit1, {x, 1, 40}];

Show[pl11, pl12]

Velmi dulezita poznamka!

pl11 = 
        ListPlot[tab1, 
          DisplayFunction -> Identity];
pl12 = 
        Plot[fit1, {x, 1, 40}, 
          DisplayFunction -> Identity];

Show[pl11, pl12, DisplayFunction -> $DisplayFunction]

Show[GraphicsArray[{pl11, pl12}]]

Priklad 2 (Nelinearni vyrovnavani vynosove krivky, nelinearni regrese)

Needs["Statistics`"]

Data ve forme {{(Subscript[T, 1]), 
          (Subscript[r, 1])}, .{(Subscript[T, 2]), (Subscript[r, 2])},  ... }, (Subscript[T, i]) doba do splatnosti [dny], (Subscript[r, i]) urokova mira [%]:

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`}};

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

V pripade, ze nebyl natazen balik Statistics:

?NonlinearFit

Naprava a po ni se znovu vratit k vypoctu.

Remove[NonlinearFit]
Needs["Statistics`"]
?NonlinearFit

p2ir3 = ListPlot[ir3];

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

Show[p1ir3, p2ir3];

Lecture 4

FAQ

Velke tabulky

velkatabulka = Table[tabulka[i] = Range[i], {i, 5}]

tabulka[3]

velkatabulka[[3]]

tabulkagrafu = 
      Table[Plot[Sin[k x], {x, 0, π}], {k, 
          3}]
tabulkagrafu[[2]]

Jednotkova matice

IdentityMatrix[3] // MatrixForm

Moznosti systemu II

Numerika

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

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

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

(tabsol = Table[{j, FindRoot[eq1, {x, j}]}, {j, -15, 1}]) // 
      TableForm

FullForm[eq1]
lhsminusrhs = Apply[List, eq1]

Transpose[
        tabsol][[2]]

Map[
        lhsminusrhs[[1]] /. # &, Transpose[tabsol][[2]]] // 
      Chop

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

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

∫Sign[x]/(x^(1/3)) \[DifferentialD]x

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

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

(∫_0^∞) 
          Exp[(-1 + I) 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

Linearni algebra

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

MatrixRank[mat2]

nu2 = NullSpace[mat2]

mat2 . Transpose[nu2] // Simplify

Numericke reseni diferencialnich rovnic II

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

Plot[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

i1 = Integrate[(1/(1 + (x^4))), x]

d1 = D[i1, x]

d1 // Simplify

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

Elementarni symbolicke vypocty

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[y2, {x, 140}];]

{Timing[d1 = D[e13, {x, 140}];], 
Timing[d2 = 
            D[e14, {x, 140}];], 
      Timing[d3 = D[y2, {x, 140}];]}

MaxMemoryUsed[]

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

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

Substituce

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

{x, y}

Konec 15. 3. 2006

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

Reseni soustav nelinearnich rovnic

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

NSolve[eq3, {x, y}]

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

Listy

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]

Sin[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 5

FAQ

Jeste dalsi moznost vymazani obsahu promennych (symbolu)

xxx = 5; x3 = 4; y = 5;
Clear["Global`x*"]
{xxx, x3, y}

xxx = 5; x3 = 4; y = 5;
Clear["Global`xx*"]
{xxx, x3, y}

Odkazy

Obecne:

http://www.wolfram.com

http://www.wolfram.cz

The Mathematica Journal
Vyvolejte 
http://www.wolfram.com/services/premiersupport/

a jako licencni cislo zadejte L4667-9660.  Tam najdete vyznaceny The Mathematica Journal. Je pravdepodobne, ze bude fungovat pouze ze site na fakulte

Klavesove zkratky (vyhodna prevence neustaleho prehmatavani mezi mysi a klavesnici): http://media.wolfram.com/brochures/mathematicashortcuts.pdf

Velmi dukladna uvodni instruktaz, viz take dalsi zajimave soubory v adresari Demos, GetingStarted:

"C:\Program Files\Wolfram \
Research\Mathematica\5.2\Documentation\English\Tour\Tour.nb"

Dalsi zajimave zdroje pouceni

"C:\Program Files\Wolfram \
Research\Mathematica\5.2\AddOns\ExtraPackages"

Help⟶ Master Index ⟶ Built - 
      in Functions ⟶ Advanced Documentation

Konecne je take vyhodne se podivat na nektere zdrojove programy ze Standard Add On Packages. Nektere jsou snadno srozumitelne, u nekterych je tomu naopak.

Lecture 6

FAQ

Matice, aplikaca souctu na radky

mat1 = Partition[Range[9], 3]

    Apply[Plus, mat1]   

    Map[Apply[Plus, #] &, mat1]   

    Map[Total, mat1]

Aplikace uzivatelsky definovane funkce

Clear[fjh];
fjh[x__] := 
      Range[Length[{x}]]*{x}
fjh[1, 2, 
      3]
fjh[a, b, c, 
      d]
Apply[
      fjh, {1, 2, 3}]
Map[Apply[fjh, #] &, 
      mat1]

Robustnejsi definice funkce fjh:

fjh2[x_List] := Range[Length[x]]*x
fjh2[{1, 2, 3}]
fjh2[{a, b, c}]
Map[fjh2, mat1]

Plus[{1, 2, 3}]

Rozdil mezi Map a Apply

Map[f, {a, b, c}]
Apply[f, g[a, b, c]]
Apply[Plus, {Sin[x], Cos[x]}]
Apply[Times, {Sin[x], Cos[x]}]
Head[mat1]

Cyklu je mozne se vyhnou vzdy:
Priklad (Sparse Array)

mat2 = SparseArray[{i_, i_} -> 1, {4, 4}]

mat2 // MatrixForm

(x^4) /. x -> 2
x

Moznosti systemu III

Vicenasobne integraly

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

TeXForm[i1]

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

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

HoldForm[Plus[1, 1]]

Integrace osklive funkce

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

i0 = Integrate[Sin[200  x], {x, 0, 10}]
i01 = 
      Integrate[Sin[200  x], {x, 0, 10}] // N
i1 = NIntegrate[Sin[200  x], {x, 0, 10}]
Options[NIntegrate]
i2 = 
      NIntegrate[Sin[200  x], {x, 0, 10}, 
        MaxRecursion -> 16]
i0 - i2

xs = {}; 
    rr = NIntegrate[Sin[200  x], {x, 0, 10}, 
        EvaluationMonitor :> AppendTo[xs, x]];

rr

xs

Spatne!!!

xs = {}; 
    rr = NIntegrate[Sin[200  x], {x, 0, 10}, 
        StepMonitor :> AppendTo[xs, x]];

$MachinePrecision

Graficke moznosti

ArrayPlot[tab3 = Table[Random[], {10}, {10}]]

ListDensityPlot[tab3]

ListDensityPlot[IdentityMatrix[10] + 100]

ArrayPlot[IdentityMatrix[10] + 1]

ArrayPlot[IdentityMatrix[10]]

Nejvetsi spolecny delitel (GCD)

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

Histogram3D

Needs["Graphics`"]

Histogram3D[array5]

Obrazky v TeXu:

\input epsf

\vglue 4 mm
\epsfysize=4cm
\centerline{\epsfbox{var13.eps}}
\vskip 2 mm
\centerline{Obr.3 (PRIBOR6M)}
\vskip 2 mm

Jeste jedno pouziti Outer

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

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

Funkce definovane uzivatelem I

Priklad 1 (Najdete prvni prvocislo nasledujici po prirozenem cislu n)

Uvodem:

{VectorQ[{a, b}], 
VectorQ[
        IdentityMatrix[2]], 
MatrixQ[{a, 
          b}]}
LowerCaseQ["A"]

Attributes[PrimeQ]

PrimeQ[Range[6]]
Map[PrimeQ, Range[6]]

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

nextpr2[n_] := 
      Module[{k = n}, 
While[! PrimeQ[k], k++]; 
        k]
nextpr[a]

nextpr3[n_] := 
      Module[{k = n}, 
While[! PrimeQ[k], k++]; k] /; 
        NumberQ[n]
SetAttributes[nextpr3, Listable]
nextpr3[Range[1000, 1050]]

Priklad 2 (Prumer a rozptyl)

Clear[mean];
mean[x_?VectorQ] := 
      Total[x]/Length[x] /; 
        Length[x] > 0
mean[{a, b, 
        c}]
mean[{a, {b}, 
        c}]
mean[{{a, b}, {c, 
          d}}]
mean[{}]

variance[x_?VectorQ] := 
      mean[(((x - mean[x]))^2)]
variance[Range[5]]

Lecture 7

FAQ

Fortransky zapis cisel a jeho zpracovani

Vystup bezproblemovy

FortranForm[10. ^70]

Vstup

$Path

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

ReadList["Fortrandata.txt", Number]

Associative Arrays JH

mat=(({{"names", "BC", "EFG", "HIJ", "KLM"}, {"a", 1, 2, 3, 4}, {"b", 5, 6, 7, 8}, {"c", 9, 10, 11, 12}, {"d", 13, 14, 15, 16}}));

Clear[data, firstrow];
firstrow = 
      mat[[1]]
 (* POZOR,  z vystupu neni patrne,  
      ze se jedna o promenne typu String ! Viz:*)

mat[[1]] // FullForm

fields = Map[data[#] &, firstrow]

matdrop = Drop[mat, 1] // Transpose

MapThread[Set, {fields, matdrop}]

fields

data["KLM"]

Select[mat, #[[2]] > 7 &]
Select[
      fields, #[[4]] == 
          16 &]

?? fields

fields // Dimensions

fields

{fields[[1, All]], fields[[All, 1]]}

firstrow

Toto je obtiznejsi:

(symbols = ToExpression[firstrow]) // 
      FullForm
nnn = Length[symbols]
MapThread[Set, {symbols, Range[nnn]}]
?symbols
symbols

BC je nyni promenna:

BC

Import dat z Excelu I

prex = Import["PrikladExcel.xls"]

Import dat z Excelu II (neuspesny)

str = OpenRead["PrikladExcel.xls"]
Read[str, {Word, Word, Word, Word}]

Import["Fortrandata.txt"]

Racionalizace, cela cast atd.

{Rationalize[4. ], Rationalize[4.8563453487, 0]}

{4.8563453487 // Floor, 4.8563453487 // Ceiling, 4.8563453487 // Round, 
      4.2 // Round}

Funkce definovane uzivatelem II

Priklad 3 (Nahodna prochazka 1, Random Walk)

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

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

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

Pro ucely mozne rekonstrukce

SeedRandom[13];
randomwalk[3]

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

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

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

Cviceni 1: Zjistete, kde nahodna prochazka skoncila. Zvyraznete tento bod v grafu.

Cviceni 2: Zobecnete na p-rozmernou nahodnoou prochazku.

Priklad 4 (Nahodna prochazka 2, demonstrace rozdilu mezi okamzitym a odlozenym prirazenim)

addone[x_] := x + Random[] - 1/2

addtwo[x_] = x + Random[] - 1/2

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

NestList[Sin, x, 3]

NestList[addone, 0, 3]

walk1 = NestList[addone, 0, {2}]

walk1 = NestList[addone, 0, 1000];
walk2 = NestList[addtwo, 0, 1000];
ListPlot[walk1, 
        PlotJoined -> True];
ListPlot[walk2, PlotJoined -> True];

Priklad 5 (Je dany list permutace prirozenych cisel? PermutationQ)

Needs["DiscreteMath`Combinatorica`"]
RandomPermutation[4]

permutationQ[list_?VectorQ] := 
      Sort[list] == Range[Length[list]]
permutationQ[{1, 1, 2}]
permutationQ[RandomPermutation[40]]
permutationQ[{a, b}]
permutationQ[{a, {b}}]

Priklad 6 (Definice logaritmu)

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

Ryzi funkce (Pure functions), Map, Apply

Ryzi funkce je funkce beze jmena.

Obycejna definice

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

Function[x, (x^2)][a]

Ekvivalentne

\(\(#\^2 &\)[a]\)

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

Priklad 4 revidovany

\(NestList[#\^2 &, x, 3]\)

pj = PlotJoined -> True

SeedRandom[13];
ListPlot[NestList[# + Random[] - 1/2 &, 0, 1000], pj]

Dalsi priklady

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

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

tabr = Table[Random[Integer], {10}]

Select[tabr, # == 1 &]

Z matice n×2 vyber radky, kde druhy element je vetsi nez 15.

mat1 = {Range[6], Range[6]^2} // Transpose

Select[mat1, # == {4, 16} &]

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

Setridit radky matice n×2 sestupne podle druheho sloupce.

(mat2 = Table[{Random[], Random[]}, {5}]) // 
      MatrixForm
Sort[mat2] // MatrixForm
Sort[
        mat2, #2[[2]] < #1[[2]] &] // 
      MatrixForm
Sort[mat2, 
        OrderedQ[{#2[[2]], #1[[2]]}] &] // 
      MatrixForm

Map a Apply

Rozdil

Clear[f, g];
Map[f, {a, b, c}]
Apply[f, {a, b, c}]
Apply[f, g[x]]

Plus[a, b, c]
Plus[{a, b, c}]
{Apply[Plus, {a, b, c}], Total[{a, b, c}]}
Plus[{a, b, c}, {d, e, f}]
Plus[{a, b, c}, d]
Plus[{a, b, c}, {d, e}]

Pro listovatelne funkce neni Apply zapotrebi:

Sin[{a, b, c}]

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

Ale

IntegerQ[cisla1 = {π, 2, 3. , 5.5}]

Attributes[IntegerQ]

Map[IntegerQ, cisla1]

mattrans = mat1 // Transpose

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

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

Faktorizace (revidovana)

{f1 = FactorInteger[c2 = 2 2 3 5 7 7 7 11 11 11 11], 
      c2}

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

Balicky programu (Packages)

Jsou zabudovane (AddOns) , externi a uzivatelske.
AddOns viz Help.

Needs["Statistics`"]

data3 = Table[Random[], {10000}];
data4 = Table[Random[], {10000}];
lr = LocationReport[data3]
dr = DispersionReport[data3]

{med, mea} = {Median, Mean} /. lr

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

TrimmedMean[data3, {0.1, 0.1}]

Correlation[data3, data4]

{Variance[l3 = {a, b, c}], VarianceMLE[l3]}

Intervaly spolehlivosti (Confidence Intervals)

MeanCI[data3]

Options[MeanCI]

MeanCI[data3, KnownVariance -> 1/12]

ceny = Range[5];
(Drop[ceny, 1] - Drop[ceny, -1])/Drop[ceny, -1]

Lecture 8

Ortonormalni baze

Needs["LinearAlgebra`Orthogonalization`"]

{w1, w2, w3} = 
      GramSchmidt[{{3, 4, 2}, {2, 5, 2}, {1, 2, 6}}]

{w1 . w2, w2 . w3, w1 . w3, w1 . w1, w2 . w2, w3 . w3}

Graf s popisy os

Plot[Sin[x], {x, 0, 2 π}, 
        Ticks -> {{0, (π/2), π, ((3 π)/2), 2 π}, 
            Automatic}];

tab1 = Table[{j*(π/2), ToString[1996 + j]}, {j, 0, 4}]

Plot[Sin[x], {x, 0, 2 π}, 
        Ticks -> {tab1, Automatic}];

Balicky programu (Packages) II

Needs["Statistics`"]

data3 = Table[Random[], {10000}];
data4 = Table[Random[], {10000}];

MeanDifferenceCI[data3, data4]
MeanDifferenceCI[data3, data4, 
      KnownVariance -> {1. /12, 1. /12}]
MeanDifferenceCI[data3, data4, KnownVariance -> {1. /12, 1. /12}, 
      ConfidenceLevel -> 0.99]

distr1 = 
      WeibullDistribution[α, β]
PDF[distr1, x]
CDF[distr1, x]
Quantile[distr1, y]

{Mean[distr2 = NormalDistribution[μ, σ]], Variance[distr2], 
      Skewness[distr2], Kurtosis[distr2], Quantile[distr2, y]}

Φ[x_] := 0.5*(1 + Erf[(x/(Sqrt[2]))])

{Φ[-10], Φ[10]}

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

CharacteristicFunction[distr2, z]

{ExpectedValue[(x^2), distr2, x], 
      ExpectedValue[#^2 &, distr2]}

SeedRandom[13];
Table[Random[distr4], {5}]
SeedRandom[13];
RandomArray[distr4, {5}]
SeedRandom[13];
RandomArray[distr4, {5, 5}] // MatrixForm

distr5 = PoissonDistribution[5];
SeedRandom[13];
t5 = RandomArray[distr5, 100];
f5 = Frequencies[t5]

Needs["Graphics`"]

MaxMemoryUsed[]

?BarChart

?Histogram

BarChart[f5]

Histogram[
      f5tr = Transpose[
            f5][[1]]]

Histogram[f5tr, HistogramCategories -> Max[f5tr], 
      ApproximateIntervals -> False]

Options[Histogram]

Vyhlazovaci techniky (Smoothing Techniques)

Klouzave prumery

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

MovingMedian[data6 = {1, 3, 2, 4, 7, 5, 8, 3, 2, 4, 5, 6, 1000}, 
      3]
MovingAverage[data6, 3]
MovingAverage[data6 // N, 3]

Vlastni (uzivatelsky) package

Viz Skeleton.m

\(\(BeginPackage["ProgrammingInMathematica`Skeleton`", \ "ProgrammingInMathematica`Package1`", \ "ProgrammingInMathematica`Package2`"]\)\(
\)
      \)

\(\(Skeleton::usage\  = \ "Skeleton.m is a package that does nothing."\)\(
\)\)
Function1::usage  =  "Function1[n] does nothing."
\(\(Function2::usage\  = \ "Function2[n, (m:17)] does even more nothing."\)\(
\)\)
Begin["`Private`"]     (*begin the private context (implementation part)*) 
Null

\(\(Protect[\ Evaluate[protected]\ ]\)\(\ \ \ \ \ \)\( (*\ 
        restore\ protection\ of\ system\ symbols\ *) \)\(
\)\)
\(\(End[\ ]\)\(\ \ \ \ \ \ \ \ \ \)\( (*\ 
        end\ the\ private\ context\ *) \)\(
\)\)
\(\(Protect[\ Function1, \ Function2\ ]\)\(\ \ \ \ \)\( (*\ 
        protect\ exported\ symbols\ *) \)\(
\)\)
EndPackage[ ]    
Null

Priklad (prumer a smerodatna odchylka)

BeginPackage["Hurt`april`"]
april::usage = "april.m contains two descriptive statistics."
avg::usage = "avg[data] returns the arithmtic mean"
\(\(\(sd::usage = "sd[data] returns the standard error of the data"\)\(
\)
    \)\)
Begin["Private`"]
Unprotect[avg, sd]
avg[data_] := Total[data]/Length[data]
sd[data_] := avg[(((data - avg[data]))^2)]
Protect[avg, sd]
End[]
EndPackage[]

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

?avg

Needs["Hurt`april`"]

Financni aplikace

Charakteristiky peneznich toku I - Kalendar Actual/Actual.

Pozn.: Specializovane kalendare v baliku Financial Essentials, Derivatives Expert

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

Needs["Miscellaneous`Calendar`"]
DayOfWeek[today]
DaysBetween[today, tom = {2007, 12, 15}]
DaysBetween[tom, today]
DaysPlus[today, 50]

Soucasna hodnota

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

Join[cff, {cf4}]

Clear[pv1];
pv1[cf_, r_] := 
      cf . (1/(((1 + r))^(Range[Length[cff]] - 1)))
pv1[Join[cff, {cf4}] // Flatten, i]

Clear[pv1];
pv1[cf_, r_] := 
      cf . (1/(((1 + r))^(Range[Length[cf]] - 1)))
pv1[Join[cff, {cf4}] // Flatten, i]

Attributes[DaysBetween]

dat7 = {{2006, 5, 12}, {2006, 12, 12}, {2007, 1, 13}, {2007, 2, 
            13}};
diff = Map[DaysBetween[today, #] &, dat7]

Clear[pv2];
pv2[cf_, today_, datums_, r_] := 
      Module[{diff, diskontnifaktory}, 
diff = 
          Map[DaysBetween[today, #] &, 
            datums]; 
diskontnifaktory = 
          (1. /(((1 + r))^(diff/360. ))); 
cf . 
          diskontnifaktory]

pv2[cff, today, dat7, r]

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

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

Profil ciste soucasne hodnoty

Plot[
      pv2[{-100, 10, 10, 110}, today, dat7, r], {r, 0.01, 
        0.3}]

Lecture 9

FAQ

Pouziti jinych baliku ve vlastnim baliku

Mean[Range[3]]

Rozdil mezi Join a AppendTo

Join[list5 = {a, b, c}, d]

{Join[list5, {d}], AppendTo[list5, d], list5}

?{a,b,c}

Insert[{1, 2, 3, 4, 5}, a, 3]

Cviceni I

Needs["DiscreteMath`Combinatorica`"]

rp = RandomPermutation[12]

fl = 
      FoldList[Max, rp[[1]], 
        rp]

Union[{4, 2, 3, 6, 6, 1, 2, 3, 4, 5, 6, 7, 8, 9}]

un = Union[fl]

Rest[un]

keblusek[li_List] := 
      Union[FoldList[Max, li[[1]], 
          li]]
keblusek[rp]

Flatten

Flatten[{x1, {x2, {x3, x4}}, x5}]

Balicky programu (Packages) III

?ExpectedValue

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:

1

BeginPackage["Hurt`fmdm`"]
\(fmdm::usage = "fmdm.m contains two descriptive statistics."
\)
\(avg::usage = "avg[data] returns the arithmetic mean of the data"
\)
\(sd::usage = "sd returns the standard deviation of the data"
\)
exval::usage = "exval calculates exval"

Begin["`Private`"]
Needs["Statistics`"];
Unprotect[avg, sd, exval]
avg[
        data_] := (Apply[Plus, data]/Length[data] // 
            N) /; (VectorQ[data] && 
            Length[data] > 0)
sd[data_] := 
        avg[(data - avg[data])^2] // Sqrt;
exval[f_] := 
      ExpectedValue[f[x], NormalDistribution[0, 1], 
        x]
Protect[avg, sd, exval]
End[]
EndPackage[]

?avg

?exval

avg[Range[3]]

sd[Range[3]]

\({exval[# &], exval[#\^2 &], exval[Sin]}\)

Select the cell with the package 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], l4}

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

plcall = 
      Plot[call[s, 50, 0.2, 20/52, 0, 0.05], {s, 40, 100}, 
        PlotRange -> All]
plput = 
      Plot[put[s, 50, 0.2, 20/52, 0, 0.05], {s, 40, 100}, 
        PlotRange -> All]
plintr = 
      Plot[s - 50, {s, 40, 100}, 
        PlotStyle -> {Hue[0.3]}]
Show[plcall, plput, plintr]

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 
pos=+1 long, pos=-1 short

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)

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

plbut1 = Plot[valuepofportfolio[s], {s, 55, 85}];

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

plbut2 = Plot[valueBS[s, 30/360], {s, 55, 85}];

Show[plbut1, plbut2, PlotRange -> All]

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]

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

Plot3D[
      rho /. {k -> 50, σ -> 0.2, T -> 20/52, t -> 0}, {s, 20, 70}, {r, 
        0.02, 0.07}]

Plot[
        Table[rho /. {k -> 50, σ -> 0.2, T -> 20/52, t -> 0, 
                s -> j}, {j, 20, 70, 10}] // Evaluate, {r, 0.02, 
          0.07}];

Vyzkousejte bez Evaluate a podivejte se do Helpu.

Markowitz portfolio I

Dalsi poznamky k importu a exportu

Import["C:\Dokumenty\VPFPM\Returns.csv"]

{"Month", "AMR", "BS", "UTX", "WX", "Z"}
{1, 11.594`, 1.515`, -1.579`, -9.36`, -0.408`}

Import["C:\Dokumenty\VPFPM\Returns.xls"]

Import["C:\Dokumenty\VPFPM\Returns.dat", "Table"]

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]

Lecture 10

FAQ

Protect a Unprotect

Sin[x_] := x

log[x_*y_] := log[x] + log[y]
log[d 2 e]
Protect[log]
\(log[x_\^y_] := y\ log[x]\)
Unprotect[log]
\(log[x_\^y_] := y\ log[x]\)
log[(x^z)]

Substituce (vybere si pouze ty, ktere potrebuje)

{x -> 2, y -> 3, z -> 4}

x /. {x -> 2, y -> 3, z -> 4}

Uziti Flatten pri reseni soustavy rovnic

sol1 = 
      Solve[{a x +  b y == c, d x  + e y == f}, {x, 
          y}]

{x, y} /. sol1

{x, y} /. sol1 // Flatten

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

Poznamka ke vzorkum (patterns)

hu1[z_List] := (z^z)

hu1[{1, 2, 3, 4}]

hu1[3]

hu2[z_Real] := Sin[z]
{hu2[3. ], hu2[π/4], hu2[π/4. ], hu2[3], 
      hu2[3 + I]}

{Head[π], Head[3], Head[3. ]}

Priklad na reseni rovnic uzitim NMinimize

eq5 = 
        Floor[(x^2)] == 196 && 
          Mod[x, 10] == 4;
FindRoot[eq5, {x, 6}]

NMinimize[{1, eq5}, x]

Markowitzovo portfolio II

Mean[{1, 2}, {3, 4}]
Map[Mean, {{1, 2}, {3, 4}}]

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]

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

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

$1*(xxxx . x5) // Expand

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

ineq1 = {y, a, b, π, z} > 2
ineq2 = {1, 2, 3} > 2
FullForm[ineq1]

Thread[ineq1]
Thread[ineq2]

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

por1

{Apply[Plus, por1], Total[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.

Example 1 (short sales allowed)

Clear[xx];
xx = Table[Unique[], {Length[expected]}];

crit1 = xx . (expected/(Sqrt[xx . cov . xx])) // Simplify

?NMax*

sol11 = NMaximize[crit1, xx]

xa = 
      xx /. sol11[[2]]
xnorm = xa/Total[xa]
Total[xnorm]

val1 = crit1 /. {xx -> xnorm // Thread}

First[val1]
val1[[1]]

Example 2 (short sales not allowed)

ineq = xx ≥ 0 // Thread
ineq = Apply[And, ineq] && Total[xx] == 1

sol2 = NMaximize[{crit1, ineq}, xx]

xx /. Last[sol2]

Example 2 (Revisited)

ineq2 = xx ≥ 0 // Thread
ineq2 = 
      Apply[And, ineq2] && 
        Total[xx] ≤ 1
sol2r = NMaximize[{crit1, ineq2}, xx]

Example 3 (Portfolio with prescribed expected return and minimum risk, short sales not allowed)

crit3 = xx . cov . xx // Simplify
cons3 = 
      Apply[And, Thread[xx ≥ 0]] && Total[xx] == 1 && 
        xx . expected == 4

sol3 = NMinimize[{crit3, cons3}, xx]

{58.14350766103593`, {$16 -> 0.`, $17 -> 0.`, $18 -> 
          0.8412876654310535`, $19 -> -3.604767387366259`*^-10, $20 -> 0.15871233492942294`}}

Total[xx3 = xx /. Last[sol3]]

xx3
expected

xx3 . expected

Unfeasible requirement

crit34 = xx . cov . xx // Simplify;
cons34 = 
        Apply[And, Thread[xx ≥ 0]] && Total[xx] == 1 && 
          xx . expected == 4.3;
sol3 = NMinimize[{crit34, cons34}, xx]

Example 4

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:

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

The criterion is given by

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

This loads the package: (Not necessary in version 5.2)

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

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

Apply[Plus, xxPS /. soldif[[2]]]
Print["Expected return of the optimal portfolio = ", (xxPS /. 
            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}, xxPS ]

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 = xxPS /. soldif[[2]];
sol = sol/Apply[Plus, sol]
sol . returns

0.007859525763768638`

Dalsi poznamky k charakteristikam peneznich toku

Clear[npv, dur, con];
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}, PlotRange -> All]

altIRR[rancf5]

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

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

cf2a = Join[{-1050}, Drop[cf2, 1]]
cf2b = Join[{-950}, Drop[cf2, 1]]

Vynosy do splatnosti (YTM)

{altIRR[cf2], altIRR[cf2a], altIRR[cf2b]}

Priklad

rancf6 = Table[Random[Real, {-1, 1}], {6}];

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

Vnitrni vynosove procento pri neekvidistantnich okamzicich plateb

Needs["Miscellaneous`"]
Clear[pv2];
pv2[cashflow_, startdate_, dates_, r_] := 
      cashflow . (1/(1 + r)^(Map[DaysBetween[startdate, #] &, dates]/
                  365))
today = Date[][[{1, 2, 
          3}]]
dates2 = 
      Map[DaysPlus[today, #] &, 
        Table[182*(j - 1), {j, Length[cf2]}]]

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

pv2[cf2, today, dates2 // N, r]

NMinimize[Abs[pv2[cf2, today, dates2, r]], {r}]

NMinimize[{1, pv2[cf2, today, dates2, r] == 0}, {r}]

FindRoot[pv2[cf2, today, dates2, r], {r, 1 + 0.1}]

Table[
      FindRoot[pv2[cf2, today, dates2, r], {r, j + 0.1}], {j, 0, 
        10}]

Konec 19.4.2006

Lecture 11

FAQ

Packages IV (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:

* Content-type: application/mathematica **************

                     CreatedBy='Mathematica 5.2'



                    Mathematica-Compatible Notebook



This notebook can be used with any Mathematica-compatible

application, such as Mathematica, MathReader or Publicon. The data

for the notebook starts with the line containing stars above.



To get the notebook into a Mathematica-compatible application, do

one of the following:



* Save the data starting with the line of stars above into a file

  with a name ending in .nb, then open the file inside the

  application;



* Copy the data starting with the line of stars above to the

  clipboard, then use the Paste menu command inside the application.



Data for notebooks contains only printable 7-bit ASCII and can be

sent directly in email or through ftp in text mode.  Newlines can be

CR, LF or CRLF (Unix, Macintosh

BeginPackage["Hurt`fmdm`"]
  
\(avg::usage = "avg[data] returns the arithmetic mean of the data"
\)
\(sd::usage = "sd returns the standard deviation of the data"
\)
exval::usage = "exval calculates exval"

Begin["`Private`"]
Needs["Statistics`"];
  
avg[
        data_] := (Apply[Plus, data]/Length[data] // 
            N) /; (VectorQ[data] && 
            Length[data] > 0)
sd[data_] := 
        avg[(data - avg[data])^2] // Sqrt;
exval[f_] := 
      ExpectedValue[f[x], NormalDistribution[0, 1], x]
  
End[]
EndPackage[]

?avg

?exval

avg[Range[3]]

sd[Range[3]]

\({exval[# &], exval[#\^2 &], exval[Sin]}\)

Select the cell with the package and save it by Menu -> Edit -> Save Selection As -> Package Format ... -> c:\Dokumenty\VPFPM\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], l4}

$Path

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

AppendTo[$Path, "C:\Dokumenty\VPFPM\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["Statistics`DataSmoothing`"]

Needs["Hurt`fmdm`"]

Needs["Statistics`fmdm`"]

?avg

avg[{a, b, c}]

Clear a Remove

x = 5

?x

Clear[x]
?x

Remove[x]
?x

Jine pouziti pri omylem volane funkci z balicku, ktery nebyl natazen.

tab11 = Table[Random[NormalDistribution[0, 1]], {100}]

?NormalDistribution

Needs["Statistics`ContinuousDistributions`"]

tab11 = Table[Random[NormalDistribution[0, 1]], {10}]

?NormalDistribution

Remove[NormalDistribution]

?NormalDistribution

tab11 = Table[Random[NormalDistribution[0, 1]], {10}]

Associative Arrays JH

mat=(({{"names", "BC", "EFG", "HIJ", "KLM"}, {"a", 1, 2, 3, 4}, {"b", 5, 6, 7, 8}, {"c", 9, 10, 11, 12}, {"d", 13, 14, 15, 16}}));

Clear[data, firstrow];
firstrow = 
      mat[[1]]
 (* POZOR,  z vystupu neni patrne,  
      ze se jedna o promenne typu String ! Viz:*)

mat[[1]] // FullForm

fields = Map[data[#] &, firstrow]

matdrop = Drop[mat, 1] // Transpose

MapThread[Set, {fields, matdrop}]

fields

data["KLM"]

Select[mat, #[[2]] > 7 &]
Select[
      fields, #[[4]] == 
          16 &]

?? fields

fields // Dimensions

fields

{fields[[1, All]], fields[[All, 1]]}

firstrow

Toto je obtiznejsi:

(symbols = ToExpression[firstrow]) // 
      FullForm
nnn = Length[symbols]
MapThread[Set, {symbols, Range[nnn]}]
?symbols
symbols

BC je nyni promenna:

BC

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]

Spatna definice:

\(hu2[x___] := 
      Which[Length[{x}] == 0, 
        Print["No arguments today: ", 
          Date[]], 
			\ \ \ \ Length[{x}] == 1, 
        Print[x], 
			\ \ \ \ True, Plus[x]]\)

{hu2[], hu2[2], hu2[r], hu2[r, g, 2005, 12, a]}

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

?MatchQ

_ zastupuje vsechno:

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

_h zastupuje vsechno s hlavickou h:

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

?Cases

Cases[{2, x, 4}, _Integer]

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}

ff[1]/ .1 -> 
      2

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

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

f[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)}

list2

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

list3

Okamzita a odlozena substituce:

w = 2
{w, w} /. w -> Random[]

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

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

Zaludnejsi priklad

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, demonstrace na faktorialu

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

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

MaxMemoryUsed[]

?Thread

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

?Select

list22 = Range[6] - 3

Select[list22, Positive]

Cases[list22, n_ /; n > 0]

Cases[list22, _?Positive]

Select[list22, _?Positive]

Cases[list22, Positive]

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

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]

Trace[Mystery[list33]]

Reverse[list33]

Priklad 1 (Pocitani minci)

coins = {p, p, q, n, d, d, p, q, q, p};

values = {p -> 1, n -> 5, d -> 10, 
          q -> 25};

v1 = coins /. values

Total[v1]

countchange[
        coins_List] := (coins /. {p -> 1, n -> 5, d -> 10, 
              q -> 25}) // Total

countchange[coins]

Priklad 2 (Opet rekordy)

rek1[x_List] := 
      Union[Rest[
          FoldList[Max, x[[1]], 
            x]]]
list44 = 
      Table[Random[Integer, {1, 20}], {20}]
rek1[list44]

rek2[x_List] := 
      x //. {a___, b_, c___, d_, 
            e___} /; ((d ≤ b): ->  {a, b, c, 
                e})

rek2[list44]

Priklad 3 (1. ukazka kodovani dle Caesara)

alphabet = CharacterRange["a", "z"]

caesarcoderules = 
      Thread[alphabet -> RotateLeft[alphabet]]

caesardecoderules = 
      Thread[alphabet -> RotateRight[alphabet]]

pozdrav1 = "ahoj"
pozdrav2 = Characters[pozdrav1]

"23ergvjgkj785^&%\n wryu"

c1 = Characters["23ergvjgkj785^&%\n wryu"]

c1 // FullForm

StringJoin["\\", "n"]

c1 // StringJoin

Table[FromCharacterCode[j], {j, 125, 255}]

\(\\[Hacek]\)

pozdrav3code = pozdrav2 /. caesarcoderules

pozdrav4code = pozdrav3code // StringJoin

dekod4 = pozdrav4code // Characters
dekod3 = dekod4 /. caesardecoderules
dekod2 = dekod3 // StringJoin

Lecture 12

FAQ

Zobrazeni vice casovych rad soucasne

ser1 = {2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 
          2};
ser2 = {1, 3, 5, 6, 7, 8, 9, 10, 4, 4, 
          4};
ser3 = {2, 3, 6, 7, 7, 9, 7, 4, 3, 2, 1};

joi1 = {ser1, ser2, ser3};
joi1tr = Transpose[joi1];
Dimensions[joi1]
uni1 = Union[joi1 // Flatten]

Map[Count[ser3, #] &, uni1]

tab1 = 
      1.  - Table[
            Map[Count[joi1tr[[i]], #] &, 
              uni1], {i, n = Length[ser1]}]/3

list12 = 
      ListPlot[{{1, 2}}, 
        PlotStyle -> {GrayLevel[
              tab1[[1, 2]]], 
            PointSize[0.02]}]

tab2 = 
        Table[ListPlot[{{i, uni1[[j]]}}, 
            PlotStyle -> {GrayLevel[
                  tab1[[i, j]]], 
                PointSize[0.02]}, DisplayFunction -> Identity], {i, 
            n}, {j, Length[uni1]}];

Show[tab2 // Flatten, 
      DisplayFunction -> $DisplayFunction]

list1 = Map[ListPlot[{{Range[n], ser1}} // Transpose, 
          PlotStyle -> {GrayLevel[
                tab1[[1, 2]]], 
              PointSize[0.02]}]

Plot[Sin[x], {x, 0, π}, 
        PlotStyle -> {GrayLevel[1.0]}];

li1 = 
      ListPlot[{Range[6], Range[6]^2} // Transpose, 
        PlotStyle -> {GrayLevel[ .1], PointSize[0.02]}]

li2 = 
      ListPlot[{Range[6], Range[6]^2 + 1} // Transpose, 
        PlotStyle -> {GrayLevel[ .5], PointSize[0.02]}]

Show[li1, li2]

Empiricka distribucni funkce

1. moznost

Attributes[UnitStep]
?UnitStep

data1 = 
        Table[Random[Integer, {1, 10}], {50}];
data2 = Table[Random[], {100}];

Clear[f1emp];
f1emp[x_, data_] := 
      Apply[Plus, UnitStep[x - data]]/Length[data]

Plot[f1emp[x, data1], {x, 0, 15}];

Plot[f1emp[x, data2], {x, -0.1, 1.1}];

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

Priklad 4 (2. ukazka kodovani dle Caesara)

Needs["DiscreteMath`Combinatorica`"]

?RandomPermutation

alphabet = CharacterRange["a", "z"]

SeedRandom[13];
ranper = RandomPermutation[na = Length[alphabet]]

ranperalpha = 
      Map[alphabet[[#]] &, 
        ranper]

percoderules = Thread[alphabet -> ranperalpha]

perdecoderules = Thread[ranperalpha -> alphabet]

pozdrav1 = "ahoj"
pozdrav2 = Characters[pozdrav1]

pozdrav3 = pozdrav2 /. percoderules

pozdrav4 = pozdrav3 /. perdecoderules

Rekurentni (rekursivni) vypocty

Fibonacciova posloupnost

Nejjednodussi je uzit zabudovanou funkci

Table[Fibonacci[j], {j, 10}]

fibA[n_] := 
      fibA[n - 1] + fibA[n - 2] /; 
        n > 0 && n ∈ Integers
fibA[1] = 1;
fibA[2] = 1;
?fibA

Map[fibA, Range[10]]

fibA[25] // Timing

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

fibB[25] // Timing

?fibB

fibB[250] // Timing

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

Map[fibC, Range[10]]

Diferencni rovnice

?RSolve

solr = 
      RSolve[{x[n] == x[n - 1] + x[n - 2], 
          x[1] == x[2] == 1}, x[n], n]

solr // FullSimplify // TraditionalForm

Table[x[n] /. (solr // Flatten) // Simplify, {n, 1, 10}]

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:

Recall the functions call and put for the European options in a previous lecture:

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

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]

Lecture 10

Lecture 10

Lecture 10

Lecture 1 - 12 (Notes)

Packages

Priklad (prumer a smerodatna odchylka)

Unprotect[avg, sd]
Remove[avg, sd]

BeginPackage["April1`"]
April::usage = "April.m contains two descriptive statistics."
avg::usage = "avg[data] returns the arithmtic mean"
sd::usage = "sd[data] returns the standard error of the data"
Begin["`Private`"]
  
avg[data_] := Total[data]/Length[data]
sd[data_] := avg[(((data - avg[data]))^2)]
  
End[]
EndPackage[]

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

Needs["April1`"]

avg[{a, b, c}]

Lecture 13

FAQ

Dimense matic

mat1[n_] := Table[Random[], {n}, {n}]
mat200 = mat1[200];
Timing[Inverse[mat200];]

mat1000 = mat1[1000];
Timing[Inverse[mat1000];]

mat5000 = mat1[5000];
Timing[Inverse[mat5000];]

Timing[
          Inverse[mat1[
                1000];]][[1]]/
      Second

ListPlot[
      Table[(Timing[
                Inverse[mat1[
                      n];]][[1]]/
            Second), {n, 100, 2000, 50}], 
      PlotJoined -> True]

MaxMemoryUsed[]

EXCEL

Datum

date1=5/10/2006

147

ghjfghrrghphpoghgdfgh121235/10/20064565/10/
                  20067895/10/2006855/10/20065  January855/10/
          200610  February855/10/200615  March855/10/
          200620  April855/10/200625  May855/10/200630  June85/10/
          200635  July85/10/200640  August85/10/
          200645  September850October855November85/10/
          200660  December85/11/2006  January85/12/
          2006  February85/15/2006  MondayMarch85/16/
          2006  TuesdayApril8885/17/2006  WednesdayMay85/18/
          2006  ThursdayJuneAlena5/19/2006  FridayJulyJarmila5/22/
          2006  SaturdayaZuzana5/23/2006  SundaybAlena5/24/
          2006  MondaycJarmila5/25/2006  
        caZuzanacbAlenaaJarmilabZuzanaccacb45

Range[5]

Lecture nest jonm@wolfram.co.uk

{Sin[4.5 + I], Sin[4.5`100 + I]}

fjo[x_] := 11*x - 2

Reduce[fjo[x] == x]

NestList[fjo, 0.2, 30]

{BaseForm[0.2, 2], BaseForm[0.2`30, 2]}

NestList[fjo, 0.2`30, 30]

FixedPoint[(11*# - 2) &, 0.2, 10]

Odkazy

Obecne:

http://www.wolfram.com

http://www.wolfram.cz

The Mathematica Journal
Vyvolejte 
http://www.wolfram.com/services/premiersupport/

a jako licencni cislo zadejte L4667-9660.  Tam najdete vyznaceny The Mathematica Journal. Je pravdepodobne, ze bude fungovat pouze ze site na fakulte

Klavesove zkratky (vyhodna prevence neustaleho prehmatavani mezi mysi a klavesnici): http://media.wolfram.com/brochures/mathematicashortcuts.pdf

Velmi dukladna uvodni instruktaz, viz take dalsi zajimave soubory v adresari Demos, GetingStarted:

"C:\Program Files\Wolfram \
Research\Mathematica\5.2\Documentation\English\Tour\Tour.nb"

Dalsi zajimave zdroje pouceni

"C:\Program Files\Wolfram \
Research\Mathematica\5.2\AddOns\ExtraPackages"

Help⟶ Master Index ⟶ Built - 
      in Functions ⟶ Advanced Documentation

Konecne je take vyhodne se podivat na nektere zdrojove programy ze Standard Add On Packages. Nektere jsou snadno srozumitelne, u nekterych je tomu naopak.