VPFPM LS 2010/2011

Úvod

Termíny zkoušek

Zkousky: 30. 5., 3. 6. (prednost maji uchazeci prihlaseni na letni termin SZZ, 29. 6., vzdy v 11:00 K11
SZZ Bc FM: 27. 6., 28. 6., vzdy v 8:00 KPMS
SZZ Bc FM: 12. 9., 13. 9., vzdy v 8:00 KPMS

Doporučený metodický postup pro studium systému Mathematica  (součást přednášky NFAP008)

Úplný začátek

Podívejte se na Help ⟶ Welcome Screen ...  ⟶ Find Your Learning Path, Get Demonstrations with Source Code, Search Complete Documentatation. Z přednášky si čiňte poznámky. Záznamy z přednášek na mé domovské stránce.

Důkladný průvodce

Help ⟶  Documentation Center. Projít letmo všechny kapitoly. V odstavci Find Your Learning Path prostudovat vsechny nabídky. Specielně:
Printed Tutorial Collection  
http://www.wolfram.com/learningcenter/tutorialcollection/

Demonstrations
http://demonstrations.wolfram.com/

Pro opravdové zájemce “The Mathematica Journal”
http://www.mathematica-journal.com/

Přednášející radí následující postup pro porozumění systému (historické i soudobé poznámky)

(Nečiní si nicméně nárok na to, že jeho doporučení jsou suboptimální, natož optimální. Veškeré konstruktivní připomínky, podněty a náměty k přednášce VPFPM sdělené přednášejícímu jsou vítány. Je možné je podávat ústně, písemně i anonymně a e-mailem. Na připomínky a dotazy podané e-mailem bude přednášející reagovat.)

(i)

Projít celou The Mathematica Book (záložka Help, Virtual Book). Tím se získá přehled, co všechno je možné pomocí systému realizovat. V této fázi není rozumné se věnovat detailům, navíc jest v prvním čtení nemožné porozumět všemu. Cílem je učinit si mlhavou představu o více než 3000 zabudovaných (build-in) funkcí, případně funkcí z Add-On Packages. Pak pomocí "Help" je v případě potřeby syntakticky správně zabudovat do uživatelského projektu.

(ii)

Jako jeden z efektivních postupů se jeví tento: Help ⟶ přibližný název funkce ⟶ přesný název funkce ⟶ (vybrat příbuznou funkci) ⟶ Further Examples ⟶ příslušnou část označit a  Copy+Paste do vlastního notebooku. U některých funkcí je k dispozici "Advanced Documentation", někdy též "Tutorials". Doporučuji její přečtení. Čtenář se zde srozumitelným způsobem seznámi se současnými metodami, které obvykle nejsou v syllabech přednášek konaných na MFF.

(iii)

Není možné se naučit, natož pochopit všechno ihned. K nejmocnějším funkcím systému Mathematica (Apply, Map, Thread, případně Fold, Nest a funkcím podobným) je nezbytné se neustále vracet a zkoušet je v konkrétních situacích.

(iv)

Jádro systému Mathematica bylo naprogramováno v jazyku C. Půjdeme-li zpátky, pak C byl naprogramován v Assembleru, a ten přímo ve strojovém kódu. Pro efektivní užití systému Mathematica však není vhodná zakořeněná znalost ideologie procedurálních jazyků jako jsou Pascal, FORTRAN, C++ a další. Znalosti objektově orientovaného programování jsou naopak vítány. V systému Mathematica je sice možné prakticky doslova přepsat program napsaný v některém z těchto procedurálních jazyků, ale tím by byla zcela negována hlavní zbraň systému Mathematica, tj. funkcionalita, přehlednost, rychlost, kontrolovatelná přesnost a verifikovatelnost při modelování. Zejména z hlediska přehlednosti konstrukce modelů (a ve vyšší fázi použití nadstavby WorkBench pro vytváření projektů) jsou funkce typu "cyklus" naprosto nevhodné. Výhodou pro porozumnění systému Mathematica  je znalost logiky jazyků Prolog, LISP (byť velmi starých, které v době svého vzniku neměly podporu v potřebné počítačové technice) .

(v)

Pro úspěšné zvládnutí přednášené látky stačí zhruba 100 až 200 hodin samostudia. Aktivní zvládnutí celého systému je prakticky nemožné. Předpoklájí se samozřejmě znalosti pravděpodobnosti, statististiky, úvodu do financí a finančního mangementu.

Instalace a licence (studentská L4697-4335)

MFF má do 27. 11. 2011 neomezenou licenci (Unlimited licence), v současnosti verzi 8. Instalační DVD k zapůjčení na 2 dny v knihovně. Po vytvoření účtu na “user.wolfram.com” je možné si také stáhnout instalační soubor přímo (Windows, Linux, Mac)

Pokyny pro registraci

Name: Jméno Příjmení (žádné pseudonymy se nepřipouštějí)
Organization: Charles University
Department: Statistics (Mathematics)

Pokyny pro vlastní instalaci (jsou v obalu instalačního DVD, zde zkopírovány v jiném českém kódování)

Instalace/Registrace MATHEMATICA 8.0.0 (WINDOWS)
1.	Po vložení DVD do mechaniky se automaticky spustí instalaèní program; zde odsouhlaste pøednastavené volby
2.	Na stránce https://user.wolfram.com si vytvoøte „Wolfram Account“ (odkaz Create Account) 
3.	Pøi pøechodu na stránku vytvoøení úètu se vás prohlížeè zeptá, zda má stahovat zabezpeèené i nezabezpeèené položky (zde POZOR!! otázka je nìkdy formulována „Chcete stahovat jen zabezpeèené položky“ pak je tøeba odpovìdìt NE; a v pøípadì, že je otázka formulována „chcete stahovat zabezpeèené i nezabezpeèené položky, pak je tøeba odpovìdìt ANO)
4.	Vyplníte registraèní informace a formuláø odešlete
5.	Obratem vám dojde ovìøovací e-mail s dalším odkazem, kterým potvrdíte správnost e-mailu
6.	Na stránce http://www.karlin.mff.cuni.cz/network/mathematica.html je uveden odkaz pro získání Aktivaèního klíèe
anebo opište následující odkaz
https://user.wolfram.com/portal/requestAK/19bce9af2a5e13230bdc385a0e41a7f8066ca8fb

7.	Obratem vám pøijde e-mail  aktivaèním klíèem s dalšími istrukcemi

Lekce 1 (21. 2. 2011)

Motivace, odkazy

Zakladni stranky:

http://www.wolfram.com

http://www.wolfram.cz

Vyhradni distributor v CR a SR: http://www.elkan.cz/

Ukazky:

http://demonstrations.wolfram.com/

Podivejte se na Latest, pak Search: Finance, Insurance, Economics, Probability, Statistics. Demonstrace je mozne stahnout a studovat. Zacnete studovat ty jednodussi. Cim hezci grafika, tim komplikovanejsi kod a pro zacatecnika hure srozumitelny.

Klavesove zkratky

http://reference.wolfram.com/mathematica/guide/SelectingAndTypingInNotebooks.html

Velmi uzitecne: F1, Ctrl+K, Ctrl+Shift+K, Ctrl+Shift+B

Klavesove zkratky (vyhodna prevence neustaleho prehmatavani mezi mysi a klavesnici, da se vytisknout na dvou strankach podelne): http://media.wolfram.com/brochures/mathematicashortcuts.pdf

Mathematica book v papirove podobe v knihovne (starsi verze)

Novinky verze 8

http://www.wolfram.com/mathematica/new-in-8/

Mathematica pro zacatecniky (Antonin Slavik, KDM)

http://www.karlin.mff.cuni.cz/~slavik/info.html

Mathematica pro pokrocile (Antonin Slavik, KDM)

tamtez

Mathematica Journal Specializovana cisla, statistika viz  Vol. 9(4)

Tutorials (Wolfram Tutorial Collection)

http://www.wolfram.com/learningcenter/tutorialcollection/

Integrated Data Sources

http://reference.wolfram.com/mathematica/guide/ComputableDataOverview.html

FinancialData["EUR/CZK"]

Documentation Center ⟶ guide/ComputableDataOverview

http://reference.wolfram.com/mathematica/guide/NewIn70ComputableData.html

http://reference.wolfram.com/mathematica/guide/FinancialAndEconomicData.html

Další odkazy

http://www.wolfram.com/products/applications/unrisk/

http://www.wolfram.com/products/applications/timeseries/

http://www.wolfram.com/products/applications/finance/features.html

http://www.karlin.mff.cuni.cz/~prokesov/Sbirka%20 uloh %20 z %20 financni %20 matematiky %20-%20 Klara %20 Jelenova.pdf

Zkousky: orientacni terminy 25. 5., 1. 6., 30. 6.

Odehrávají se v laboratoři K11, problémově orientované otázky, řešení je v elektronické podobě odesíláno examinátorovi.

4 pilire

Numericke vypocty

Vypocet pomoci Shift + Enter nebo Enter na numericke klavesnici

1+1
2+3

(2^1000)
N[(2^1000)]

N[(2^100000)]

π

npi=N[π,15]

N[π,{15,10}]

Sdruzovani vice vypoctu do jedne bunky

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

{π,N[π],"log 2 numericky = ",Log[2.]," log 2 presne = ",Log[2],N[Log[2]]}

3.141592653589793`

Symbolicke vypocty

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

e2=Factor[e1]

TeXForm[e2/e1]

(((x-1))^5)//Expand
Expand@(((x-1))^5)

TraditionalForm@Expand[(((x-1))^5)]
TraditionalForm[Expand[(((x-1))^5)]]
Expand[(((x-1))^5)]//TraditionalForm

Pozor na zapis funci

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

Clear[a,b,c,x]
reseni=Solve[a (x^2)+b x +c==0,x]

?x

x=4
Solve[a (x^2)+b x +c== 0,x]
?x

Clear[x]
WolframAlpha[solve (a x^2+b x+c==0),TimeConstraint->60]

a (x^2)+b x +c=0
a (x^2)+b x +c//FullForm

?dbknsdfashgjfsfge

?reseni

Vypreparovani korenu

{x,x^2,y,z}/.{x->a,y->(y^2),t->0}

koreny=x/.reseni

Clear[a,b,c,x]

WolframAlpha["solve a x^2+b x+c=0",IncludePods->"Result",AppearanceElements->{"Pods"}]

i1=(∫_0^∞)(E^(a (x^2)/2))\[DifferentialD]x

Assuming[a>0,(∫_0^∞)(E^(-a (x^2)/2))\[DifferentialD]x]

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

Refine[((Sqrt[(π/2)])/(Sqrt[a])),a>0]

i1/.If[a_,b_,c___]->b

Grafika

Plot[Sin[x+x^2],{x,0,3π},PlotRange->{-5,5},Ticks->{{{0.01,"0"},Pi,2Pi,3Pi},{-5,1,4,5}},FormatType->TraditionalForm,AxesLabel->{x,"Sin[x+(x^2)]\n"}]
Plot[Sin[x+x^2],{x,0,π},PlotRange->{-5,5},Ticks->{{0,(Pi/2),Pi},{-Range[5],0,Range[5]}//Flatten}]
Plot3D[Sin[x+y^2],{x,0,π},{y,0,π}]

Options[Plot3D]//Length

Programovani

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

hu[x_]:=x^2

"="je okamzite prirazeni, ":=" je odlozene prirazeni (globalni pravidlo)

{hu[2],hu[Log[2]],hu[Log[2.]],hu[e1]}

Expand[hu[e1]]

{(Sin^2)[1.],((Sin[1.])^2)}

Poznamka k sazbe: (x^2) se vysazi jako x Ctrl ^ 2, (Subscript[x, 3]) jako xCtrl _ 3, (Sqrt[x]) Ctrl2 x nebo z palety (Menu File ⟶ Palettes nebo pomoci Esc sekvenci

Π ⟶⟶⟵↔↘↗⥮\[ControlKey]\[EscapeKey]⌘

hu[x_,y_]:=x^y

{hu[2.],hu[2,3],hu[2,3,4]}

?hu

hu[2]=0
{hu[3],hu[2],hu[2.]}

?hu

{Unset[hu[2]],hu[2]}

?hu

hu[x_]:=(x^2)

hu[2]

hu[x_]:=(x/2)

hu[y]

?hu

Clear[hu] (* Totez je hu=. *)

?hu

Remove[hu]
?hu

?Global`*

Odstraneni casti definice (viz Unset)

hu[2]=0;
hu[x_]:=(x^2)
hu[x_,y_]:=(x^y)

{hu[2],hu[3]}
hu[2]=.
{hu[2],hu[3]}

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

hil1[3]

h3m=hil1[3]//MatrixForm

Poznamka: f[x] je totez, co x//f resp. f@x

Poznamka: h3m neni totez co hil1[3]!!!

FullForm[h3m]

FullForm[hil1[3]]

S h3m se neda pocitat:

{Det[h3m],Det[hil1[3]]}

Jak odstranit MatrixForm?
1. Nedefinovat symbol jako MatrixForm, ale MatrixForm pouzit pouze pro znazorneni:

(h4m=hil1[3])//MatrixForm

S h4m vsak mohu dale pocitat:

Det[h4m]

2. Aplikovat funkci Identity. Ta vymeni hlavicku vyrazu, tj. MatrixForm, za Identity.

Apply[Identity,h3m]//FullForm

Inverse[hil1[3]]//MatrixForm
Inverse[h3m]

h20=hil1[20];
First[h20inv=Inverse[h20]]
dh20=Det[h20]
N[dh20]
Det[h20inv]
Det[h20]*Det[h20inv]
First[h20.h20inv]
(* 2. radek matice*)
h20.h20inv[[2]] 
First[Rest[h20.h20inv]]

?Global`*

Priklad (nasobeni vektoru a matic)

xx={x1,x2,x3};yy={y1,y2,y3};
xx.yy
xx*yy
Transpose[xx]
{xx}
Transpose[{xx}]
{xx}.{yy}
Transpose[{xx}].{yy}//MatrixForm
{xx}.Transpose[{yy}]
({xx}.Transpose[{yy}])[[1,1]]

Outer[Times,xx,yy]//MatrixForm

Detaily

Integrated Data Sources

FinancialData["EUR/CZK"]

(*Dynamic[Refresh[FinancialData["USD/CZK","LatestTrade"],UpdateInterval->1]]
Dynamic[Refresh[FinancialData["EUR/CZK","LatestTrade"],UpdateInterval->1]] *)

Definice funkce s jednim a vice argumenty

Pri vyvolani procedury se nejdrive vyhodnoti argument. Nevyhovuje-li predepsanemu vzoru (v nasem pripade dvojite __ znamena jeden a vice aergumentu), procedura se vubec nespusti. Viz priklad nize hj2[].

Clear[hj2];
hj2[x__]:=
Which[Length[{x}]==0,Print[x,"     No arguments"],
	  Length[{x}]==1,Print[x,"    Just one argument"],
          Length[{x}]>1,Print[{x},"   More than one argument so that ","\nThe sum of the elements = ",Total[{x}]]		]
{hj2[],hj2[3],hj2[{1,2,3}],hj2[1,2,34,5,6]}

Konec 21. 2. 2011

Funkce hj (libovolny pocet argumentu)

Clear[hj];
hj[x___]:=
Which[Length[{x}]==0,Print[x,"     No arguments"],
	  Length[{x}]==1,Print[x,"    Just one argument"],
          Length[{x}]>1,Print[{x},"   More than one argument so that ","\nThe sum of the elements = ",Total[{x}]]		]
{hj[],hj[3],hj[{1,2,3}],hj[1,2,34,5,6]}

Clear[hj1];
hj1[x___]:=
Which[Length[{x}]==0,{x,"     No arguments"},
	  Length[{x}]==1,{x,"    Just one argument"},
          Length[{x}]>1,{{x},"   More than one argument so that ","\nThe sum of the elements = ",Total[{x}]}		]
{hj1[],hj1[3],hj1[{1,2,3}],hj1[1,2,34,5,6]}

hj[{1,2}⋂{3,4}]

Lekce 2 (28. 2. 2011)

FAQ

Poznamky k maticim

Remove [aa];Clear[bb,cc,dd];
aa=Table[a[i,j],{i,2},{j,3}]
(bb=Table[b[i,j],{i,3},{j,2}])//MatrixForm
cc={c[1],c[2],c[3]};
dd={d[1],d[2]};

Za cviceni vysvetlete vystupy nasledujicich prikazu:

a[1,1]=4444;
aa
aa[[1,2]]=55555;
aa
{a[1,1],a[1,2]}
Clear[aa]
aa
{a[1,1],a[1,2]}

Maticove nasobeni

aa.bb

bb.aa

* u matic znamena nasobeni po slozkach (Hadamarduv soucin):

aa*bb

aa*Transpose[bb]

{cc cc,(cc^2)}

Pouziti cisla π

{N[π],N[π,100000]}

((N[π,1000000])/(N[π,100000]))

{5 π,(4/7)π,(4/7.)π,"\n"N[(4/7)π,100],"\n"N[(4/(7.`100))π,100],"\n"(4/(7.`100))π}

1.7951958020513104`

Print["\n"N[(4/7)π,100],"\n"N[(4/(7.`100))π,100]]

Vnitrni struktura systemu Mathematica a jeho jazyk  I

Vyrazy (Expressions)

Hlavicky

{Head[7],Head[(1/7)],Head[7.0],Head[2+I],Head[π],Head[Pi]}

Atomy = Symbol ∨ Cislo ∨ Retezec (String)

h=Head;
{h[f],h[Integrate],h[π],h[EulerGamma],h[hj2]}

h["Toto je retezec."]

Struktura vyrazu h[(Subscript[e, 1]),...,(Subscript[e, n])], kde h je hlavicka a (Subscript[e, 1]),...,(Subscript[e, n]) jsou vyrazy pripadne atomy.

Head[a+b+c]

FullForm[a+b+c]

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

Length[{a,b,c}]
{Length[{{a,b,c}}],Dimensions[{{a,b,c}}]}
Length[a+b+c]

expr=Sin[x](a (x^2)+b x +c)

Head[expr]
FullForm[expr]
TreeForm[FullForm[expr]]

FullForm[(5/7)]

{FullForm[(5./7)],FullForm[HoldForm[(5./7)]]}

Casti vyrazu( [[ se dosahne escape sekvenci Esc[[Esc)

soucet=a+b;
FullForm[soucet]
{Part[soucet,1],Part[soucet,2],Part[soucet,3],Part[soucet,0]}
{soucet[[1]],soucet[[2]],soucet[[3]],soucet[[0]]}

TreeForm[soucet]

Clear[aa,a];
aa=Table[a[i],{i,-4,4}]
{aa[[0]],aa[[5]]}

{1,4,7}[[2]]

"Co dostanu ted?"//Length
FullForm["Co dostanu ted?"]
Head["Co dostanu ted?"]

1[[1]]

{"Co dostanu ted?"[[1]],"Co dostanu ted?"[[0]]}

{π[[1]],π[[0]]}

expr

expr[[2,1]]

Table[expr[[j]],{j,0,3}]

{expr[[1,3]],expr[[1,3,2]],expr[[1,3,2,2]]}

(c+b x+a (x^2)) Sin[x]//TreeForm

Definice

Okamzite a odlozene prirazeni (globalni pravidlo, Immediate and delayed assignment): "=" a ":="

SeedRandom[13]
r1=Random[]
SeedRandom[13]
r2:=Random[]

{r1,r1,r1}
{r2,r2,r2}
SeedRandom[13]
{Random[],Random[],Random[]}
{r2,r2,r2}

SeedRandom[13]
Table[Random[],{6}]

?r1

?r2

Nahodna cisla z jinych rozdeleni

Histogram[Table[Random[NormalDistribution[0,1]],{1000}]]

Nahodna cisla z diskretniho rovnomerneho rozdeleni

Histogram[Table[RandomVariate[DiscreteUniformDistribution[{100,999}]],{9000}]]

Funkce definovane po castech

Clear[g];
g[x_]:=(x^3)/;(*za podminky *)x≤0
g[x_]:=x/;0<x≤1
g[x_]:=Sin[x]/;x>1
Plot[g[x],{x,-1,2}]

Funkce definovane po castech, Exclusions, periodicke funkce

Clear[g1];
g1[x_]:=Piecewise[{
{(x^3),x≤0},
{x,0<x<1},
{Sin[x],x≥ 1}
				}]
Plot[g1[x],{x,-1,2},Exclusions->{x==1}]

Clear[g2];
g2[x_]:=Which[x<0,(x^3),
			0≤x<1,x,
			x≥1,Sin[x]
			]
Plot[g2[x],{x,-1,2}]

Clear[g3];
g3[x_]:=Which[x≤3&&x≥2,Print["Haf"],
			  True,Print["Haf, Haf"]]
{g3[0],g3[2.5]}

g3[a]

Definice periodickych funkci I

Clear[g4];
g4[x_]:=Mod[x,2]-1/;-1≤Mod[x,2]-1<0
g4[x_]:=(((Mod[x,2]-1))^2)/;0≤Mod[x,2]-1<1
Plot[g4[x],{x,-10,10}]

?g2

D[g2[x],x]

Plot[D[g2[x],x],{x,-1,2}]

Zde je videt dalsi zaludnost systemu: bez aplikace funkce Evaluate se funkce nezobrazi a objevi se chybova zprava. Viz Help k funkci Evaluate. Souvisi to s atributy funkce Plot. Prilis to nemusite studovat, spokojte se s tim, ze kdyz nektera funkce se nechce zobrazit, aplikujte na ni Evaluate.

Attributes[Plot]

Nyni jiz funguje :

Plot[[D[g2[x],x]]//Evaluate,{x,-1,2}]

Pouziti Evaluate

der2[x_]:=D[g2[x],x]
?der2
Plot[der2[x],{x,-1,2}] (* s odlozenym prirazenim nefunguje. *)

der2i[x_]=D[g2[x],x]
?der2i
Plot[der2i[x],{x,-1,2}] (* s okamzitym prirazenim funguje. *)

Dalsi mozny, ale velmi spatny zpusob je pomoci funkce If.

g33[x_]:=If[x<0,(x^3),If[0≤x<1,x,If[x≥1,Sin[x]]]]
Plot[g2[x],{x,-1,2}]

g44[x_]:=If[x<0,(x^3),If[0≤x<1,x,If[True,Sin[x]]]]

Odstraneni svislych carek v grefech

Plot[Tan[x],{x,-5,5},Exclusions->Tan[x]==ComplexInfinity]

Plot[Tan[x],{x,-5,5},Exclusions->Cos[x]==0]

Plot[Tan[x],{x,-5,5},Exclusions->True]

Konec 28. 2. 2011

Lekce 3 (7. 3. 2011)

FAQ

Definice periodickych funkci II

Necht f(x) je funkce definovana na [a,b), a g(x) jeji periodicke rozsireni.

Clear[gperiod];
gperiod[f_,x_,a_,b_]:=f[a+(b-a)Mod[((x-a)/(b-a)),1]]
Clear[p1,p2];
p1[x_]:=Which[-1≤x≤0,x,
			     0<x<1,(x^2)]
Plot[gperiod[p1,x,-1,1],{x,-5,5}]
p2[x_]:=Sin[x]
Plot[gperiod[p2,x,0,π],{x,-15,15}]
Plot[gperiod[p2,x,0,2π],{x,-14,15}]

Plot[Sin[x]-gperiod[p2,x,0,2π],{x,-500,500},PlotRange->All]

{Max[Abs[Table[Sin[x]-gperiod[p2,x,0,2π],{x,5000000000000π,5000000000005π,0.01}]]],Min[Abs[Table[Sin[x]-gperiod[p2,x,0,2π],{x,5000000000000π,5000000000005π,0.01}]]]}

Plot[gperiod[p2,x,0,2π],{x,5000000000000π,5000000000003π}]

{Sin[5000000000000π],Sin[5000000000000.π],Sin[5000000000000`30.π]}

Table[{j,Mod[0.5+(10^j),1]},{j,20}]
Table[{j,Mod[0.5`30+(10^j),1]},{j,20}]

Button["Dotazy",Print[Date[],"\n",Date[][[{1,2,3}]],Text[Style["\nPrávě teď očekávám dotazy.\n",Red,30]]]]

Speak["Are there any questions?"]
Speak["Jsou nejake dotazy?"]

ListAnimate[Table[Text[Style["???  Zase je všechno jasné?",Blue,i]],{i,10,50}]]

Button["OK",Print[Text[Style["Tak to je dobře ...",Red,30]]]]

Kazda vestavena funkce ma atribut Protected

{N=2,Sin=sin}

Attributes[sin]

Attributes[Sin]

Unprotect[Sin];
Sin[x_]:=5+Cos[x]
Protect[Sin];
Sin[0]
Plot[Sin[x]//Evaluate,{x,-π,π}]

Sin[x_]=5+Cos[x]

Plot[Sin[x],{x,-π,π}]

?Sin

Baliky programu (Packages)

Needs["Calendar`"]

1+1

?Calendar`*

{DaysBetween[Date[],{1986,12,7}],DaysBetween[{1986,12,7},Date[]]}

?Global`*

{DaysPlus[Date[],100],DaysPlus[Date[],-100]}
{Take[Date[],3],Drop[Date[],-3],Take[Date[],{4,6}]}

Vnitrni struktura systemu Mathematica a jeho jazyk  III

Predikaty a booleovske operace

Predikaty (predicates)

?*Q

{PrimeQ[144],NextPrime[10000000],PrimeQ[NextPrime[10000000]]}

{OddQ[21],EvenQ[21]}

IntegerQ[5/9]

{NumericQ[π],NumericQ[∞]}

{NumberQ[5/9],NumberQ[π],NumberQ[1+4.5I]}

AtomQ["asbhjo \\[dfuhdb\\. n"]
"asbhjo \\[dfuhdb\\. n"

{ListQ[{a,{b,c}}],VectorQ[{a,{b,c}}],MatrixQ[{{a,b},{b,c}}]}

MemberQ[{a,b,c},b]

IntervalMemberQ[Interval[{3,4}],π]

?Interval*

Promyslete:

MatchQ[12345,_Integer]

{MatchQ[Plus[a,],Plus[_,___]],MatchQ[Plus[a],Plus[___]],MatchQ[Plus[],Plus[___]]}

MatchQ[(x-1)(1 + 2 x + 3 x^2),Plus[_,__]]
(x-1)(1 + 2 x + 3 x^2)//FullForm

MatchQ[Expand[(x-1)(1 + 2 x + 3 x^2)],Plus[_,__]]
Expand[(x-1)(1 + 2 x + 3 x^2)]//TraditionalForm
Expand[(x-1)(1 + 2 x + 3 x^2)]//FullForm

?GoldenRatio

Uziti Simplify

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

D[int1,x]
D[int1,x]//Simplify

1+1/GoldenRatio-GoldenRatio
1+1/GoldenRatio-GoldenRatio==0
1+1/GoldenRatio-GoldenRatio=== 0
1+1/GoldenRatio-GoldenRatio//Simplify

MatchQ[1+1/GoldenRatio-GoldenRatio,0]

MatchQ[Simplify[1+1/GoldenRatio-GoldenRatio],0]

Vzory (patterns)

x_, x__, x___, x_Head

Clear[f];f[x_]:=(x^2)
{f[a],f[{a,b}],f[a,b]}

Clear[f1];f1[x_List]:=(x^2)
hodnoty1={f1[a],f1[{a}],f1[{a,b}]}

Chceme vypreparovat (a^2):

Cases[hodnoty1,{(a^2),___}]

{hodnoty1[[2]],hodnoty1[[2]][[1]],Part[hodnoty1,2]//First,Most[hodnoty1]//Last//First}
{Take[hodnoty1,{2}],Take[hodnoty1,{2}]//First//First}

Clear[f2];f2[x_:1]:=x
{f2[3],f2[]}

Clear[f3];
f3[f_,x0_:1]:=FindRoot[f[y]==y,{y,x0}]
{f3[Sin,2],f3[Sin],f3[Cos]}

Clear[f4];
f4[x_/;x<0]:=-1
f4[x_/;x==0]:=0
f4[x_/;x>0]:=1
Table[f4[i],{i,-4,4}]

Booleovske funkce a logicke operace

7<5

Equal[3,7-4,6/2]

3==7-4==6/2
HoldForm[3==7-4==6/2]//FullForm

expr=((x^2)-1== (((x^4)-1)/((x^2)+1)))

Tan[z]==((Sin[z])/(Cos[z]))

Nekdy je nutne pouzit funkci Simplify.

(((x^4)-1)/((x^2)+1))//Simplify
expr//Simplify

Reduce[(x^2)-1== (((x^4)-1)/((x^2)+1)),x]

{6≠12/2,Unequal[6,12/2]}

Doslovne (literal) porovnani (===)

Remove[a,b];
{a==b,a===b,a===a}

Logicke spojky

{p||q,p&&q,!p}

{3<5||5<3,3<5&&5<3,!(3<5)}
{3<5∨5<3,3<5∧5<3,¬(3<5)}

Truth Table (Example of Distribute and a Table) (Exercise)

TruthTable[vars_List,expr_]:=Distribute[Table[{True,False},{Length[vars]}],List,List,List,{##,Function[vars,expr][##]}&]
DisplayTruthTable[vars_List,expr_]:=With[{t=ReplaceAll[TruthTable[vars,expr],{True->"T",False->"F"}]},Grid[Prepend[t,Append[vars,expr]],GridFrame->True,RowLines->Prepend[Table[0,{Length[t]-1}],2],ColumnLines->Append[Table[0,{Length[vars]-1}],2]]]

Sestrojte pravdivostni tabulku vyrazu (A∨B)⟹C

DisplayTruthTable[{A,B,C},Implies[Or[A,B],C]]

DisplayTruthTable[{A,B,C,D},Implies[Or[A,B],Or[C,D]]]

DisplayTruthTable[{A,B,C,D},Implies[(A∨B),(C∨D)]]

DisplayTruthTable[{A,B},A∧B]

DisplayTruthTable[{A,B},Xor[A,B]]

Atributy

Attributes[Plus]

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

Plus[Plus[a,b]]//FullForm

Attributes[Sin]

Range[10]

Sin[{1,2,3}]
Sin[(π/4)(Range[10]-1)]

Remove[druhemocniny];
druhemocniny[x_]:=(x^2)+Length[x]
druhemocniny[Range[3]]
Map[druhemocniny,Range[3]]
{Length[3],druhemocniny[3]}
Attributes[druhemocniny]
SetAttributes[druhemocniny,Listable]
Attributes[druhemocniny]
druhemocniny[Range[3]]

Clear[f6];
f6[x_]:=Sin[x].Cos[x]

f6[{1,2}]

Remove[plus];
plus[a,plus[b,c]]
SetAttributes[plus,Flat]
plus[a,plus[c,b]]
SetAttributes[plus,Orderless]
plus[a,plus[c,b]]

?IntegerQ

?Positive

Attributes[IntegerQ]

Funkce koncici na Q daji odpoved True nebo False.

{IntegerQ[-3],IntegerQ[a]}
IntegerQ[{-3,a}]

Attributes[Positive]

Positive[{5,a,-5}]

IntegerQ[{1,2.5,3}]

Tuto funkci ucinime listovatelnou:

Clear[integerQ];
integerQ[x_]:=IntegerQ[x]
integerQ[{1,2.5,3}]
SetAttributes[integerQ,Listable]
integerQ[{1,2.5,3}]

Pozor, Clear nezrusi atributy u symbolu integerQ!!! Nejjednodussi je uplne odstraneni symbolu pomoci Remove.

Clear[integerQ];
?integerQ

Attributes[IntegerQ]

Ryzi, anonymni funkce (Pure functions)

Ryzi funkce ma formu Function[promenna, predpis]

Function[x,(x^2)]

f=Function[x,(x^2)]
{f[x],f[3],f[I]}

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

v uspornem zapise  # znaci promennou ryzi funkce a & znaci ryzi funkci. Ekvivalentni zapis druhe mocniny

((#)^2)&

{((#)^2)&[x],((#)^2)&[3],((#)^2)&[I]}

{Function[{x,y},x+y][cc,bbbb],Function[#1+#2][cc,bbbb],
(#1+#2)&[cc,bbbb]}

Function[x,x.x][{aa,bb,cc}]

Function[x,x.x][qq]
Function[x,x.x][{qq}]
Function[x,x.x][(({{a, b, c}, {d, e, f}, {g, h, i}}))]//MatrixForm
Function[x,x.x][(({{a, b, c}, {d, e, f[g]}, {g, h, i}}))]

Definice ryzich funkci vice promennych

Function[{x,y},x+(y^2)][a,b]
(#1+((#2)^2))&[a,b]

Co se stane, kdyz nasobime matice, ktere nemaji stejne stycne rozmery?

(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}})).(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))
(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))

Off[Dot::"dotsh"]
(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}})).(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))
On[Dot::"dotsh"]

Seznamy (Lists) I

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

{b,c,a,aa}⋃{d,e,f,e,a}
{b,c,a,aa}⋂{d,e,f,e,a}
Complement[{a,b,c,d,e},{a,c},{d}]

Konec 7. 3. 2011

Lekce 4 (14. 3. 2011)

FAQ

Rozdil mezi DaysBetween a DateDifference

Needs["Calendar`"]

DaysBetween["Jan. 1, 2000","July 4, 2004"]
DateDifference["Jan. 1, 2000","July 4, 2004","Week"]

Seznamy (Lists) II: Pokracovani

r8=RandomInteger[{-4,4},15]
Sort[r8]
{Select[r8,Positive],Select[r8,#>0&]}

Table[(i^2),{i,1,5}] (* spatne *)
Table[(i^2),{i,5}] (* spatne *)
((Range[5])^2)(* spravne *)

Range[4,8]
malaabeceda=CharacterRange["a","z"]
CharacterRange["a","z"]//FullForm
{a,b,c}
{a,b,c}//FullForm
{malaabeceda[[1]],a,malaabeceda[[1]]===a}

t1=Table[3k,{k,1,10,2}]
Length[t1]
Range[3,27,6]

Table[{i,Expand[(((1+α))^i)]},{i,3}]//TableForm

{Range[3],Expand[(((1+α))^(Range[3]))]}//Transpose//TableForm

Table[i+j,{i,5},{j,i}]//TableForm

{(∑_(i=1)^n)(∑_(j=1)^i)1,(∑_(i=1)^n)i,Total[Range[100]],Apply[Plus,Range[100]]}
{(∑_(i=1)^n)(i^6),(∑_(n=1)^∞)(1/n),(∑_(n=1)^∞)(1/(n^2)),(∑_(n=1)^∞)(1/(n^4))}

Sum[1,{i,n},{j,i}]

Integrate[x y,{x,0,1},{y,0,x}]

Length[mat24={Range[4],((Range[4])^2)}]
Length[mat42={Range[4],((Range[4])^2)}//Transpose]
{Dimensions[mat24],Dimensions[mat42]}

Dimensions[r10={Range[4],((Range[4])^2)}]
Dimensions[{Range[3],((Range[4])^2)}]
Dimensions[{Range[3],((Range[4])^2)},Infinity]
Dimensions[{{{{a,b}}}}]

FullForm[1+(x^2)]

a+(1+x^2)^2//TreeForm

Position[{1+x^2,1+x,5,x^4,a+(1+x^2)^2},x^_] (* x sammotne neni ve forme x^ !!!*)

r10

{r10[[0]],r10[[2]],r10[[2,1]],r10[[All,2]],r10[[2,All]],r10[[6]]}
{Sin[x][[0]],Sin[x][[1]]}

Remove[mat];
matmatrixform=(mat=Table[(Subscript[a, i,j])=4,{i,3},{j,4}])//MatrixForm

%.Transpose[%]

matmatrixform.Transpose[matmatrixform]

Musi byt v zavorce, jinak by mat melo hlavicku MatrixForm, a ne List!!!

mat.Transpose[mat]

{(Subscript[a, 1,1])//FullForm,(Subscript[a, 1,1])=.,(Subscript[a, 1,1])}

Timing[Table[Sin[π k/4],{k,0,150000}];]
Timing[Sin[(Range[150001]-1)π/4];]

Vektory a matice

vektor={a,b};
xx={(Subscript[x, 1]),(Subscript[x, 2])};
yy={(Subscript[y, 1]),(Subscript[y, 2])};
matice={{x,y},{z,v}};
vektor.matice
matice.vektor
{xx.yy,xx*yy,xx yy}

{Transpose[xx],Transpose[{xx}]}
Grid[Transpose[{xx}].{yy},Frame->All]
Transpose[{xx}].{yy}//Grid[#,Frame->All]&

Transpose[{xx}].yy
yy.Transpose[{xx}]

Grid[Outer[Times,xx,yy],Frame->All]

Vektorovy soucin

Cross[{x1,x2,x3},{y1,y2,y3}]//TableForm

(aa=Array[a,{2,3}])//MatrixForm

aa.Transpose[aa]//MatrixForm

IdentityMatrix[3]//MatrixForm

(dm=DiagonalMatrix[{x1,x2,x3}])//MatrixForm

dm+1

Diagonal[dm+1]

Priklad

Vytvoreni positivne definitni matice. Vytvorime nahodnou ctvercovou matici C. Ta je s. j. regularni Potom C.(C^T) je s. j. positivne definitni..

c[n_]:=RandomReal[{0,1},{n,n}]
c[2]
Clear[posdef];posdef[n_]:=Module[{pom},pom=c[n];pom.Transpose[pom]]
Det[posdef[50]]
Det[(p=posdef[50]).Inverse[p]]

Procedura s vystupem positivne def. matice a jejim maximalnim a minimalnim vlastnim cislem

Clear[posdefeig];
posdefeig[n_]:=Module[{pom=c[n],eigvalues},
{pom=pom.Transpose[pom],Max[eigvalues=Eigenvalues[pom]],Min[eigvalues]}]
SeedRandom[131313456]
posdefeig[2]
posdefeig[2]
posdefeig[2][[2]]

prkhpfadoknhojsfpnb=%%%

?pom

Eigenvalues[p]

cisla=HilbertMatrix[2]//Eigenvalues
vektory=HilbertMatrix[2]//Eigenvectors
eig1=HilbertMatrix[5]//Eigenvalues
N[eig1]

HilbertMatrix[5]//MatrixForm

Map[Norm,vektory]//N
First[vektory].Last[vektory]

cisla50=Eigenvalues[p];
vektory50=Eigenvectors[p];

Map[Norm,vektory50]

Ukazka spatneho zadani. Nutno pouzit Abort, pripadne Quit Kernel:

Outer[Dot,vektory50,vektory50]

p5=posdef[5];
cisla5=Eigenvalues[p5];
vektory5=Eigenvectors[p5];
tup5=Tuples[vektory5,2];Length[tup5]
tup5[[{1,2}]](* {(Subscript[x, 1]),(Subscript[x, 1])} a {(Subscript[x, 1,])(Subscript[x, 2])} *)
Map[Dot[#[[1]],#[[2]]]&,tup5]//Chop

Transpose[vektory5].DiagonalMatrix[cisla5].vektory5-p5

Transpose[vektory5].DiagonalMatrix[cisla5].vektory5-p5//Chop

Odmocninova matice

Spektralni rozklad symetricke realne matice

Transpose[vektory5].DiagonalMatrix[(Sqrt[cisla5])].vektory5-MatrixPower[p5,(1/2)]//Chop

Timing[HilbertMatrix[30]//Eigenvalues//N;]
Timing[HilbertMatrix[30]//N//Eigenvalues;]

Ukazka programu viz napr. C:\Program Files\Wolfram Research\Mathematica\7.0\AddOns\Packages\Calendar

Testovani a hledani v listech

{Position[li={a,b,c,d,e},c],Position[{a,b,c,d,e},cc]}

Count[li2={a,b,a,c,a,d,c,b,b,b,a,e,e},a]

Union[li2]

Cetnosti prvku v listu li2

{Union[li2],Map[Count[li2,#]&,Union[li2]]}//Transpose//TableForm

Avsak  existuje zabudovana funkce Tally:

Tally[li2]

{MemberQ[li2,a],MemberQ[li2,f]}

$ExportFormats

Seznamy (Lists) II: Casti seznamu

Obecne

r11=Range[7];
Take[r11,2]
{Take[r11,-2],Drop[r11,2]}
Take[r11,{3,5}]
(* Pracuje s libovolnou hlavickou, nejenom s Head, jako i mnoho jinych funkci *)
Take[a+b+c+d+e+f,3]

Take[(({{11, 12, 13, 14, 15}, {21, 22, 23, 24, 25}, {31, 32, 33, 34, 35}, {41, 42, 43, 44, 45}, {51, 52, 53, 54, 55}})),{2,4},{3,5}]//MatrixForm

Take[(({{11, 12, 13, 14, 15}, {21, 22, 23, 24, 25}, {31, 32, 33, 34, 35}, {41, 42, 43, 44, 45}, {51, 52, 53, 54, 55}})),{1,-1,2},{1,-1,2}]//MatrixForm

Drop[r11,2]
Drop[r11,-2]

Viz Help k Take.

Priklad (Vypocet vynosu z vektoru cen aktiv)

Clear[ceny,s];
ceny=Table[s[i],{i,10}]

(vynosy(* returns *)=((Rest[ceny]-Most[ceny])/(Most[ceny])))//TraditionalForm

Konec 14. 3. 2011

Lekce 5 (21. 3. 2011)

FAQ

Formaty pro export a pro import

$ExportFormats(* //Length *)

$ImportFormats (* //Length *)

tutorial/ImportingAndExportingData

Pro “R” se doporucuje “Table”, “CSV”, “Text”. Otazka je, proc to delat, kdyz v systemu  Mathematica jsou statisticke funkce zabudovany a je jich vice nez v “R”?

Odstraneni problemu s vypisovanim cesty k souboru

$Path

Directory[]
SetDirectory[NotebookDirectory[]]
Directory[]

Tabulky s ruznym poctem poli v radku

Grid

Grid[{{1,2,3},{4,SpanFromLeft,6},{7,8,9}},Frame->All]

Seznamy (Lists) III: Retezce

Vybrane operace (viz take tutorial/CharactersInStrings, tutorial/OperationsOnStrings)

s1=StringJoin["Toto" ," ","jsou" ," retezce","."]
(* alternativne *)
"Toto"<>" "<>"jsou"<>" retezce"<>"."

Jednotlive znaky a znovu dohromady:

ch1=Characters[s1]
StringJoin[ch1]
StringLength[%]
Length[%%]
(* Vyzkousejte FullForm *)

StringReverse[{"ALENA","ROMA"}]

Viz take Reverse:

Reverse[{a,b,c}]
Reverse[Range[7]]
Reverse[{{a,b},{c,d}}]

CharacterRange["1","4"]

Prostudovat Outer:

indices=Outer[StringJoin,CharacterRange["1","4"],CharacterRange["1","8"]]//FullForm

List[List["11","12","13","14","15","16","17","18"],List["21","22","23","24","25","26","27","28"],List["31","32","33","34","35","36","37","38"],List["41","42","43","44","45","46","47","48"]]

%//ToExpression

Priklad (Aktualita: matice prechodu ratingovych trid, Credit Rating Transition Probabilities)

Zpracovani dat ze souboru PDF!!!!

Source

Dimitrios Kavvathas 
Goldman Sachs Group, Inc.
Estimating Credit Rating Transition Probabilities for Corporate Bonds,
AFA 2001 New Orleans Meetings

http://papers.ssrn.com/sol3/papers.cfm?abstract_id=252517

Pozor, po nasledujicim prikazu bude mit vytvoreny notebook 166M!!!

Import["CreditTransitionProbabilities.pdf"]

Processing

mat1="AAA AA A BBB BB B CCC D
AAA 94.46% 4.77% 0.60% 0.10% 0.07% 0.01% 0.00% 0.00%
AA 0.49% 91.79% 6.94% 0.59% 0.07% 0.09% 0.02% 0.01%
A 0.10% 1.87% 91.99% 5.22% 0.55% 0.25% 0.01% 0.01%
BBB 0.04% 0.32% 5.47% 88.19% 4.96% 0.83% 0.07% 0.11%
BB 0.04% 0.14% 0.84% 6.75% 83.88% 7.13% 0.72% 0.50%
B 0.01% 0.07% 0.27% 0.60% 5.62% 84.79% 4.81% 3.84%
CCC 0.12% 0.01% 0.51% 0.46% 1.70% 8.91% 54.01% 34.28%";

mat2=StringReplace[mat1,"%"->" "]

(m3=StringReplace[mat2,x:NumberString..:>ToExpression[x]])//FullForm
h1=First[m3]
(h2=StringReplace[h1,"\n"~~__->""])//FullForm
(h3=StringSplit[h2])//FullForm

m4=Select[m3,NumberQ]

m5=Apply[List,m4]/100

m6=Append[m5,{Table[0.,{7}],1.}]//Flatten;

m7=Map[PaddedForm[#,{4,4}]&,m6]

m8=Partition[m7,8]

mh9=Append[h3,m7]//Flatten

mh10=Append[{"   ",h3}//Flatten,Transpose[Partition[mh9,8]]//Flatten]//Flatten

mh11=Partition[mh10,9]//Transpose

Text@Grid[mh11,Frame->All,Alignment->Left]

(Text@Grid[mh11,Frame->All,Alignment->Left])//TeXForm

Ruzne poznamky

Zpetna carka za cislem: jedna se o cislo realne, nepresne:

{2.,FullForm[2.],2.`30}

Rationalize[2.`]

Attributes[StringReplace]
StringReplace[{"AAA","BAB"},"B"->""]
StringReplace["AAA","BB","B"->""]
StringReplace[StringExpression["AAA","BB"],"B"->""]//FullForm

Map[ToExpression,{"AOOOO","3.4"}]//FullForm

mat2=StringReplace[mat1,"%"->" "];
(mat3petr=StringSplit[mat2])//FullForm
mat4petr=Map[ToExpression,mat3petr]
Partition[{" ",mat4petr}//Flatten,9]

Ryzi, anonymni funkce (Pure functions)

Dalsi ukazka aplikace ryzi funkce

data=Range[20];
(* Nasledujicimi dvema prikazy vytvorime velmi odlehla pozorovani (outliers) *)
data[[5]]=1000;
data[[15]]=-1000;
ListPlot[data,PlotStyle->PointSize[0.02]]

Vnitrni algoritmus pro grafiku extremy nezobrazi. Muzeme vsak pouzit nasledujici Option, pak se nam vsak s nejvetsi pravdepodobnosti ztrati to, co jsme chteli videt.

ListPlot[data,PlotStyle->PointSize[0.02],PlotRange->All]

Pouzitim ryzi fce v zabudovane funkci Select vyradime odlehla pozorovani. Za odlehla pozorovani budeme uvazovat ta, ktera jsou v avbsolutni hodnote ≥ 30. Podivejte se do Help na ListPlot, pro zduvodneni nasledujiciho vystupu.

withoutoutliers=Select[data,Function[x,Abs[x]<30]]
withoutoutliers=Select[data,Abs[#]<30&]
ListPlot[withoutoutliers,PlotStyle->PointSize[0.02]]

Nahradte vsechny prvky v absolutni hodnote vetsi nez 29 nulou.

Map[Abs[#]>30&,data]
Map[Boole[(Abs[#]>30)]&,data]
data*(1-Map[Boole[(Abs[#]>30)]&,data])

Pridame x-souradnice

data2={Range[20],((Range[20])^2)}//Transpose
ListPlot[data2,PlotStyle->PointSize[0.01]]
data2[[5,2]]=-10000;
ListPlot[data2]

withoutoutliers2=Select[data2,Abs[#[[2]]]<1000&]
ListPlot[withoutoutliers2]

Map, Apply, Thread (provleceni?, vlakno?)  (Obecny princip)

Map preskoci hlavicku druheho argumentu a aplikuje prvni argument na to, co je za hlavickou.
Apply vymeni hlavicku druheho vyrazu za prvni.

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

Clear[fjh];
fjh[x_]:=NumberQ[x]
fjh[{1,2,3}]
Map[fjh,{1,2,3}]

Apply[Times,a+b+c]
Apply[Times,3+4+5]

Thread[f[{a,b,c}]]

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

Thread[f[{a,b,c},{x,y,z}]]

Convert equations for lists to lists of equations:

Thread[{a,b,c}=={x,y,z}]

Apply a function to both sides of an equation:

Thread[Log[x==y],Equal]

Map a Apply (specialni pripady)

Clear[f,g,a,b,c];
li={a,b,c};
soucin=a b c;
{Head[li],Head[soucin]}
{FullForm[li],FullForm[soucin]}

Map preskoci hlavicku druheho argumentu a applikuje prvni argument na to, co je za hlavickou druheho argumentu

Map[Plus,soucin]

Apply[Plus,soucin]

🐺

🐺//FullForm

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

Plus[{a,b,c}]

{Total[{a,b,c}],Apply[Times,{a,b,c}]}

Apply je funkce, ktera vymeni hlavicku druheho argumentu za prvni argument

Apply[Times,li]
Apply[List,soucin]
Apply[Plus,soucin]

Priklad vypoctu aritmetickeho prumeru s pojistkou na to, aby argument funkce byl vektor delky alespon jedna. ( _? je test predlohy, PatternTest.) (Existuje zabudovana funkce Mean.)

prumer[x_]:=Apply[Plus,x]/Length[x]/;VectorQ[x]&&Length[x]>0
{prumer[{a,b,c}],prumer[{{a},b,c}],prumer[a,b,c],prumer[2]}

{{{a}},b,c}//Total

{{{a}},b,c}//Total//Flatten
{{{a}},b,c}//Flatten//Total

{{{a}},b,c}//Total//Flatten//First
{{{a}},b,c}//Total//Flatten//Flatten

{{{a}},b,c}//Flatten

Pouziti zabudovane funkce Total (pozor na formu vstupu!!!)

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

Clear[pr];
pr[data_]:=((Total[dt=data])/(Length[dt]))
pr[{{{a,d}},b,c}]

Clear[pr];
pr[data_]:=((Total[dt=data//Flatten])/(Length[dt]))
{pr[{{{a,d}},b,c}],pr[],pr[{}]}

Priklad (vypocet rozptylu). Zabudovana funkce Variance

rozptyl[x_?VectorQ]:=prumer[(((x-prumer[x]))^2)]
{rozptyl[Range[5]],rozptyl[{Range[5]}]}
ch5=CharacterRange["a","e"]
{prumer[ch5],prumer[{a,b,c,d,e}]}
(* Analyzujte pomoci FullForm *)
{Variance[ch5],Variance[{a,b,c,d,e}]}

Dalsi pouziti Map a Apply

Velmi dulezity priklad

f[x_,y_]:=x (y^2)
{f[2,1],f[{2,1}],Apply[f,{2,1}]}

ff[{x_,y_}]:=x (y^2)
{ff[2,1],ff[{2,1}]}

fff[x_,y_]:=x (y^2)
fff[{x_,y_}]:=x (y^2)
?fff

{fff[2,1],fff[{2,1}]}

{fff[{x,y},{z,w}],fff[{{x,y},{z,w}}],fff[{{x,y},{z,v,w}}]}

fff[🐺,🐺]

Priklad (Sin posloupnosti)

{Sin[],Sin[2.,1],Sin[{2.,1}]}

Clear[sin];
sin[y__]:=Sin[{y}]//Flatten
{sin[],sin[1.],sin[1.,2],sin[1.,2,3],{sin[1.,2],sin[1.,2,3]},
{sin[1.,2],sin[{{1.,2,3}}]},sin[{1.,2,3}]}

Clear[sin1];
sin1[y__]:=Sin[{y}]
{sin1[],sin1[1.],sin1[1.,2],sin1[1.,2,3],{sin1[1.,2],sin1[1.,2,3]},sin1[{1.,2,3}]}

Konec 21. 3. 2011

Lekce 6 (28. 3. 2011)

FAQ

Normal

Tayloruv rozvoj

expan1=Series[Exp[-x],{x,0,3}]
Normal[expan1]
expan2=Series[(n-1)!(*=Gamma[n]*),{n,∞,3}]
Normal[expan2]
Normal[expan2]//Simplify

Seznamy (Lists) IV

Pridavani, odstranovani a modifikovani prvku seznamu

Prepend[m={a,b,c},aa]
Append[m,aa]
AppendTo[m,aa]
m
AppendTo[m,aa]

Nasledujici Riffle je michani (karet). Neplest s Rifle.

Riffle[m={a,b,c},x]
Riffle[m={a,b,c},{x,y,z}]

ReplacePart[m,2->x]
mat🐺=ReplacePart[mat={{a,b,c},{d,e,f}},{1,2}->🐺]

{mat🐺[[2]]=🐺🐺🐺,mat🐺}

mat[[All,1]]={x,y};
mat

Kombinace listu

Join[{a,b,c},{x,y},{t,u},{b,c},a]
Join[{a,b,c},{x,y},{t,u},{b,a,c},{a}]
Join[{a,b,c},{x,y},{t,u},{{b,a},c},{a}]

Union[{a,b,c},{x,y},{t,u},{b,c},{a}]
DeleteDuplicates[{{a,b,c},{x,y},{t,u},{b,c},{a}}//Flatten]

Deleni a vyplnovani seznamu (partitioning and padding)

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

Deleni s prekrytim (offset)

Partition[{a,b,c,d,e,f,g},3,1]//MatrixForm
Partition[{a,b,c,d,e,f,g},3,3,1]
Partition[{a,b,c,d,e,f,g},3,3,1,x]
Partition[{a,b,c,d,e,f,g},3,3,1,{x,y,z,v,w}]
Partition[{a,b,c,d,e,f,g},3,1,{1,1},x]

Poznamky k indexovani

mat={{x,b,c},{y,e,f},{h,i,j}};
{mat[[1,2]],mat[[1]][[2]],mat[[{1,2}]],mat[[{1,2},1]],mat[[2;;3]]}

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

Nuance, ale tez jine moznosti viz Help:

Extract[mat,{1,3}]
Extract[mat,{{1,3},{1,2}}]
mat[[{1,3},{1,2}]]
mat[[{{1,3},{1,2}}]]

Konstrukce listu

NestList, FoldList, Range, Table, Array

NestList[f,x,5]

Nekolik moznosti, jak vytvorit posloupnost prvnich 6 mocnin cisla 2 (vcetne nulte):

NestList[2#&,1,5]
NestList[Function[x,2x],1,5]
(2^(Range[6]-1))
Array[(2^(#-1))&,6]

Select[{1,3,7,4,10,2},#<5&]

TakeWhile[{1,3,7,4,10,2},#<5&]

Pick[{a,b,c,d,e},{1,1,0,1,0},1]
Pick[{a,b,c,d,e},{1,1,0,1,0},0]

Priklad na postupna maxima (rekordy)

{First[r5=Range[5]],Last[r5],Most[r5],Rest[r5]}

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

FoldList[f,x,{a,b,c,d}]
FoldList[Plus,0,{a,b,c,d}]
r6=FoldList[Max,0,{5,3,6,3,2,77,4,888}]
r6modif=FoldList[Max,First[seznam={5,3,6,3,2,77,4,888}],seznam]
r7=Rest[r6]
Union[r7]

Clear[rekordy];
rekordy[x_?VectorQ]:=Union[Rest[FoldList[Max,First[x],x]]]
rekordy[{5,3,6,3,2,77,4,888}]
rekordy[-{5,3,6,3,2,77,4,888}]

Rozdil mezi NestList a FoldList

NestList nemuze pracovat se suoborem napr. {a,b,c,d,e}.

NestList[f,x,5]
FoldList[f,x,{a,b,c,d,e}]

Modelovani ceny akcie

Clear[cena];
cena[s0_?NumericQ,delka_?IntegerQ]:=FoldList[Plus,s0,RandomReal[{-s0/20,s0/20},delka]]

cenanest[s0_,delka_]:=NestList[#+RandomReal[{-s0/20,s0/20}]&,s0,delka]
cenanest2[s0_,delka_]:=NestList[Function[x,x+RandomReal[{-s0/20,s0/20}]],s0,delka]

SeedRandom[131313];
cena[100,10]
SeedRandom[131313];
cenanest[100,10]
SeedRandom[131313];
cenanest2[100,10]

Cviceni. Pro stejny ucel pouzijte funkci Accumulate

ListPlot[Table[cena[100,1000],{5}],Joined->True,PlotRange->All]

Cviceni. Ve vyse uvedenem priklade jsme pouzili maximalni rust  resp. minimalni pokles o 5 %. Vytvorte proceduru, kde dalsim parametrem bude _procento znamenajici  maximalni rust  resp. minimalni pokles.

Map[Length,{49+Range[1001],cena[100,1000]}]

ListPlot[Table[{49+Range[1001],cena[100,1000]}//Transpose,{5}],Joined->True,PlotRange->All,AxesOrigin->{0,Automatic}]

Rust nebo pokles z predchozi (a ne pocatecni) ceny

cenanestdynamicky[s0_,delka_]:=NestList[#+RandomReal[{-#/20,#/20}]&,s0,delka]

SeedRandom[131313];
ListPlot[Table[cenanestdynamicky[100,1000],{5}],Joined->True,PlotRange->{0,160}]
SeedRandom[131313];
ListPlot[Table[cena[100,1000],{5}],Joined->True,PlotRange->{0,160}]

Cviceni a priklad (opakovani: funkce, ktera neni listovatelna)

NumericQ[1,π,Sin[(Sqrt[E])],a,Exp[x]]

NumericQ[{1,π,Sin[(Sqrt[E])],a,Exp[x]}]

Vsechna nasledujicic vyjadreni jsou ekvivalentni:

Map[NumericQ,{1,π,Sin[(Sqrt[E])],a,Exp[x]}]
Map[NumericQ[#]&,{1,π,Sin[(Sqrt[E])],a,Exp[x]}]
Map[Function[x,NumericQ[x]],{1,π,Sin[(Sqrt[E])],a,Exp[x]}]

Pridani listovatelnosti k NumericQ

Clear[numericQ];
numericQ[x_]:=NumericQ[x];
SetAttributes[numericQ,Listable]
numericQ[{1,π,Sin[(Sqrt[E])],a,Exp[x]}]

Sparse arrays (Ridka pole, matice)

sa1=SparseArray[{3->x,4->y},15]

Normal[sa1]

sa1.sa1
sa1*sa1

Normal[sa1*sa1]
Normal[sa1*sa1]//FullForm
Normal[sa1*sa1]//Head
HoldForm[Normal[sa1*sa1]]//FullForm
Head@(4 5)

SparseArray[{i_->(i^2)},5]//Normal

SparseArray[{2->a,5->b}]//Normal

SparseArray[{{1,2}->a,{3,2}->b,{3,3}->c}]//Normal

(tab1=Table[{i,j},{i,4},{j,i}])//Length
tab1
Flatten[tab1,1]
(SparseArray[(indexy1=Flatten[tab1,1])->Range[Length[indexy1]]]//Normal)//MatrixForm

SparseArray[{Band[{1,1}]->x,Band[{2,1}]->y},{5,5}]//MatrixForm

SparseArray[{i_,i_}->1,{4,4}]//MatrixForm

SparseArray[{{i_,i_}->1,{_,3}->a,{2,_}->b},{4,4}]//MatrixForm

SparseArray[{i_,i_}->RandomReal[],{4,4}]//MatrixForm
SparseArray[{i_,i_}:> RandomReal[],{4,4}]//MatrixForm

(* Vektor, ktery bude mit na sudych pozicich p a na ostatnich druhou mocninu *)SparseArray[{_?EvenQ->p,i_-> (i^2)},{15}]//Normal

(* V prvnim nebo tretim radku a ve druhem nebo ctvrtem sloupci a *)
SparseArray[{1|3,2|4}->a,{4,4}]//MatrixForm

SparseArray[i_/;3<i<7->p,{10}]//Normal

SparseArray[{i_,j_}/;Abs[i-j]==1->i+j,{5,5}]//MatrixForm

SparseArray[{i_,j_}/;Abs[i-j]<2->i+j,{5,5}]//MatrixForm

SparseArray[{Band[{1,1}]->2Range[5],Band[{1,2}]->2Range[4]+1,Band[{2,1}]->2Range[4]+1}]//MatrixForm

(sprec=SparseArray[{Band[{1,1}]->2Range[5],Band[{1,2}]->2Range[4]+1,Band[{2,1}]->2Range[10]+1}])//MatrixForm
Dimensions[sprec]

(sprec=SparseArray[{Band[{1,1}]->2Range[5],Band[{1,2}]->2Range[4]+1,Band[{2,1}]->2Range[10]+1},{10,10}])//MatrixForm
Dimensions[sprec]

Nespecifikovana pole jsou misto nul nahrazeny karetnimi symboly

SparseArray[{i_,j_}/;Abs[i-j]<2->i+j,{5,5},♡♢♣♠]//MatrixForm

♠

Obrázky

A collection of random squares:

Graphics[{EdgeForm[Black],Table[{Hue[RandomReal[]],Rotate[Rectangle[RandomReal[2,2]],RandomReal[2Pi]]},{20}]}]

A color wheel:

Animate[Graphics[{EdgeForm[Opacity[.5]],Table[{Hue[(t+k)/20,1,.9],Rotate[Rectangle[{1,-Sin[Pi/20]},{1+2Sin[Pi/20],Sin[Pi/20]}],2Pi t/20,{0,0}]},{t,20}]}],{k,1,20},AnimationRunning->False]

Digital petals:

Graphics[Table[{EdgeForm[Opacity[.6]],Hue[(-11+q+10 r)/72],With[{t=2Pi q/12},Rotate[Rectangle[(8-r)({Cos[t],Sin[t]}-1/4),(8-r)({Cos[t],Sin[t]}+1/4)],t+Pi/4]]},{r,6},{q,12}]]

Principy programovani

Definice (funkce logaritmus)

Definice je prepisovaci pravidlo, ktere je automaticky pouzito na libovolny vyraz. Priklad na definici logaritmu:

Clear[log];
log[x_*y_]:=log[x]+log[y]
{log[2x (y^2)],log[2x y y]}
{(y^2)//FullForm, (y y)//FullForm,HoldForm[y y]//FullForm}

Prestoze z matematickeho hlediska je (y^2) (Power) totez, co y*y (Times), potrebujeme dalsi prepisovaci pravidlo, protoze porovnani vzoru (pattern matching) je v zasade syntakticke, a ne semanticke!

log[((a_)^(n_))]:=n*log[a]
log[2x (y^2)]

log[x*2*y^3*z]

Clear[klaralog];
klaralog[((a_)^(n_))*b_]:=n klaralog[a]+klaralog[b]
klaralog[b_*((a_)^(n_))]:=n klaralog[a]+klaralog[b]
klaralog[x*2*y^3*(z^1)]

Vysvetleni

{FullForm[(y^1)],FullForm[y],FullForm[(y^2)],Power[y,1]}

Definice (linearni funkce na linearnim prostoru, M220, p. 63)

Clear[f,x,y,a];
f[c_?NumericQ*x_]:=c*f[x]
f[x_+y_]:=f[x]+f[y]
f[π+(Sqrt[2])x+3 a (x^2)]

Ale zcela v poradku to neni:

{f[2π],f[π 2]}

Pridame

f[c_?NumericQ]:=c f[1]/;c≠1
{f[2π],f[π 2]}
f[π+(Sqrt[2])x+3 a (x^2)]

Ruzne programovaci styly I

Nasledujici funkce zjistuji cetnosti v posloupnosti prirozenych nahodnych cisel 1,...,10. (Diskretni rovnomerne rozdeleni R(0,...,10)

Proceduralni styl

SeedRandom[130];
(* data=Table[Random[Integer,{1,10}],{100}]; *)
data=RandomInteger[{1,10},{100}];

proc[x_]:=Module[{bins=Table[0,{10}]},
Do[bins[[x[[i]]]]++,{i,1,Length[x]}];bins]
pr=proc[data]
Total[pr]

Funkcionalni styl

Clear[fun];
fun[x_?VectorQ]:=fun[Map[{#,1}&,x]]
fun[{1,2,3,2}]

Poznamka: _ prave jeden argument, __ jeden a vice argumentu, ___ nula a vice argumentu, x___

fun[{a___,{b_,n_},c___,{b_,m_},d___}]:=fun[{a,{b,n+m},c,d}]

SeedRandom[131313];
rc1=RandomChoice[Range[5],15]
fun[{1,2,3,2}]
fun[rc1]

fun[a_List]:=Map[Last,a];
data
fun[data]
fun[data//Sort]

fun[a_List]:=a//Sort;
cetnosti=fun[data]

?fun

Pochazeji tato cisla skutecne z diskretniho rovnomerneho rozdeleni?

Konec 28. 3. 2011

Lekce 7 (4. 4. 2011)

FAQ

Zjistit cetnosti pozorovani, ktere jsou napr. kladne

SeedRandom[130];
(* data=Table[Random[Integer,{-10,10}],{100}]; *)
dataneg1=RandomInteger[{-10,10},{100}];
dataneg2=Select[dataneg1,#>0&];
tallyneg1=Tally[dataneg2]
Total[Map[Last,tallyneg1]]

Ruzne programovaci styly II

SeedRandom[130];
(* data=Table[Random[Integer,{1,10}],{100}]; *)
data=RandomInteger[{1,10},{100}];

Mathematica styl

Table[Count[data,i],{i,10}]

Map[Count[data,#]&,Union[data]]

Pouziti funkce Tally

Tally[{c,b,a,a,a,b,a,c,b,a}]
Tally[data]//Sort

Use test argument to count expressions with the same Head:

Tally[{{a,b},{w,x,y,z},E,a,vvv,11111,{w,x,y,z},2,3,4,5,E},Head[#1]===Head[#2]&]

Map[Head,{{a,b},{w,x,y,z},E,a,vvv,11111,{w,x,y,z},2,3,4,5,E}]

ExampleData[]

ExampleData["Text"]

ExampleData[{"Text","UNHumanRightsChinese"}]

text=StringTake[ExampleData[{"Text","Hamlet"}],10000]

StringTake[ExampleData[{"Text","Hamlet"}],1000]

thamlet=Tally[Characters[text]]
Sort[Map[Reverse,thamlet]]//Reverse

SeedRandom[130];
soubor=RandomVariate[NormalDistribution[],{5,3}]
Sort[soubor,#1[[2]]<#2[[2]]&]
(* na vysvetlenou *)
SeedRandom[130];
soubor=RandomVariate[NormalDistribution[],{5}]
Sort[soubor]
Sort[soubor,#1>#2&]

Substitucni pravidla (lokalni pravidla, replacement rules, rule-based programming)

Lokalni pravidla se aplikuji zleva doprava, kazde jenom jednou

{1,2,3,2,3,1,2,4}/. 3->"a"/.x_Integer->0

(* Pocet vyskytu cisla 3 v uvedenem seznamu *)v1={1,2,3,2,3,1,2,4}/. 3->"a"/.x_Integer->0/."a"->1
Total[v1]

pravidla[data_List,n_]:=
Total[data/.n->"a"/.x_Integer->0/."a"->1]
pravidla[{1,2,3,2,3,1,2,4},2]
pravidla[{1,2,3,2,3,1,2,4},3]

Map[pravidla[data,#]&,Range[10]]
Map[{#,pravidla[data,#]}&,Union[data]]

char2=Characters["Map[pravidla[data,#]&,Range[10]] Map[pravidla[data,#]&,Union[data]]'"]

pravidla2[data_List,n_]:=
Count[data,n]
pravidla2[{"a","p","[","p","r","a","v"},"p"]

(cetnosti2={Union[char2],Map[pravidla2[char2,#]&,Union[char2]]}//Transpose)//TableForm

Pozor, pri importu souboru (napr. z Excelu) mohou vznikonout zaludnosti tohoto typu (opticky se mezera v obou pripadech tvari stejne).

{" "//FullForm," "//FullForm}

" "

Nejcetnejsi znak:

sort2=Sort[cetnosti2,#1[[2]]>#2[[2]]&]

Vybrat znaky s cetnosti 2:

vsechnycetnosti=Union[Map[Last,sort2]]//Reverse
Select[sort2,#[[2]]==2&]//Sort

Znaky se stejnymi cetnostmi jsou v jednom listu (ukazka vnorene ryzi (anonymni) funkce):

Map[Function[x,Select[sort2,#[[2]]==x&]//Sort],vsechnycetnosti]

Prehledne usporadani

(tab2=Map[Function[x,Select[sort2,#[[2]]==x&]//Sort],vsechnycetnosti])//TableForm

Partition[Flatten[tab2],2]//TableForm

Map[Reverse,Partition[Flatten[tab2],2]]//TableForm

Poznamka (rozdil v pouziti Reverse, v prvnim pripade na radky, v druhem pomoci Map na sloupce)

mat123=(({{a, b}, {c, d}, {e, f}}));
Reverse[mat123]//MatrixForm
Map[Reverse,mat123]//MatrixForm

Zakladni financni funkce I

PV: ekvidistatni casove intervaly mezi platbami

Uvodni poznamky

Poznamka 1

x=.;
NestList[Function[x,(x/(1+r))],c,4]
FoldList[Function[{x,r},(x/(1+r))],c,{(Subscript[r, 1]),(Subscript[r, 2]),(Subscript[r, 3])}]
1/FoldList[#1 * (1+#2)&,1,{(Subscript[i, 1]),(Subscript[i, 2]),(Subscript[i, 3]),(Subscript[i, 4])}]
1/FoldList[Function[{x1,x2},x1(1+x2)],1,{(Subscript[i, 1]),(Subscript[i, 2]),(Subscript[i, 3]),(Subscript[i, 4])}]

Poznamka 2 (vyhodnocovani substitucnich (lokalnich) pravidel: hodnoty promennych se nemeni):

{a,b,c,a}/.{a->s,s->d} (* ReplaceAll *)
{a,b,c,a}//.{a->s,s->d} (* ReplaceRepeated *)

Clear[x];
x/.x->x+1
x//.x->x+1
?x

Options[pv]=.;
interests=interestrates/.Flatten[{{interestrates->"variable"},Options[pv]}]
interests=interestrates/.Flatten[{Options[pv]}]
Options[pv]={interestrates->"constant"};
interests=interestrates/.Flatten[{Options[pv]}]
interests=interestrates/.Flatten[{{interestrates->"variable"},Options[pv]}]

Soucasna hodnota I

Predpoklada se, ze cf je penezni tok, cf={(Subscript[cf, 0]),(Subscript[cf, 1]),...,(Subscript[cf, T])}, i hodnotici urokova mira, ii={(Subscript[i, 1]),...,(Subscript[i, T])} promenne urokove miry pro obdobi 1,...,T, rr={(Subscript[R, 1]),...,(Subscript[R, T])} okamzite (spotove) urokove miry pro investice se splatnostmi 1,...,T, (vynosova krivka).

Clear[pv];
Options[pv]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
pv[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.Flatten[{opts,Options[pv]}];
Which[interests=="constant",cf.(1/(((1+i))^(Range[n]-1))),
	interests== "variable",cf.(1/(FoldList[#1 * (1+#2)&,1,i])),
	interests=="yieldcurve",First[cf]+Rest[cf].(1/(((1+i))^(Range[n-1])))]
	 ]

Options[pv]

cf=.;i=.;ii=.;rr=.;i1=.;i2=.;i3=.;
cf={cf0,cf1,cf2,cf3};
ii={i1,i2,i3};
rr={r1,r2,r3};
{pv[1,ii],pv[{},ii],pv[,ii]}

pv[cf,i]//TraditionalForm
pv[cf,ii,interestrates->"variable"]
pv[cf,rr,interestrates->"yieldcurve"]
pv[cf,i,in] (* in nemuze byt uvazovano jako option, vrati tedy pv[... *)

Nejsou sladeny argumenty a volby:

pv[cf,i,interestrates->"yieldcurve"]//TraditionalForm
(* Pozor, i je skalar!!! *)

Soucasna hodnota II (alternativni postup pri zpracovani voleb)

Clear[pv2,interestrates];
Options[pv2]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
pv2[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.{opts}/.Options[pv2];
Which[interests=="constant",cf.(1/(((1+i))^(Range[n]-1))),
	interests== "variable",cf.(1/(FoldList[#1 * (1+#2)&,1,i])),
	interests=="yieldcurve",First[cf]+Rest[cf].(1/(((1+i))^(Range[n-1])))]
	 ]

pv2[cf,i]//TraditionalForm
pv2[cf,ii,interestrates->"variable"]
pv2[cf,rr,interestrates->"yieldcurve"]
pv2[cf,i,in] (* in nemuze byt uvazovano jako option, vrati tedy pv[... *)
pv2[cf,ii,inter->"variable"]

IRR: ekvidistatni casove intervaly mezi platbami

Nejjednodussi pristup (vyhledani vsech korenu) spociva v pouziti funkce NSolve. Odpovidajici algebraicka rovnice je v sice v promenne v, ale Mathematica rozpozna, ze se jedna o algebraickou rovnici i kdyz je promennou urokova mira i.

Nejprve ukazka projektu, ktery asi zadny potencialni investor nebude povazovat za realny:

i=.
cfnerealny={-1000,3600,-4310,1716};
Plot[pv[cfnerealny,i],{i,0,0.4}]

Vnitrni mira vynosnosti se ziska takto:

i=.
irrsolv=NSolve[pv[cfnerealny,i]==0,i]

i/.irrsolv

Zde je podstatne, ze se pouzije numericka varianta Solve, tj. NSolve. Vsimneme si, co udela obecna procedura pro vyhledani korenu FindRoot (zde je nutne zadat pocatecni reseni x0)

i0=.;i=.;
Table[{i0,FindRoot[pv[cfnerealny,i],{i,i0}]},{i0,0,0.5,0.05}]//TableForm

Nyni priklad nahodneho penezniho toku. Uvadime tez prislusne statisticke charakeristiky

i=.;
SeedRandom[313131];
{Mean[cf13=RandomReal[{-25,33},{30}]],Variance[cf13],"   Smerodatna odchylka = ",StandardDeviation[cf13],Quantile[cf13,0.05],Quantile[cf13,0.95]}
{Min[cf13],Max[cf13]}
ListPlot[cf13,Joined->True]
Plot[pv[cf13,i],{i,-0.18,.5},PlotRange->All]
NSolve[pv[cf13,i]==0,i]

Samozrejme, ze nas nezajimaji imaginarni, a z hlediska ekonomickeho ani zaporna vnitrni vynosova procenta. Separace v nasledujici obecne procedure:

Clear[irr];
(* Predpoklada vyse definovanou funkci pv *)
irr[cf_]:=Module[{solutions,i},
solutions=i/.NSolve[pv[cf,i]==0,i];
Select[solutions,(Im[#]==0&&Re[#]>0)&]
				]

irr[cf13]

Alternativne 1 (Kovacova)

Clear[irrkovacova];
(* Predpoklada vyse definovanou funkci pv *)
irrkovacova[cf_]:=Module[{solutions,i},
Select[solutions=i/.NSolve[pv[cf,i]==0,i,Reals],Positive]
				]
irrkovacova[cf13]

Alternativne 2

solutions=i/.NSolve[pv[cf13,i]==0,i]
Cases[solutions,_?((Im[#]==0&&Re[#]>0)&)]
Cases[solutions,i_Real]
Cases[solutions,i_Real/;i≥0]

Alternativne 3

Cases[i/.NSolve[pv[cf13,i]==0,i],i_Real]
Cases[i/.NSolve[pv[cf13,i]==0,i],i_Real/;i≥0]

Obecne poznamky k presnosti vypoctu, ruzne I

Obecna poznamka k presnosti I - ilustrace na vypoctu IRR

Clear[irrsimple,i];
(* Predpoklada vyse definovanou funkci pv *)
irrsimple[cf_]:=Module[{solutions,i},
solutions=i/.NSolve[cf.(1/(((1+i))^(Range[Length[cf]]-1)))==0,i];
Select[solutions,(Im[#]==0&&Re[#]>0)&]
				]
irrsimple[cfnerealny={-1000,3600,-4310,1716}]

Zkopirovana hodnota 0.3:

0.3000000000000002`

irrsimple[cfnerealny2={-1000.05,3600,-4310,1716}]

Zkopirovana prvni hodnota

0.294064372011615`

Avsak presnost vypoctu na 30 mist dosahneme

irrsimple[cfnerealny3={-1000.05`30,3600,-4310,1716}]

Zkopirovana druha hodnota

0.294064372011614836801956765663300974059618719095`29.69897000433602

Obecna poznamka k presnosti II - vypocet funkce sin komplexniho argumentu

Sin[4.5`100+I]

Obecna poznamka k presnosti III - hledani pevneho bodu zobrazeni (porovnani s programem Excel)

3 Try in MS Excel: 0.2->C2. Then copy formula 11*C2-2 15 times. This results (dependent on the used computer) in:

0.20.20.20.20.20.20.20.20.20.20.20.2000050.2000510.2005570.2061320.2674570.9420288.36231289.98543987.839710864.24

4 The same in Mathematica is solved as follows:

fjo[x_]:=11*x-2
NestList[fjo,0.2,30]

Obviously, in rationals it works:

NestList[fjo,(1/5),30]

Look at

BaseForm[1/5,2]
BaseForm[0.2,2]
BaseForm[0.2`30,2]

But (e.g.):

fjo30=NestList[fjo,0.2`30,30]

Konec 4. 4. 2011

Lekce 8 (11. 4. 2011)

FAQ

Otazka 1

Vev vektoru z nul a jednicek najit nejdelsi retezec za sebou jdoucich jednicek.

v10[n_]:=Module[{},SeedRandom[11042011];Table[RandomInteger[],{n}]]
v10[40]

sp1=Split[v10[40]]
(* sel1=Select[sp1,First[#]==1&] *)
cetnosti=Map[Count[#,1]&,sp1]
Max[cetnosti]

sp1=Split[v10[2000000]];
(* sel1=Select[sp1,First[#]==1&] *)
cetnosti=Map[Count[#,1]&,sp1];
Max[cetnosti]

Obecne poznamky k presnosti vypoctu, ruzne  II

Nalezeni bodu zlomu pro linearni funkce

Priklad 1

Clear[g1,g2,g3];
g1[x_]:=2x+3;g2[x_]:=3x+4;g3[x_]:=-x+5;
Plot[{g1[x],g2[x],g3[x]},{x,-3,6}]
gmin[x_]:=Min[g1[x],g2[x],g3[x]]
gminexpand[x_]=PiecewiseExpand[gmin[x]]
Plot[gmin[x],{x,-3,6}]

regions=Map[Reduce[(gminexpand[x]==#),x]&,{g1[x],g2[x],g3[x]}]
regions//FullForm
(regions2=regions//.{Inequality->List,LessEqual->List,GreaterEqual->List}//Flatten)//FullForm
crossings=Select[regions2,NumericQ]
Union[crossings]

Priklad 2 (nahodne linearni funkce)

Clear[garray, x]; (* generuje pole n linearnich funkci *)
garray[n_]:=RandomInteger[{-150,150},n] x+RandomInteger[{-150,50},n]
garray[3]
(* obecneji *)
garrayanonym[n_]:=Function[x,RandomInteger[{-150,150},n] x+RandomInteger[{-150,50},n]]
garrayanonym[3]
garrayanonym[3][z]

SeedRandom[1313131]
g23=garray[23]
Plot[g23,{x,-30,30}]
Plot[g23,{x,-3,3}]

garrayminimum[x_]=Min[g23]
garrayminimumpiecewise[x_]=PiecewiseExpand[garrayminimum[x]]
Plot[garrayminimumpiecewise[x],{x,-10,10},PlotPoints->500]

Najit indexy funkci v g23, ktere se vyskytuji v garrayminimumpiecewise[x]

Position[g23,25-132 x]

garrayminimumpiecewise[x]
garrayminimumpiecewise[x]//FullForm

garrayminimumpiecewise[x][[1]]
garrayminimumpiecewise[x][[1]]//Length

poslednifunkce=garrayminimumpiecewise[x]//Last

prvnichpetfunkci=Map[First,garrayminimumpiecewise[x][[1]][[All]]]

vsechnyfunkce=Join[prvnichpetfunkci,{poslednifunkce}]//Expand

Map[Position[g23,#]&,vsechnyfunkce]

(regionsarray23=Map[Reduce[(garrayminimumpiecewise[x]==#),x]&,g23])//FullForm

(regions2array23=regionsarray23//.{Inequality->List,LessEqual->List,GreaterEqual->List}//Flatten)//FullForm

crossingsarray23=Select[regions2array23,NumericQ]
Union[crossingsarray23]

Uzitecne formy funkci Expand, FindSequenceFunction

PiecewiseExpand

Expand[(((x+y))^3)]

{Expand[Sqrt[x y]],PowerExpand[Sqrt[x y]],Expand[(((a b))^2)]}

{Fibonacci[n],Fibonacci[Range[8]]}

Plot[Fibonacci[n],{n,0,4}]
FunctionExpand[Fibonacci[n]]
FunctionExpand[Fibonacci[n],n∈Integers]
Simplify[FunctionExpand[Fibonacci[n]],n∈Integers]

Plot[n!,{n,-4.5,3}]

FindSequenceFunction[{1,1,2,3,5,8,13,21}]
FindSequenceFunction[{1,1,2,3,5,8,13,21},n]

determinantsHilbert=Table[Det[HilbertMatrix[n]],{n,10}];
FindSequenceFunction[determinantsHilbert,n]

Tuples

(tuples6=Tuples[{1,0},6])//Length
Map[{#,FindSequenceFunction[#,n]}&,tuples6]//MatrixForm

Table[DifferenceRoot[Function[{\[FormalY],\[FormalN]},{-1+\[FormalY][\[FormalN]]+\[FormalY][1+\[FormalN]]==0,\[FormalY][1]==0,\[FormalY][2]==0}]][n]/.n->j,{j,30}]

A jeste pridavek: najit vsechny vektory, kde jsou dve jednicky za sebou

Cases[Tuples[{1,0},5],{___,1,1,___}]
Count[Tuples[{1,0},5],{___,1,1,___}]

Soucty pri podminkach kladenych na indexy, funkce Boole

rr=.;k=.;aa=.;
rr=RandomChoice[Range[10],10]
k=4;
aa=Array[a,{10}]
Sum[a[i]*Boole[rr[[i]]≥k],{i,Length[aa]}]

Boole[rr≥k]
Thread[rr≥k]
Boole[Thread[rr≥k]]
Total[aa*Boole[Thread[rr≥k]]]

{Count[Thread[rr≥k],True],Total[Boole[Thread[rr≥k]]]}

aa[[2]]

Jine pouziti funkce Boole

ineq=y^2-4x^2+4x^4≤-z^2;
RegionPlot3D[ineq,{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->50]
Integrate[Boole[ineq],{x,-Infinity,Infinity},{y,-Infinity,Infinity},{z,-Infinity,Infinity}]

{(∫_(-∞)^∞)(∫_(-∞)^∞)Boole[(x^2)+(y^2)≤1]\[DifferentialD]x\[DifferentialD]y,(∫_(-∞)^∞)(∫_(-∞)^∞)Boole[(x^4)+2(y^2)≤1]\[DifferentialD]x\[DifferentialD]y}

(∫_(-∞)^∞)(∫_(-∞)^∞)Boole[(x^4)+2(y^2)≤1]\[DifferentialD]x\[DifferentialD]y//N

Nahodny vyber s navracenim a bez navraceni

SeedRandom[13];
rwith=RandomChoice[Range[12],10]
Length[rwith]==Length[Union[rwith]]
SeedRandom[13];
rwithout=RandomSample[Range[12],10]
Length[rwithout]==Length[Union[rwithout]]

Zjistit, zda dana mnozina je permutace

Zabudovana funkce PermutationQ

Needs["Combinatorica`"]
{PermutationQ[{3,1,2}],PermutationQ[{3,1,1}]}

perQ[x_?VectorQ]:=Sort[x]==Range[Length[x]]
{perQ[{2,1,3}],perQ[{2,2,1}],perQ[{{1},2,3}]}

per2Q[x_,y_]:=Sort[x]===Union[y]
{per2Q[{a,c,b},{c,b,a}],per2Q[{a,c,b},{c,c,b,a}],per2Q[{a,a,b},{a,b,c}]}

Opakovane pouziti Thread

aa=Array[a,{10}]
Thread[aa==(aa==3)]
Thread[Thread[aa==3]]

Thread[aa==(aa==3)]
Thread[Thread[aa==(aa==3)]]

Zvyrazneni bodu a spojovani v DateListPlot

Konstrukce, ktera je obecne uzitecna.

Show[DateListPlot[{1,1,2,3,5,8,11},{2000,8},PlotStyle->PointSize[0.04]],DateListPlot[{1,1,2,3,5,8,11},{2000,8},Joined->True]]

Primocarejsi konstrukce

DateListPlot[{1,1,2,3,5,8,11},{2000,8},PlotStyle->PointSize[0.04],Joined->True,Mesh->All]

Zakladni financni funkce II

Durace (ekvidistantni casove intervaly, bez overovani vstupu)

Clear[pv];
Options[pv]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
pv[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.Flatten[{opts,Options[pv]}];
Which[interests=="constant",cf.(1/(((1+i))^(Range[n]-1))),
	interests== "variable",cf.(1/(FoldList[#1 * (1+#2)&,1,i])),
	interests=="yieldcurve",First[cf]+Rest[cf].(1/(((1+i))^(Range[n-1])))]
	 ]

Clear[dur];
Options[dur]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
dur[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.Flatten[{opts,Options[pv]}];
Which[interests=="constant",((cf*(Range[n]-1)).(1/(((1+i))^(Range[n]-1))))/pv[cf,i,opts],
	interests== "variable",((cf*(Range[n]-1)).(1/(FoldList[#1 * (1+#2)&,1,i])))/pv[cf,i,opts],
	interests=="yieldcurve",((Rest[cf]*Range[n-1]).(1/(((1+i))^(Range[n-1]))))/pv[cf,i,opts]]
	 ]

cf=.;ii=.;rr=.;i1=.;i2=.;i3=.;
cf={cf0,cf1,cf2,cf3};
ii={i1,i2,i3};
rr={r1,r2,r3};
{dur[1,ii],dur[{},ii],dur[,ii]}

dur[cf,i]//TraditionalForm
dur[cf,ii,interestrates->"variable"]
dur[cf,rr,interestrates->"yieldcurve"]
dur[cf,i,in] (* in nemuze byt uvazovano jako option, vrati tedy dur[... *)

Cviceni: Jaky je profil durace, zahrneme-li do grafu IRR?

Zakladni financni funkce III: Fair values of bonds and durations of bonds

Present value of a coupon bond

Fair = dirty = hruba = spravedliva,
Pure = cista 
c...kuponova sazba!!!
n...pocet obdobi do splatnosti

(∑_(j=1)^n)(v^j)//TraditionalForm

Clear[dirtyvalue,purevalue,v];
dirtyvalue[c_,f_,n_,i_]:=c f (v^(n-Floor[n])) (∑_(j=0)^(Floor[n]))(v^j)+f (v^n)/.v->(1/(1+i))
purevalue[c_,f_,n_,i_]:=f ((v^n)+((c v ((v^n)-1))/(v-1)))/.v->(1/(1+i))

Kuponova obligace s kuponovou sazbou c=0.1, splatnost T=10, nominalni hodnota F=1, hodnotici urokova sazba i=0.15, tj. diskontni obligace.

Znazorneno ve skutecnem (a ne invertovanem) case

p12=Plot[{dirtyvalue[0.1,1,10-n,0.15],purevalue[0.1`,1,10-n,0.15]},{n,0,10},PlotPoints->500,Exclusions->None,AxesLabel->{t,PV}]

Kuponova obligace s kuponovou sazbou c=0.1, splatnost T=10, nominalni hodnota F=1, hodnotici urokova sazba i=0.1.

p3=Plot[dirtyvalue[0.1`,1,10-n,0.1],{n,0,10},Exclusions->None,PlotStyle->{{Dashing[{0.01`,0.01`}],Black}},AxesLabel->{t,PV}]

Kuponova obligace s kuponovou sazbou c=0.1, splatnost T=10, nominalni hodnota F=1, hodnotici urokova sazba i=0.05, tj. premiova obligace.

p45=Plot[{dirtyvalue[0.1,1,10-n,0.05],purevalue[0.1,1,10-n,0.05]},{n,0,10},PlotPoints->500,Exclusions->None,AxesLabel->{t,PV}]

Show[p12,p3,p45,PlotRange->All]

Manipulate[Plot[{dirtyvalue[c,f,tT-t,i],purevalue[c,f,tT-t,i],purevalue[c,f,tT-t,c]},{t,0,tT},PlotPoints->500,Exclusions->None,AxesLabel->{t,PV},PlotRange->{{0,tT},{0,2f}}],
{{tT,10,"T"},5,30,1},
{{c,0.1,"c"},0.01,0.20,0.01},
{{f,1,"F"},1,1000,50},
{{i,0.15,"i"},0.01,0.20,0.01},SaveDefinitions->True]

Duration of a coupon bond

((c f(∑_(j=1)^n)j (v^j)+f n (v^n))/(c f(∑_(j=1)^n) (v^j)+f  (v^n)))/.v->(1/(1+i))//PowerExpand//FullSimplify

Poznamka
Co je mozne, spocitame predem, tim usetrime cas.

{Mean[LogNormalDistribution[μ,σ]],Variance[LogNormalDistribution[μ,σ]]}//Timing
PDF[LogNormalDistribution[μ,σ],x]

{ex=Integrate[x ((E^(-((((-μ+Log[x]))^2)/(2 (σ^2)))))/((Sqrt[2 π]) x σ)),{x,0,∞}],Integrate[(x^(2 ))((E^(-((((-μ+Log[x]))^2)/(2 (σ^2)))))/((Sqrt[2 π]) x σ)),{x,0,∞}]-(ex^2)}//Timing

{ex=Integrate[x ((E^(-((((-μ+Log[x]))^2)/(2 (σ^2)))))/((Sqrt[2 π]) x σ)),{x,0,∞},Assumptions->{σ>0}]//Simplify,Integrate[(x^(2 ))((E^(-((((-μ+Log[x]))^2)/(2 (σ^2)))))/((Sqrt[2 π]) x σ)),{x,0,∞},Assumptions->{σ>0}]-(ex^2)//Simplify}//Timing

Clear[pureduration,dirtyduration,v];
pureduration[c_,f_,n_,i_]:=1+(1/i)-((1+i+c n-i n)/(i+c (-1+(((1+i))^n))))
dirtyduration[c_,f_,n_,i_]:=(c f (∑_(j=0)^(Floor[n])) (n-Floor[n]+j) (v^(n-Floor[n]+j))+f n (v^n))/(c f (v^(n-Floor[n])) (∑_(j=0)^(Floor[n]))(v^j)+f (v^n))/.v->(1/(1+i))

Plot[{dirtyduration[0.1,1,10-n,0.15],pureduration[0.1,1,10-n,0.15]},{n,0,10},PlotPoints->500,Exclusions->None]

Plot[{dirtyduration[0.1,1,10-n,0.05],pureduration[0.1,1,10-n,0.05]},{n,0,10},PlotPoints->500,Exclusions->None]

Financni funkce ve verzi 8

Reference

Financni funkce

Annuity

AnnuityDue

Cashflow

TimeValue

EffectiveInterest

FinancialBond

FinancialDerivative

TradingChart

InteractiveTradingChart

Dalsi financni temata

In Documentation Center find “financial”.

See also

Financial & Economic Data
Mathematica has immediate built-in access to current and historical financial and economic data.
FinancialData — current and historical stock, fund, index, currency, etc. data
“Price”  ▪  “CumulativeReturn”  ▪  “Volatility50Day”  ▪  “MarketCap”  ▪  “PERatio”  ▪  “EBITDA”  ▪  “DividendYield”  ▪  “SICCode”  ▪  ...
	
CountryData — socioeconomic data and time series for all countries
“GDP”  ▪  “ExportValue”  ▪  “FixedInvestment”  ▪  “ElectricityProduction”  ▪  “OilReserves”  ▪  ...
“CurrencyCode”  ▪  “CurrentUnit”  ▪  “ExchangeRate”  ▪  ...
	
DateListPlot — plot a time series
MovingAverage  ▪  ExponentialMovingAverage  ▪  LinearModelFit

See also

Finance
Mathematica has fully integrated support for many of the tools used in classical and modern finance. These capabilities include financial instrument valuation, advanced time value of money computations, and advanced financial charting with a library of technical indicators. Mathematica also provides immediate access to a large array of financial and economic data, and contains financial import and export tools for working with external data.
Time Value Computations
TimeValue — time value of money computations
EffectiveInterest — effective interest rate conversions
Annuity  ▪  AnnuityDue  ▪  Cashflow
	
Financial Instruments
FinancialDerivative — financial derivatives valuation
FinancialBond — financial bond valuation
	
Financial Visualization »
CandlestickChart — candlestick and OHLC chart of open-high-low-close prices
InteractiveTradingChart — interactive trading chart including financial indicators
TradingChart  ▪  DateListPlot  ▪  RenkoChart  ▪  ...
	
Financial Indicators »
"Open"  ▪  "High"  ▪  "Low"  ▪  "Close"  ▪  "Volume"  ▪  "TypicalPrice"  ▪  "TrueRange"  ▪  "SimpleMovingAverage"  ▪  "ExponentialMovingAverage"  ▪  "BollingerBands"  ▪  "KeltnerChannels"  ▪  "RelativeStrengthIndex"  ▪  "FastStochastic"  ▪  "SlowStochastic"  ▪  "FullStochastic"  ▪  "RateOfChange"  ▪  "UltimateOscillator"  ▪  ...
	
Financial and Economic Data »
FinancialData — current and historical stock, fund, index, currency, etc. data
CountryData — socioeconomic data and time series for all countries
DateListPlot  ▪  MovingAverage  ▪  ExponentialMovingAverage  ▪  LinearModelFit
	
Importing & Exporting External Data »
Import — import from files and URLs
Export — export data, graphics, and expressions to files
"Table"  ▪  "CSV"  ▪  "TSV"  ▪  "XLS"  ▪  "ODS"  ▪  "SXC"  ▪  "MDB"  ▪  "DBF"  ▪  "DIF"  ▪  ...

Financni vizualizace

http://www.wolfram.com/mathematica/new-in-8/financial-visualization/

Ukazky financnich vypoctu na serveru Wolfram

http://www.wolfram.com/mathematica/new-in-8/built-in-financial-computations/

http://www.wolfram.com/solutions/industry/financial-risk-management/

http://www.wolfram.com/solutions/industry/financial-risk-management/

Cviceni a ukol

Prohlednete si
http://www.wolfram.com/solutions/FinancialRiskManagement/

Financni funkce verze 8, predbezna informace

TimeValue

Clear[i,v,T,PMT,FV,PV];
{(* FV *) TimeValue[PV,i,T]//TraditionalForm,
(* PV *) TimeValue[FV,i,-T]//TraditionalForm}

Cashflow

Clear[CF];
Print["PV(CF,i) = ",
TimeValue[ccff=Cashflow[ccfftab=Table[(Subscript[CF, i]),{i,0,5}]],i,0]]
Print["FV(CF,i) = ",TimeValue[ccff=Cashflow[ccfftab=Table[(Subscript[CF, i]),{i,0,5}]],i,Length[ccfftab]-1]]

Promenne urokove miry po obdobich

ii=Table[(Subscript[i, j]),{j,1,5}];
Print["PV(CF,i) = ",TimeValue[ccff,ii,0]]
Print["FV(CF,ii) = ",TimeValue[ccff,ii,5]]

Urokove miry ze spotove vynosove krivky

rr=Table[(Subscript[r, j]),{j,1,5}];
rulesrr=Thread[Range[Length[rr]]->rr]
Print["PV(CF,rr) = ",TimeValue[ccff,rulesrr,0]]

Annuity, AnnuityDue

Soucasna hodnota polhutni anuity (Annuity Immediate)

{TimeValue[Annuity[PMT,T],i,0]//TraditionalForm,
(pvAnnuityImmediate=TimeValue[Annuity[PMT,T],i,0]//Simplify),(pvAnnuityImmediate/.{i->((1-v)/v)})//FullSimplify[#,Assumptions->{v>0}]&}//TraditionalForm

Budouci hodnota polhutni anuity

TimeValue[Annuity[PMT,T],i,T]//TraditionalForm

Soucasna hodnota predlhutni anuity (Annuity Due)

{pvAnnuityDue=TimeValue[AnnuityDue[PMT,T],i,0]//FullSimplify,pvAnnuityDue/.{i->((1-v)/v)}//Simplify[#,Assumptions->{v>0}]&}//TraditionalForm

(pvAnnuityDue/pvAnnuityImmediate)//Simplify

Budouci hodnota predlhutni anuity (Annuity Due)

TimeValue[AnnuityDue[PMT,T],i,T]//TraditionalForm

Zakladni financni funkce IV

PV, FV, D, C, IRR, nestejne intervaly mezi platbami

Calendar I

Calendar conventions play a crucial role in time value of money. There is a variety of calendar conventions for  financial instruments. Within Mathematica framework, they are included in Applications Packages like Financial Essentials, UnRisk, or Derivatives Expert. They must be bought separately, however. In Mathematica Add-On Packages the only available calendar is Actual/Actual, i.e., actual number of days in a month and actual number of days in a year. For the illustration of the financial calculations we will make use of this one. Nevertheless, looking on the source code (available in the Package "`Calendar"), you can simply develop the codes for any calendar convention.

Remark

The capabilities of Financial Essentials will be included in version 8. Also the statistical functions in the mentioned version will be comparable with R.

Calendar II

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

Date[]

We may extract just the date:

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

Needs["Calendar`"]

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

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

DaysPlus[today,50]

Present Value, calendar Actual/360

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

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

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

(* Budouci hodnoty  pocitane pomoci funkce pv2 *)dates1={{2004,3,28},{2003,1,13},{2003,1,1}};
Map[DaysBetween[#,today]&,dates1]
dates2=Map[DaysPlus[today,#]&,Map[DaysBetween[#,today]&,dates1]]
cf={1000,-900,-50};
{pv2[cf,dates1,0.1,today],pv2[cf,dates2,0.1,today]}
cfplus={1000,2000,3000};

Kolik vynese koruna ulozena v roce zalozeni University Karlovy na 4% p.a.? Na zacatku roku, v den zalozeni, na konci roku? Jake jsou rozdily? Poznamenejme, ze vklad pri uvedene urokove mire se zdojnasobi asi za 18 let.

koruna1348={pv2[{1},{{1348,1,1}},0.04,today],pv2[{1},{{1348,4,7}},0.04,today],pv2[{1},{{1348,12,31}},0.04,today]}
Most[koruna1348]-Rest[koruna1348]

Konec 11. 4. 2011

Lekce 9 (18. 4. 2011)

FAQ

Body zlomu (reseni Daniel Krycha)

Clear[a,bodyZlomu]
bodyZlomu[fce_]:=Reduce[Limit[D[fce[a],a],a-> x, Direction-> -1]≠  Limit[D[fce[a],a],a-> x, Direction-> 1] ,x]

(* funkce bodyZlomu najde body, kde se meni derivace *)
(* Direction meni jednostranne limity *)

Clear[garray, x]; (* generuje pole n linearnich funkci *)
garray[n_]:=RandomInteger[{-150,150},n] x+RandomInteger[{-150,50},n]
SeedRandom[1313131]
g23=garray[23];
garrayminimum[x_]=Min[g23]

(* zkopirovano z prednasky *)

bodyZlomuPrep=bodyZlomu[garrayminimum]

(* odtud je preparace snazsi *)

Map[bodyZlomuPrep[[#,2]]&,Range[Length[bodyZlomuPrep]]]

Vypreparovani reseni jest jeste snazsi (JH):

bZP2 = {bodyZlomuPrep // ToRules}
x /. bZP2

Bohuzel, ToRules funguje pouze na ==, a ne na nerovnosti.

Casove omezeni vypoctu

TimeConstraint

IRR (Internal Rate of Return)

cf3={-1000,Table[100,{10}],1100}//Flatten;dates3=Map[DaysPlus[today,#]&,Table[365*(j-1),{j,Length[cf3]}]]

Plot[pv2[cf3,dates3,r,today],{r,0.09,0.11}]

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

Cviceni
Vyzkousejte pro skutecne neekvidistantni casove intervaly.

Poznamky

Poznamky k analytickemu zjednoduseni

i1=Integrate[v1=(Sqrt[x]) ArcTan[x],x]
i1=Integrate[v1=(Sqrt[x]) ArcTan[x],x]//Simplify
i1=Integrate[v1=(Sqrt[x]) ArcTan[x],x]//FullSimplify
d1=D[i1,x]
d1==v1
d1===v1
d1//Simplify
(d1//Simplify)==v1

Clear[k,n];
soucet=(∑_(k=1)^n)(1/(k^6))
Limit[soucet,n->∞]
(∑_(k=1)^∞)(1/(k^6))

Poznamka (Neurcite soucty)

sumk/.k->n-sumk
{Sum[(1/(k^6)),k],sumk=Sum[k,k],sumk/.k->n-sumk/.k->0}

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

Series[n!,{n,∞,3}]//PowerExpand//FullSimplify

Series[n!,{n,∞,3}]//PowerExpand//FullSimplify//Normal

((Cos[x])^4)//FullSimplify
TrigReduce[((Cos[x])^4)]

?Trig*

((Cos[x])^4)//TraditionalForm
TrigReduce[%]

Distribucni a kvantilova funkce normalniho rozdeleni

CDF[NormalDistribution[0,1],x]
Quantile[NormalDistribution[0,1],y]//InputForm

Dalsi poznamky ke grafice

Options

Options[Plot]

optionsPlotbutnotListPlot=Complement[Options[Plot],Options[ListPlot]]
optionsListPlotbutnotPlot=Complement[Options[ListPlot],Options[Plot]]

a=.;x=.;
FullForm[x->a]
x->a//First

Map[First,optionsPlotbutnotListPlot]

Vyuziti vlastnich definic

Poznamka 3. Jak zaridit, abychom stale nemuseli psat napr. Joined -> True?

j=Joined->True;
ListPlot[Sin[Range[100]],Joined->False]
ListPlot[Drop[Sin[Range[100]],20],Joined->False]

ListPlot[Sin[Range[100]],j]

(funkcnihodnotysesouradneicemi={Range[100],Sin[Range[100]]}//Transpose)//Short

Vytvorime funkci, ktera zobrazi pouze hodnoty od m-te do n-te hodnoty

takeListPlot[souradnice_,funkce_,m_,n_]:=Module[{pom},
		pom=Take[souradnice,{m,n}];
		pom=Map[{#,funkce[#]}&,pom];
		 ListPlot[pom]]
takeListPlot[Range[100],Sin,50,60]

listPlotJH[x_]:=ListPlot[x,Joined->True]
listPlotJH[{1,2,3}]
ListPlot[{1,2,3}]

jj=Sequence[Joined->True,AxesLabel->{x,y},AxesStyle->Directive[Orange,12]]
pp=ListPlot[Sin[Range[10]],jj,AxesOrigin->{1,0}]

Show: ridi se direktivami prvniho grafu

p=Plot[Sin[x/2],{x,0,10},AxesOrigin->{0,1}]
Show[pp,p]
Show[p,pp]

Prekryvani

{disk1=Graphics[Disk[]],disk2=Graphics[{Orange,Disk[{0,0},1,{Pi/4,3Pi/4}]}]}

Show[disk1,disk2]
Show[disk2,disk1]

Plot[Table[BesselJ[n,x],{n,4}],{x,0,10},Filling->Axis]
Plot[Table[BesselJ[n,x],{n,4}]//Evaluate,{x,0,10},Filling->Axis]

Zaludnost (jedna z mnohych)

U vetsiny, zvlaste grafickych funkci, je siroka moznost voleb. To, co je mozne volit se ziska uzitecnym prikazem (uvadime na konkretnim priklade):

Toto ocividne nedava kruznici:

?AspectRatio

Plot[{-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])},{x,-1,1}]

Options[Plot]

Musime zmenit AspectRatio. Takova ponekud nasilna zmena vsak neni obvykla.

Plot[{-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])},{x,-1,1},AspectRatio->Automatic (* 1 *)]

Zvyrazneni bodu a spojovani v DateListPlot II (uzitecna myslenka i v jinych pripadech)

data1={{{2006,10,1,0,0,0},10},{{2006,10,15,0,0,0},12},{{2006,10,30,0,0,0},15}, {{2006,11,20,0,0,0},20}};
data2={{{2006,10,5,0,0,0},15},{{2006,10,20,0,0,0},8},{{2006,11,10,0,0,0},5}, {{2006,11,15,0,0,0},1}};
DateListPlot[{data1,data2},Joined->{True,False},PlotStyle->PointSize[0.04]]
(* Pouziti na jednu mnozinu dat *)
DateListPlot[{data1,data1},Joined->{True,False},PlotStyle->PointSize[0.04]]

Complement[Options[ListPlot],Options[DateListPlot]]

<<PlotLegends`
DateListPlot[{data1,data1},Joined->{True,False},PlotStyle->PointSize[0.04],Frame->False,AxesLabel->{"Sine","Cosine"}]

Diferencialni rovnice

Diferencialni rovnice symbolicky

eqdif=y'[x]+x y[x]==x Exp[-(x^2)]
sol1=DSolve[eqdif,y[x],x]
funkce=y[x]/.sol1//First

Zabudovani pocatecni podminky pridanim dalsi rovnice

dsol=DSolve[{eqdif,y[0]==0},y[x],x]
DSolve[{eqdif,y[0]==0},y,x]
(y/.%//First)[z]

{yx=y[x]/.dsol[[1]],(y[x]/.dsol)//First}

?yx

y1[x_]:=yx
y1[z]

y2[x_]=yx
y2[z]

yx/.x->z

y3[x_]=y[x]/.dsol[[1]]
y3[z]

Obycejna diferencialni rovnice numericky

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

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

Short[ndsol//FullForm,10]

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

Reseni soustavy dif. rovnice

DSolve[{y'[x]==x+1,z'[x]==y[x]+x,y[0]==0,z[0]==0},{y[x],z[x]},x]
{y[x],z[x]}/.%

Parcialni diferencialni rovnice numericky (sledujeme zaroven potrebu pameti)

Nejdrive se musi spocitat derivace, pak je mozne dosadit konkratni hodnotu argumentu.

{D[u[t,x],0],D[u[t,x],t] /. t->0}

U funkci jedne promenne to neni zpotrebi (ale jenom u zapisu y'[x] !!!

{FullForm[y'[x]],FullForm[y'[0]]}

mmu = MaxMemoryUsed[];ndsolvepartial=NDSolve[{D[u[t,x],t,t]==D[u[t,x],x,x]+Sin[u[t,x]],u[0,x]==Exp[-x^2],(D[u[t,x],t] /. t->0)==0,u[t,-10]==u[t,10]},u,{t,0,10},{x,-10,10}]

D[u[t,x],t,t]==D[u[t,x],x,x]+Sin[u[t,x]]//TraditionalForm

MaxMemoryUsed[]- mmu

Graf muzete libovolne natacet!!!

Plot3D[u[t,x]/.ndsolvepartial,{t,0,10},{x,-10,10}]

Rekurence, diferencni rovnice (recurrence equations)

Pro reseni diferencnich rovnic (nebo jejich soustav) je zabudovana mohutna funkce RSolve. Syntaxe zadani je prakticky stejna jako u obycejnych rovnic, diferencialnich rovnic, &c. Jak jest obecne z kalkulu znamo, spojite pripady jsou vetsinou daleko snadneji zvladnutelne nez pripady diskretni. Mathematica poskytuje neuveritelne rozsahly aparat pro symbolicke i numericke reseni diferencnich rovnic. Ilustraci vsak zacneme na banalnim prikladu vypoctu Fibonacciovych cisel. Pomineme na chvili skutecnost, ze funkce pro vypocet Fibonacciovych cisel jest funkce zabudovana, a z metodickych duvodu uvedeme nekolik algoritmu jejich vypoctu. Take se pritom poucime, jak efektivne vyuzit prikazy pro rekursivni procedury.

Fibonacci A

Spatna definice

Remove[fibS];
fibS[1]=1;
fibS[2]=1;
fibS[n_]:=fibS[n-1]+fibS[n-2]

Table[fibS[j],{j,6}]
Map[fibS,Range[6]]

Toto nefunguje

?$*

fibS[3.5]

fibS[270]

$RecursionLimit=5000;

fibS[270]

fibS[-3]

Trivialni definice:

Remove[fibA];
fibA[1]=1;
fibA[2]=1;
fibA[n_?IntegerQ]:=fibA[n-1]+fibA[n-2]/;Positive[n]

Tady uzivame test, zda n je cele cislo a navic, jestli je kladne.

{fibA[3],fibA[4],fibA[-6],fibA[3.4],fibA[10]}

Casova narocnost je velka (za cviceni zjistete řád a proč)

fibA[29]//Timing

{6.546999999999999`,514229}

(fibA[24]//Timing)[[1]]

ListPlot[Table[(fibA[j]//Timing)[[1]],{j,24}],Joined->True,PlotRange->All]

Fibonacci B (tzv. Dynamicke programovani)

Vyhodou je, kdyz si jiz vypocitana cisla muzeme zapamatovat. To Mathematica umoznuje nasledujici konstrukci:

$RecursionLimit=256;

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

?fibB

Vsimneme si, jak se po vyvolani definice funkci A a B lisi:

{fibA[6],fibB[6]}

?fibA

?fibB

Vsimneme si tez casu potrebneho pro vypocet:

fibB[240]//Timing

?fibB

Ale i tato konstrukce selze pri vyvolani

fibB[2000]

fibB[350]

To lze odstranit zmenou hodnoty systemove promenne, uvidime, ze to bude nezbytne pri hodnoceni americkych opci:

$RecursionLimit

$RecursionLimit=5000;
fibB[500]//Timing

Na tomto miste je take mozne demonstrovat "listovatelnost" funkci, ktera jiz byla demonstrovana pri vyvolavani BuiltIn funkci. Pro uzivatelem definovane funkce pouzijeme nasledujici postup:

Nejprve zkusme aplikovat funkci fibA na list:

fibA[Range[6]]
Attributes[fibA]

Map[fibA,Range[6]]

Po zmene vlastnosti vsak jiz dosahneme kyzeneho vysledku:

SetAttributes[fibA,Listable];
fibA[Range[6]]
Attributes[fibA]

Fibonacci M

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

fibM[3]

Timing[fibM[160000];]

Fibonacci W

fibW[n_]:=Module[{f1=1,f2=1},Do[{f1,f2}={f2,f1+f2},{n-2}];f2]

Map[fibW,Range[6]]

Timing[fibW[160000];]
fibW[160000]//N

Fibonacci C (Reseni diferencnich rovnic)

Vcetne pocatecnich podminek se jedna o soustavu 3 rovnic.

Clear[n];
reseniC=RSolve[{fibC[n]==fibC[n-1]+fibC[n-2],fibC[1]==fibC[2]==1},fibC[n],n]

To je obecne reseni. Neni vsak patrne, jak vypada v konkretnim pripade. Uzijeme

reseniC2=reseniC//FunctionExpand

Uvedomme si, ze cos(n π)=(((-1))^n) pro n cele. Kyzeny tvar dosahneme pouzitim predpokladu

Assuming[n∈Integers,Simplify[reseniC2]]

Konkretni hodnoty dostaneme napr. takto:

reseniC/.n->200

Fibonacci D (zabudovana funkce)

Fibonacci[200]//Timing

Fibonacci F (zabudovana funkce, IPM, p. 218)

More efficient algorithm is based on the formula:

(Subscript[F, 2n])=2(Subscript[F, n-1])(Subscript[F, n])+(Subscript[F, n]^2),	 n≥1

(Subscript[F, 2n+1])=(Subscript[F, n+1]^2)+(Subscript[F, n]^2),                         n≥1

Clear[fibF];
fibF[n_?(IntegerQ[#]&&Positive[#]&)]:=fibF[n]=
	Which[n==1,1,
     	      n==2,1,
		EvenQ[n],2 fibF[n/2-1] fibF[n/2]+((fibF[n/2])^2),
		OddQ[n],((fibF[(n+1)/2])^2)+((fibF[(n-1)/2])^2)
		]
{fibF[2.5],fibF[3],fibF[-1]}

Map[fibF,100 Range[10]]

Timing[fibF[160000];]

Faktorial (zalozene na lokalnich pravidlech)

Substituce jednoducha (ReplaceAll)

Clear[fac];
fac[4]/.{fac[0]->1,fac[n_]:>n*fac[n-1]}

?fac

Nahrazovani se deje pouze pouze jednim pruchodem listu pravidel. Nasledujici "//." (ReplaceRepeated) aplikuje pravidla tak dlouho, az se vyraz nemeni.

fac[4]//.{fac[0]->1,fac[n_]:>n*fac[n-1]}

Trace[fac[4]//.{fac[0]->1,fac[n_]:>n*fac[n-1]}]

?fac

⚠ Varovny priklad

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

Duvod spociva v tom, ze prvni pravidlo se neaplikuje nikdy, a tudiz druhe se aplikuje do nekonecna. Odpomoc spociva v pojistce, dotazu, na typ argumentu. V pripade, ze skutecny parametr nevyhovuje, vypocet se neprovede.

fact[-3]//.{fact[1]->1,fact[n_Integer?Positive]:>n*fact[n-1]}

fact[5]//.{fact[1]->1,fact[n_Integer?Positive]:>n*fact[n-1]}

Diference vyssich radu

Nejdrive cviceni (diference pro vektory dat, srovnej vypocet vynosu z cen)

(********* differences ***********)
ListDifference[data_,n_Integer,s_?(#>5&)]:=Nest[diff[#,s]&,data,n]
diff[data_,s_]:=Drop[data,s]-Drop[data,-s]
diff[((Range[10])^2),2]

data=((Range[10])^3)
ListDifference[data,1,6]

Clear[Δ,f,j];
Δ[f_,1,n_]:=f[n]-f[n-1]
Δ[f_,j_,n_]:=Δ[f,j,n]=Δ[f,j-1,n]-Δ[f,j-1,n-1]
Map[{#,"     ",Δ[f,#,n]}&,Range[4]]//TableForm

Zabudovana funkce (F1)

Differences

Srovnej: Derivace vyssich radu a diference vyssich radu

D[f[x],{x,4}]//FullForm

Derivative[-1][f][x]
Derivative[-1][Function[x,(x^n)]][x]//Evaluate

Clear[x];
g[x_]:=(x^5)
g''[x]
(∂_({x,4})) g[x]
Derivative[-1][g][z]

Dalsi priklad (diferencni rovnice, pokracovani)

Ryzi funkce je nepojmenovana funkce urcena k okamzitemu pouziti

(#1+Sin[#2]&)[x,y]

Function[{x,y},x+Sin[y]][x,y]
Function[{x,y},x+Sin[y]][z,v]

Zde dostaneme reseni pomoci ryzi funkce (Function).

solpure=RSolve[{y[x+2]==ay[x+1]+y[x],y[0]==0,y[1]==1},y,x]

y/.solpure[[1]]

(y/.solpure[[1]])[z]

Chceme-li reseni ve forme funkce x, pouzijeme kostrukci (predposledni argument zmenime na y[x])

solr=RSolve[{y[x+2]==ay[x+1]+y[x],y[0]==0,y[1]==1},y[x],x]

Pripomenme, jak vypreparovat reseni:

y[x]/.solr[[1]]

a definovat tuto funkci jako dale pouzitelnou (zkuste, co takova definice udela bez pouziti Evaluate).

yy[x_]:=Evaluate[(y[x]/.solr[[1]])]
yy[z]

(y[x]/.solr[[1]])/.x->z

MaxMemoryUsed[]

Statisticke funkce I

guide/StatisticalVisualization

TUTORIALS
Basic Statistics
Continuous Distributions
Discrete Distributions
Descriptive Statistics
Convolutions and Correlations
Curve Fitting
Statistical Model Analysis
Importing and Exporting Data
Pseudorandom Numbers
Random Number Generation

Statisticke funkce I odkazy verze 8

Generally

http://www.wolfram.com/mathematica/new-in-8/

New in 8.0: Alphabetical Listing

http://reference.wolfram.com/mathematica/guide/NewIn80AlphabeticalListing.html

More specific

Documentation Center: guide/ProbabilityAndStatistics

Probability and Statistics Solvers and Properties

http://www.wolfram.com/mathematica/new-in-8/probability-and-statistics-solvers-and-properties/

Parametric Probability Distributions

http://www.wolfram.com/mathematica/new-in-8/parametric-probability-distributions/

Nonparametric, Derived, and Formula Distributions

http://www.wolfram.com/mathematica/new-in-8/nonparametric-derived-and-formula-distributions/

Parameter Estimation and Testing

http://www.wolfram.com/mathematica/new-in-8/parameter-estimation-and-testing/

Parametric Statistical Distributions

http://reference.wolfram.com/mathematica/guide/ParametricStatisticalDistributions.html

Statistical Visualization

http://www.wolfram.com/mathematica/new-in-8/statistical-visualization/

Financial Visualization

http://www.wolfram.com/mathematica/new-in-8/financial-visualization/

Built-in Financial Computations

http://www.wolfram.com/mathematica/new-in-8/built-in-financial-computations/

http://www.wolfram.com/solutions/industry/financial-risk-management/

http://www.wolfram.com/solutions/industry/financial-risk-management/

Financial Functions

Annuity

AnnuityDue

Cashflow

TimeValue

EffectiveInterest

FinancialBond

FinancialDerivative

TradingChart

InteractiveTradingChart

Other Financial Topics

In Documentation Center find “financial”.

See also

Financial & Economic Data
Mathematica has immediate built-in access to current and historical financial and economic data.
FinancialData — current and historical stock, fund, index, currency, etc. data
“Price”  ▪  “CumulativeReturn”  ▪  “Volatility50Day”  ▪  “MarketCap”  ▪  “PERatio”  ▪  “EBITDA”  ▪  “DividendYield”  ▪  “SICCode”  ▪  ...
	
CountryData — socioeconomic data and time series for all countries
“GDP”  ▪  “ExportValue”  ▪  “FixedInvestment”  ▪  “ElectricityProduction”  ▪  “OilReserves”  ▪  ...
“CurrencyCode”  ▪  “CurrentUnit”  ▪  “ExchangeRate”  ▪  ...
	
DateListPlot — plot a time series
MovingAverage  ▪  ExponentialMovingAverage  ▪  LinearModelFit

See also

Finance
Mathematica has fully integrated support for many of the tools used in classical and modern finance. These capabilities include financial instrument valuation, advanced time value of money computations, and advanced financial charting with a library of technical indicators. Mathematica also provides immediate access to a large array of financial and economic data, and contains financial import and export tools for working with external data.
Time Value Computations
TimeValue — time value of money computations
EffectiveInterest — effective interest rate conversions
Annuity  ▪  AnnuityDue  ▪  Cashflow
	
Financial Instruments
FinancialDerivative — financial derivatives valuation
FinancialBond — financial bond valuation
	
Financial Visualization »
CandlestickChart — candlestick and OHLC chart of open-high-low-close prices
InteractiveTradingChart — interactive trading chart including financial indicators
TradingChart  ▪  DateListPlot  ▪  RenkoChart  ▪  ...
	
Financial Indicators »
"Open"  ▪  "High"  ▪  "Low"  ▪  "Close"  ▪  "Volume"  ▪  "TypicalPrice"  ▪  "TrueRange"  ▪  "SimpleMovingAverage"  ▪  "ExponentialMovingAverage"  ▪  "BollingerBands"  ▪  "KeltnerChannels"  ▪  "RelativeStrengthIndex"  ▪  "FastStochastic"  ▪  "SlowStochastic"  ▪  "FullStochastic"  ▪  "RateOfChange"  ▪  "UltimateOscillator"  ▪  ...
	
Financial and Economic Data »
FinancialData — current and historical stock, fund, index, currency, etc. data
CountryData — socioeconomic data and time series for all countries
DateListPlot  ▪  MovingAverage  ▪  ExponentialMovingAverage  ▪  LinearModelFit
	
Importing & Exporting External Data »
Import — import from files and URLs
Export — export data, graphics, and expressions to files
"Table"  ▪  "CSV"  ▪  "TSV"  ▪  "XLS"  ▪  "ODS"  ▪  "SXC"  ▪  "MDB"  ▪  "DBF"  ▪  "DIF"  ▪  ...

Statistika: ukazky I

Dokumentace viz heslo "Statistics" (Help -> Documentation Center). Vetsina funkci, ktere byly puvodne pouze v Standard Add Ons Packages jsou nyni funkce zabudovane. Ilustrace na Poissonove rozdeleni. Pro dalsi kliknuti na modre polozky je nutne Menu: "Evaluation   ->   Dynamic Evaluation Enabled"!

?Poi*

PDF (Power Density Function) = hustota 
CDF (Cumulative Distribution Function) = distribucni funkce

λ=2;k=7;
λ=.;k=.;
PDF[distr=PoissonDistribution[λ],k]
distrGauss=NormalDistribution[μ,σ]

{CDF[distr,k],CDF[distrGauss,x],Quantile[distrGauss,q]}

Cviceni: Jaky je vztah mezi Erfc a Φ?

Momenty

{Mean[distr],Variance[distr],Skewness[distr],Kurtosis[distr],Kurtosis[distrGauss]}

ExpectedValue[(x^4),distr,x]
ExpectedValue[((#)^4)&,distr]
ExpectedValue[Function[c,(c^4)],distr]
ExpectedValue[((#)^(4/3))&,distr]

Table[{k,ExpectedValue[(((x-λ))^k),distr,x]},{k,6}]//TableForm

Moment (* pouze vyberovy *)

Obecne podrobnosti: Help -> DiscreteDistributions a ContinuousDistributions

Nahodna cisla z Poissonova rozdeleni s parametrem λ

RandomInteger[distr,{4}]

Neni mozne generovat nahodna cisla s obecnym parametrem. Obecny parametr nutne speciikovat:

RandomInteger[distr/.λ->3,{4}]

To sice funguje, ale lepsi je prehlednejsi

tabPois1=RandomInteger[PoissonDistribution[3],{4}]

Jak videt, tato 4 nahodna cisla jsou jina nez v predchozim pripade. Po kazdem novem vyvolani dostaneme jedinecnou posloupnost. Pomoci nasledujici konstrukce je mozne zajistit reprodukovatelnost nahodnych posloupnosti. Vrele se doporucuje!!! Ten, kdo vystupy sleduje si muze nahodne experimenty overit.

SeedRandom[13];
tabPois2=RandomInteger[PoissonDistribution[3],{4}]
SeedRandom[13];
tabPois3=RandomInteger[PoissonDistribution[3],{5}]
tabPois4=Table[Random[PoissonDistribution[3]],{5}]

Matice nahodnych cisel

SeedRandom[1113];
(tabPois3=RandomInteger[PoissonDistribution[3],{3,5}])//MatrixForm

SeedRandom[1113];
(tabPois4=RandomInteger[PoissonDistribution[3],{5,3}])//MatrixForm

MaxMemoryUsed[]

SeedRandom[1113];tabPois5=RandomInteger[PoissonDistribution[3],{100000}];
MaxMemoryUsed[]

Na programu v blizke budoucnosti

V pristi prednasce se bude hovorit o opcich, binomickem modelu hodnoceni opci a teorii optimalniho portfolia. Pripomente si prislusne partie z prednasek Financni management.

Pomocne odkazy

Vyhledejte aktualni demonstrace tykajici se temat finance, portfolio, opce, kreditni riziko na serveru:
http://demonstrations.wolfram.com/

Obecne odkazy

http://reference.wolfram.com/mathematica/tutorial/StatisticalModelAnalysis.html

http://reference.wolfram.com/mathematica/ref/LinearModelFit.html

http://demonstrations.wolfram.com/ExpectedReturnsOfTheDowIndustrialsBetaModel/

http://demonstrations.wolfram.com/ExpectedReturnsOfTheDowIndustrialsFamaFrenchModel/

http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html

tutorial/UnconstrainedOptimizationOverview

Konec 18. 4. 2011

Lekce 10 (2. 5. 2011, předposlední)

FAQ

Prispevek p. Dufek

Solve ve verzi 8:

Solve ve verzi 8:

Solve[x^3+x^2+x+1==0,x]
Solve[x^3+x^2+x+1==0,x,Reals]
Solve[{x^3+x^2+x+1==0,x>0},x]

Pozn.: Pokud přidáme druhou rovnici jako podmínku, že x>0 automaticky budou vynechána všechna komplexní řešení.

NSolve ve verzi 8:

NSolve[x^3+x^2+x+1==0,x]
NSolve[x^3+x^2+x+1==0,x,Reals]
NSolve[{x^3+x^2+x+1==0,x>0},x]

Reduce

U Reduce se nám nelíbí, že místo {} dává False:

Remove[rovnice];
rovnice[x_]:=x^3+x^2+x+1
Reduce[rovnice[x]==0,x]
Reduce[rovnice[x]==0,x,Reals]
Reduce[{rovnice[x]==0,x>0},x]

Aplikace Solve

rovnice[x_]:=x^3+x^2+x+1;
Solve[{rovnice[x]==0,x>0},x]
reseni=If[Length[%]==0,{},x/.%]
reseni1=If[%%=={},{},x/.%%]
reseni2=x/.%%%

Solve ve verzi 7:

Výstupy z verze 7:

rovnice[x_]:=x^3+x^2+x+1;
Solve[rovnice[x]==0,x]
Solve[rovnice[x]==0,x,Reals]
Solve[{rovnice[x]==0,x>0}x]

NSolve[rovnice[x]==0,x]
NSolve[rovnice[x]==0,x,Reals]
NSolve[{rovnice[x]==0,x>0}x]

Reduce[rovnice[x]==0,x]
Reduce[rovnice[x]==0,x,Reals]
Reduce[{rovnice[x]==0,x>0},x]

Aplikace na výpočet IRR

Aplikace předchozí konstrukce na výpočet irr

ccff={-1100,100,100,100,1000};

pv[ccff_,i_]:=Total[ccff*(1/(1+i))^(Range[Length[ccff]]-1)]

Remove[irr];
irr[ccff_]:=Module[{reseni},reseni=NSolve[{pv[ccff,irr]==0,irr>0},irr];
                                 If[Length[reseni]==0,{},irr/.reseni]
                        ]

irr[ccff]

Plot[pv[ccff,i],{i,-0.5,2}]

Range

U funkce Range je možné zadat minimální a maximální prvek seznamu.

Range[11]-1

Range[0,10,1]

Otázkou zůstává co je přehlednější.

Pokud maximální prvek není dělitelný krokem, vezme se jako největší prvek nejbližší nižší násobek:

Range[0,7,2]

Mean = průměr, Variance = výběrový rozptyl

Průměr fce Mean:

prumer[x_]:=Apply[Plus,x]/Length[x]/;VectorQ[x]&&Length[x]>0
{prumer[{a,b,c}],prumer[{{a},b,c}],prumer[a,b,c],prumer[2],prumer[{2}]}
Mean[{a,b,c}]

rozptyl[x]=prumer[(x-prumer[x])^2] neni shodný se zabudovanou funkcí Variance:

rozptyl[x_?VectorQ]:=Mean[(((x-Mean[x]))^2)];

SeedRandom[13];
vyber=RandomInteger[{-5,5},10];
rozptyl[vyber]
Variance[vyber]

Funkce Variance (pro náhodný výběr) je definovaná následovně: Variance[X]:=(1/(Length[X]-1))(∑_(i=1)^(Length[X]))((((Subscript[X, i])-Mean[X]))^2).

Remove[vyberovyRozptyl];

vyberovyRozptyl[X_]:=(1/(Length[X]-1))*Total[(X-Mean[X])^2]/;VectorQ[X]&&Length[X]>1

vyberovyRozptyl[vyber]

Test odolnosti vůči špatnému zadání:

{vyberovyRozptyl[{1,{1}}],vyberovyRozptyl[{1}]}

Vypnuti chybovych hlasek

Solve[(x^2)=0,x]
Off[Set::write]
Solve[(x^2)=0,x]

Off[Solve::naqs]
Solve[(x^2)=0,x]
On[Set::write,Solve::naqs]
Solve[(x^2)=0,x]

Mean a ExpectedValue

ExpectedValue

Statisticke funkce II

Statistika: ukazky II

Ilustracemi tykajici se pameti

MaxMemoryUsed[]

{ByteCount[0],ByteCount[0.]}
Table[ByteCount[(4^j)],{j,0,200,10}]

Compatibility/tutorial/Statistics/DescriptiveStatistics

SeedRandom[1113];tabPois5=RandomInteger[PoissonDistribution[3],{100000}];
MaxMemoryUsed[]

Show[Histogram[tabPois5,{1},"ProbabilityDensity",PlotRange->All],Plot[PDF[PoissonDistribution[3],x]//Evaluate,{x,0,15},PlotStyle->Thick,PlotRange->All]]

Clear[LocationReport];LocationReport[data_]:={Mean->Mean[data],HarmonicMean->HarmonicMean[data],Median->Median[data]}

LocationReport[tabPois5]
%//N

Cviceni: Vytvorte DispersionReport (Min, Max, Rozptyl, prumer absolutnich odchylek, dolni a horni kvartil, mezikvartilove rozpeti, rozpeti)

DispersionReport[tabPois5]//N

Dalsi charakteristiky

{Min[tabPois5],m5=Max[tabPois5]}
Print["P{Po[3]>15} = ",1-CDF[PoissonDistribution[3.],15]]

fr5=Tally[tabPois5] (* zuctovani, tabulka *)

fr5//Sort

Sort[fr5,#1[[2]]<#2[[2]]&]

Map[Reverse,Sort[Map[Reverse,fr5]]]

BarChart[fr5]

Histogram[tabPois5,PlotRange->All]

Veta (Metoda inverzni transformace, bude o ni jeste pojednano v metodach MC) :

Nechť náhodná veličina X má distribuční funkci F a kvantilovou funkci (F^(-1))(y)=inf{x|F(x)≥y}, 0<y<1. Nechť U∼R(0,1). Potom (F^(-1))(U) má rozdělení s ditribuční funkcí F, tj., má rozdělení jako X. Platí pro libovolné rozdělení, spojité, diskrétní a libovolné jiné!!!

Priklad (vlastni distribucni funkce a generovani nahodnych cisel)

dFuniform[x_]:=Piecewise[{{0,x≤0},{x,0<x<1},{1,x≥1}}]
Plot[dFuniform[x],{x,-1,2},AxesOrigin->{0,-1}]

Solve[dFuniform[x]==y,x]

Cviceni: Upravte do funkcni podoby.

Jiny priklad (exponencialni rozdeleni se stredni hodnotou 1)

dFexp[x_]:=1-(E^(-x))
Plot[dFexp[x],{x,0,10}]

sol1=Solve[dFexp[x]==y,x]
invF[y_]=First[x/.sol1]
invF[z]

nahexp=invF[RandomReal[UniformDistribution[{0,1}],{100000}]];
{Mean[nahexp],Variance[nahexp]}

Reduce[dFexp[x]==y,x,Reals]

Vyrovnani, regrese

fr5={{2,22494},{7,2152},{6,4991},{3,22455},{1,14834},{4,16850},{5,10019},{0,5065},{8,781},{9,249},{10,78},{12,5},{11,25},{15,1},{13,1}};

ListPlot[fr5]

x^(Range[6]-1)
fit5=Fit[fr5,x^(Range[6]-1),x]

Dalsi varovani pred vyrovnavanim (fitting) polynomy:

p1=ListPlot[fr5,PlotStyle->PointSize[0.02]];
p2=Plot[fit5,{x,0,25},PlotRange->All];
Show[p1,p2]

FindMinimum[fit5,{x,9}]

FindMinimum[fit5,{x,9.52618271607493`}]

regr5=LinearModelFit[fr5,x^(Range[6]-1),x]

regr5//Normal

Help (F1)

regr5[{"ParameterTable","RSquared"}]

regr5["ANOVATable"]

RSquared is the ratio of the model sum of squares to the total sum of squares.

AdjustedRSquared is equal to 1-(n-1)/(n-p)(1-RSquared)
where n=Length[data]andp is the number of parameters in the model.

Goodness-of-fit measures are used to assess how well a model fits or to compare models. The coefficient of determination "RSquared" is the ratio of the model sum of squares to the total sum of squares. "AdjustedRSquared" penalizes for the number of parameters in the model and is given by . 
"AIC"  (Akaike Information Criterion) and "BIC" (Bayesian Information Criterion) are likelihood-based goodness-of-fit measures. Both are equal to -2 times the log-likelihood for the model plus k p where p is the number of parameters to be estimated including the estimated variance. For "AIC" k is 2, and for "BIC" k is log (n).

{regr5["RSquared"],regr5["AdjustedRSquared"],regr5["AIC"]}

Letmo prostudujte (pro ucely pozdejsiho aktivniho uziti):

Rozdeleni matematicke statistiky

Nejdříve vytvoříme obecnou proceduru pro ilustraci všech tří základních rozděleních. Výstupem je list o dvou položkách, první je graf hustoty a druhá je tiskne základní charakteristiky rozdělení.

Clear[ilustracerozdeleni];
Options[ilustracerozdeleni]={graf->False};
ilustracerozdeleni[rozdeleni_,opts___]:=Module[{q01=Quantile[rozdeleni,0.001],q99=Quantile[rozdeleni,0.999],mean=Mean[rozdeleni],std=StandardDeviation[rozdeleni],igraf},
igraf=graf/.Flatten[{opts,Options[ilustracerozdeleni]}];
Which[igraf,
Plot[{PDF[NormalDistribution[mean,std],x],PDF[rozdeleni,x]}//Evaluate,{x,q01,q99},AxesOrigin->{0,0}],
True,Print["Mean = ",mean,"  Variance = ",Variance[rozdeleni],"\nSkewness = ",Skewness[rozdeleni]//PowerExpand//Simplify,"  Kurtosis = ",Kurtosis[rozdeleni]//Simplify]]
]

(χ^2)-kvadrát

ilustracerozdeleni[ChiSquareDistribution[n]]

Manipulate[ilustracerozdeleni[ChiSquareDistribution[n],graf->True],{n,2,300}]

{PDF[ChiSquareDistribution[2],x],PDF[ExponentialDistribution[1/2],x]}

Studentovo t-rozdělení o n stupních volnosti

ilustracerozdeleni[StudentTDistribution[n]]

Manipulate[ilustracerozdeleni[StudentTDistribution[n],graf->True],{n,3,100}]

{PDF[StudentTDistribution[1],x],PDF[CauchyDistribution[0,1],x]}

Fisher-Snedecorovo F-rozdělení

ilustracerozdeleni[FRatioDistribution[n,m]]

Manipulate[ilustracerozdeleni[FRatioDistribution[n,m],graf->True],{n,1,100},{m,5,100}]

Vybrane statisticke metody Mathematica v. 7: Vsechny nove rysy (zastarale)

(Tento odkaz se spatne hleda.Ostatni odkazy nejsou uplne.)

http://reference.wolfram.com/mathematica/guide/SummaryOfNewFeaturesIn70.html

Statisticke funkce

Toto je vypis z vyse uvedene stranky:

Statistics & Data Analysis
FindCurvePath for reconstruction of sampled curves
Median, Quartiles, QuartileDeviation, ... added support for built - in distributions
LinearModelFit, NonlinearModelFit to construct regression models from data
LogitModelFit, ProbitModelFit, GeneralizedLinearModelFit added
Interpolation now allows immediate calculation of interpolated values
InverseGammaDistribution, InverseChiSquareDistribution added
LevyDistribution added
StudentTDistribution extended to have multiple parameters
ExpectedValue extended to support sampled data
Standardize to adjust constant shift and variance of the data
FittedModel for symbolic representation of fitted models
DesignMatrix to construct the design matrix for a linear model
Entropy function to measure of randomness in data added
ContraharmonicMean added

FindCurvePath

({{, FindCurvePath[{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},…}] 
gives an ordering of the {(Subscript[x, i]),(Subscript[y, i])} that corresponds to one or more smooth curves.}})

list=Table[{Cos[t],Sin[t]},{t,RandomReal[{0,2π},(* 1500 *)15]}];
ListLinePlot[list,AspectRatio->Automatic]
curve=FindCurvePath[list]
ListLinePlot[list[[curve[[1]]]],AspectRatio->Automatic]
ListCurvePathPlot[list]

Nejkratší cesta

({{, FindShortestTour[{(Subscript[e, 1]),(Subscript[e, 2]),…}] 
attempts to find an ordering of the (Subscript[e, i]) that minimizes the total distance on a tour that visits all the (Subscript[e, i]) once.}})

Find the length and ordering of a shortest tour through six points in the plane:

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

Median, Quartiles, QuartileDeviation, ...

Nyni funguje i pro zabudovana rozdeleni

Clear[λ];
{Median[ExponentialDistribution[λ]],Quartiles[ExponentialDistribution[λ]]}

Quartiles[distrGauss=NormalDistribution[μ,σ]]//N

NonlinearModelFit

data={{0,1},{1,0},{3,2},{5,4},{6,4},{7,5}};
nlm=NonlinearModelFit[data,Log[a+b x^2],{a,b},x]
Normal[nlm]

Show[ListPlot[data],Plot[nlm[x],{x,0,7}],Frame->True]

nlm["BestFitParameters"]

ProbitModelFit

Pouziva se pro odhad pravdepodobnosti nesplaceni (default risk). V probitove analyze se predpoklada, ze pravdepodobnost P jefunkci nezavisle promennych (Subscript[x, 1]),...(Subscript[x, k]):

P=Φ((Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, 1])+...+(Subscript[β, k])(Subscript[x, k])),

kde (Subscript[β, 0]),...,(Subscript[β, k]) jsou nezname parametry a Φ je distribucni funkce N(0,1). Pozorovany jsou hodnoty Y nabyvajici hodnot 0,1 a jim odpovidajici hodnoty x. V nasledujicim priklade je k=1 a data jsou ve standardnim tvaru {x,y}.

data={{1,0},{1,1},{1,1},{2,0},{2,0},{2,1},{2,1},{2,1},{2,1},{3,1},{3,1},{3,1},{3,1},{3,1}};
probit=ProbitModelFit[data,x,x]
Normal[probit]
Show[ListPlot[data],Plot[probit[x],{x,1,3}],PlotRange->All,AxesOrigin->{.9,-.1}]
pr1=probit["PredictedResponse"]

Podobna je logitova analyza. Misto Φ se vsak pouziva distribucni funkce logistickeho rozdeleni.

LogitModelFit

Pouziva se take pro odhad pravdepodobnosti nesplaceni (default risk)

{PDF[LogisticDistribution[μ,β],x],CDF[LogisticDistribution[μ,β],x]}

logit=LogitModelFit[data,x,x]
Normal[logit]
Show[ListPlot[data],Plot[logit[x],{x,1,3}],PlotRange->All,AxesOrigin->{.9,-.1}]

lo1=logit["PredictedResponse"]

Porovnani predikce pomoci probitu a logitu:

Max[Abs[pr1-lo1]]

LevyDistribution

Rozdeleni s tezkym pravym chvostem (heavy right tail). Pozor, nosic rozdeleni je (μ,∞)!!!

PDF[LevyDistribution[μ,σ],x]
Mean[LevyDistribution[μ,σ]]
Plot[PDF[LevyDistribution[0,1],x],{x,0,5}]
{Median[LevyDistribution[μ,σ]],Median[LevyDistribution[μ,σ]]//N}

Integrate[Exp[-a (x^2)],{x,-∞,∞}]
Integrate[Exp[-a (x^2)],{x,-∞,∞},Assumptions->a>0]
Assuming[a>0,Integrate[Exp[-a (x^2)],{x,-∞,∞}]]

Assuming[σ>0&&μ∈Reals,(∫_μ^∞)PDF[LevyDistribution[μ,σ],x]\[DifferentialD]x]

SeedRandom[13];rannorm=RandomReal[NormalDistribution[0,1],{100000}];
SeedRandom[13];ranlev=RandomReal[LevyDistribution[0,1],{100000}];
SeedRandom[13];rancau=RandomReal[CauchyDistribution[0,1],{100000}];

Map[Max,{rannorm,ranlev,rancau}]

Map[Min,{rannorm,ranlev,rancau}]

Median[LevyDistribution[0.,1]]

Map[Mean,{rannorm,ranlev,rancau}]
Map[Median,{rannorm,ranlev,rancau}]

Select[rancau,Abs[#]>5000&]

Mean[CauchyDistribution[0,1]]

StudentTDistribution

Nyni take s parametry polohy (location) μ a meritka (scale) σ.

PDF[StudentTDistribution[μ,σ,ν],x]
{Mean[StudentTDistribution[μ,σ,ν]],Variance[StudentTDistribution[μ,σ,ν]]}

Standardize

Prekvapujici moznosti

Bezne

strs1=Standardize[rs1=RandomReal[{0,10},20]];
{Mean[strs1],Variance[strs1]}

Pouze nulovy prumer:

strs2=Standardize[rs1,Mean,1&];
{Mean[strs2],Variance[strs2]}

strs3=Standardize[rs1,0&,StandardDeviation];
{Mean[strs3],Variance[strs3]}

Standardize[rs1,Median,InterquartileRange]

Standardize[list,(Subscript[f, 1]),(Subscript[f, 2])] is effectively (list-(Subscript[f, 1])[list])/(Subscript[f, 2])[list].

DesignMatrix (Matice planu)

Matice planu pro kvadratickou regresi s absolutnim clenem a bez nej.

data=Table[{i,i^(3/2)+RandomReal[]},{i,5}]
DesignMatrix[data,{x,x^2},x]//MatrixForm
DesignMatrix[data,{x,x^2},x,IncludeConstantBasis->False]//MatrixForm

ContraharmonicMean

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

(* the order p Lehmer contraharmonic mean *)ContraharmonicMean[{a,b,c,d},p]

Statistical Packages

Directory[]
SetDirectory["C:\\Program Files\\Wolfram Research\\Mathematica\\8.0\\AddOns\\Packages"];
Directory[]
FileNames[]
SetDirectory[NotebookDirectory[]]

Needs["MultivariateStatistics`"]
Names["MultivariateStatistics`*"]

MaxMemoryUsed[]
Needs["NumericalDifferentialEquationAnalysis`"]
Needs["PrimalityProving`"]
MaxMemoryUsed[]

Metody Monte Carlo I

Pokus: dokaze clovek generovat nahodne cislice?

vera=;

michal=;

michal5=IntegerDigits[michal]
Length[michal5]

(* Inovace *)
vera5=IntegerDigits[vera]

Clear[datareport];
datareport[data_,noclass___]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=If[noclass===Null,10,noclass];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data,k]]

datareport[michal5//N]
(* Vsimnete si rozdilu mezi //N a bez numericke hodnoty! *)

𝒟unif=UniformDistribution[{0,9}]
{Mean[𝒟unif],Variance[𝒟unif],StandardDeviation[𝒟unif]}//N

Pokus: dokaze clovek generovat nahodne cislice? (Vyhodnoceni dat Vera)

Puvodne zvoleny postup. Ten sice funguje, ale je zbytecne slozity. Ten vyse uvedeny potrebuje jedinou instrukci.

vera1=ToString[vera];
vera2=Characters[vera1];
vera3=ToExpression[vera2];
(* Pozor, nutne iniciovat proceduru datareport! *)
datareport[vera3]

Test hypotezy, zda data maji diskretni rovnomerne rozdeleni 
DiscreteUniformDistribution[{0, 9}]. Pouziji vlastni balik obsahujici test dobre shody. Za cviceni si vytvorte analogickou proceduru.

cetnosti1=Tally[michal5]
cetnosti2=Sort[cetnosti1]
cetnosti3=Map[Last,cetnosti2]
pocetpozorovani=Length[michal5]
pocetruznychcisel=Union[Map[First,cetnosti1]]//Length
chikvadrat=Total[((((cetnosti3-0.1 pocetpozorovani))^2)/(0.1 pocetpozorovani))]

Print["Hodnota (χ^2) na ",pocetruznychcisel-1," stupnich volnosti = ",chikvadrat,"\nOdpovidajici p - hodnota = ",1-CDF[ChiSquareDistribution[pocetruznychcisel-1],chikvadrat]]

cetnosti2=Tally[michal5]
pocetruznychcisel=Length[cetnosti2]
pocetpozorovani=Total[cetnosti2]
((((cetnosti2-pocetpozorovani/pocetruznychcisel))^2)/(pocetpozorovani/pocetruznychcisel))
chi2=Total[((((cetnosti2-pocetpozorovani/pocetruznychcisel))^2)/(pocetpozorovani/pocetruznychcisel))]//N

Metody Monte Carlo II

Nahodna cisla

tutorial/RandomNumberGeneration

Generator typu (Subscript[X, n+1])=a (Subscript[X, n]) mod M se nazyva multiplikativni,
(Subscript[X, n+1])=a (Subscript[X, n])+b  mod M smisena kongruencni metoda. Pro ilustraci uzijeme zabudovanou funkci Mod:

Mod[5,3]

Priklad 1 (Perioda) Generovani posloupnosti 17 nahodnych cisel     (Subscript[X, n+1])=5(Subscript[X, n])+1 mod 16. Perioda je maximalni rovna 16

gen16=NestList[Mod[5 # +1,16]&,0,16]

Priklad 2    (Vetsi problem, rozdeleni dvojic. Stejny generator (Subscript[X, n+1])=5(Subscript[X, n])+1 mod 16.)   Hypoteticky kazda dvojice ma pst  (1/(16^2))=(1/256). Vzhledem k rekurenci vsak neni mozne generovat vice nez 16-1 dvojic; gnerovana cisla lezi na dvou systemech primek:

gen16shift=RotateLeft[gen16]

ListPlot[Transpose[{gen16,gen16shift}],PlotStyle-> PointSize[0.02]]

Priklad 3 (Generator kapesniho kalkulatoru TI-59). Nasledujici funkce generuje nahodna cisla smisenou kongruencni metodou delky n+1. Perioda by mela byt 199016.

Clear[gen59];
gen59[n_]:=NestList[Mod[24298 # +9991,199017]&,0,n]

gen59[10]

Po sobe jdouci dvojice

gen592[n_]:=Module[{shift,g},
shift=RotateLeft[g=gen59[n]];
Transpose[{g,shift}]
]

gen592[50000]//ListPlot

ListPlot[gen59compl=Select[gen592[199016],(#[[1]]<1000)&]]

Generovani nahodnych cisel z velmi zvlastnich rozdeleni (prevzato z MathGroup Archive: December 1998 [00243], Paul Abbott)

Pripominka : Metoda inversni transformace

Clear[pdf];
pdf[A_,s_][r_]:=(r/s^2*BesselI[0,(r*A)/s^2])/E^((A^2+r^2)/(2*s^2));

After some 93 (version 7) 25 (version 8) seconds we get that pdf is really a PDF.

Integrate[pdf[A,s][r],{r,0,∞},Assumptions->Re[(s^2)]>0&&(Re[(A/(s^2))]>0||(Re[(A/(s^2))]==0&&Im[(A/(s^2))]>0))]//Simplify

{"Mean =  ",e11=Integrate[x pdf[1,1][x],{x,0,∞}]//N,"   Variance =   ",Integrate[(x^2) pdf[1,1][x],{x,0,∞}]-(e11^2)//N}

Plot[pdf[1,1][x],{x,0,5}]

Use NDSolve (and dynamic programming) instead of NIntegrate tp determine the CDF:

Clear[cdf,r];
cdf[A_,s_]:=cdf[A,s]=Module[{y},
First[y/.NDSolve[{y'[r]==pdf[A,s][r],y[0]==0},y,{r,0,10}]]]

Porovnejte s vypoctem CDF pomoci integrace.

From the plot,

cdfplot=Plot[cdf[1,1][x],{x,0,10},PlotRange->All]

Short[InputForm[cdfplot],10]

extract the plotpoints (which takes advantage of the adaptive code built into Plot):

Clear[x];
pts=Cases[cdfplot,Line[{x__}]->x,Infinity];

ListPlot[pts]

and then numerically compute the inverse function using interpolation:

inverseCDF[1,1]=Interpolation[Map[Reverse,pts]];

Plot[inverseCDF[1,1][x],{x,0.0000001,0.9999999999}]

InverseCDF

Table[inverseCDF[1,1][Random[]],{10}]

Testovani nahodnosti (znamenka diferenci, extremalni body, korelogram), specielne generatoru nahodnych cisel

Integrate[Boole[x1<x2<x3],{x1,0,1},{x2,0,1},{x3,0,1}]

Generatory nahodnych cisel se testuji baterii statistickych testu, z nichz nektere uvadime. Podrobnosti HURT Simulacni metody, jeste podrobneji KNUTH The Art of Computer Programming.

1 Test dobre shody viz vyse a nize

2 Test podle znamenek diferenci
Vhodny proti alternative trendu. Pro spojite iid X,Y je P(X≥Y)=(1/2). Pro diskretni priblizne take. Vezmeme-li  X,Y z rozdeleni UniformDistribution[{0, m-1}], je

P(X≥Y)=(∑_(i=0)^(m-1))(∑_(j=0, i≥j)^(m-1))(1/(m^2))=(1/(m^2))(∑_(i=0)^(m-1))(∑_(j=0)^i)1=

spocitano jiz v systemu Mathematica

pxy=(1/(m^2))(∑_(i=0)^(m-1))(∑_(j=0)^i)1//Simplify
eC=(n-1)pxy

Necht testovana posloupnost je (Subscript[X, 1]),...,(Subscript[X, n]). Polozme (Subscript[Y, i])=I{(Subscript[X, i])≤(Subscript[X, i+1])}, i=1,...,n-1,

C=(∑_(i=1)^(n-1))(Subscript[Y, i]). Spocteme var C.

varY1=pxy-(pxy^2)//Simplify

Plati cov((Subscript[Y, 1]),(Subscript[Y, 2]))=E (Subscript[Y, 1])(Subscript[Y, 2]) - (((E (Subscript[Y, 1])))^2)  a E (Subscript[Y, 1])(Subscript[Y, 2])= P((Subscript[X, 1])≤(Subscript[X, 2])≤(Subscript[X, 3]))

pxyz=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=j)^(m-1))1
covY1Y2=pxyz-(pxy^2)//Simplify

varC=(n-1)varY1+2(n-2)covY1Y2//Simplify

Nyni staci aplikovat CLV pro konecne zavisle nahodne veliciny.

rozdily=Rest[michal5]-Most[michal5];

Count[{1,2,3,4,0,0,-6,-7,I},_?NonNegative]

cc=Count[rozdily,_?NonNegative]

kriter=((cc-eC)/(Sqrt[varC]))/.{n->Length[michal5],m->9.}

P(|N(0,1)|>k)=1-Φ(k)+Φ(-k)=2(1-Φ(k)). Ale P-hodnotu muzeme jednodusse spocitat jako

1- Integrate[PDF[NormalDistribution[0,1],x],{x,-k,k}]/.k->kriter

Partition[{1,2,3,4,0,0,-6},3,1]
Select[Partition[{1,2,3,4,0,0,-6},3,1],#[[1]]≤#[[2]]≥#[[3]]&]

3 Test extremalnich hodnot (Lokalni extremy)
Pro spojite iid X,Y, Z je P(X≤ Y≥Z)=(1/3). (Staci vzit vsechny permutace trojic (tech je 6) a zjistit, kolik z nich ma lokalni maximum.) Podobne pst lokalniho minima je (1/3). Pst lokalniho extremu v trojici je tedy (2/3). Vezmeme-li  X,Y, Z z rozdeleni UniformDistribution[{0, m-1}] a dostaneme

P(X≤Y≥Z)=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=0,j≥i)^(m-1))(∑_(l=0,l≤j)^(m-1))1=

pstmax=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=0)^j)1

pstextr=2pstmax (* maxima + minima *)

V nahodne posloupnosti delky n je tedy prumerne

eP=(n-2)(2/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=0)^j)1

lokalnich extremu. Rozptyl P je priblizne

varP=((16n-29)/90);

Needs["Combinatorica`"]

Select[Permutations[Range[3]],(#[[1]]≤#[[2]]≥#[[3]]&)]

Select[{{1,2,3},{1,2,1},{5,6,7},{3,4,1}},(#[[1]]≤#[[2]]≥#[[3]]&)]
Length[%]

Count[{{1,2,3},{1,2,1},{5,6,7},{3,4,1},{7,3,8}},_?(#[[1]]≤#[[2]]≥#[[3]]&)]

Count[{{1,2,3},{1,2,1},{5,6,7},{3,4,1},{7,3,8}},_?(Or[#[[1]]≤#[[2]]≥#[[3]],#[[1]]≥ #[[2]]≤ #[[3]]]&)]

Clear[extremalnibody];extremalnibody[data_,m_:10]:=Module[{par,pp,n,eP,varP},
	n=Length[data];
	par=Partition[data,3,1];			pp=Count[par,_?(Or[#[[1]]≤#[[2]]≥#[[3]],#[[1]]≥ #[[2]]≤ #[[3]]]&)];
	eP=(((1+m) (1+2 m) (-2+n))/(3 (m^2)));
varP=((16n-29)/90);
Print["n = ",n,"  P = ",pp,"  P-hodnota = ",
2 CDF[NormalDistribution[0.,1],-Abs[((pp-eP)/(Sqrt[varP]))]]]
]

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

extremalnibody[michal5]

Konec 2. 5. 2011

Lekce 11 (9. 5. 2011, poslední)

FAQ

Definice vlastnich rozdeleni

Minimum z nezavislych Poissonovych rozdeleni a momenty

t2=TransformedDistribution[Min[u,v],{u\[Distributed]PoissonDistribution[(Subscript[μ, 1])],v\[Distributed]PoissonDistribution[(Subscript[μ, 2])]}]

60 seconds

fk=PDF[t2,k]

fk//FunctionExpand//Simplify

Nespocita:

Mean[t2]

Correlation[michal5,RotateLeft[michal5]]//N

Funkce nedefinovane pro nektere argumenty

(0^0)

FullForm[HoldForm[(0^0)]]

Unprotect[Power]
Power[0,0]=1
Protect[Power]

(0^0)

Volba vstupu (napr. Wolfram Alpha, IRR)

WolframAlpha["IRR",IncludePods->"Input",AppearanceElements->{"Pods"},InputAssumptions->{"*C.IRR-_*Formula-"}]

Hodnocení derivátů typu Cap (strop) a Floor  (dno)

Důvod, proč ve funkci FinancialDerivative není volba Cap resp. Floor (přesněji Caplet a Floorlet) vidím v tom, že tyto pojmy nejsou jednoznačně definovány a možností je celá řada. Obvykle jsou však kombinací elementárních finančních instrumentů (obligací, opcí, apod.)

Uvedeme ilustraci Cap na evropskou CALL opci s dobou vypršení T a realizační cenou K. Chrání investora v krátké pozici před neomezenou ztrátou. V případě, že cena (Subscript[S, T]) podkladového aktiva je v okamžiku vypršení vyšší než B,  vyplácí se B.

vyplataCAPnaCALL[St_,K_,B_]:=Min[Max[St-K,0],B]
Plot[vyplataCAPnaCALL[st,1,2],{st,0,5},PlotStyle->Thick]

Dlouhá pozice v tomto finančním instrumentu je ekvivalentní držení portfolia dvou CALL opcí, obě se stejnou dobou vypršení, jedna v dlouhé pozici s realizační cenou K, druhá v krátké pozici s realizační cenou K+B.

vyplataportfolia[St_,K_,B_]:=Max[St-K,0]-Max[St-K-B,0]
Plot[vyplataportfolia[st,1,2],{st,0,5},PlotStyle->Thick]

FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50, "Dividend"->0.05}]

Cena tohoto instrumentu podle Black-Scholesova vzorce je tedy (musíme přidat bezrizikovou úrokovou míru r,0, volatitlitu σ a dobu do vypršení T,t.

Clear[cenaCAPnaCALL];
cenaCAPnaCALL[St_,K_,B_,r0_,σ_,Tt_]:=FinancialDerivative[{"European","Call"}, {"StrikePrice"-> K, "Expiration"->Tt},  {"InterestRate"-> r0, "Volatility" -> σ, "CurrentPrice"-> St}]-FinancialDerivative[{"European","Call"}, {"StrikePrice"-> K+B, "Expiration"->Tt},  {"InterestRate"-> r0, "Volatility" -> σ, "CurrentPrice"-> St}]

Plot[cenaCAPnaCALL[st,1,2,0.05,0.4,1],{st,0,5},PlotRange->{0,2},PlotStyle->Thick]

Vyzkoušejte si manipulovat s některými parametry:

Manipulate[Plot[cenaCAPnaCALL[st,1,2,0.05,σ,Tt],{st,0,5},PlotRange->{0,2},PlotStyle->Thick],{{σ,0.3},0.1,0.5},{{Tt,1},0.01,2},SaveDefinitions->True]

Metody Monte Carlo III (Metody redukce rozptylu)

Korelogram

Korelogram (Data michal5 v Metody Monte Carlo I)

ListPlot[Table[{k,Correlation[michal5,RotateLeft[michal5,k]]},{k,0,15}]//N,Joined->True,PlotRange->All]

Testovaci funkce

Podrobneji viz napr. skripta Hurt, J.: Simulacni metody

f1[u_]:=u
f2[u_]:=(E^u)-1
du=UniformDistribution[{0,1}];

{intf1=(∫_0^1)f1[u]\[DifferentialD]u,intf2=(∫_0^1)f2[u]\[DifferentialD]u}//N

Plot[{f1[u],f2[u]},{u,0,1}]

Obycejna metoda Monte Carlo

Needs["HypothesisTesting`"]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},100000];]
(f1observ=f1[randomuniform])//Length
Mean[f1observ]
Variance[f1observ]
MeanCI[f1observ]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},100000];]
(f1observ=f1[randomuniform])//Length
Mean[f1observ]
Variance[f1observ]
MeanCI[f1observ]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},100000];]
(f2observ=f2[randomuniform])//Length
{Integrate[f2[x],{x,0.,1}],Mean[f2observ]}
Variance[f2observ]
MeanCI[f2observ]

Korelacni koeficient mezi X a (X^2), X\[Distributed]R(0,1)

Clear[correlationcoefficient];
correlationcoefficient[f_,g_,distribution_]:=
	Module[{e1,e2,v1,v2,cov12},
		e1=ExpectedValue[f[u],distribution,u]//N;
		e2=ExpectedValue[g[u],distribution,u];
		v1=ExpectedValue[((f[u])^2),distribution,u]-(e1^2);
		v2=ExpectedValue[((g[u])^2),distribution,u]-(e2^2);
		cov12=ExpectedValue[f[u]g[u],distribution,u]-e1 e2;
		(cov12/(Sqrt[v1 v2]))//Simplify
		]

correlationcoefficient[Function[u,u],Function[u,(u^2)],UniformDistribution[{0,1}]]

e1=ExpectedValue[u,du,u]
e2=ExpectedValue[(u^2),du,u]
cov12=ExpectedValue[(u^3),du,u]-e1 e2
v1=Variance[du]
v2=ExpectedValue[(u^4),du,u]-(e2^2)
corr=(cov12/(Sqrt[v1 v2]))//Simplify
corr//N

e1=ExpectedValue[f1[u],du,u]
e2=ExpectedValue[f2[u],du,u]
cov12=ExpectedValue[f1[u] f2[u],du,u]-e1 e2
v1=Variance[du]
v2=ExpectedValue[((f2[u])^2),du,u]-(e2^2)
corr=(cov12/(Sqrt[v1 v2]))//Simplify
corr//N

Korelovany vyber f(u)=(e^u)-1, g(u)=u

Mean[fg=f2[randomuniform]-f1[randomuniform]]
Length[fg]
Variance[fg]
MeanCI[fg]+intf1

Antiteticke promenne

Zakladni myslenka: Pouzit dva nebo vice zaporne korelovanych odhadu. Nejjednodussi pripad:

Priklad 1

correlationcoefficient[Function[u,u],Function[u,1-u],UniformDistribution[{0,1}]]

correlationcoefficient[Function[u,u],Function[u,1-u],ExponentialDistribution[1]]

Ukazka, ze neplati pro "rychle se menici" funkce:

correlationcoefficient[Function[u,Sin[10u]],Function[u,Sin[10(1-u)]],UniformDistribution[{0,1}]]

correlationcoefficient[Function[u,Sin[2u]],Function[u,Sin[2-2u]],UniformDistribution[{0,1}]]

Needs["HypothesisTesting`"]
f1[u_]:=u
f2[u_]:=(E^u)-1
du=UniformDistribution[{0,1}];
SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},10000];]
(f2observAnt1=f2[randomuniform])//Length
(f2observAnt2=f2[1-randomuniform])//Length
f2observAnt=(1/2)(f2observAnt1+f2observAnt2);
Mean[f2observAnt]
Variance[f2observAnt]
MeanCI[f2observAnt]

Regrese

Idea: Optimalne zkombinovat dva nebo vice odhadu pomoci v Gauss-Markovove modelu

Clear[yy];
yy[f_,u_]:={0.5(f[u]+f[1-u]),
0.25(f[0.5u]+f[0.5(1-u)]+f[0.5(1+u)]+f[1-0.5u])}

simulations[f_,n_]:=Map[yy[f,#]&,RandomReal[{0,1},n]]

(∫_0^1)Sin[5(x^2)]\[DifferentialD]x//N
s3=simulations[Sin[5((#)^2)]&,30000];
me=Mean[s3]
cov=Covariance[s3]
Det[cov]

βest=(1/({1,1}.(inv=Inverse[cov]).{1,1})){1,1}.inv.me

Total[inv]

βest=(1/(Total[Flatten[inv=Inverse[cov]]]))Total[inv].me

Vysvetleni:

{1,1}.(({{a, b}, {c, d}}))

(({{a, b}, {c, d}}))//Total
((({{a, b}, {c, d}}))//Flatten)//Total

Clear[regr];
regr[f_,n_,seed___]:=Module[
{sim,me,cov,inv,den1,β},
If[seed>0,SeedRandom[seed]];
sim=Map[{0.5(f[#]+f[1-#]),
0.25(f[0.5#]+f[0.5(1-#)]+f[0.5(1+#)]+f[1-0.5#])}&,
RandomReal[{0,1},n]];
me=Mean[sim];
cov=Covariance[sim];
inv=Inverse[cov];
den1=(1/(Total[Flatten[inv]]));
β=den1 Total[inv].me;
Print["β est = ",β,"  with  variance = ",(den1/n),"  and standard error = ",(Sqrt[(den1/n)])]
]

regr[Sin[5((#)^2)]&,100000]

regr[(E^(#))-1&,10000,13]
(∫_0^1)((E^x)-1)\[DifferentialD]x//N

Balicek (Package): struktura a tvorba

Pripona souboru Package je *.m. Muze byt matouci, MATLAB ma to same. Struktura viz priklad

Na siti vyhledejte Skeleton.m. Odkazy viz Maeder: Programming in Mathematica

Kde muze Mathematica najit balicek?

Tri moznosti: (i) tam, kam vede cesta, (ii) pridat cestu (iii) tam, odkud je vyvolana seance (spusten notebook). (iii) je nejjednodussi.

$Path (* (i) sem vedou cesty *)

AppendTo[$Path,"C:\\Documents and Settings\\hurt"] (* (ii) Pridatcestu *)

?*Directory*

Nejjednodussi:

(* iii *)
Directory[]
SetDirectory[NotebookDirectory[]]
Directory[]

Take je mozne mit package jako soucast notebooku, ve kterem jej potrebuji.

Priklad IPM (An Introduction to Programming with Mathematica)

Autori:
Paul R. Wellin
Organization: Wolfram Research, Inc.

Richard J. Gaylord 
Organization: University of Illinois at Urbana-Champaign
Department: Department of Material Science and Engineering

Samuel N. Kamin

http://library.wolfram.com/infocenter/Books/5169/#downloads

Names["ANOVA`*"]

SetDirectory[NotebookDirectory[]]
<<RandomWalks`
Names["IPM3`RandomWalks`*"]

?ShowWalk

walk1D[n_] := NestList[#1 + (-1)^Random[Integer] & , 0, n]

walk1D[10]

RandomWalk[10,Dimension->2]

Package RandomWalks

?RandomWalk

BeginPackage["IPM3`RandomWalks`"]
Needs["Graphics`ArgColors`"]
Needs["Graphics`Polyhedra`"]
RandomWalk::"usage"="RandomWalk[n] generates an n-step walk in two dimensions. The default behavior gives a lattice walk with steps in one of the four compass directions. Setting the option LatticeWalk->False will generate an off-lattice walk. Setting the option Dimensions->3 will generate a 3-dimensional walk; Dimensions->1 gives a 1-dimensional walk.";

LatticeWalk::"usage"="LatticeWalk->val is an option to RandomWalk that determines whether the random walk will be a lattice walk or an off-lattice walk. Possible values are True and False.";

ShowWalk::"usage"="ShowWalk[walk] displays a two-dimensional random walk connecting each site with a line. Any Graphics options can be passed to ShowWalk. E.g., ShowWalk[walk,  Background->GrayLevel[0]] to produce a black background.";

AnimateWalk::"usage"="AnimateWalk[walk, opts] creates an animation of a two-dimensional random walk. A red ball will be seen to move to the current position in the walk to aid in visualizing the animation.";
RandomWalk::rwn="Argument `1` is not a positive integer.";
Options[RandomWalk]={LatticeWalk->True,Dimension->2}
Begin["`Private`"]
walk1D[n_] := NestList[#1 + (-1)^Random[Integer] & , 0, n]

walk1DOffLattice[n_]:=
  FoldList[Plus,0,Table[Random[Real,{-1,1}],{n}]]

walk2D[n_]:=
  Module[{NSEW={{0,1},{1,0},{0,-1},{-1,0}}},
    FoldList[Plus,{0,0},
      NSEW[[
        Table[Random[Integer,{1,4}],{n}]]]]]

walk2DOffLattice[n_]:=
  FoldList[Plus,{0,0},
    Map[{Cos[#],Sin[#]}&,Table[Random[Real,{0,2 π}],{n}]]]

walk3D[n_] := Module[{NSEW3 = Sqrt[2]*Vertices[Cube]}, 
   FoldList[Plus, {0, 0, 0}, 
    NSEW3[[Table[Random[Integer, {1, 8}], {n}]]]]]

walk3DOffLattice[n_] := Module[{x, y, z}, 
   x[ϕ_, θ_] := Sqrt[1 - Cos[ϕ]^2]*Cos[θ]; 
    y[ϕ_, θ_] := Sqrt[1 - Cos[ϕ]^2]*Sin[θ]; 
    z[ϕ_] := Cos[ϕ]; FoldList[Plus, {0, 0, 0}, 
     Apply[{x[#1, #2], y[#1, #2], z[#1]} & , 
      Table[Random[Real, {0, 2*Pi}], {n}, {2}], {1}]]]

RandomWalk[n_,opts___?OptionQ]:=Module[{dim,latticeQ},
    If[Not[IntegerQ[n]&&n>0],
      Message[RandomWalk::rwn,
        n],{latticeQ,
          dim}={LatticeWalk,Dimension}/.Flatten[{opts,Options[RandomWalk]}];
      Which[
        dim==1,If[latticeQ,walk1D[n],walk1DOffLattice[n]],
        dim==2,If[latticeQ,walk2D[n],walk2DOffLattice[n]],
        dim==3,If[latticeQ,walk3D[n],walk3DOffLattice[n]]
        ]]]

ShowWalk[coords_,opts___]:=Which[
    Length[Dimensions[coords]]==1,
    ListPlot[coords,opts,PlotJoined->True],
    Dimensions[coords][[2]]==2,
    Show[Graphics[Line[coords],opts,AspectRatio->Automatic]],
    Dimensions[coords][[2]]==3,
    Show[Graphics3D[Line[coords],opts,AspectRatio->Automatic]]]
AnimateWalk[coords_,opts___]:=
  Scan[Show[
        Graphics[{{RGBColor[1,0,0],PointSize[0.025],
              Point[coords[[#1]]]},
            Line[Take[coords,#1]]}],opts,AspectRatio->Automatic,
        PlotRange->({Min[#1]-0.2,Max[#1]+0.2}&)/@Transpose[coords]]&,
    Range[2,Length[coords]]]
End[]
EndPackage[]

Vicerozmerna nahodna prochazka (m-rozmerna)

Table[(((-1))^(RandomInteger[{0,1},4])),{5}]//TableForm

Clear[walkJH];
walkJH[m_(* dimenze *),n_(* delka *)]:=NestList[#1+(((-1))^(RandomInteger[{0,1},m]))&,Table[0,{m}],n]

walkJH[1,10]

ListPlot[walkJH[1,10000]//Flatten]

ListPlot[walkJH[2,10000],Joined->True]

ListPlot3D[walkJH[3,10000]]

Financni funkce verze 8

Financni funkce (pripomenuti)

Annuity

AnnuityDue

Cashflow

TimeValue

EffectiveInterest

FinancialBond

FinancialDerivative

TradingChart

InteractiveTradingChart

TimeValue (nektere detaily)

Cashflow umoznuje zadavat kalendarni data:

K datu {2009,1,1}

(************************************* 
Hodnoceni k datu {2009,1,1}, hodnotici urokova mira .08
*****************************)
TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],.08,{2009,1,1}]

K datu {2008,1,1}

(************************************* 
Hodnoceni k datu {2008,1,1}, hodnotici urokova mira .08
*****************************)
TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],.08,{2008,1,1}]

Plot[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2000,1,1}],{r,0,0.12},AxesLabel->{"r","PV"}]

K datu {2010,1,1}

Plot[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2010,1,1}],{r,0,0.12},AxesLabel->{"r","PV"}]

{FindRoot[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2000,1,1}],{r,0.1}],FindRoot[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2010,1,1}],{r,0.1}]}

EffectiveInterest

EffectiveInterest[r,f]//TraditionalForm

Najit nominalni urokovou miru odpovidajici efektivni pri ctvrtletnim uroceni:

FindRoot[EffectiveInterest[r,1/4]==.05,{r,.2}]

Prevod okamzitych (spotovych) urokovych mer na forwardove spolu s intervaly, na kterych plati:

EffectiveInterest[{1/12->.01,3/12->.015,6/12->.02,1->.025,2->.03,3->.032,5->.035,7->.04,10->.045}]
ListPlot[%]

FinancialBond

Issue price of a 30-year annual coupon bond of $1000 with a 6% yield:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->30},{"InterestRate"->0.06,"Settlement"->0}]

Price of a 10-year semiannual coupon bond with a 5% yield 9 months after the issue date:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->10,"CouponInterval"->(1/2)},{"InterestRate"->.05,"Settlement"->(9/12)}]

Price of a 4% quarterly coupon bond maturing December 31, 2030 and settling September 5, 2013:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.04,"Maturity"->{2030,12,31},"CouponInterval"->(1/4)},{"InterestRate"->.06,"Settlement"->{2013,9,5}}]

Accrued interest for a semiannual bond maturing on December 31, 2030 and settling on November 12, 2010, using the "Actual/360" day count convention:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.07,"Maturity"->{2030,12,31},"CouponInterval"->(1/2)},{"InterestRate"->0.06,"Settlement"->{2010,11,12},"DayCountBasis"->"Actual/360"},"AccruedInterest"]

YTM: Implied yield to maturity for a quarterly coupon bond settling on January 6, 2018 and valued at $900:

FindRoot[FinancialBond[{"FaceValue"->1000,"Coupon"->0.05,"Maturity"->{2018,6,31},"CouponInterval"->(1/4)},{"InterestRate"->y,"Settlement"->{2012,1,6}}]==900,{y,.1}]

Poznamka: Properties:

({{"Value", price adjusted for accrued interest}, {"FullValue", unadjusted price}, {"AccruedInterest", accrued interest at settlement}, {"Duration", Macaulay duration}, {"ModifiedDuration", modified duration}, {"Convexity", convexity}, {"CouponPeriodDays", days in coupon period containing the settlement}, {"CouponToSettlementDays", days from previous coupon to settlement}, {"SettlementToCouponDays", days from settlement to next coupon}, {"NextCouponDate", next coupon date}, {"PreviousCouponDate", previous coupon date}, {"RemainingCoupons", remaining coupon payments}})

FinancialDerivative

FinancialDerivative[]//Length

Plot[FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> st, "Dividend"->0.05}],{st,30,60}]

Plot[Table[FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->t},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> st, "Dividend"->0.05}],{t,{0,0.01,0.1,0.5,1}}],{st,30,60}]

Manipulate[Plot[FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->t},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> st, "Dividend"->0.05}],{st,40,60},PlotRange->{0,15}],
{{t,0,"Doba do vyprseni"},0,1}]

Hodnota evropske CALL a citlivosti (Value and all Greeks of a European call option):

FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50, "Dividend"->0.05}, {"Value", "Greeks"}]

Podle (podle mne pochybne konstrukce, viz DHS)  teorie by mela byt hodnota americke CALL stejna jako evropske. Ejhle, neni!!!

FinancialDerivative[{"American","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50, "Dividend"->0.05}, {"Value", "Greeks"}]

TradingChart

InteractiveTradingChart

Statistical functions in version 8

References

Generally

http://www.wolfram.com/mathematica/new-in-8/

New in 8.0: Alphabetical Listing

http://reference.wolfram.com/mathematica/guide/NewIn80AlphabeticalListing.html

More specific

Documentation Center: guide/ProbabilityAndStatistics

Probability and Statistics Solvers and Properties

http://www.wolfram.com/mathematica/new-in-8/probability-and-statistics-solvers-and-properties/

Parametric Probability Distributions

http://www.wolfram.com/mathematica/new-in-8/parametric-probability-distributions/

Nonparametric, Derived, and Formula Distributions

http://www.wolfram.com/mathematica/new-in-8/nonparametric-derived-and-formula-distributions/

Parameter Estimation and Testing

http://www.wolfram.com/mathematica/new-in-8/parameter-estimation-and-testing/

Parametric Statistical Distributions

http://reference.wolfram.com/mathematica/guide/ParametricStatisticalDistributions.html

Statistical Visualization

http://www.wolfram.com/mathematica/new-in-8/statistical-visualization/

Financial Visualization

http://www.wolfram.com/mathematica/new-in-8/financial-visualization/

Built-in Financial Computations

http://www.wolfram.com/mathematica/new-in-8/built-in-financial-computations/

http://www.wolfram.com/solutions/industry/financial-risk-management/

http://www.wolfram.com/solutions/industry/financial-risk-management/

Financial Functions

Annuity

AnnuityDue

Cashflow

TimeValue

EffectiveInterest

FinancialBond

FinancialDerivative

TradingChart

InteractiveTradingChart

Other Financial Topics

In Documentation Center find “financial”.

See also

Financial & Economic Data
Mathematica has immediate built-in access to current and historical financial and economic data.
FinancialData — current and historical stock, fund, index, currency, etc. data
“Price”  ▪  “CumulativeReturn”  ▪  “Volatility50Day”  ▪  “MarketCap”  ▪  “PERatio”  ▪  “EBITDA”  ▪  “DividendYield”  ▪  “SICCode”  ▪  ...
	
CountryData — socioeconomic data and time series for all countries
“GDP”  ▪  “ExportValue”  ▪  “FixedInvestment”  ▪  “ElectricityProduction”  ▪  “OilReserves”  ▪  ...
“CurrencyCode”  ▪  “CurrentUnit”  ▪  “ExchangeRate”  ▪  ...
	
DateListPlot — plot a time series
MovingAverage  ▪  ExponentialMovingAverage  ▪  LinearModelFit

See also

Finance
Mathematica has fully integrated support for many of the tools used in classical and modern finance. These capabilities include financial instrument valuation, advanced time value of money computations, and advanced financial charting with a library of technical indicators. Mathematica also provides immediate access to a large array of financial and economic data, and contains financial import and export tools for working with external data.
Time Value Computations
TimeValue — time value of money computations
EffectiveInterest — effective interest rate conversions
Annuity  ▪  AnnuityDue  ▪  Cashflow
	
Financial Instruments
FinancialDerivative — financial derivatives valuation
FinancialBond — financial bond valuation
	
Financial Visualization »
CandlestickChart — candlestick and OHLC chart of open-high-low-close prices
InteractiveTradingChart — interactive trading chart including financial indicators
TradingChart  ▪  DateListPlot  ▪  RenkoChart  ▪  ...
	
Financial Indicators »
"Open"  ▪  "High"  ▪  "Low"  ▪  "Close"  ▪  "Volume"  ▪  "TypicalPrice"  ▪  "TrueRange"  ▪  "SimpleMovingAverage"  ▪  "ExponentialMovingAverage"  ▪  "BollingerBands"  ▪  "KeltnerChannels"  ▪  "RelativeStrengthIndex"  ▪  "FastStochastic"  ▪  "SlowStochastic"  ▪  "FullStochastic"  ▪  "RateOfChange"  ▪  "UltimateOscillator"  ▪  ...
	
Financial and Economic Data »
FinancialData — current and historical stock, fund, index, currency, etc. data
CountryData — socioeconomic data and time series for all countries
DateListPlot  ▪  MovingAverage  ▪  ExponentialMovingAverage  ▪  LinearModelFit
	
Importing & Exporting External Data »
Import — import from files and URLs
Export — export data, graphics, and expressions to files
"Table"  ▪  "CSV"  ▪  "TSV"  ▪  "XLS"  ▪  "ODS"  ▪  "SXC"  ▪  "MDB"  ▪  "DBF"  ▪  "DIF"  ▪  ...

Probability

Available distributions

distributions=?*Distribution

Probability, Nprobability, Conditional probability, ProbabilityDistribution

Probability[event, distribution]

𝒟01=NormalDistribution[0,1];
{Probability[x≤3,x\[Distributed]𝒟01],NProbability[x≤3,x\[Distributed]𝒟01]}

{Probability[x>4\[Conditioned]x>3,x\[Distributed]𝒟01],NProbability[x>4\[Conditioned]x>3,x\[Distributed]𝒟01]}

{Probability[x>4\[Conditioned]x>3,x\[Distributed]ExponentialDistribution[λ]]//FullSimplify,Probability[x>w+h\[Conditioned]x>w,x\[Distributed]ExponentialDistribution[λ],Assumptions->w>0&&h>0]}

Examples of multivariate distributions

?DirichletDistribution

PDF[DirichletDistribution[{α1,α2,α3,α4}],{x1,x2,x3}]

Coincidence with B-distribution:

{PDF[DirichletDistribution[{α1,α2}],{x}]//FullSimplify,PDF[BetaDistribution[α1,α2],x]}

{Probability[x+(y^2)+(z^2)<(3/10)\[Conditioned](x==1/7∧y==1/4),{x,y,z}\[Distributed]DirichletDistribution[{1,3,4,7}]],NProbability[x+(y^2)+(z^2)<(3/10)\[Conditioned](x==1/7∧y==1/4),{x,y,z}\[Distributed]DirichletDistribution[{1,3,4,7}]]}

Empirical distribution

r=RandomVariate[NormalDistribution[],10000];
{Probability[x^2+3x<11,x\[Distributed]r],Probability[x^2+3x<11,x\[Distributed]NormalDistribution[]]}
N[%]

Independent random variables

Probability[1<x<5&&-3<y<8,{x\[Distributed]ExponentialDistribution[1],y\[Distributed]NormalDistribution[]}]

Nonliear conditions in events

Probability[(x^2<5&& (2x+3)<5)||(x^3+x>11),x\[Distributed]ExponentialDistribution[2]]
N[%]

Product distribution (the join distribution of independent components)

𝒟product=ProductDistribution[ExponentialDistribution[1],NormalDistribution[]];
PDF[𝒟product,{x,y}]

Probability[1<x<5&&-3<y<8,{x,y}\[Distributed]𝒟product]

Tables of distributions

ℋ=Table[NormalDistribution[μ,(Sqrt[μ+1])],{μ,Range[5]-1}];
Plot[Map[PDF[#,x]&,ℋ],{x,-3,10}]

Derived distributions

TransformedDistribution

𝒟trans=TransformedDistribution[(x^2)+(y^2),{x,y}\[Distributed]ProductDistribution[NormalDistribution[],2]]
𝒟trans2=TransformedDistribution[(x^2)+(y^2),{x\[Distributed]NormalDistribution[],y\[Distributed]NormalDistribution[]}]
{PDF[𝒟trans,x]//Evaluate,PDF[𝒟trans2,x],PDF[ChiSquareDistribution[2],x]}

ChiSquareDistribution[2]

Copula distribution

Motivation: The class of multivariate distributions based on univariate generalizations is rather poor. Copulas provide fairly wide range of distributions which are flexible enough. Here we deal only with the case of two-dimensional copulas. Suppose that C((Subscript[u, 1]),(Subscript[u, 2]))  is a distribution function of the pair of random variables taking on the values from the interval [0,1] (kernel of the copula), F and G distribution functions of random variables X and Y, respectively. Then C(F(x),G(y)) is the distribution function of the copula generated by the kernel C and marginal distributions F and G. In Mathematica more than 10 kernels.

Product copula

CDF[𝒞1=CopulaDistribution["Product",{UniformDistribution[],UniformDistribution[]}],{x,y}]

Plot3D[Evaluate@PDF[𝒞1, {x,y}], {x,-0.2,1.2}, {y,-0.2,1.2},PlotRange->All]

𝒟cop1=CopulaDistribution["Product",{GeometricDistribution[.3],PoissonDistribution[4]}];

DiscretePlot3D[PDF[𝒟cop1,{x,y}],{x,0,9},{y,0,9},ExtentSize->Full]

data=RandomVariate[𝒟cop1,(10^5)];
Histogram3D[data,{{0,9,1},{0,9,1}},"PDF"]

ProbabilityDistribution

Enables to define users probability distribution by a density

𝒟prob=ProbabilityDistribution[(Sqrt[2]/Pi)(1/(1+x^4)),{x,-Infinity,Infinity}];

{Mean[𝒟prob],Variance[𝒟prob]}

Discrete distribution

𝒟discr=ProbabilityDistribution[2^(-k-1),{k,0,∞,1}];
Sum[PDF[𝒟discr,k],{k,0,∞}]

Estimation

EstimatedDistribution

data=RandomVariate[GammaDistribution[1,2],(10^2)];
edist=EstimatedDistribution[data,GammaDistribution[α,β]]
Plot[{PDF[GammaDistribution[1,2],x],PDF[edist,x]},{x,0,10}]

EstimatedDistribution[data,GammaDistribution[α,β],ParameterEstimator->"MethodOfMoments"]

data=RandomVariate[BetaDistribution[2,3],1000];EstimatedDistribution[data,BetaDistribution[α,β]]EstimatedDistribution[data,BetaDistribution[α,β],ParameterEstimator->{"MaximumLikelihood",Method->"FindMaximum"}]

See 15 very interesting applications at the documantation of EstimatedDistribution.

LogLikelihood

LogLikelihood[NormalDistribution[μ,σ],Array[x,5]]

Empirical distribution

data=RandomVariate[NormalDistribution[],50];
𝒟emp=EmpiricalDistribution[data];
{DiscretePlot[PDF[𝒟emp,x],{x,data}],DiscretePlot[CDF[𝒟emp,x],{x,-4,4,.01}]}
{Moment[𝒟emp,2],Quantile[𝒟emp,0.95]}

Kernel estimator (jadrove odhady)

data=RandomVariate[NormalDistribution[],10^3];
𝒟=SmoothKernelDistribution[data];
Plot[PDF[𝒟,x],{x,-3,3}]

Choice of bandwidth:

𝒟1=SmoothKernelDistribution[data,.05];
𝒟2=SmoothKernelDistribution[data,1.];
Table[Plot[PDF[d,x],{x,-3,3}],{d,{𝒟1,𝒟2}}]

Best strategy: use Manipulate

Manipulate[Plot[Evaluate@PDF[SmoothKernelDistribution[data,bw],x],{x,-4,4},Filling->Axis,PlotLabel->{"band width",bw},Ticks->None],{bw,Table[0.025 j,{j,20}]},SaveDefinitions->True]

Choice of the kernel

Table[Plot[Evaluate@PDF[SmoothKernelDistribution[data,Automatic,i],x],{x,-4,4},Filling->Axis,PlotLabel->i,Ticks->None],{i,{"Biweight","Cosine","Epanechnikov","Gaussian","Rectangular","SemiCircle","Triangular","Triweight"}}]

Testing

DistributionFitTest (test-of-goodness of fit, test dobre shody)

Test for normality

data=RandomVariate[NormalDistribution[],2500];
DistributionFitTest[data]

ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]

Comparison of data and fitted

Show[Histogram[data,Automatic,"ProbabilityDensity"],Plot[PDF[ℋ["FittedDistribution"],x],{x,-5,5}]]

Location test

LocationTest

Variance test

VarianceTest

Various problems

Tukey Λ distribution

PPCCPlot[data_,{λ_,a_,b_}]:=Module[{len=10,qq,tlq,dq,cor},qq=(Range[len]-1/2)/len;dq=Quantile[data,qq];
cor[la_?NumberQ]:=Correlation[dq,Quantile[TukeyLambdaDistribution[la],qq]];
Plot[cor[λ],{λ,a,b}]]

pC=Interpolation[Transpose[{n,Table[Probability[x≤0.05,x\[Distributed]i],{i,ℋ}]}],InterpolationOrder->1];

Probability[x≤3,x\[Distributed]NormalDistribution[]]

Create Population Pyramids

{age,male,female}=Transpose[WolframAlpha["us population pyramid",{{"Result",1},"ComputableData"}]];

PairedBarChart[male,female,ChartStyle->{{ColorData["Pastel",1],ColorData["Pastel",0.2]},None,None},ChartLabels->{Placed[{"Men","Women"},Top],None,Placed[age,Axis,Row[#[[1]],"-"]&]},BarSpacing->{2 10^6,0,0},PlotLabel->"United States Population by Age, 2008",ImageSize->550,BaseStyle->{FontFamily->"Helvetica"}]

Tests of goodness of fit

𝒟=ProductDistribution[{StudentTDistribution[0,1, 9],2}];
sample=RandomVariate[𝒟,10^3];
edist=EstimatedDistribution[sample,MultivariateTDistribution[{{1,0},{0,1}},ν]]
DistributionFitTest[sample,edist]
DistributionFitTest[sample,edist,Automatic]
DistributionFitTest[sample,edist,All]

DistributionFitTest[sample,edist,Automatic]
DistributionFitTest[sample,edist,All]

data=RandomVariate[NormalDistribution[],2500];
DistributionFitTest[data]
ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]
Show[Histogram[data,Automatic,"ProbabilityDensity"],Plot[PDF[ℋ["FittedDistribution"],x],{x,-5,5}]]
data=RandomVariate[𝒹=ParetoDistribution[1,2],100];
DistributionFitTest[data,𝒹,{"TestDataTable","AndersonDarling"}]
ProbabilityPlot[data,𝒹]

Kernel estimators

SmoothHistogram[{RandomReal[NormalDistribution[1,0.1],100]},ColorFunction->Function[{x,y},Hue[y]],Filling->Axis,FillingStyle->Automatic]

data=RandomVariate[NormalDistribution[],1000];{Histogram[data,{0.25},"PDF"],SmoothHistogram[data,PlotStyle->Thick]}
Show[%]

SmoothKernelDistribution

As the bandwidth approaches infinity, the estimate approaches the shape of the kernel:

data={-.5,0,.5};

𝒦={"Gaussian","Epanechnikov","Biweight","Triweight","Rectangular","Triangular","Cosine","SemiCircle"};

Table[SmoothHistogram[data,{200,k},Ticks->None,PlotLabel->k],{k,𝒦}]

Convert between Formal Moments

Manipulate[TraditionalForm[f[n]]->TraditionalForm@MomentConvert[f[n],g],{{f,Moment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{g,CentralMoment,"target"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{n,3},1,20,1}]

Pictures

Two densities with the same moments

Assemble graphics from diverse sources (combining two pictures)

http://www.wolfram.com/mathematica/new-in-8/enhanced-illustration-and-drawing-tools/assemble-graphics-from-diverse-sources.html

Sinus random points

x=RandomReal[10,1000];
y=Sin[x] RandomReal[1,1000];
data=Transpose[{x,y}];
ListPlot[data]

Probability

Probability, NProbability, Probability distributions

How to enter: random variable x distributed according to a distribution 𝒟: 
x \[Distributed] 𝒟?
xEscdistEsc EscscDEsc

x\[Distributed]𝒟 (* or *)
Distributed[x,𝒟]

Keyboard characters

tutorial/StructuralElementsAndKeyboardCharacters

Various problems

Tukey Λ distribution

PPCCPlot[data_,{λ_,a_,b_}]:=Module[{len=10,qq,tlq,dq,cor},qq=(Range[len]-1/2)/len;dq=Quantile[data,qq];
cor[la_?NumberQ]:=Correlation[dq,Quantile[TukeyLambdaDistribution[la],qq]];
Plot[cor[λ],{λ,a,b}]]

pC=Interpolation[Transpose[{n,Table[Probability[x≤0.05,x\[Distributed]i],{i,ℋ}]}],InterpolationOrder->1];

Probability[x≤3,x\[Distributed]NormalDistribution[]]

Create Population Pyramids

{age,male,female}=Transpose[WolframAlpha["us population pyramid",{{"Result",1},"ComputableData"}]];

PairedBarChart[male,female,ChartStyle->{{ColorData["Pastel",1],ColorData["Pastel",0.2]},None,None},ChartLabels->{Placed[{"Men","Women"},Top],None,Placed[age,Axis,Row[#[[1]],"-"]&]},BarSpacing->{2 10^6,0,0},PlotLabel->"United States Population by Age, 2008",ImageSize->550,BaseStyle->{FontFamily->"Helvetica"}]

Tests of goodness of fit

𝒟=ProductDistribution[{StudentTDistribution[0,1, 9],2}];
sample=RandomVariate[𝒟,10^3];
edist=EstimatedDistribution[sample,MultivariateTDistribution[{{1,0},{0,1}},ν]]
DistributionFitTest[sample,edist]
DistributionFitTest[sample,edist,Automatic]
DistributionFitTest[sample,edist,All]

DistributionFitTest[sample,edist,Automatic]
DistributionFitTest[sample,edist,All]

data=RandomVariate[NormalDistribution[],2500];
DistributionFitTest[data]
ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]
Show[Histogram[data,Automatic,"ProbabilityDensity"],Plot[PDF[ℋ["FittedDistribution"],x],{x,-5,5}]]
data=RandomVariate[𝒹=ParetoDistribution[1,2],100];
DistributionFitTest[data,𝒹,{"TestDataTable","AndersonDarling"}]
ProbabilityPlot[data,𝒹]

Kernel estimators

SmoothHistogram[{RandomReal[NormalDistribution[1,0.1],100]},ColorFunction->Function[{x,y},Hue[y]],Filling->Axis,FillingStyle->Automatic]

data=RandomVariate[NormalDistribution[],1000];{Histogram[data,{0.25},"PDF"],SmoothHistogram[data,PlotStyle->Thick]}
Show[%]

SmoothKernelDistribution

As the bandwidth approaches infinity, the estimate approaches the shape of the kernel:

data={-.5,0,.5};

𝒦={"Gaussian","Epanechnikov","Biweight","Triweight","Rectangular","Triangular","Cosine","SemiCircle"};

Table[SmoothHistogram[data,{200,k},Ticks->None,PlotLabel->k],{k,𝒦}]

Convert between Formal Moments

Manipulate[TraditionalForm[f[n]]->TraditionalForm@MomentConvert[f[n],g],{{f,Moment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{g,CentralMoment,"target"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{n,3},1,20,1}]

Pictures

Two densities with the same moments

Assemble graphics from diverse sources (combining two pictures)

http://www.wolfram.com/mathematica/new-in-8/enhanced-illustration-and-drawing-tools/assemble-graphics-from-diverse-sources.html

Sinus random points

x=RandomReal[10,1000];
y=Sin[x] RandomReal[1,1000];
data=Transpose[{x,y}];
ListPlot[data]

Probability

Probability, NProbability, Probability distributions

How to enter: random variable x distributed according to a distribution 𝒟: 
x \[Distributed] 𝒟?
xEscdistEsc EscscDEsc

x\[Distributed]𝒟 (* or *)
Distributed[x,𝒟]

Keyboard characters

tutorial/StructuralElementsAndKeyboardCharacters

Zaverecna vsehochut (avsak velmi podstatna!!!)

Mentor, přítel Odysseův, rádce jeho syna Telemacha praví:

Poznámka I (obecná): Budete-li kdekoli ve svých pracích užívat poznatky získané z této a jiných přednášek, citujte je. Vyhnete se možným problémům. Kompilace cizích výsledků není trestné. Jejich uvedení bez citace jest postižitelné.

Poznámka II (obecná, pramenící z nedávné zkušeností ): Absolventi MFF měli dlouho vysoký kredit. Teď jej částečně ztrácejí nejen svými neznalostmi, ale i bezcharakterním jednáním. Dbejte na etiku. Při změně zaměstnání nevynášejte know-how získané u předchozího zaměstnavatele.

Nelinearni regrese a vynosove krivky

There is a number of parametric models used for the either fitting or interpolation of empirical data coming from markets. We will deal with just two of them which, under the author's opinion are most vital.

Realna data, dvojice ve tvaru {splatnost ve dnech, urokova mira}:

plot3=ListPlot[(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}}),PlotStyle->{Red,PointSize[0.01]}]

Nelson-Siegel

Clear[nelsonsiegelspot]; (* Spot curve or zero-coupon *)
nelsonsiegelspot[b0_,b1_,b2_,c_,t_]:=b0+(((b1+b2) (1-(E^(-(t/c)))))/(t/c))-b2 (E^(-(t/c)))

The parameters have the following meaning: the function starts at level b0+b1 at time 0 and has asymptote at level b0; a local extreme is determined by parameter c.

Clear[t];
Limit[nelsonsiegelspot[b0,b1,b2,c,t],t->0]
Limit[nelsonsiegelspot[b0,b1,b2,c,t],t->∞,Assumptions->{b0>0∧b1>0∧b2>0∧c>0}]

Manipulating spot:

Manipulate[Plot[nelsonsiegelspot[b0,b1,b2,c,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3}],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{c,400,"c"},100,1000,10}
]

(* Forward curve *)
Clear[nelsonsiegelforward];
nelsonsiegelforward[b0_,b1_,b2_,c_,t_]:=b0+b1 (E^(-t/c))+b2 (t/c) (E^(-t/c))

Manipulating forward:

Manipulate[Plot[nelsonsiegelforward[b0,b1,b2,c,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3},PlotStyle->Orange],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{c,400,"c"},100,1000,10}
]

Fitting Nelson-Siegel: on the graph, estimated spot is blue and estimated forward is orange.

fitnelsonsiegel3= FindFit[ir3,nelsonsiegelspot[b0,b1,b2,c,t],{b0,b1,b2,c},t]
curvenel3=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel3;
curvenelfc3[t_]:=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel3;
curvenel3forward=nelsonsiegelforward[b0,b1,b2,c,t]/.fitnelsonsiegel3;
Show[plot3,Plot[curvenel3forward,{t,0,10000},PlotRange->All,PlotStyle->Orange],Plot[curvenel3,{t,0,10000},PlotRange->All],PlotRange->All]

Swensson

By adding the 5th parameter in Nelson-Siegel model we get Swensson model. Note that possibly two extremes are determined by parameters c1 and c2.

Clear[swenssonforward];
swenssonforward[b0_,b1_,b2_,b3_,c1_,c2_,t_]:=b0+b1 (E^(-t/c1))+b2 (t/c1) (E^(-t/c1))+b3 (t/c2)(E^(-t/c2))
Clear[swenssonspot];
swenssonspot[b0_,b1_,b2_,b3_,c1_,c2_,t_]=(1/t)(∫_0^t)swenssonforward[b0,b1,b2,b3,c1,c2,τ]\[DifferentialD]τ//Simplify

Limit[swenssonspot[b0,b1,b2,b3,c1,c2,t],t->0,Assumptions->{b0>0&&b1>0&&b2>0&&b3>0&&c1>0&&c2>0}]
Limit[swenssonspot[b0,b1,b2,b3,c1,c2,t],t->∞,Assumptions->{b0>0&&b1>0&&b2>0&&b3>0&&c1>0&&c2>0}]

Manipulate[Plot[swenssonspot[b0,b1,b2,b3,c1,c2,t],{t,0,5000},PlotRange->{{0,5000},{3,6}},AxesOrigin->{0,3}],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{b3,-3,"b3"},-5,-1,0.2},
{{c1,200,"c1"},100,1000,10},
{{c2,5000,"c2"},100,10000,10}
]

Manipulate[Plot[swenssonforward[b0,b1,b2,b3,c1,c2,t],{t,0,5000},PlotRange->{{0,5000},{3,6}},AxesOrigin->{0,3},PlotStyle->Orange],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{b3,-3,"b3"},-5,-1,0.2},
{{c1,200,"c1"},100,1000,10},
{{c2,5000,"c2"},100,10000,10}
]

fitswensson3= FindFit[ir3,swenssonspot[b0,b1,b2,b3,c1,c2,t],{b0,b1,b2,b3,c1,c2},t]
(curveswe3=swenssonspot[b0,b1,b2,b3,c1,c2,t]/.fitswensson3)//Simplify;
curveswefc3[t_]:=swenssonspot[b0,b1,b2,b3,c1,c2,t]/.fitswensson3
(curveswe3forward=swenssonforward[b0,b1,b2,b3,c1,c2,t]/.fitswensson3)//Simplify;
Show[plot3,Plot[curveswe3,{t,0,10000},PlotRange->All],Plot[curveswe3forward,{t,0,10000},PlotRange->All,PlotStyle->Orange],PlotRange->All]

Interpolation by polynomials

Polynomials of  the degree n-1

Stupid but sometimes observable method: for data size n, use polynomial of the degree n-1. Since the maximum allowed order for interpolation is 21 we use it for the demonstration of the catastrophic behaviour of the interpolated function

intpol3=Interpolation[ir3,InterpolationOrder->Min[Length[ir3],21]];
plotpol3=Plot[intpol3[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All];
Show[plot3,plotpol3,PlotRange->All]

Manipulate[intpol3=Interpolation[ir3,InterpolationOrder->order];
Show[plotpol3=Plot[intpol3[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All],plot3,PlotRange->All],
{{order,4,"order"},1,21,1}]

Dalsi priklad na vynosove krivky (PRIBOR)

PRIBOR (Prague Interbank Offered Rate) je prumerna urokova sazba, za kterou si banky navzajem poskytuji uvery na mezibankovnim trhu. PRIBID (Prague Interbank Bid Rate) je prumerna urokova sazba, za kterou si ceske banky mohou ukladat penize u jinych ceskych bank.

http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.jsp

http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=05.05.2009

Soubor ulozit jako Pribor20090505.txt

Stazeni 7 dnu stare tabulky (pozor na pracovni dny!)

"http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=11.05.2009"

Needs["Calendar`"]
Names["Calendar`*"]

DaysPlus[Date[],-6]

dennipribor1="http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=06.05.2009"
dennipribor2=StringDrop[dennipribor1,-10]
(* Posun o 7 dnu zpet *)
(actualdate1=DateString[DaysPlus[Date[],-7],{"Day",".","Month",".","Year"}])//FullForm
cestaaktualnipriborCNB=StringJoin[dennipribor2,actualdate1]

priboraktualni=Import[cestaaktualnipriborCNB,"Text"]

(* Desetinne carky nahradit teckami, oddelovace | carkami *)pribortxt2=StringReplace[priboraktualni,","->"."];
pribortxt3=StringReplace[pribortxt2,"|"->","];
(* Rozdelit podle \n *)
(pribortxt4=StringSplit[pribortxt3,"\n"])//FullForm
(* Prvni dva radky k nicemu nejsou *)
(pribortxt5=Drop[pribortxt4,2])//FullForm

Map[StringReplace[#," "->""]&,pribortxt5]//InputForm

"d"~~__->"d"

(* Rozdelime radky na retezce *)(pribortxt6=Map[StringSplit[#,","]&,pribortxt5])

Tato konstrukce je obtiznejsi. Misto jiraskovskeho popisu dodaneho poskytovatelem dat potrebuje strucne a jednoznacne vyjadreni. Povsimnete si, ze mezera pred "d", "m", "r" je pohlcena, a vse, co nasleduje za temito znaky je vynechano.

(pribortxt7=Map[StringReplace[#,{" d"~~__->"d"," m"~~__->"m"," r"~~__->"r"}]&,pribortxt6])//FullForm

Nyni vsechno vyjadrime ve dnech.

pribortxt7[[All,1]]//FullForm

(rules7=Thread[pribortxt7[[All,1]]-> Map[ToString,{1,7,14,30,60,90,180,270,360}]])

Zatim ale vsechno je ve forme retezcu.

(pribortxt8=Map[StringReplace[#,rules7]&,pribortxt7])//FullForm

Pripomenme: 1. sloupec termin(splatnost), 2. sloupec PRIBID, 3. sloupec PRIBOR. Konecna podoba dat, se kterymi jiz muzeme snadno pracovat.

(pribortxt9=pribortxt8//ToExpression)//TableForm

Zajimat nas bude jenom PRIBOR.

pribortxt10=Map[Drop[#,{2}]&,pribortxt9]

plot10=ListPlot[pribortxt10,PlotStyle->{Red,PointSize[0.015]}]

Zpracovani souboru pomoci ReadList (mozno vynechat)

"05.05.2009\ntermín|PRIBID|PRIBOR\n1 den|1,75|2,04\n7 dní|1,57|1,85\n14 dní|1,64|1,92\n1 m\:011bsíc|1,82|2,16\n2 m\:011bsíce|2,04|2,39\n3 m\:011bsíce|2,13|2,49\n6 m\:011bsíc\:016f|2,23|2,60\n9 m\:011bsíc\:016f|2,33|2,70\n1 rok|2,43|2,80"

SetDirectory[NotebookDirectory[]];
stream=OpenRead["Pribor20090505.txt"]
Read[stream,String]
Read[stream,String]
ReadList[stream,{Word,Word,Word},WordSeparators->"|",RecordSeparators->"\n"]
Close[stream]

Dale se postupuje stejne jako v pripade Import.

Odhad parametru vynosove krivky

Clear[nelsonsiegelspot]; (* Spot curve or zero-coupon *)
nelsonsiegelspot[b0_,b1_,b2_,c_,t_]:=b0+(((b1+b2) (1-(E^(-(t/c)))))/(t/c))-b2 (E^(-(t/c)))

fitnelsonsiegel10= FindFit[pribortxt10,nelsonsiegelspot[b0,b1,b2,c,t],{b0,b1,b2,c},t]
curvenel10=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel10;
curvenelfc10[t_]:=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel10;
Show[plot10,Plot[curvenel10,{t,0,360},PlotRange->All],PlotRange->All]

Cviceni: Opticky Nelson-Siegelova krivka prilis nesedi. Vyzkousejte Swenssonovu.

Automaticke stazeni aktualni tabulky PRIBOR

Stazeni 7 dnu stare tabulky (pozor na pracovni dny!)

"http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=11.05.2009"

Needs["Calendar`"]
Names["Calendar`*"]

DaysPlus[Date[],-6]

dennipribor1="http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=06.05.2009"
dennipribor2=StringDrop[dennipribor1,-10]
(* Posun o 7 dnu zpet *)
(actualdate1=DateString[DaysPlus[Date[],-7],{"Day",".","Month",".","Year"}])//FullForm
cestaaktualnipriborCNB=StringJoin[dennipribor2,actualdate1]

priboraktualni=Import[cestaaktualnipriborCNB,"Text"]

Zpracovani stazene tabulky

(* Desetinne carky nahradit teckami, oddelovace | carkami *)pribortxt2=StringReplace[priboraktualni,","->"."];
pribortxt3=StringReplace[pribortxt2,"|"->","];
(* Rozdelit podle \n *)
(pribortxt4=StringSplit[pribortxt3,"\n"])//FullForm
(* Prvni dva radky k nicemu nejsou *)
(pribortxt5=Drop[pribortxt4,2])//FullForm

Aktualni data jsou vsak k dispozici az po 13. hodine. Proto pro jistotu aktualni data budeme stahovat az po 14. hodine, jinak stahneme data z predchoziho dne. Pozor, nejsou osetreny neobchodni dny!!! Musi byt aktivovan kalendar.

DaysPlus[Date[],-1]

dennipribor1="http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=11.05.2009"
dennipribor2=StringDrop[dennipribor1,-10]
actualdate2=DateString[If[Date[][[4]]<14,DaysPlus[Date[],-1],Date[]],{"Day",".","Month",".","Year"}]
cestaaktualnipriborCNB2=StringJoin[dennipribor2,actualdate2]

pribotaktualni2=Import[cestaaktualnipriborCNB2,"Text"]

Import[cestaaktualnipriborCNB2,"Table"]

$ImportFormats

Options and 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 LN(ln (Subscript[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^(2r Δ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}],Joined->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]

Demonstrations

http://demonstrations.wolfram.com/search.html?query=Black-Scholes

Sharpe ratio and Markowitz portfolio

Data (Mesicni vynosy 5 akcii)

Nemate-li ve svem adresari soubor "Returns.dat", pouzijte "tvrdy vstup":

d1={{5},{"Month","AMR","BS","UTX","WX","Z"},{1,11.594`,1.515`,-1.579`,-9.36`,-0.408`},{2,22.078`,5.299`,9.947`,7.578`,0.`},{3,-8.511`,-8.`,9.446`,-16.327`,0.694`},{4,-3.488`,0.395`,5.164`,-8.571`,-3.862`},{5,-9.639`,-9.764`,-3.571`,-18.243`,-10.949`},{6,-12.`,6.195`,-0.943`,-11.57`,-8.197`},{7,-10.606`,-0.417`,-7.619`,-8.441`,-6.786`},{8,-1.695`,-6.276`,11.856`,-8.22`,-13.725`},{9,-12.069`,-9.091`,-3.286`,-10.227`,-2.273`},{10,27.451`,15.`,11.186`,-3.797`,2.791`},{11,-21.538`,-4.465`,6.55`,-4.021`,-8.14`},{12,-19.608`,-6.132`,8.75`,12.676`,-5.063`},{13,48.78`,27.156`,4.215`,12.5`,35.2`},{14,19.672`,-2.372`,7.355`,15.271`,12.121`},{15,1.37`,15.19`,13.889`,18.`,-0.901`},{16,4.054`,15.751`,17.073`,-0.847`,18.548`},{17,-15.584`,-9.81`,1.042`,23.884`,-10.958`},{18,3.077`,4.27`,20.833`,4.196`,14.912`},{19,-4.478`,-6.485`,-7.543`,-4.027`,-5.058`},{20,4.688`,11.679`,3.263`,-6.333`,4.918`},{21,-11.94`,-7.285`,3.189`,-18.939`,-3.125`},{22,-1.695`,-3.214`,-7.506`,-2.856`,27.742`},{23,15.517`,0.369`,-3.58`,5.882`,16.026`},{24,5.97`,-1.866`,-7.25`,-0.926`,-2.762`},{25,11.268`,26.616`,15.903`,16.822`,13.295`}};

(* Needs["Statistics`"] *)
(* Needs["Statistics`MultiDescriptiveStatistics`"] *)
SetDirectory[NotebookDirectory[]]
d1=Import["Returns.dat","Table"]
d1[[{1,2,3,4}]]
d2=Rest[d1];
d3=Rest[d2];
datareturns=Map[Rest,d3];
expected=Mean[datareturns]
(cov=Covariance[datareturns])//MatrixForm
(Correlation[datareturns])//MatrixForm

Cviceni: Prostudujte nasledujici alternativni pristup k extrakci potrebnych dat.

ret1=Import["E:\Dokumenty\VPFPM\VPFPM2007 2008\Returns.dat","Table"];
Short[ret1,5]
ret2=Rest[Rest[ret1]];
Short[ret2,5]
returns=Map[Rest,ret2];
Short[returns,5]

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.

2+2

Example 1 (short sales allowed)

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

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

?Maxim*

?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]
Total[xx/.sol2r[[2]]]

Example 3a (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]

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

xx3
expected

xx3.expected

Example 3b (inequality instead of equality)

crit3b=xx.cov.xx//Simplify;
cons3b=Apply[And,Thread[xx≥0]]&&Total[xx]≤ 1&&xx.expected==4;
sol3b=NMinimize[{crit3b,cons3b},xx]
Total[xx3b=xx/.Last[sol3b]]
xx3b.expected

Example 3c (inequalities instead of equalities)

crit3c=xx.cov.xx//Simplify;
cons3c=Apply[And,Thread[xx≥0]]&&Total[xx]≤ 1&&xx.expected≥4;
sol3c=NMinimize[{crit3c,cons3c},xx]
Total[xx3c=xx/.Last[sol3c]]
xx3c
xx3c.expected

Total[xx3d={1.289839352143948`*^-32,-4.880587180699747`*^-33,0.8613009345035485574`,3.132574895506314`*^-34,0.13860717880556636`}]
xx3d.expected
xx3d.cov.xx3d

Unfeasible requirement
(Neither the asset with maximum expected return reaches 4.3.)

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)

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](*,Method->"DifferentialEvolution"]*)

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`

Lagrange multipliers

Example 3 (Markowitz, short sales allowed, μ prescribed expected return)

Clear[n,ones,xx,λ1,λ2,μ];
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

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

expected
μ=3;
(* !!!! *)
eq1=Map[D[lagrange,#]&,xxx={xx,λ1,λ2}//Flatten]==0
Take[Flatten[xxx/.Solve[eq1,xxx]],Length[xx]]

Optimalizace

V optimalizacnich ulohach je mozne volit normy (Subscript[L, 2]), (Subscript[L, 1]), (Subscript[L, ∞]).

tutorial/ConstrainedOptimizationLinearProgramming

LinearProgramming

FindMinimum, FindMaximum, NMinimize, NMaximize, Minimize, and Maximize. LinearProgramming

Prace s databazi

Approach in elder versions

Simulation of a simple insurance database with fields {names, sex (Male, Female), age, insured amount}

c1=CharacterRange["A","Z"]
ToCharacterCode[c1]//Flatten

ToCharacterCode[{"M","F"}]//Flatten

genernames:=FromCharacterCode[RandomInteger[{65,90},5]]

sex2=Table[FromCharacterCode[RandomInteger[] (77-70)+70],{30}];
names2=Table[genernames,{30}];
ages2=RandomInteger[{20,85},30];
insured2=Table[RandomInteger[{20,100}] 1000,{30}];
database=Transpose[{names2,sex2,ages2,insured2}];

database

Find all men elder than 60, insured over 50000. Both of the following constructions give the same result.

Cases[database,_?(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]

sele1=Select[database,(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]

Sort these selecting records according to the insured amount:

Sort[sele1,#1[[4]]>#2[[4]]&]

Export the database in the Excel format:

SetDirectory[NotebookDirectory[]]
Export["database.xls",database,"XLS"]

Modern approach

SetDirectory[NotebookDirectory[]]

Simulation of a simple insurance database with fields {names, sex (Male, Female), age, insured amount}

n=30; (* The size of the database *)

genername:=StringJoin[RandomChoice[CharacterRange["A","Z"],RandomInteger[{3,10}]]]
names=Table[genername,{n}];

sex=RandomChoice[{"M","F"},n]

ages=RandomInteger[{20,85},n]

insuredamount=Table[RandomInteger[{20,100}] 1000,{n}]

head={"Name","Sex ","Age ","Insured amount"}

database=Transpose[{names,sex,ages,insuredamount}]
(database=PrependTo[database,head])//TableForm

Find all men elder than 60, insured over 50000. Both of the following constructions give the same result.

Join[{head},Cases[database,_?(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]]//TableForm

The same

sele=Select[database,(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]

Sort these selecting records according to the insured amount:

Sort[sele,#1[[4]]>#2[[4]]&]//TableForm

Export the database in the Excel format:

Directory[]

Export["database.xls",database]

Databaze: Asociovana pole (Associated Arrays) (precist)

Clear["names"];

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]]
fields=Map[data[#]&,firstrow]
matdrop=Drop[mat,1]//Transpose
MapThread[Set,{fields,matdrop}]

data["KLM"]

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

Znat az sem.

??fields

fields//Dimensions

fields

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

firstrow

{"names","BC","EFG","HIJ","KLM"}

FullForm[firstrow]
ToExpression[firstrow]

List["names","BC","EFG","HIJ","KLM"]

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

BC

Priklady naivniho sifrovani (precist)

Julius Caesar

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

a===First[abeceda]

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

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

text="dominika"
roztext=Characters[text]
codetext=roztext/.caesarcoderules
decodetext=codetext/.caesardecoderules
StringJoin[decodetext]

Lepsi zpusob sifrovani

Needs["Combinatorica`"]

SeedRandom[13];
permcoderules=Thread[abeceda->(perm=RandomPermutation[abeceda])]

permc=roztext/.permcoderules

permc/.Thread[perm->abeceda]
%//StringJoin

Kernel estimate

Jadrovy odhad - neparametricky odhad hustoty
x ... hodnota, ve ktere odhad pocitame
k ... jadro
h ... sirka okenka - nutno vyzkouset experimentalne; cim vetsi, tim hladsi krivku dostaneme

Clear[kernel];
kernel[x_,k_,h_,observ_]:=Module[{n=Length[observ],au},
		au=Map[k[(x-#)/h]&,observ];
		Total[au]/(n*h)
		             ]

k[x_]:=Exp[-x^2/2]/Sqrt[2*Pi] (* jadro = hustota N(0,1) *)

observ=Table[Random[]+2*Random[],{50}];

Plot[kernel[x,k,0.2,observ],{x,-1,4}]

Za cviceni sestrojte proceduru pomoci funkce Manipulate, kde budete manipulovat sirkou okenka.

Poznamka k praci s jednotkami, import atypickych formatu a jejich zpracovani

Import kursovniho listku CNB a jeho zpracovani

SetDirectory[NotebookDirectory[]]

Aktualni kursy

http ://www.cnb.cz/cs/financni_trhy/devizovy_trh/kurzy_devizoveho _trhu/denni_kurz.jsp

-> Stazeni v textovem formatu (nazev souboru denni_kurz.txt)

Rozsireni $Path:

(* AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008"] *)
SetDirectory[NotebookDirectory[]]

Cteni souboru

kursytxt=Import["denni_kurz20090429.txt"];
kursytxt//FullForm

Nejprve zamenit desetinnou carku (vymozenost z roku 1941 za Protektoratu Boehmen und Maehren, kdy CSN prevzala normu DIN, kterou se dosud v ceske matematicke literature psala desetinna tecka)

kursydesetinnatecka=StringReplace[kursytxt,","->"."]

Dale oddelovac "|" nahradit carkami:

(kursyoddelovace=StringReplace[kursydesetinnatecka,"|"->","])//FullForm

Jeednotlive zaznamy (radky) jsou zde oddeleny (jako v UNIX) "\n". Vyse uvedeny jeden kompaktni retezec "kursyoddelovace" tedy rozdelime na retezce podle radku:

(kursyporadcich=StringSplit[kursyoddelovace,"\n"])//FullForm

Prvni dva radky jsou pro dalsi praci nezadouci: odstranime je:

(kursyporadcichbezhlavicky=Drop[kursyporadcich,2])//FullForm

Vytvorime matici pocet men × 5:

(kursymatice5=Map[StringSplit[#,","]&,kursyporadcichbezhlavicky])//MatrixForm

Kazdy zaznam jej jeden zaznam tvoreny 5 polozkami. Vytvorime tedy novy soubor ve tvaru matice pocet men × 2:

Potrebujeme jen posledni dva sloupce: (Pozor, i cislo ma tvar String!!!)

(kursymatice2string=kursymatice5[[All,{4,5}]])//FullForm

(kursymatice2=kursymatice2string//ToExpression)//FullForm

kursymatice2

((kursymatice2[[1]])/(kursymatice2[[2]]))

(kursyAUD=Map[((#)/(kursymatice2[[1]]))&,kursymatice2])//TableForm

HoldForm[(AUD/AUD)]

Select[kursyAUD,#[[2]]>.15&]

Select[kursyAUD,0.63092≤ #[[2]]≤ 0.63093&]

Select[kursyAUD,#[[1]]==RUB/AUD&]

CZK=1;
EUR=25.30 CZK
15 EUR

Soucasne nahrazeni desetinne carky desetinnou teckou a ponechani carky v gramaticke vete.

Vychazi z pravidel ceskeho pravopisu: v cisle se za desetinnou carkou nepise mezera, zatimco v textu ano. Provizorne v textu nahradime retezec "carka, mezera" retezcem "XXX", v cislech nahradime carku teckou a pote v textech opet nahradime "XXX" carkou a mezerou.

retez={"245,345","Nevim, jestli to dokazu.", "5,34"}
r1=StringReplace[retez,", "->"XXX"]
(r2=StringReplace[r1,","->"."])
(r3=StringReplace[r2,"XXX"->", "])//FullForm

Ted si jeste postavime ulohu: Text ponechat textem a z cisel ve forme retezec vytvorime skutecne cisla.

ToExpression[r3]//FullForm

Takto jednodusse to vsak nejde. Ty elementy listu, ktere vrati $Failed, ponechame v nezmenene podobe.

Map[If[ToExpression[#]===$Failed,#,ToExpression[#]//Evaluate]&,r3]

FinancialData a CountryData

FinancialData

FinancialData["GE"]

DateListPlot[fd=FinancialData["GE","Jan. 1, 2000"],Joined->True]

fd//First
{Min[fd[[All,2]]],Max[fd[[All,2]]]}

f1=FinancialData["EUR/USD",{2009,2,3}]

Country Data

CountryData["SA","Flag"]

CountryData["CZ","Population"]

ListPlot[Map[Last,CountryData["CZ",{{"GDP"},{1995,2008}}]]]

DateListPlot[CountryData["CZ",{{"GDP"},{1995,2009}}],PlotStyle->PointSize[0.02]]

Konvoluce Cauchyovych rozdeleni

Clear[μ,σ,x];
PDF[CauchyDistribution[μ,σ],x]

Clear[x,y];
PDF[CauchyDistribution[0,1],x]

Clear[pdf

Konvoluce dvou a ctyr

Clear[pdf2];
pdf2[x_]=Integrate[PDF[CauchyDistribution[0,1],y] PDF[CauchyDistribution[0,1],x-y],{y,-∞,∞},Assumptions->x∈Reals]
pdf4=Integrate[pdf2[y] pdf2[x-y],{y,-∞,∞},Assumptions->x∈Reals]

Rozdeleni aritmetickeho prumeru nahodneho vyberu z Cauchyova rozdeleni

Rozdeleni souctu n Caychy iid ma hustotu

Clear[n,pdfn];
pdfn[x_,n_]=PDF[CauchyDistribution[0,n],x]//Simplify
pdfn[x,4]

Prumer ma tedy hustotu:

Clear[pdfnprumer];
pdfnprumer[x_,n_]=n pdfn[n x,n]
Clear[pdfnprumer];
pdfnprumer[x_,n_]=n pdfn[n x,n]//Simplify

(∫_(-∞)^∞)pdfnprumer[x,n]\[DifferentialD]x

Plot[Table[pdfnprumer[z,n],{n,1,10}]//Evaluate,{z,-5,5},PlotRange->All,AxesOrigin->{0,0}]

Rests

KJ questions

Q1 Jsou v systemu Mathematica zabudovany funkce IRR, "vnitřní míra výnosnosti" (nebo i např. jiné finanční funkce jako je durace apod.), neboť se mi nic takového nepodařilo najít, ale mám pocit, že to existuje. Našla jsem pár demonstrací a tam byla durace i IRR počítána "klasicky", dle známé rovnice..... Samozřejmě, že se to dá obejít, pokud si vyvtoříme funkci NPV(i) a pak hledáme řešení rovné 0, ale umí to Mathematica nějak sama pouhým zadáním jednotlivých cash flow (a popřípadě časovým určením)?

A1 Nejsou , ale

 (i) velmi jednoduche zpusoby programovani charakteristik peneznich toku budou uvedeny na prednasce, 

(ii) financni funkce jsou v komercnich balicich (ty jsou sice drahe, ale porad podstatne levnejsi a vykonejsi nez treba v systemu MATLAB) napr.

Unrisk
http://www.unrisk.com/
dokumentace viz
http://www.unrisk.com/download/pe/UnRiskPE30Docu.pdf

Financial Essentials (v pristi verzi bude soucasti systemu)
http://www.wolfram.com/products/applications/finance/

Derivatives Expert IV
http://www.wolfram.com/products/applications/derivatives/
dokumentace viz
http://www.wolfram.com/products/applications/derivatives/manual.pdf

V teto souvislosti jest vhodne tez upozornit na balik 

Time Series
http://www.wolfram.com/products/applications/timeseries/
inspirujici dokumentace viz
http://media.wolfram.com/documents/TimeSeriesDocumentation.pdf

Q2 Existuje nějaká klávesová zkratka na přepínání mezi jednotlivými notebooky?  
A2  Mozna, ze ano, ale nejjednodussi je: Menu ⟶ Window ⟶ spodni okno vybrat zadouci notebook. Ja sam pouziji organizaci na desktopu  s vhodnymi prekryvy, musim take pocitat s velmi castym volanim Help resp. F1.

Jak zasejes, tak sklidis. (Sow, Reap)

NestList[((#)^2)&,2,6]
Nest[((#)^2)&,2,6]

Reap[Nest[(If[#>1000,Sow[#]];((#)^2))&,2,6]]

Apply[Part[mat,Sequence[#]]&,{1,3}]

Poznamky

Monitor, EvaluationMonitor, Slavik

Priklad z FinancialData

Intervalova aritmetika

Interval[{3,1}]
Interval[{3,1}]//Min

intervalComplement[interval_Interval,universe_:Interval[{-∞,∞}]]:=
IntervalIntersection[Interval[{-∞,Min[interval]},{Max[interval],∞}],universe]

intervalComplement[Interval[{-3,3}]]
intervalComplement[Interval[{-3,3}],Interval[{-3,8}]]

Interval[{-∞,Min[Interval[{-3,3}]]},{Max[Interval[{-3,3}]],∞}]

Min[{-3,3}]

Komentare, Miscellaneous (nepovinne, neusporadane nicmene se zajimavymi podnety)

FAQ

Nula na osach v grafu 2D

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

Show[Plot[Sin[(x^2)],{x,0,π},Ticks->{{0,π/4,π/2,3π/4},Automatic}],Graphics[Text[Style["0",10],{0.06,-0.0.06}]]]

Funkce nedefinovane pro nektere argumenty

(0^0)

FullForm[HoldForm[(0^0)]]

Unprotect[Power]
Power[0,0]=1
Protect[Power]

(0^0)

Random matrices

ranmatposdef[n_]:= mat=RandomVariate[UniformDistribution[],{n,n}]

r30=ranmatposdef[30]//Eigenvalues

Table[{n,Max[Abs[ranmatposdef[n]//Eigenvalues]]},{n,30,200,10}]//TableForm

ranmat[n_]:= mat=RandomVariate[NormalDistribution[],{n,n}]

Table[{n,Max[Abs[ranmat[n]//Eigenvalues]]},{n,30,200,10}]//TableForm

$LicenseExpirationDate

Balicek (Package): zabudovane balicky (Add-On Packages, StandardExtraPackages)

Cim novejsi verze, tim vice funkci je primo v jadru (kernel). Z metodickych duvodu doporucuji proto nemazat starsi verze: zdrojove kody jsou zdrojem zabavy a pouceni. Nektere jsou vsak mimo jadro a musi se vyvolat pomoci funkci Get (Unix verse <<), resp. Needs. Velmi citlive na syntaxi!!! V nejnovejsi versi viz (za cviceni zkuste vytvorit odkaz na Documentation Center s timto heslem)

guide/StandardExtraPackages

Manipulace s $Path a adresari

{$UserDocumentsDirectory,$Path} (* Poznamenejme, ze $ je povazovano za pismeno, v Mathematice oznacuje systemove promenne a vnitrni promenne modulu. *)

Directory[]
SetDirectory["C:\\Program Files\\Wolfram Research\\Mathematica\\6.0\\AddOns\\Packages"];
Directory[]

FileNames[]

SetDirectory[$UserDocumentsDirectory]

SetDirectory[NotebookDirectory[]]

SetDirectory["C:\Program Files\Wolfram Research\Mathematica\7.0\AddOns\ExtraPackages"]
FileNames[]

SetDirectory["C:\Program Files\Wolfram Research\Mathematica\7.0\AddOns\Packages"]
FileNames[]

Ukazka balicku `Units

Help: UnitsPackage

Needs["Units`"]
units=Names["Units`*"]

Copied from Help (Units of volume): results in Input form. UnitsJH added.

unitsJH=({{Bag, Barrel}, {BoardFoot, Bucket}, {Bushel, Butt}, {Cord, Cup}, {Drop, Ephah}, {Fifth, Firkin}, {FluidDram, FluidOunce}, {Gallon, Gill}, {Hogshead, Jeroboam}, {Jigger, Last}, {Liter, Magnum}, {Minim, Noggin}, {Omer, Pint}, {Pony, Puncheon}, {Quart, RegisterTon}, {Seam, Shot}, {Stere, Tablespoon}, {Teaspoon, Tun}, {UKGallon, UKPint}, {WineBottle, \[Null]}})

There is the comma after "WineBottle". This command will eliminate it and, morevoer flatten the buttons in an one-dimensional list:

unitsJH=unitsJH//Flatten//Most

Methodology

unitsJH[[1]]//Head
unitsJH[[1]]//FullForm
f1=Apply[List,unitsJH[[1]]]
f1[[1]]//FullForm

unitsJH2=Map[First,Map[Apply[List,#]&,unitsJH]]

Conversion

Convert[Bag,Liter]
Convert[Liter,Pony]
Convert[Pony,Liter]

Transpose[{unitsJH2,Map[Convert[#,Liter]&,unitsJH2]}]//TableForm

Vetsina nasledujiciho vystupu je nesmyslna. Proto provedem selekci nasledujicim prikazem.

(konverzevsechjednoteknalitr={units//ToExpression,Map[Convert[#,Liter]&,units//ToExpression]}//Transpose)//Short

Select[konverzevsechjednoteknalitr,#[[1]]=!=#[[2]]&]//TableForm

Ukazka jinych balicku

Needs["ANOVA`"]
Names["ANOVA`*"]

?StudentNewmanKeuls

Needs["RegressionCommon`"]
Names["RegressionCommon`*"]

?AdjustedRSquared

FAQ

Alternativni vypreparovani korenu rovnic

(x^(Range[5])).RandomReal[{0,1},5]
sol=NSolve[(x^(Range[5])).RandomReal[{0,1},5]==0,x]

Cases[sol,Rule[x__,y__]->y,Infinity]

Vlastni funkce v modulu

mojefce1[x_]:=Sin[(x^2)]
mojefce2[x_]:=Cos[x]

procedura1[fce_]:=Module[{f1=mojefce1},
{Plot[2f1[x],{x,0,π}],Plot[fce[x],{x,0,π}]}]

procedura1[Cos]

Poznamka k vyvolavani funkci z balicku a informace o techto funkcich v Help

DaysPlus[Date[],10]

Needs["Calendar`"]

DaysBetween

FAQ

Zvyrazneni Pouze jednoho bodu v ListPlot

p1=ListPlot[{{5,Sin[5]}},PlotStyle->PointSize[0.03]];
p2=ListPlot[Table[{t,Sin[t]},{t,0,2π,2π/10}]];
Show[p1,p2]

Uvedeni veci na spravnou miru: spolecne grafy, direktivy podle prvniho grafu
(Nagy Gergely, {2010,04,19})

"V tom prikladu s disky: i tam se Show ridi se direktivami prvniho grafu:: Show[h1 = Graphics[{Orange, EdgeForm[Directive[Thick, Blue]], Disk[{-0.5, 0}, 0.7]}, Axes -> True, AxesOrigin -> {-0.5, 0}], h2 = Graphics[{Darker[Green], Disk[{0.5, 0}, 0.7]}, Axes -> True, AxesOrigin -> {0.5, 0}]] Show[h2, h1]"

Show[h1 = Graphics[{Orange, EdgeForm[Directive[Thick, Blue]], Disk[{-0.5, 0}, 0.7]}, Axes -> True, AxesOrigin -> {-0.5, 0}], h2 = Graphics[{Darker[Green], Disk[{0.5, 0}, 0.7]}, Axes -> True, AxesOrigin -> {0.5, 0}]] 
Show[h2, h1]

FAQ I (Zvetseni fontu pro popisy os v grafech)

Zvetseni fontu pro popisy os v grafech pro ucely exportu v EPS a zabudovani do textu

data={{{2006,10,1},8},{{2006,10,10},10},{{2006,10,20},12},{{2006,10,24},14},{{2006,11,5},15}, {{2006,11,15},20}};

(* Pozor na poradi! Zde  pouzit cesky zpusob, DD/MM!!! *)DateListPlot[data,DateTicksFormat->{"Day","/","MonthShort"},Frame->False,TicksStyle->Directive[Orange,16],PlotStyle-> PointSize[0.03]]
(* Vysledny graf oznacit kliknutim a za pravy roh roztahnout do citelne podoby *)

(* Pozor na poradi! Zde  pouzit anglosasky zpusob,MM/DD!!! *)DateListPlot[data,DateTicksFormat->{"MonthShort","/","Day"},Frame->False,TicksStyle->Directive[Orange,16],AxesLabel->{Style["   Datum",Italic,20],Style["Cena\n",Italic,20]},AxesStyle->Directive[16],Joined->True]
(* Vysledny graf oznacit kliknutim a za pravy roh roztahnout do zadane podoby *)

Plot[Fibonacci[x],{x,1,5},AxesOrigin->{0,0}]

FAQ II (Rozdeleni rezidui v regresi v pripade nejistoty o normalnim rozdeleni poruch)

Dotaz se týká problému, jak se změní šíře intervalů spolehlivosti a hladiny testů při porušení předpokladů o normalitě a případně homoskedasticitě. Tyto problémy souvisí s robustností statistických procedur. Elementární příklady (např. robustnost t-testu) viz Hurt, J.: Simulační metody, skripta SPN Praha 1982, str. 140.

Robustness of  the test in regression models violating assumptions of normality

Model Y=Xβ+ε, ε∼(Subscript[N, n])(0,(σ^2)(Subscript[I, n])) . We will simulate Y=α+X,β+ ε; without loss of generality α=0 (but unknown!!! and it should be estimated!!!), β=0,1,4,10, n=10,20,50,100, N ... number of simulated models. Usually N,*n is chosen to be a constant. Here we choose N*n=100000 so that N=10000,5000,2000,1000.

For ε we choose N(0,1) and Exp(1)-1, both with vanishing mean and unit variance. Exponential distribution is an exteme violation of normality. (Other possible choices for a more serious study are double exponential, t(1), t(2).) But for both distributions we will use the test procedures developed for a normal distribution.

For the sake of simplicity we will consider model (Subscript[Y, i])= α + β(10 i)/n+(Subscript[ε, i]), i=1,...,n.

Clear[simulationofregressionmodelnormal];simulationofregressionmodelnormal[α_,β_,n_,nN_,randomseed_:13]:=
Module[{yy=α+β(10/n) Range[n],xx=(10/n) Range[n],randomnumbers},
SeedRandom[randomseed];randomnumbers=RandomReal[NormalDistribution[0,1],{nN,n}];
Map[Function[εε,{xx,yy+εε}//Transpose],randomnumbers]
	    ]

Clear[simulationofregressionmodelexponential];simulationofregressionmodelexponential[α_,β_,n_,nN_,randomseed_:13]:=
Module[{yy=α+β(10/n) Range[n],xx=(10/n) Range[n],randomnumbers},
SeedRandom[randomseed];randomnumbers=RandomReal[ExponentialDistribution[1],{nN,n}]-1;
Map[Function[εε,{xx,yy+εε}//Transpose],randomnumbers]
	    ]

Clear[simulationofregressionmodelcauchy];simulationofregressionmodelcauchy[α_,β_,n_,nN_,randomseed_:13]:=
Module[{yy=α+β(10/n) Range[n],xx=(10/n) Range[n],randomnumbers},
SeedRandom[randomseed];randomnumbers=RandomReal[CauchyDistribution[0,1],{nN,n}];
Map[Function[εε,{xx,yy+εε}//Transpose],randomnumbers]
	    ]

The preliminary illustrations

ListPlot[simulationofregressionmodelnormal[0,4,50,3],Joined->True]

Simulation of 3 normal regression models with n=10 observations

(s15=simulationofregressionmodelnormal[0,1,10,3])//TableForm

ListPlot[s15,Joined->True]

The same for the exponetial disturbances

(s15e=simulationofregressionmodelexponential[0,1,10,3])//TableForm

Linear model fit based on normal distribution of errors with simulated normal errors

lmf15=Map[LinearModelFit[#,x,x]&,s15]

Linear model fit based on normal distribution of errors with simulated exponential errors

lmf15e=Map[LinearModelFit[#,x,x]&,s15e]

95 per cent confidence intervals for α and  β in 3 simulated models: one row = one simulated model = two intervals:

Normal errors

ci15=Map[#["ParameterConfidenceIntervals"]&,lmf15]

The average bounds

Mean[ci15]

Exponetial erors

ci15e=Map[#["ParameterConfidenceIntervals"]&,lmf15e]

Mean[ci15e]

The confidence intervals for α and β in separate files, transformed into Mathematica intrervals.

alpha15=Map[First,ci15]
beta15=Map[Last,ci15]
int15=Map[Interval,beta15]
(* Is 1.2 in the interval? *)
Map[IntervalMemberQ[#,1.2]&,int15]
(* Transform to 0, 1 by Boole *)
Boole[%]
(* Calculate the relative frequency *)
Mean[%]

Normal disturbances, n=100,N=1000

sn100N1000 = simulationofregressionmodelnormal[0, 1, 100, 1000];
lmfn100N1000 =Map[LinearModelFit[#,x,x]&,sn100N1000 ]; cin100N1000=Map[#["ParameterConfidenceIntervals"]&,lmfn100N1000];
Mean[ cin100N1000]

alphan100N1000=Map[First,cin100N1000];
betan100N1000=Map[Last,cin100N1000];
intervalalphan100N1000=Map[Interval,alphan100N1000];
intervalbetan100N1000=Map[Interval,betan100N1000];
(* Hits *)
alphahitn100N1000=Map[IntervalMemberQ[#,0]&,intervalalphan100N1000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn100N1000//N]}
betahitn100N1000=Map[IntervalMemberQ[#,1]&,intervalbetan100N1000]//Boole;
{"Observed significance value for β = ",Mean[betahitn100N1000//N]}

Exponential disturbances, n=100,N=1000

sn100N1000 = simulationofregressionmodelexponential[0, 1, 100, 1000];
lmfn100N1000 =Map[LinearModelFit[#,x,x]&,sn100N1000 ]; cin100N1000=Map[#["ParameterConfidenceIntervals"]&,lmfn100N1000];
Mean[ cin100N1000]

alphan100N1000=Map[First,cin100N1000];
betan100N1000=Map[Last,cin100N1000];
intervalalphan100N1000=Map[Interval,alphan100N1000];
intervalbetan100N1000=Map[Interval,betan100N1000];
(* Hits *)
alphahitn100N1000=Map[IntervalMemberQ[#,0]&,intervalalphan100N1000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn100N1000//N]}
betahitn100N1000=Map[IntervalMemberQ[#,1]&,intervalbetan100N1000]//Boole;
{"Observed significance value for β = ",Mean[betahitn100N1000//N]}

Cauchy disturbances, n=100,N=1000

sn100N1000 = simulationofregressionmodelcauchy[0, 1, 100, 1000];
lmfn100N1000 =Map[LinearModelFit[#,x,x]&,sn100N1000 ]; cin100N1000=Map[#["ParameterConfidenceIntervals"]&,lmfn100N1000];
Mean[ cin100N1000]

alphan100N1000=Map[First,cin100N1000];
betan100N1000=Map[Last,cin100N1000];
intervalalphan100N1000=Map[Interval,alphan100N1000];
intervalbetan100N1000=Map[Interval,betan100N1000];
(* Hits *)
alphahitn100N1000=Map[IntervalMemberQ[#,0]&,intervalalphan100N1000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn100N1000//N]}
betahitn100N1000=Map[IntervalMemberQ[#,1]&,intervalbetan100N1000]//Boole;
{"Observed significance value for β = ",Mean[betahitn100N1000//N]}

Normal disturbances, n=20,N=5000

sn20N5000 = simulationofregressionmodelnormal[0, 1, 20, 5000];
lmfn20N5000 =Map[LinearModelFit[#,x,x]&,sn20N5000 ]; cin20N5000=Map[#["ParameterConfidenceIntervals"]&,lmfn20N5000];
Mean[ cin20N5000]

alphan20N5000=Map[First,cin20N5000];
betan20N5000=Map[Last,cin20N5000];
intervalalphan20N5000=Map[Interval,alphan20N5000];
intervalbetan20N5000=Map[Interval,betan20N5000];
(* Hits *)
alphahitn20N5000=Map[IntervalMemberQ[#,0]&,intervalalphan20N5000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn20N5000//N]}
betahitn20N5000=Map[IntervalMemberQ[#,1]&,intervalbetan20N5000]//Boole;
{"Observed significance value for β = ",Mean[betahitn20N5000//N]}

Exponential disturbances, n=20,N=5000

sn20N5000 = simulationofregressionmodelexponential[0, 1, 20, 5000];
lmfn20N5000 =Map[LinearModelFit[#,x,x]&,sn20N5000 ]; cin20N5000=Map[#["ParameterConfidenceIntervals"]&,lmfn20N5000];
Mean[ cin20N5000]

alphan20N5000=Map[First,cin20N5000];
betan20N5000=Map[Last,cin20N5000];
intervalalphan20N5000=Map[Interval,alphan20N5000];
intervalbetan20N5000=Map[Interval,betan20N5000];
(* Hits *)
alphahitn20N5000=Map[IntervalMemberQ[#,0]&,intervalalphan20N5000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn20N5000//N]}
betahitn20N5000=Map[IntervalMemberQ[#,1]&,intervalbetan20N5000]//Boole;
{"Observed significance value for β = ",Mean[betahitn20N5000//N]}

Cauchy disturbances, n=20,N=5000

sn20N5000 = simulationofregressionmodelcauchy[0, 1, 20, 5000];
lmfn20N5000 =Map[LinearModelFit[#,x,x]&,sn20N5000 ]; cin20N5000=Map[#["ParameterConfidenceIntervals"]&,lmfn20N5000];
Mean[ cin20N5000]

alphan20N5000=Map[First,cin20N5000];
betan20N5000=Map[Last,cin20N5000];
intervalalphan20N5000=Map[Interval,alphan20N5000];
intervalbetan20N5000=Map[Interval,betan20N5000];
(* Hits *)
alphahitn20N5000=Map[IntervalMemberQ[#,0]&,intervalalphan20N5000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn20N5000//N]}
betahitn20N5000=Map[IntervalMemberQ[#,1]&,intervalbetan20N5000]//Boole;
{"Observed significance value for β = ",Mean[betahitn20N5000//N]}

Exponential disturbances, n=10,N=10000

sn10N10000 = simulationofregressionmodelexponential[0, 1, 10, 10000];
lmfn10N10000 =Map[LinearModelFit[#,x,x]&,sn10N10000 ]; cin10N10000=Map[#["ParameterConfidenceIntervals"]&,lmfn10N10000];
Mean[ cin10N10000]

alphan10N10000=Map[First,cin10N10000];
betan10N10000=Map[Last,cin10N10000];
intervalalphan10N10000=Map[Interval,alphan10N10000];
intervalbetan10N10000=Map[Interval,betan10N10000];
(* Hits *)
alphahitn10N10000=Map[IntervalMemberQ[#,0]&,intervalalphan10N10000]//Boole;
{"Observed significance value for α = ",Mean[alphahitn10N10000//N]}
betahitn10N10000=Map[IntervalMemberQ[#,1]&,intervalbetan10N10000]//Boole;
{"Observed significance value for β = ",Mean[betahitn10N10000//N]}

FAQ III

(χ^2)-test dobre shody pro tento konkretni pripad (starsi pristup)

cetnosti2=Map[Last,cetnosti]
pocetruznychcisel=Length[cetnosti2]
pocetpozorovani=Total[cetnosti2]
((((cetnosti2-pocetpozorovani/pocetruznychcisel))^2)/(pocetpozorovani/pocetruznychcisel))
chi2=Total[((((cetnosti2-pocetpozorovani/pocetruznychcisel))^2)/(pocetpozorovani/pocetruznychcisel))]//N

Print["Hodnota (χ^2) na ",pocetruznychcisel-1," stupnich volnosti = ",chi2,"\nOdpovidajici p - hodnota = ",1-CDF[ChiSquareDistribution[pocetruznychcisel-1],chi2]]

Complement[data,{1}]//Length

DeleteCases[data,1]//Length
fun[DeleteCases[data,1]]

1

Ruzne

Pocasi

WeatherData["Prague","Temperature"]

Text[Grid[{CityData[#,"Name"],Row[{Round@WeatherData[#,"Temperature"],"°C"}]}&/@Take[CityData[],6],Frame->All]]

WeatherData["Prague", #]& /@ {"Temperature", "Pressure", "WindSpeed"}

DateListPlot[WeatherData["Prague","MeanTemperature",{{2008,1,1},{2009,5,2},"Day"}],Joined->True]

DateListPlot[WeatherData["Prague","Temperature",{{2008,1,1},{2009,4,26},"Day"}],Joined->True]

stations=WeatherData[{"Prague",5},"NearestStation"]

WeatherData[{"Prague",5},"StationDistance"]

WeatherData[{"Prague",5},"StationDistance","Units"]

Slider (=posuvnik) bez poklepnuti mysi

Slider[0.5,AutoAction->True]

Speak

Speak["zmrzlina, zmrzlina, zmrzlina, zmrzlina. Once upon a midnight dreary, while I pondered, weak and weary, Over many a quaint and curious volume of forgotten lore-- While I nodded, nearly napping, suddenly there came a tapping, As of some one gently rapping--rapping at my chamber door."]

Konvoluce

Convolve[Exp[-x] UnitStep[x],Sin[x],x,y]

h=Exp[-x] UnitStep[x];
Convolve[h,UnitStep[x],x,y]

pdf=PDF[NormalDistribution[μ,σ],x];
Convolve[pdf,pdf,x,y]//Simplify[#,Assumptions->{σ>0}]&

Poloha pocitace nebo IP

FindGeoLocation[]

FAQ

Graficke zobrazeni minima nekolika funkci

Clear[f1,f2,fmin];
f1[x_]:=(x^2)
f2[x_]:=Sin[2x]
fmin[x_]:=Min[f1[x],f2[x]]
Plot[fmin[x],{x,-5,5}]

Nalezeni bodu zlomu

FullSimplify[Min[f1[x],f2[x]],x∈Reals]

FullSimplify[f1[x]-fmin[x]]

Reduce[{f1[x]==fmin[x],f2[x]==fmin[x]},x]

fminexpand[x_]=PiecewiseExpand[fmin[x]]

Plot[fminexpand[x],{x,-5,5}]

i1=Reduce[{f1[x]≤fminexpand[x],-5≤x<5},x]
i2=Reduce[{f2[x]≤fminexpand[x],-5≤x<5},x]

forder[x_]:=Ordering[{f1[x],f2[x]}]

Position[forder[2],2]

Reseni pro linearni funkce viz Lekce 9 -> FAQ -> Nalezeni bodu zlomu pro linearni funkce

Nalezeni bodu zlomu pro linearni funkce

Clear[g1,g2,g3];
g1[x_]:=2x+3;g2[x_]:=3x+4;g3[x_]:=-x+5;
Plot[{g1[x],g2[x],g3[x]},{x,-3,6}]
gmin[x_]:=Min[g1[x],g2[x],g3[x]]
gminexpand[x_]=PiecewiseExpand[gmin[x]]
Plot[gmin[x],{x,-3,6}]

(regions=Map[Reduce[(gminexpand[x]==#),x]&,{g1[x],g2[x],g3[x]}])//FullForm
(regions2=regions//.{Inequality->List,LessEqual->List,GreaterEqual->List}//Flatten)//FullForm
crossings=Select[regions2,NumericQ]
Union[crossings]

Problem z bakalarske prace (Celkove doby setrvani ve stavech markovskeho retezce)

Data ve tvaru {{(Subscript[dobasetrvani, 1]),(Subscript[stav, 1])},...,{(Subscript[dobastervani, n]),(Subscript[stav, n])}}

SeedRandom[13];
markcasy=RandomReal[ExponentialDistribution[1],1000];
markstavy=RandomInteger[{1,10},{1000}];
moznestavy=Union[markstavy]
markdata={markcasy,markstavy}//Transpose;

Vybereme vsechna data, kdy je retezec ve stavu 1

(s1=Select[markdata,#[[2]]==1&])//Length

Z toho secteme prvni slozky (nejdrive vybereme prvni slozky)

Map[First,s1] (* ekvivalentne *)
Transpose[s1][[1]]

Total[Map[First,s1]]

Short[(ss=Map[Function[x,Select[markdata,#[[2]]==x&]],moznestavy]),10]

Total[Map[First,ss[[1]]]]

Map[{#,Total[Map[First,ss[[#]]]]}&,moznestavy]//TableForm

Reseni uloh celocoselneho linearniho programovani

FindMinimum[{x+y,x+2y≥3&&x≥0&&y≥0&&y∈Integers},{x,y}]

Prohlednete si funkci FindInstance.

Viz tez tutorial/ConstrainedOptimizationExact #683534434.

Manualy k versi 7.0

http://url.wolfram.com/6eVW8E

FAQ

Mnohorozmerne normalni rozdeleni

Needs["MultivariateStatistics`"];
pdf[xx_,μμ_,Σ_]:=PDF[MultinormalDistribution[μμ,Σ],xx]
pdf[Array[x,5],Array[μ,5],Array[σ,{5,5}]+Transpose[Array[σ,{5,5}]]]

Nacitani dat z Excelu

$ImportFormats

SetDirectory[NotebookDirectory[]]

Import["pokusnadataExcel.xls"]//TableForm

({{({{"a"}, {"b"}, {"a"}, {"b"}, {"a"}, {"b"}}), ({{"b"}, {1.`}, {2.`}, {3.`}, {4.`}, {5.`}}), ({{"a"}, {2.`}, {3.`}, {4.`}, {5.`}, {6.`}}), ({{"b"}, {3.`}, {4.`}, {5.`}, {6.`}, {7.`}}), ({{"a"}, {4.`}, {5.`}, {6.`}, {7.`}, {8.`}}), ({{"b"}, {5.`}, {6.`}, {7.`}, {8.`}, {9.`}}), ({{"a"}, {6.`}, {7.`}, {8.`}, {9.`}, {10.`}}), ({{"b"}, {7.`}, {8.`}, {9.`}, {10.`}, {11.`}}), ({{"a"}, {8.`}, {9.`}, {10.`}, {11.`}, {12.`}}), ({{"b"}, {9.`}, {10.`}, {11.`}, {12.`}, {13.`}}), ({{"a"}, {10.`}, {11.`}, {12.`}, {13.`}, {14.`}}), ({{"b"}, {11.`}, {12.`}, {13.`}, {14.`}, {15.`}})}})

FAQ

Mnohorozmerne normalni rozdeleni

Needs["MultivariateStatistics`"];
pdf[xx_,μμ_,Σ_]:=PDF[MultinormalDistribution[μμ,Σ],xx]
pdf[Array[x,5],Array[μ,5],Array[σ,{5,5}]+Transpose[Array[σ,{5,5}]]]

Nacitani dat z Excelu

$ImportFormats

SetDirectory[NotebookDirectory[]]

Import["pokusnadataExcel.xls"]//TableForm

({{({{"a"}, {"b"}, {"a"}, {"b"}, {"a"}, {"b"}}), ({{"b"}, {1.`}, {2.`}, {3.`}, {4.`}, {5.`}}), ({{"a"}, {2.`}, {3.`}, {4.`}, {5.`}, {6.`}}), ({{"b"}, {3.`}, {4.`}, {5.`}, {6.`}, {7.`}}), ({{"a"}, {4.`}, {5.`}, {6.`}, {7.`}, {8.`}}), ({{"b"}, {5.`}, {6.`}, {7.`}, {8.`}, {9.`}}), ({{"a"}, {6.`}, {7.`}, {8.`}, {9.`}, {10.`}}), ({{"b"}, {7.`}, {8.`}, {9.`}, {10.`}, {11.`}}), ({{"a"}, {8.`}, {9.`}, {10.`}, {11.`}, {12.`}}), ({{"b"}, {9.`}, {10.`}, {11.`}, {12.`}, {13.`}}), ({{"a"}, {10.`}, {11.`}, {12.`}, {13.`}, {14.`}}), ({{"b"}, {11.`}, {12.`}, {13.`}, {14.`}, {15.`}})}})

Petrikuv problem

tabDP=Table[NMaxValue[{Exp[-n(n-1)/(2m)]-(∏_(k=1)^n)(1-(k-1)/m),m≥n},m],{n,2,100}]

ListPlot[tabDP]

Ocividne spatne!!! Chyba ve funkci NMaxvalue.

(∏_(k=1)^n)(1-(k-1)/m)

Clear[m,n,approxDP];approxDP[m_,n_]:=Exp[-n(n-1.`30)/(2m)]-(((-(1/m)))^n) Pochhammer[-m,n]

t40=Table[{m,approxDP[m,40]},{m,40,1000}];
maxt40=Max[t40[[All,2]]]
Position[t40[[All,2]],maxt40]
ListPlot[t40]

Comments

Random matrices

ranmatposdef[n_]:= mat=RandomVariate[UniformDistribution[],{n,n}]

r30=ranmatposdef[30]//Eigenvalues

Table[{n,Max[Abs[ranmatposdef[n]//Eigenvalues]]},{n,30,200,10}]//TableForm

ranmat[n_]:= mat=RandomVariate[NormalDistribution[],{n,n}]

Table[{n,Max[Abs[ranmat[n]//Eigenvalues]]},{n,30,200,10}]//TableForm

$LicenseExpirationDate

Estimation

EstimatedDistribution

data=RandomVariate[GammaDistribution[1,2],(10^2)];
edist=EstimatedDistribution[data,GammaDistribution[α,β]]
Plot[{PDF[GammaDistribution[1,2],x],PDF[edist,x]},{x,0,10}]

EstimatedDistribution[data,GammaDistribution[α,β],ParameterEstimator->"MethodOfMoments"]

data=RandomVariate[BetaDistribution[2,3],1000];EstimatedDistribution[data,BetaDistribution[α,β]]EstimatedDistribution[data,BetaDistribution[α,β],ParameterEstimator->{"MaximumLikelihood",Method->"FindMaximum"}]

See 15 very interesting applications at the documantation of EstimatedDistribution.

LogLikelihood

LogLikelihood[NormalDistribution[μ,σ],Array[x,5]]

Empirical distribution

data=RandomVariate[NormalDistribution[],50];
𝒟emp=EmpiricalDistribution[data];
{DiscretePlot[PDF[𝒟emp,x],{x,data}],DiscretePlot[CDF[𝒟emp,x],{x,-4,4,.01}]}
{Moment[𝒟emp,2],Quantile[𝒟emp,0.95]}

Kernel estimator (jadrove odhady)

data=RandomVariate[NormalDistribution[],10^3];
𝒟=SmoothKernelDistribution[data];
Plot[PDF[𝒟,x],{x,-3,3}]

Choice of bandwidth:

𝒟1=SmoothKernelDistribution[data,.05];
𝒟2=SmoothKernelDistribution[data,1.];
Table[Plot[PDF[d,x],{x,-3,3}],{d,{𝒟1,𝒟2}}]

Best strategy: use Manipulate

Manipulate[Plot[Evaluate@PDF[SmoothKernelDistribution[data,bw],x],{x,-4,4},Filling->Axis,PlotLabel->{"band width",bw},Ticks->None],{bw,Table[0.025 j,{j,20}]},SaveDefinitions->True]

Choice of the kernel

Table[Plot[Evaluate@PDF[SmoothKernelDistribution[data,Automatic,i],x],{x,-4,4},Filling->Axis,PlotLabel->i,Ticks->None],{i,{"Biweight","Cosine","Epanechnikov","Gaussian","Rectangular","SemiCircle","Triangular","Triweight"}}]

Testing

DistributionFitTest (test-of-goodness of fit, test dobre shody)

Test for normality

data=RandomVariate[NormalDistribution[],2500];
DistributionFitTest[data]

ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]

Comparison of data and fitted

Show[Histogram[data,Automatic,"ProbabilityDensity"],Plot[PDF[ℋ["FittedDistribution"],x],{x,-5,5}]]

Location test

LocationTest

Variance test

VarianceTest