VPFPM ZS2019/20 - cvičení

NMFM308 Výpočetní prostředky finanční a pojistné matematiky

NMFM308 Computational Tools for Finance and Insurance

$LicenseExpirationDate

Externí zdroje

http://www.karlin.mff.cuni.cz/~pesta/NMFM301/rozdeleni.pdf

http://www.karlin.mff.cuni.cz/~pesta/NMFM301/statistika_fm.pdf

Poznamka. Asymptoticky rozptyl V () vyberoveho kvantilu se spatne odhaduje,
protoze nemame k disposici univerzalne pouzitelny a spolehlivy odhad hustoty.

https://reference.wolfram.com/language/guide/CombinatorialFunctions.html

https://www.wolfram.com/language/fast-introduction-for-math-students/en/probability/

https://reference.wolfram.com/language/guide/ProbabilityAndStatistics.html

https://www.wolfram.com/language/fast-introduction-for-math-students/en/statistics/

2019-10-10 Cvičení 2

FAQ

Pravděpodobnostní rozdělení 1

Binomické rozdělení

?*Distribution

Clear[n,p,bin]
bin[n_,p_]:=BinomialDistribution[n,p]

Row@{"Prumer = ",Mean[bin[n,p]],"   Rozptyl = ",Variance[bin[n,p]]}

(* hustota PDF probability density function *)
PDF[bin[n,p],j]

(* distribucni funkce CDF cumulative distribution function *)
CDF[bin[n,p],j]

DiscretePlot[Table[PDF[BinomialDistribution[40,p],k],{p,{0.1,0.5,0.7}}]//Evaluate,{k,0,36},PlotRange->All,PlotMarkers->Automatic]

Normální rozdělení

Stejné charakteristiky, navíc

Skewness (* šikmost *) Kurtosis (* špičatost *)

(χ^2)- rozdělení a jeho grafické porovnání s normálním

Clear[chi]
chi[n_]:=ChiSquareDistribution[n]

Through[{f,g,h}[x]]

charakteristiky=Through[(momenty={Mean,Variance,StandardDeviation,Skewness,Kurtosis})[chi[n]]]
(momentyString=Map[ToString,momenty])//FullForm
Text@Row@{"(χ^2)-rozdělení o n stupních volnosti"}
Grid[{momentyString,charakteristiky}//Transpose,Frame->All]

Clear[nd]
nd[n_]=NormalDistribution[n,(Sqrt[2n])]
p2=Plot[PDF[nd[n],x]/.n->7,{x,Quantile[nd[n],0.03],Quantile[nd[n],0.97]}/.n->7,PlotStyle->{Thick,Red}];
p1=Plot[PDF[chi[n],x]/.n->7,{x,Quantile[nd[n],0.03],Quantile[nd[n],0.97]}/.n->7,PlotRange->All];
GraphicsRow@{Show[p2,p1],Show[p1,p2]}

With[{n=7},
Clear[chi,nd];
chi[n_]:=ChiSquareDistribution[n];
nd[n_]=NormalDistribution[n,(Sqrt[2n])];
p2=Plot[PDF[nd[n],x],{x,Quantile[nd[n],0.03],Quantile[nd[n],0.97]},PlotStyle->{Thick,Red}];
p1=Plot[PDF[chi[n],x],{x,Quantile[nd[n],0.03],Quantile[nd[n],0.97]},PlotRange->All,PlotLabel->"Hustoty normálního \na (χ^2)-rozdělení o "<>ToString[n]<>" stupních volnosti"];
GraphicsRow@{Show[p2,p1],Show[p1,p2]}]

Manipulate[p2=Plot[PDF[nd[n],x],{x,Quantile[nd[n],0.03],Quantile[nd[n],0.97]},PlotStyle->{Thick,Red}];
p1=Plot[PDF[chi[n],x],{x,Quantile[nd[n],0.03],Quantile[nd[n],0.97]},PlotRange->All,AxesOrigin->{Automatic,0},PlotLabel->"Hustoty normálního \na (χ^2)-rozdělení o "<>ToString[n]<>" stupních volnosti"];
Show[p1,p2],{n,7,351,2,Appearance->"Labeled"},SaveDefinitions->True]

2019-10-17 Cvičení 3

FAQ

Pravděpodobnostní rozdělení 2

Metody Monte Carlo (Simulační metody): Princip

MC metody: odvětví experimentální matematiky, které se zabývá experimenty s náhodnými čísly. Numerická technika provádění hromadných experimentů s modely pomocí počítače umožňující generovat velké množství náhodných výběrů (umělých náhodných výběrů). Většinou jsou efektivnější v situacích, kdy sama podstata problému je náhodná. Při řešení deterministických problémů tomu jsou většinou úspěšnější klasické metody numerické matematiky.

Probability

https://reference.wolfram.com/language/guide/ProbabilityAndStatistics.html

Probability
NProbability
ProbabilityDistribution
ProbabilityPlot

Type Esc dist Esc for the \[Distributed] symbol.

Probability[x==1,x\[Distributed]BinomialDistribution[1,1/2]]

Expectation[2 x+3,x\[Distributed]NormalDistribution[]]

?*Probability*

Domácí cvičení : Zjistěte, k čemu je dobrá funkce ProbabilityPlot

Základní funkce:

Spočtěte pst jevu, že X>x, x∈ℝ, X\[Distributed]Geo(p).

Probability[X>x,X\[Distributed]GeometricDistribution[p]]

Spočtěte distribuční funkci (X^2)+(Y^2), kde X,Y iid N(0,1).

Clear[x]
Probability[(X^2)+(Y^2)<x,{X\[Distributed]NormalDistribution[],Y\[Distributed]NormalDistribution[]}]

TransformedDistribution[(X^2)+(Y^2),{X\[Distributed]NormalDistribution[],Y\[Distributed]NormalDistribution[]}]
CDF[%,x]

Clear[x,n]
CDF[ChiSquareDistribution[n],x]
Table[{n,%},{n,4}]//Simplify//TableForm

{PDF[GammaDistribution[(n/2),2],x],PDF[ChiSquareDistribution[n],x]}
{({{\[Piecewise], ({{(((2^(-n/2)) (E^(-x/2)) (x^(-1+(n/2))))/(Gamma[(n/2)])), x>0}, {0, True}})}}),({{\[Piecewise], ({{(((2^(-n/2)) (E^(-x/2)) (x^(-1+(n/2))))/(Gamma[(n/2)])), x>0}, {0, True}})}})}

Clear[x]
Probability[(X^2)+(Y^2)<x,{X,Y}\[Distributed]ProductDistribution[{NormalDistribution[],2}]]

Conditional Probability (Podmíněná pst)

A, B náhodné jevy, n. v. X\[Distributed]distribution. Pak P(X∈A|X∈B) je

Probability[X∈A\[Conditioned]X∈B,X\[Distributed]distribution]

Znak \[Conditioned]: EsccondEsc nebo z palety.

Exponenciální rozdělení

Spočtěte pst jevu x<X<x+h za podmínky X>x, kde h>0.

Assuming[h>0∧x>0,Probability[x<X<x+h\[Conditioned]X>x,X\[Distributed]ExponentialDistribution[λ]]]

Nemá paměť.

Geometrické rozdělení

Spočtěte pst jevu x<X<x+h za podmínky X>x, kde h>0.

Clear[pst]
pst[x_,h_]=Probability[x<X<x+h\[Conditioned]X>x,X\[Distributed]GeometricDistribution[p]]

Assuming[k>0&&k∈Integers&&x>0&&x∈Integers,pst[k,h]]//Simplify

pst[4,h]//Simplify

Také nemá paměť.

Stejná data narození (birthday problem)

Najděte pravděpodobnost P(n,r), že aspoň 2 lidé ze skupiny r lidí přicházejících do místnosti se narodili ve stejný den. Pro jednoduchost se předpokládá n=365  dnů v roce. Najděte r takové, že  P(n,r)>(1/2).

Stejná data narození (Řešení)

Nejprve spočítáme pravděpodobnost, že všichni příchozí se narodí v různých dnech. Ten, co přichází první se mohl narodit v libovolném z n dnů, druhý v jednom ze zbývajích n-1 dnů, atd. až ten r-tý ve zbývajících n-r+1 dnech. Je tedy n(n-1)...(n-r+1) možností. Všech možností je (n^r). Pak

P(n,r)=1-((n!)/((n-r)!))(1/(n^r)).

Clear[p,n,r];
p[n_,r_]:=1-((n!)/((n^r)(n-r)!))
Table[{r,p[365.,r]},{r,15,25}]//TableForm
p[365,23]//N

Stejná data narození  (kde hledat)

https://reference.wolfram.com/language/guide/ProbabilityAndStatistics.html

Probability
ProbabilityDistribution
ProbabilityPlot

Type ESC dist ESC for the \[Distributed] symbol.

Probability[x==1,x\[Distributed]BinomialDistribution[1,1/2]]

Expectation[2 x+3,x\[Distributed]NormalDistribution[]]

2019-10-24 Cvičení 4

FAQ

Poznámka k vyhledávání funkcí

Z minulého cvičení: Všechna zabudovaná rozdělení mají tvar *Distribution, takže jedna z možností je:

?Geo*Distribution

*Vyjádření distribuční funkce N(0,1) pomocí Φ(x)

CDF[NormalDistribution[],x]
%//FunctionExpand

(1/2) (1+Erf[(x/(Sqrt[2]))])==Φ[x]

⟹  Erfc[y]=2Φ[-(Sqrt[2 y])]

(* Lokální pravidlo (substituce) *)ruleNormal=Erfc[y_]-> 2Φ[-(Sqrt[2])y]

{CDF[NormalDistribution[],x]/.ruleNormal,CDF[NormalDistribution[μ,σ],x]/.ruleNormal}//Simplify//TraditionalForm
(* Místo dvojího použití /. ((ReplaceAll) použijeme //. (ReplaceRepeated - používá lokální pravidlo tak dlouho, dokud se daný výraz nemění): *)
{CDF[NormalDistribution[],x],CDF[NormalDistribution[μ,σ],x]}//.ruleNormal//Simplify//TraditionalForm

%//TeXForm

Pravděpodobnostní rozdělení 3

Vlastní rozdělení

ProbabilityDistribution

Syntax: ProbabilityDistribution[hustota, nosič]; místo hustoty může být třeba distr. fce, nutné ale vyznačit. Použití je pak stejné jako u zabudovaných rozdělení.

N(0,1)

Clear[standardNormal]
standardNormal=ProbabilityDistribution[(1/(Sqrt[2π]))(E^(-(x^2)/2)),{x,-∞,∞}]
charakteristiky={Mean,Variance,Skewness,Kurtosis};
Row@{"Charakteristiky:   ",Through[charakteristiky[standardNormal]]}

CDF[standardNormal,x]
(* ** Pro statistika běžné vyjádření pomocí distr. fce N(0,1) *)
ruleNormalErfc={Erfc[y_]-> 2Φ[-(Sqrt[2])y]};
ruleNormalErf={Erf[y_]-> 2Φ[(Sqrt[2])y]-1};
(%%%/.{ruleNormalErfc,ruleNormalErf})//Simplify

N(μ,(σ^2))

Clear[normalDistribution]
normalDistribution[μ_,σ_]:=ProbabilityDistribution[(1/(σ(Sqrt[2π])))Exp[- ((((x-μ))^2)/(2(σ^2)))],{x,-∞,∞},Assumptions->σ>0]

Integrate[PDF[normalDistribution[μ,σ],x],{x,-∞,∞},Assumptions->σ>0]

Row@{"Charakteristiky:   ",Through[charakteristiky[normalDistribution[μ,σ]]]}

Diskrétní rozdělení

(* Diskrétní analogie Laplaceova (= dvojitě exponenciálního rozdělení) *)Table[{k,Exp[-Abs[k]]},{k,-5,5}]//ListLinePlot

(* Normalizační konstanta *)
constDiscreteLaplace=1/Total[Table[Exp[-Abs[k]],{k,-5,5}]]//Simplify

discreteLaplaceDistribution=
ProbabilityDistribution[constDiscreteLaplace Exp[-Abs[k]],{k,-5,5}]

Row@{"Charakteristiky:   ",Through[charakteristiky[discreteLaplaceDistribution]]}//Simplify
%//N

EmpiricalDistribution

Charakterizováno empirickou distr. fcí náhodného výběru: (Subscript[F, n])(x)=(1/n)(∑_(i=1)^n)I{(Subscript[X, i])≤x}. Vrací objekt DataDistribution.

𝒟 se dostane jako EscscDEsc

(* Random sample of size 50 *)
SeedRandom[13]
randomSampleStandardNormal50=RandomVariate[NormalDistribution[],50];
𝒟50=EmpiricalDistribution[randomSampleStandardNormal50]

Through[charakteristiky[𝒟50]]
{DiscretePlot[PDF[𝒟50,x],{x,randomSampleStandardNormal50}],DiscretePlot[CDF[𝒟50,x],{x,-4,4,.01}]}
Histogram[randomSampleStandardNormal50]

(* Random sample of size 5000 *)
SeedRandom[13]
randomSampleStandardNormal5000=RandomVariate[NormalDistribution[],5000];
𝒟5000=EmpiricalDistribution[randomSampleStandardNormal5000]

1./(Sqrt[5000])

Through[charakteristiky[𝒟5000]]
GraphicsRow@{DiscretePlot[PDF[𝒟5000,x],{x,randomSampleStandardNormal5000}],DiscretePlot[CDF[𝒟5000,x],{x,-4,4,.01}]}
Histogram[randomSampleStandardNormal5000]

Charakteristiky rozdělení

Expectation, Moment, CentralMoment
NExpectation, NProbability
Quantile

{Expectation[x^2+7x+8,x\[Distributed]PoissonDistribution[μ]],
NExpectation[x^2+7x+8,x\[Distributed]PoissonDistribution[10]],
%%/.μ->10}

{Moment[PoissonDistribution[μ],r],CentralMoment[PoissonDistribution[μ],r],CentralMoment[PoissonDistribution[μ],13]}

Úloha 1

Definujte rozdělení n. v. X jejíž hustota je

const (1/(1+(x^4))),   x∈ ℝ.

Graficky znázorněte hustotu. Spočtěte obecné momenty, E,X,var X,(Subscript[γ, 1]),(Subscript[γ, 2]) (tj. šikmost a špičatost).

Řešení 1

Výpočet konstanty:

const=1/(∫_(-∞)^∞)(1/(1+(x^4)))\[DifferentialD]x
Row@{"Hustota:   ",hustota4[x_]=const (1/(1+(x^4)))}

distribution4=ProbabilityDistribution[hustota4[x],{x,-∞,∞}]

Možné je i použití volby Normalize:

distribution4UsingNormalize=ProbabilityDistribution[(1/(1+(x^4))),{x,-∞,∞},Method->"Normalize"]
PDF[%,x]

Plot[PDF[distribution4,x],{x,-4,4},PlotLabel->"Density of distribution4",PlotStyle->Thick]

Clear[r]
{Moment[distribution4,r],Moment[distribution4,3],
Moment[distribution4,2]}

Through[charakteristiky[distribution4]]

Úloha 2

Definujte rozdělení diskrétní n. v. X s nosičem {1,2,...}, kde

P(X=k)=const (1/(k^2)),   k=1,2,...

Graficky znázorněte hustotu vzhledem k čítací míře. Spočtěte obecné momenty, E,X,var X,(Subscript[γ, 1]),(Subscript[γ, 2]) (tj. šikmost a špičatost).

Řešení 2

const2=1/Sum[(1/(k^2)),{k,1,∞}]

Clear[distribution2]
distribution2=ProbabilityDistribution[(1/(k^2)), {k,1, ∞,1},Method->"Normalize"]
Sum[PDF[distribution2,k],{k,1,∞}]

DiscretePlot[PDF[distribution2,k],{k,1,15},PlotRange->All]

Through[charakteristiky[distribution2]]

StudentTDistribution

2019-10-31 Cvičení 5

FAQ

Náhodný výběr

Náhodný výběr (random sample) (Subscript[X, 1]),...,(Subscript[X, n]) jsou iid, (kopie referenční n. v. X). Ve Wolfram language většinou ve formě seznamu {(Subscript[X, 1]),...(Subscript[X, n])}.

Náhodný výběr z množiny prvků

SeedRandom[13]
seznam1=RandomReal[{1,5},20]
%//Union
Histogram[seznam1,5,LabelingFunction->Above]

SeedRandom[13]
seznam2=RandomInteger[{1,5},20]
%//Union
Histogram[seznam2,LabelingFunction->Above]

S navracením (s opakováním): funkce RandomChoice

Odpovídá tomu, co si my představujeme pod náhodným výběrem.

SeedRandom[13]
randomChoice1=RandomChoice[seznam1,20]
%//Union
Histogram[randomChoice1,5,LabelingFunction->Above]

SeedRandom[13]
randomChoice11000=RandomChoice[seznam1,1000];
%//Union//Length
Histogram[randomChoice11000,5,LabelingFunction->Above]

SeedRandom[13]
randomChoice2=RandomChoice[seznam2,20]
%//Union
Histogram[randomChoice2,LabelingFunction->Above]

SeedRandom[13]
randomChoice21000=RandomChoice[seznam2,1000];
%//Union
Histogram[randomChoice21000,LabelingFunction->Above]

Bez navracení (bez opakování): funkce RandomSample

Pro statistiky matoucí: neodpovídá tomu, co si my představujeme pod náhodným výběrem.

V případě RandomSample[množina,Length[množina]] dostaneme náhodnou permutaci.

SeedRandom[13]
randomSample1=RandomSample[seznam1,20]
%//Union//Length

RandomSample[seznam1,21]

seznam2
SeedRandom[13]
RandomSample[seznam2,20]

Tally[seznam2]//Sort
SeedRandom[13]
Tally[RandomSample[seznam2,20]]//Sort

SmoothHistogram (Vyhlazený histogram) a příbuzné funkce

Motivace: Histogram je možné jednoduše transformovat na hustotu (vzhledem k Lebesgueově míře), ale sloupkový diagram není pro spojitá rozdělení to nejlepší vizuální znázornění. Odpomoc spočívá ve vyhlazení (smoothing). Nejjednodušší bez potřeby dalších znalostí je Vyhlazený histogram.

?*Smooth*

Příklad

SeedRandom[13]
randomSample=RandomVariate[ExponentialDistribution[1/2],5000];
p1=SmoothHistogram[randomSample,PlotRange->All]
p2=Histogram[randomSample,Automatic,"PDF"];
Show[p1,p2,PlotRange->All]

Vyhlazené jádrové rozdělení (SmoothKernelDistribution)

(Subscript[X, 1]),...,(Subscript[X, n ]) náhodný výběr, k symetrická okolo nuly hustota zvaná jádro (kernel). Úkolem je najít neparametrický odhad hustoty f (neznáme její analytický tvar).

Jádrový odhad je  (OverHat[f])(x)=(1/(n h))(∑_(i=1)^n)k(((x-(Subscript[X, i]))/h)) pro vyhlazovací jádro k(x) a šířku okénka (bandwidth) h.

The probability density function for (RefLink[SmoothKernelDistribution, paclet:ref/SmoothKernelDistribution]) for a value x is given by a linearly interpolated version of (1/(n h))(∑_(i=1)^n)k(((x-(Subscript[x, i]))/h)) for a smoothing kernel k(x) and bandwidth parameter h.

Odpovídající zabudovaná funkce je:

SmoothKernelDistribution

(*Jádrový odhad rozdělení*)
SeedRandom[13];
data=RandomReal[NormalDistribution[10,3 ],500000];
distribuce=SmoothKernelDistribution[data]
GraphicsRow[{Plot[PDF[distribuce,x],{x,-5,20},PlotLabel->"Hustota"],
Plot[CDF[distribuce,x],{x,-5,20},PlotLabel->"Distr. fce"]}]

Apple

SetDirectory[NotebookDirectory[]]
cenyApple=Get["aapl20191030.save"]
%//DateListPlot

TimeSeries je objekt. Z něj je možné získat následující data (mimochodem totéž platí i pro ostatní objekty podobného typu) pomocí  “Properties”:

cenyApple["Properties"]

(* Ceny Apple *)
cenyAppleHodnoty=cenyApple["Values"]
ListLinePlot[%,PlotLabel->"Ceny Apple"]

cenyAppleHodnoty//Histogram
cenyAppleHodnoty//SmoothHistogram

***VW

(* Výnosy Apple *)returnsVW=FinancialData["F:VOW","Return",{{2014,9,22},{2015,(* 8 *)11,30}},"Value"];
grafyJH[returnsVW]

smooothHistogramVW=SmoothHistogram[returnsVW]
Text@Row@{"Šířka okénka:   ",SmoothKernelDistribution[returnsVW]["Bandwidth"]}

Manipulate[SmoothHistogram[returnsVW,bw,"PDF", Filling->Axis,PlotLabel->"Odhadnutá hustota výnosů VW",PlotRange->All],{{bw,.005},0.0000001,0.1,Appearance->"Labeled"},SaveDefinitions->True]

***Apple Details

sp50020190826=FinancialData["SP500",{{2009,1,1},{2019,8,26}},Method->"Legacy"]
Length@%

aapl20191030=FinancialData["AAPL",{{2009,1,1},{2019,8,26}},Method->"Legacy"]
Length@%

Save["aapl20191030.save",aapl20191030]

Export["aapl20191030.xls",Get["aapl20191030.save"]["Path"]]

Export["aapl20191030.xlsx",aapl20191030]

aapl20191030["Properties"]
pathApple20191030=aapl20191030["DatePath"]
Export["pathApple20191030.xls",pathApple20191030]

SetDirectory[NotebookDirectory[]]
Export["cenyApple20191030.XLSX",aapl20191030]

SmoothHistogram[hodnotyAAPL]

Text@Row@{"Šířka okénka:   ",SmoothKernelDistribution[hodnotyAAPL]["Bandwidth"]}

Manipulate[SmoothHistogram[hodnotyAAPL,bw,"PDF", Filling->Axis,PlotLabel->"Odhadnutá hustota výnosů Apple",PlotRange->All],{{bw,.005},0.0000001,0.1,Appearance->"Labeled"},SaveDefinitions->True]

**20191030_DJI_Yahoo_Finance

https://finance.yahoo.com/quote/%5EDJI/history?p=%5 EDJI

https://finance.yahoo.com/quote/%5EDJI/history?period1=1230764400&period2=1572390000&interval=1d&filter=history&frequency=1d

DJIYahooFinanceHead20191030="Date Open High Low Close* Adj Close** Volume"

SetDirectory[NotebookDirectory[]]
Import["20191030_DJI_Yahoo_Finance.txt","Table"]

(recordFirst={"Oct","29,",2019,27061.07`,27165.94`,27039.76`,27071.46`,27071.46`,261420000})//FullForm
StringReplace[recordFirst[[2]],","->""]//FullForm

recordFirst[[2]]=StringReplace[recordFirst[[2]],","->""];
recordFirst[[3]]=ToString[recordFirst[[3]]];
recordFirst//FullForm

recordFirst[[1;;3]]//FullForm

Riffle[recordFirst[[1;;3]],"/"]
%//StringJoin
DateList[%][[1;;3]]

Pravděpodobnostní rozdělení 4

Úloha 3

Definujte funkci firstFourMoments  se vstupním argumentem distribution (rozdělení), která bude mít jako výstup tabulku ve tvaru

Vyzkoušejte pro N(μ,(σ^2)) a Studentovo rozdělení s parametry μ (parametr polohy, location) a σ (parametr měřítka, scale) a 4.1 stupni volnosti.

Řešení 3.1

ClearAll[firstFourMoments]
firstFourMoments[distribution_]:=Module[{charakteristiky={Mean,Variance,Skewness,Kurtosis},dist,aux},
dist=distribution;
aux={charakteristiky,Through[charakteristiky[dist]]}//Transpose;
Grid[Prepend[aux,{"Distribution",dist}],Frame->All]
	]

firstFourMoments[NormalDistribution[μ,σ]]
firstFourMoments[StudentTDistribution[μ,σ,4.1]]
firstFourMoments[StudentTDistribution[k,μ,σ]]

Řešení 3.2

charakteristiky={Mean,Variance,Skewness,Kurtosis};Clear[firstFourMoments2]
firstFourMoments2[distribution_]:=
Grid[Prepend[{charakteristiky,Through[charakteristiky[distribution]]}//Transpose,{"Distribution",distribution}],Frame->All]//TraditionalForm

firstFourMoments2[NormalDistribution[μ,σ]]
firstFourMoments2[StudentTDistribution[μ,σ,4.1]]
firstFourMoments2[StudentTDistribution[k,μ,σ]]

Úloha 4

Z pracovního adresáře importujte soubory vwceny.xls a vwvynosy.xls s cenami společnosti VolskWagen. Spočtěte charakteristiky jako v úloze 3 a pro oba soubory sestrojte histogram a vyhlazený histogram.

Řešení 4

SetDirectory[NotebookDirectory[]]
FileNames[]

vwCeny=Import["vwceny.xls"][[1]];
vwCeny2=vwCeny[[All,2]];
rozdelenivwCeny=EmpiricalDistribution[vwCeny2];
firstFourMoments[rozdelenivwCeny]
GraphicsRow@{Histogram[vwCeny2],SmoothHistogram[vwCeny2]}

vwVynosy=Import["vwvynosy.xls"][[1]];
vwVynosy2=vwVynosy[[All,2]];
rozdelenivwVynosy=EmpiricalDistribution[vwVynosy2];
firstFourMoments[rozdelenivwVynosy]
GraphicsRow@{Histogram[vwVynosy2],SmoothHistogram[vwVynosy2,PlotRange->All]}

2019-11-07 Cvičení 6

FAQ

Pravděpodobnostní rozdělení 5

Test

Notes (tj. tento notebook s níže uvedeným názvem) s řešením zašlete e-mailem na mou adresu hurt@karlin.mff.cuni.cz.
Název notesu: 20191107_Jmeno_Prijmeni
Předmět zprávy: 20191107_Jmeno_Prijmeni
Nepoužívejte diakritiku.
Termín: 2019-11-07 12:30

Úloha 5

Pomocí funkce Manipulate demonstrujte vliv šířky okénka (bandwidth u funkce SmoothHistogram) na výsledný tvar hustoty pro rozdělení výnosů z dat cen společnosti Apple (soubor aapl20191105ceny.xlsx v pracovním adresáři).

Řešení 5

Viz také definici jádrového odhadu.

SetDirectory[NotebookDirectory[]];appleCeny=Import["aapl20191105ceny.xlsx"];
%//Dimensions
(* 1. prvek tohoto seznamu je matice, ceny jsou ve sloupci 2 *)
aapl20191105cenyWithoutDates=appleCeny[[1,All,2]];
(* Výnosy: *)
aapl20191105vynosyWithoutDates=With[{a=aapl20191105cenyWithoutDates},((Rest[a]-Most[a])/(Most[a]))];

Náhodně zvolíme šířku okénka (2. argument funkce SmoothHistogram). Je-li funkce příliš hladká, šířku snížíme, je-li příliš “divoká”, šířku zvýšíme.

SmoothHistogram[aapl20191105vynosyWithoutDates,.001,PlotRange->All]

Oproti zadání ještě připojíme graf hustoty normálního rozdělení s parametry μ= Mean[aapl20191105vynosyWithoutDates] a σ= StandardDeviation[aapl20191105vynosyWithoutDates].

Manipulate[Show[SmoothHistogram[aapl20191105vynosyWithoutDates,bw,PlotRange->All],Plot[PDF[NormalDistribution[Mean[aapl20191105vynosyWithoutDates],StandardDeviation[aapl20191105vynosyWithoutDates]],x],{x,-0.1,0.1},PlotStyle->{Red,Thick}]],{{bw,0.001},0.0001,0.01,Appearance-> "Labeled"},SaveDefinitions->True]

Úloha 6
Pravděpodobnosti chvostů rozdělení

Jedná se vlastně o demonstraci pravidel 2σ a 3σ. Nechť X je n. v. s rozdělením distr, E,X=μ, var X=(σ^2). Definujte funkci probabilityTales s jedním argumentem distr (rozdělení), která vytvoří tabulku typu (SD směrodatná odchylka, Standard deviation):

Ilustrujte  pro rozdělení normální, Studentovo s parametry μ, σ s počtem stupňů volnosti 4.1 a pro rozdělení výnosů společnosti VW (soubor vwvynosy.xls v pracovním adresáři).

Řešení 6

Pozor, v pravidlu 2σ resp. 3σ figuruje směrodatná odchylka n. v. X. Parametr σ v Laplaceově, Studentově a i jiných rozděleních nemá význam směrodatné odchylky!!!
Abychom dostali numerické hodnoty, přidáme desetinnou tečku třeba u výrazu 2σ ⟶ 2.σ

ClearAll[probabilityTales]
probabilityTales[distribution_]:=Module[{mean=Mean[distribution],sd=StandardDeviation[distribution],γ1=Skewness[distribution],γ2=Kurtosis[distribution],prob2σ,prob3σ},
Grid[{{"Distribution",distribution},
		{"Mean",mean},{"SD",sd},{"Skewness",γ1},{"Kurtosis",γ2},{"P(|X-μ|>2σ)",Probability[Abs[x-mean]>2.sd,x\[Distributed]distribution]},{"P(|X-μ|>3σ)",Probability[Abs[x-mean]>3.sd,x\[Distributed]distribution]}},Frame->All]
]

probabilityTales[LaplaceDistribution[μ,σ]]

probabilityTales[StudentTDistribution[μ,σ,4.1]]

SetDirectory[NotebookDirectory[]];
vwVynosy=Import["vwvynosy.xls"];
%//Dimensions
vwVynosy=vwVynosy[[1,All,2]];
vwVynosDistribution=EmpiricalDistribution[vwVynosy]
probabilityTales[vwVynosDistribution]

probabilityTales[StudentTDistribution[μ,σ,4.1]/.{μ->.0003969,σ->.01945}]

Úloha 7
Rozdělení výběrového průměru pro různá rozdělení

Je-li (OverBar[X]) výběrový průměr náhodného výběru z libovolného rozdělení s konečným rozptylem (σ^2) a střední hodnotou μ, pak E,(OverBar[X])=μ a var (OverBar[X])=(σ^2)/n. Platí CLV.

Úloha: Generovat N náhodných výběrů o rozsahu n, z každého výběru spočítat výběrový průměr. Dostaneme tedy N hodnot, vlastně náhodný výběr výběrových průměrů. Graficky znázornit empirickou hustotu (histogram), “vyhlazený” histogram (SmoothHistogram) a do stejného grafu odpovídající hustotu normálního rozdělení, pozor na stejné parametry! Použijte rozdělení z předchozího odstavce. Efektivní postup spočívá ve vytvoření obecné procedury (Module), kde vstup je N, n, rozdělení  a volitelné semínko pro generování náhodných čísel (Optional, ), celkem 4 vstupní parametry.

Vyzkoušejte na 3 rozděleních N(4,4), exponenciálním s parametrem λ=5 a rozdělení výnosů společnostiApple.

Očekávaný výstup:

rozdeleniPrumeru[10000,5,ExponentialDistribution[3]]

Řešení 7

Obecná procedura:

Clear[rozdeleniPrumeru]
rozdeleniPrumeru[pocetVyberuN_,rozsahVyberun_,rozdeleni_,seminko_:13]:=
Module[{n=rozsahVyberun,nahodneVybery,xx,histogram,smoothHistogram,normalHustota,mean,sd},
SeedRandom[seminko];
mean=Mean[rozdeleni];sd=StandardDeviation[rozdeleni];
nahodneVybery=Table[xx=RandomVariate[rozdeleni,{n}];Mean[xx],{pocetVyberuN}];
histogram=Histogram[nahodneVybery,Automatic,"PDF",PlotLabel->rozdeleni];
smoothHistogram=SmoothHistogram[nahodneVybery,PlotLabel->rozdeleni,PlotStyle->Thick];
normalHustota=Plot[PDF[NormalDistribution[mean,sd/(Sqrt[n])],x],{x,mean-4sd/(Sqrt[n]),mean+4sd/(Sqrt[n])},PlotStyle->{Thick,Red}];
Show[{histogram,smoothHistogram,normalHustota}]
]

rozdeleniPrumeru[10000,5,ExponentialDistribution[5]]

rozdeleniPrumeru[10000,4,NormalDistribution[4,2]]

GraphicsRow@{rozdeleniPrumeru[10000,10,StudentTDistribution[3]],rozdeleniPrumeru[10000,10,StudentTDistribution[30]]}

distr=EmpiricalDistribution[aapl20191105vynosyWithoutDates];rozdeleniPrumeru[2000,10,distr]

2019-11-14 Cvičení 7

FAQ

Test 2019-11-07

Notes (tj. tento notebook s níže uvedeným názvem) s řešením zašlete e-mailem na mou adresu hurt@karlin.mff.cuni.cz.
Název notesu: 20191107_Jmeno_Prijmeni
Předmět zprávy: 20191107_Jmeno_Prijmeni
Nepoužívejte diakritiku.
Termín: 2019-11-07 12:30

Úloha 5

Pomocí funkce Manipulate demonstrujte vliv šířky okénka (bandwidth u funkce SmoothHistogram) na výsledný tvar hustoty pro rozdělení výnosů z dat cen společnosti Apple (soubor aapl20191105ceny.xlsx v pracovním adresáři).

Řešení 5

Viz také definici jádrového odhadu.

SetDirectory[NotebookDirectory[]];appleCeny=Import["aapl20191105ceny.xlsx"];
%//Dimensions
(* 1. prvek tohoto seznamu je matice, ceny jsou ve sloupci 2 *)
aapl20191105cenyWithoutDates=appleCeny[[1,All,2]];
(* Výnosy: *)
aapl20191105vynosyWithoutDates=With[{a=aapl20191105cenyWithoutDates},((Rest[a]-Most[a])/(Most[a]))];

Náhodně zvolíme šířku okénka (2. argument funkce SmoothHistogram). Je-li funkce příliš hladká, šířku snížíme, je-li příliš “divoká”, šířku zvýšíme.

SmoothHistogram[aapl20191105vynosyWithoutDates,.001,PlotRange->All]

Oproti zadání ještě připojíme graf hustoty normálního rozdělení s parametry μ= Mean[aapl20191105vynosyWithoutDates] a σ= StandardDeviation[aapl20191105vynosyWithoutDates].

Manipulate[Show[SmoothHistogram[aapl20191105vynosyWithoutDates,bw,PlotRange->All],Plot[PDF[NormalDistribution[Mean[aapl20191105vynosyWithoutDates],StandardDeviation[aapl20191105vynosyWithoutDates]],x],{x,-0.1,0.1},PlotStyle->{Red,Thick}]],{{bw,0.001},0.0001,0.01,Appearance-> "Labeled"},SaveDefinitions->True]

Úloha 6
Pravděpodobnosti chvostů rozdělení

Jedná se vlastně o demonstraci pravidel 2σ a 3σ. Nechť X je n. v. s rozdělením distr, E,X=μ, var X=(σ^2). Definujte funkci probabilityTales s jedním argumentem distr (rozdělení), která vytvoří tabulku typu (SD směrodatná odchylka, Standard deviation):

Ilustrujte  pro rozdělení normální, Studentovo s parametry μ, σ s počtem stupňů volnosti 4.1 a pro rozdělení výnosů společnosti VW (soubor vwvynosy.xls v pracovním adresáři).

Úloha 8

Pozor, v pravidlu 2σ resp. 3σ figuruje směrodatná odchylka n. v. X. Parametr σ v Laplaceově, Studentově a i jiných rozděleních nemá význam směrodatné odchylky!!!
Abychom dostali numerické hodnoty, přidáme desetinnou tečku třeba u výrazu 2σ ⟶ 2.σ

ClearAll[probabilityTales]
probabilityTales[distribution_]:=Module[{mean=Mean[distribution],sd=StandardDeviation[distribution],γ1=Skewness[distribution],γ2=Kurtosis[distribution],prob2σ,prob3σ},
Grid[{{"Distribution",distribution},
		{"Mean",mean},{"SD",sd},{"Skewness",γ1},{"Kurtosis",γ2},{"P(|X-μ|>2σ)",Probability[Abs[x-mean]>2.sd,x\[Distributed]distribution]},{"P(|X-μ|>3σ)",Probability[Abs[x-mean]>3.sd,x\[Distributed]distribution]}},Frame->All]
]

probabilityTales[NormalDistribution[μ,σ]]

probabilityTales[LaplaceDistribution[μ,σ]]

probabilityTales[StudentTDistribution[μ,σ,4.1]]

SetDirectory[NotebookDirectory[]];
vwVynosy=Import["vwvynosy.xls"];
%//Dimensions
vwVynosy=vwVynosy[[1,All,2]];
vwVynosDistribution=EmpiricalDistribution[vwVynosy]
probabilityTales[vwVynosDistribution]

probabilityTales[StudentTDistribution[μ,σ,4.1]/.{μ->.0003969,σ->((.01945)/(1.3972762620115438`))}]

Úloha 7
Rozdělení výběrového průměru pro různá rozdělení

Je-li (OverBar[X]) výběrový průměr náhodného výběru z libovolného rozdělení s konečným rozptylem (σ^2) a střední hodnotou μ, pak E,(OverBar[X])=μ a var (OverBar[X])=(σ^2)/n. Platí CLV.

Úloha: Generovat N náhodných výběrů o rozsahu n, z každého výběru spočítat výběrový průměr. Dostaneme tedy N hodnot, vlastně náhodný výběr výběrových průměrů. Graficky znázornit empirickou hustotu (histogram), “vyhlazený” histogram (SmoothHistogram) a do stejného grafu odpovídající hustotu normálního rozdělení, pozor na stejné parametry! Použijte rozdělení z předchozího odstavce. Efektivní postup spočívá ve vytvoření obecné procedury (Module), kde vstup je N, n, rozdělení  a volitelné semínko pro generování náhodných čísel (Optional, ), celkem 4 vstupní parametry.

Vyzkoušejte na 3 rozděleních N(4,4), exponenciálním s parametrem λ=5 a rozdělení výnosů společnostiApple.

Očekávaný výstup:

rozdeleniPrumeru[10000,5,ExponentialDistribution[3]]

Řešení 7

Obecná procedura:

Clear[rozdeleniPrumeru]
rozdeleniPrumeru[pocetVyberuN_,rozsahVyberun_,rozdeleni_,seminko_:13]:=
Module[{n=rozsahVyberun,nahodneVybery,xx,histogram,smoothHistogram,normalHustota,mean,sd},
SeedRandom[seminko];
mean=Mean[rozdeleni];sd=StandardDeviation[rozdeleni];
nahodneVybery=Table[xx=RandomVariate[rozdeleni,{n}];Mean[xx],{pocetVyberuN}];
histogram=Histogram[nahodneVybery,Automatic,"PDF",PlotLabel->rozdeleni];
smoothHistogram=SmoothHistogram[nahodneVybery,PlotLabel->rozdeleni,PlotStyle->Thick];
normalHustota=Plot[PDF[NormalDistribution[mean,sd/(Sqrt[n])],x],{x,mean-4sd/(Sqrt[n]),mean+4sd/(Sqrt[n])},PlotStyle->{Thick,Red}];
Show[{histogram,smoothHistogram,normalHustota}]
]

rozdeleniPrumeru[(* počet náhodných výběrů N = *)100000,(* rozsah každého výběru n = 50 *)500,ExponentialDistribution[5],131]

rozdeleniPrumeru[10000,4,NormalDistribution[4,2]]

GraphicsRow@{rozdeleniPrumeru[10000,10,StudentTDistribution[3]],rozdeleniPrumeru[10000,10,StudentTDistribution[30]]}

distr=EmpiricalDistribution[aapl20191105vynosyWithoutDates];rozdeleniPrumeru[2000,10,distr]

Test 2019-11-07 Komentář k zaslaným řešením

Obecné

1. Čtěte pozorně zadání

...
Předmět zprávy: 20191107_Jmeno_Prijmeni (porušeno v jednom případě)
Nepoužívejte diakritiku.  (porušeno v jednom případě)
...

... pro rozdělení výnosů z dat cen společnosti Apple (soubor aapl20191105ceny.xlsx v pracovním adresáři) ... : v uvedeném soubory jsou ceny a z nich se měly nejdříve vypočítat výnosy! Kvůli tomu jsou výnosy v zadání zdůrazněny tučně! (porušeno vícekrát).

2. Rozlišení formátu buněk

Reseni

má formát vstupní buňky, tj. Input. Správně má být Section, Subsection, ... . Viz Menu -> Format -> Style -> Subsection nebo z klávesnice Alt+5.

⚠ Dbejte na formátování, je součástí hodnocení při zkoušce.

Užitím systému Mathematica jsou sepsány desítky knih, na MFF UK minimálně jedna bakalářská práce: Adam Kubetta, Simulace ve financích, MFF UK 2009.

3. Označení proměnných

Důsledně dodržujte pravidlo: každý uživatelský identifikátor začíná malým písmenem. Profesionálové používají značení s malým i jako prvním znakem. Místo “NN” resp. “Pocet” tedy použít “iNN” resp. “iPocet”. V zadání úlohy 7 tedy použít např. “in”, “iN”.

4. Buďte pozorní

Ze souboru “aapl20191105ceny.xlsx” správně vyberete 1. element takže dostanete matici {1527,2}, první sloupec data, druhý sloupec ceny (znovy zdůrazňuji: ne výnosy). Tuto matici si označíte jako “ceny”. Pak zkusíte EmpiricalDistribution[ceny], ale to se má aplikovat pouze na druhý sloupec a ne na data!

Úloha 5

1.

V zadání je požadavek na zpracování výnosů , a ne cen!!!

Úloha 6
Pravděpodobnosti chvostů rozdělení

Úloha 7
Rozdělení výběrového průměru pro různá rozdělení

1.

S ohledem na to, že ve vstupu nejsou specifikovány numerické hodnoty, které jsou vyžadovány (např. u konstrukce Table a všude, kde se vyskytuje iterátor), je zapotřebí místo okamžitého přiřazení použít odložené přiřazení:

(* Špatně *)
nahodnyVyber[num_Integer, n_, distribution_, seed_]= abc
(* Správně *)
nahodnyVyber[num_Integer, n_, distribution_, seed_]:= abc

Metody odhadu parametrů

Zabudované funkce

?*ParameterEstimator*

?*Estimated*

?*Likelihood*

2019-11-21 Cvičení 8

Doporučuji včasnou účast (7:20) na přednášce 2019-11-26. Téma: Asociovaná pole.

FAQ

Testování statistických hypotéz 1

Je rozhodovací proces. Formuluje se (nulová) hypotéza (Subscript[H, 0]), z dostupných dat se sestrojí testová statistika a na jejím základě se provede test hypotézy. Výsledkem té rozhodovací procedury je výrok: “(Subscript[H, 0]) zamítáme” nebo “(Subscript[H, 0]) nemůžeme zamítnout”. Ve skutečnosti (Subscript[H, 0]) platí nebo (Subscript[H, 0]) neplatí. Jsou tedy čtyři možnosti shrnuté v následující tabulce:

data={{"(Subscript[H, 0]) nezamítáme","OK","chyba 2. druhu"},{"(Subscript[H, 0]) zamítáme","chyba 1. druhu","OK"}};Text@Grid[Prepend[data,{Column[{"Skutečnost  →","Rozhodnutí  ↓ "}],Column[{"(Subscript[H, 0]) platí",""}],Column[{"(Subscript[H, 0]) neplatí",""}]}],Background->{None,{Lighter[Yellow,.9],{White,Lighter[Blend[{Blue,Green}],.8]}}},Dividers->{{Darker[Gray,.6],{Lighter[Gray,.5]},Darker[Gray,.6]},{Darker[Gray,.6],Darker[Gray,.6],{False},Darker[Gray,.6]}},Alignment->{{Left,Center,{Center}}},ItemSize->{{10,9,9,5}},Frame->Darker[Gray,.6],ItemStyle->14,Spacings->{Automatic,.8}]

Pravděpodobnost chyby 1. druhu můžeme kontrolovat, předepíšeme její pravděpodobnost zvanou hladina testu α: P((Subscript[H, 0]) zamítáme|(Subscript[H, 0])platí)=α. Pravděpodobnost chyby 2. druhu je dána použitým testem (statistik by měl volit ten nejvhodnější se zřetelem na rozdělení záladního souboru, alternativu apod. Od toho se odvíjí další důležitá charakteristika zvaná síla testu β=1-P((Subscript[H, 0]) nezamítáme|(Subscript[H, 0])neplatí), tedy β=1-P(chyba 2. druhu).

Metody odhadu parametrů 1

Zabudované funkce

EstimatedDistribution

EstimatedDistribution

?*ParameterEstimator*

Likelihood, LogLikelihood

LogLikelihood

Úlohy

Úloha 8

Uvažujte soubor výnosů firmy VW ze  “vwvynosy.xls” v pracovním adresáři. 
(i) Odhadněte parametry rozdělení jestliže se domníváte, že data pochází z N(μ,(σ^2)). 

(ii) Odhadněte parametry rozdělení jestliže se domníváte, že data pochází ze Studentova rozdělení s parametry μ,σ,ν, všechny parametry neznámé.

Řešení 8

SetDirectory[NotebookDirectory[]];
vwData=Import["vwvynosy.xls"];
vwReturns=vwData[[1,All,2]];
{Min[vwReturns],Max[vwReturns]}
Through[{Mean,StandardDeviation,Skewness,Kurtosis}[vwReturns]]

(i) Normální rozdělení

edNormal=EstimatedDistribution[vwReturns,NormalDistribution[μ,σ]]

(ii) Studentovo rozdělení

edStudent=EstimatedDistribution[vwReturns,StudentTDistribution[μ,σ,ν]]
Through[{Mean,StandardDeviation,Skewness,Kurtosis}[edStudent]]
Show[Plot[PDF[SmoothKernelDistribution[vwReturns],x],{x,-.16,.1},PlotStyle->Thick],Plot[PDF[edNormal,x],{x,-.16,.1},PlotStyle->{Thick,Red}]]

Úloha 9

Testujte hypotézu, že data z předchozí úlohy pochází 
(i) z normálního rozdělení
(ii) ze Studentova rozdělení.
Použijte Kolmogorov-Smirnovův test, zabudovaná funkce je KolmogorovSmirnovTest. Výsledky slovně interpretujte.

Řešení 9

(i) Normální rozdělení

KolmogorovSmirnovTest[vwReturns,NormalDistribution[μ,σ]]

testKSNormal=KolmogorovSmirnovTest[vwReturns,edNormal]

Detaily

testKSNormalDetails=KolmogorovSmirnovTest[vwReturns,edNormal,"HypothesisTestData"]

propertiestestKSNormalDetails=testKSNormalDetails["Properties"]

testKSNormalDetails["TestConclusion"]//Quiet//TraditionalForm

(ii) Studentovo rozdělení

KolmogorovSmirnovTest[vwReturns,StudentTDistribution[μ,σ,ν]]

testKSStudent=KolmogorovSmirnovTest[vwReturns,edStudent]

KolmogorovSmirnovTest[vwReturns,edStudent,"TestConclusion"]//Quiet//TraditionalForm

Úloha 10

1. Pro data úlohy 8 odhadněte momentovou metodou parametry
(i) N(μ,(σ^2))
(ii) StudentTDistribution[μ,σ,ν].

2. Pro data úlohy 8 odhadněte metodou maximální věrohodnosti parametry
(i) N(μ,(σ^2))
(ii) StudentTDistribution[μ,σ,ν].

Řešení 10.1

(* i *)

estimatedNormalMoments=EstimatedDistribution[vwReturns,NormalDistribution[μ,σ],ParameterEstimator->"MethodOfMoments"]

(* ii *)

estimatedStudentMoments=EstimatedDistribution[vwReturns,StudentTDistribution[μ,σ,ν],ParameterEstimator->"MethodOfMoments"]
Through[{Mean,StandardDeviation,Skewness,Kurtosis}[estimatedStudentMoments]]

V úloze 8 jsme  dospěli k odhadu

StudentTDistribution[0.0003182870367486177`,0.01466459251902522`,4.629194067408561`]
Through[{Mean,StandardDeviation,Skewness,Kurtosis}[%]]

Řešení 10.2

(* i *)

estimatedNormalML=EstimatedDistribution[vwReturns,NormalDistribution[μ,σ],ParameterEstimator->"MaximumLikelihood"]

(* ii *)
estimatedStudentML=EstimatedDistribution[vwReturns,StudentTDistribution[μ,σ,ν],ParameterEstimator->"MaximumLikelihood"]
Through[{Mean,StandardDeviation,Skewness,Kurtosis}[estimatedStudentMoments]]

2019-11-28 Cvičení 9

FAQ

Testování statistických hypotéz 2

Připomeňme hladinu testu (hladinu významnosti, significance level) α: P((Subscript[H, 0]) zamítáme|(Subscript[H, 0])platí)=α. Frekvenční interpretace: Máme 100 výběrů z rozdělení pro které platí nulová hypotéza (Subscript[H, 0]). Provedeme  100 testů. V nich by se měla hypotéza zamítnout v přibližně 100α případech. V případě α=0.05 by se tak nulová hypotéza zamítnout přibližně pětkrát.

Úloha 11

Navrhněte simulační proceduru pro ověření výše uvedeného tvrzení pro počet výběrů (počet testů) N, rozsahu výběru n, z rozdělení N(0,1), (Subscript[H, 0]):μ=0. Pro testování použijte zabudovanou funkci TTest.

Řešení 11

Clear[hladinaTestu]
hladinaTestu[iN_,in_,seed_:131]:=Module[{vybery,pValues,countRejected,p},
	SeedRandom[seed];
	vybery=RandomVariate[NormalDistribution[],{iN,in}];
	pValues=Map[TTest[#]&,vybery];
	Row@{"Počet zamítnutí = ",countRejected=Count[pValues,_?(Function[p,p≤0.05])]," z ",iN," výběrů\n","Relativní četnost zamítnutí = ",(countRejected/iN)//N}
			]

(* Alternativně *)
Length[Select[pvalues,#≤0.05&]]

hladinaTestu[1000,50]//Quiet
hladinaTestu[1000,50,13131]//Quiet
hladinaTestu[1000,50,131313]//Quiet

Zobecnění

Clear[hladinaTestuObecne]
hladinaTestuObecne[iN_,in_,distribution_,α_(* hladina testu *),seed_:131]:=Module[{vybery,pValues,countRejected,p},
	SeedRandom[seed];
	vybery=RandomVariate[distribution,{iN,in}];
	pValues=Map[TTest[#,SignificanceLevel->α]&,vybery];
	Row@{"Počet zamítnutí = ",countRejected=Count[pValues,_?(Function[p,p≤α])]," z ",iN," výběrů\n","Relativní četnost zamítnutí = ",(countRejected/iN)//N}
			]

TTest[RandomVariate[NormalDistribution[],{2,100}]]

hladinaTestuObecne[1000,50,NormalDistribution[],0.05]//Quiet
hladinaTestuObecne[1000,50,NormalDistribution[],0.01]//Quiet
hladinaTestuObecne[1000,50,LaplaceDistribution[0,1],0.05]//Quiet

(* Výběry o rozsahu 15 *)
hladinaTestuObecne[1000,15,NormalDistribution[],0.05]//Quiet
hladinaTestuObecne[1000,15,NormalDistribution[],0.01]//Quiet
hladinaTestuObecne[1000,15,LaplaceDistribution[0,1],0.05]//Quiet

Metody odhadu parametrů 2

Směs rozdělení (Mixture distribution)

Nechť (Subscript[p, 1]),...,(Subscript[p, K]) je rozdělení pravděpodobností, tj. (Subscript[p, k])≥0, (∑_(k=1)^K)(Subscript[p, k])=1. Nechť (Subscript[F, 1]),...,(Subscript[F, K]) jsou distribuční funkce. Pak

F=(Subscript[p, 1])(Subscript[F, 1])+...+(Subscript[p, K])(Subscript[F, K])

je distribuční funkce a nazývá se distribuční funkcí směsi (Subscript[F, 1]),...,(Subscript[F, K]) podle (Subscript[p, 1]),...,(Subscript[p, K]). Odpovídající rozdělení se nazývá směsí rozdělení (mixture distribution). Podobně pro hustoty:

f=(Subscript[p, 1])(Subscript[f, 1])+...+(Subscript[p, K])(Subscript[f, K]).

Takové rozdělení je flexibilní, modeluje i těžké chvosty (heavy tails, leptokurtosis), problémy mohou vzniknout při odhadu parametrů.

Velmi flexibilní je např. směs tří normálních se stejným rozptylem a vahami (pravděpodobnostmi) ϵ,1-2ϵ,ϵ, 0<ϵ<1/2:

ϵ,N((Subscript[μ, 1]),(σ^2))+(1-2ϵ)N((Subscript[μ, 2]),(σ^2))+ϵ,N((Subscript[μ, 3]),(σ^2))

Uživatelská funkce založená na zabudované funkci pro tento případ:

𝒟

Clear[𝒟mixtureOfNormal]
𝒟mixtureOfNormal[ϵ_,μ1_,μ2_,μ3_,σ_]:= MixtureDistribution[{ϵ,1-2ϵ,ϵ},{NormalDistribution[μ1,σ],NormalDistribution[μ2,σ],NormalDistribution[μ3,σ]}]

Příklad trimodálního rozdělení

𝒟mixture=𝒟mixtureOfNormal[1/20,-3,0,4,2/3]
Plot[PDF[𝒟mixture,x],{x,-5,5},PlotLabel->"The density of mixture\n",PlotRange->All,PlotStyle->Thick]

Manipulate[Plot[PDF[MixtureDistribution[{ϵ,1-2ϵ,ϵ},{NormalDistribution[μ1,σ],NormalDistribution[μ2,σ],NormalDistribution[μ3,σ]}],x],{x,-5,5},PlotLabel->"The density of mixture",PlotRange->All],
{{ϵ,0.05},0,0.2,Appearance->"Labeled"},
{{μ1,-2},-3,0,Appearance->"Labeled"},
{{μ2,0},-1,2,Appearance->"Labeled"},
{{μ3,2},0,3,Appearance->"Labeled"},
{{σ,0.5},0.1,2,Appearance->"Labeled"},SaveDefinitions->True]

Manipulate[Plot[Kurtosis[MixtureDistribution[{ϵ,1-2ϵ,ϵ},{NormalDistribution[μ1,σ],NormalDistribution[μ2,σ],NormalDistribution[μ3,σ]}]],{σ,0.1,5},PlotLabel->"Kurtosis",PlotRange->All,AxesLabel->{"σ",Automatic}],
{{ϵ,0.05},0,0.2,Appearance->"Labeled"},
{{μ1,-2},-3,0,Appearance->"Labeled"},
{{μ2,0},-1,2,Appearance->"Labeled"},
{{μ3,2},0,3,Appearance->"Labeled"},
SaveDefinitions->True]

2019-12-05 Cvičení 10

Úloha 12

Odhadněte parametry směsi rozdělení distribution pro data vwvynosy.xls z pracovního adresáře.

Řešení 12

Clear[𝒟mixtureOfNormal]
𝒟mixtureOfNormal[ϵ_,μ1_,μ2_,μ3_,σ_]:= MixtureDistribution[{ϵ,1-2ϵ,ϵ},{NormalDistribution[μ1,σ],NormalDistribution[μ2,σ],NormalDistribution[μ3,σ]}]

SetDirectory[NotebookDirectory[]];
vwvynosy=Import["vwvynosy.xls"][[1,All,2]];

Clear[ϵ,μ1,μ2,μ3,σ]
𝒟estimatedVWvynosy=EstimatedDistribution[vwvynosy,𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]]

Through[{Mean,StandardDeviation,Skewness,Kurtosis}[vwvynosy]]
Through[{Mean,StandardDeviation,Skewness,Kurtosis}[𝒟estimatedVWvynosy]]
Through[{Min,Max}[vwvynosy]]

Show[Plot[PDF[𝒟estimatedVWvynosy,x],{x,-0.15,0.15},PlotStyle->Thick],
SmoothHistogram[vwvynosy,Automatic,"PDF",PlotStyle->{Red,Thick}]]

Plot[PDF[𝒟estimatedVWvynosy,x],{x,-0.2,-.02},PlotStyle->Thick]

Míry rizika

Test

Notes (tj. tento notebook s níže uvedeným názvem) s řešením zašlete e-mailem na mou adresu hurt@karlin.mff.cuni.cz.
Název notesu: 20191205_Jmeno_Prijmeni
Předmět zprávy: 20191205_Jmeno_Prijmeni
Nepoužívejte diakritiku.
Termín: 2019-12-05 12:30

Úvod

Připomeneme, že za míru výnosnosti, krátce výnos považujeme za náhodnou veličinu a v daném kontextu označujeme ρ, (Subscript[ρ, t]), (Subscript[R, t]) nebo ρ(t), potřebujeme-li vyjádřit závislost na čase (pozor na časové jednotky!). Náhodná veličina

L=-ρ ,              (Subscript[L, t])=(((Subscript[P, t])-(Subscript[P, t+1]))/(Subscript[P, t]))

se nazývá ztráta (loss) nebo také ztrátová náhodná veličina (loss random variable).

V praxi, je-li výnos ρ>0, jde o  zisk, je-li výnos ρ<0, jde o ztrátu, která se však také uvádí s kladným znaménkem, tj. uvádí se L.

Rozptylové míry rizika

Směrodatná odchylka, volatilita (Standard Deviation, SD, volatility)

Je nejjednodušší mírou rizika. V teorii portfolia se přímo rizikem nazývá.

σ=(Sqrt[var ρ])=(Sqrt[var L]).

Nevýhoda: citlivá na “velké” odchylky, chyby pozorování.

Střední absolutní odchylka (mean absolute deviation, MAD)

MAD(ρ)=MAD(L)=E|ρ-E,ρ|=E|L-E,L|.

Robustní vůči chybám a velkým odchylkám.

Example

The assets A and B possess returns dependent on the state of economy in the following way:

Text@Grid[{{"State of Economy","Probability","Return A [%]","Return B [%]"},
{"Boom",0.3,100,20},
{"Normal",0.4`,15,15},
{"Recession",0.3`,-70,10}},Frame->All]

𝒟assetA=EmpiricalDistribution[{0.3,0.4,0.3}->{100,15,-70}];
𝒟assetB=EmpiricalDistribution[{0.3,0.4,0.3}->{20,15,10}];
Text@Grid[{{" ","Asset A","Asset B"},
{"Mean",Mean[𝒟assetA],Mean[𝒟assetB]},
{"Variance",Variance[𝒟assetA],Variance[𝒟assetB]},
{"S. D.",StandardDeviation[𝒟assetA],StandardDeviation[𝒟assetB]},
{"MAD",Expectation[Abs[x-Mean[𝒟assetA]],x\[Distributed]𝒟assetA],Expectation[Abs[x-Mean[𝒟assetB]],x\[Distributed]𝒟assetB]}},Frame->All]

Hodnota v riziku (Value at Risk, VaR)

Value at Risk at (confidence) level 1-α (VaR)

Hodnota v riziku na úrovni spolehlivosti 1-α:

Value at risk at (confidence) level 1-α is simply the α-quantile of L, sometimes denoted by (Subscript[q, α]).

(Subscript[VaR, α])(L):=(F^(-1))(α)=inf{x:F(x)≥α}=:(Subscript[q, α]).

The usual interpretation is that loss will exceed (Subscript[VaR, α])(L) with a (small) probability 1-α with typical values of α=0.95, 0.99 or in some extreme cases even α=0.999. In practice the time horizons are one day, one week, or one month but often longer time horizons are of interest. These need a special care, however. For the sake of simplicity of some calculations we often suppose (1/2)<α<1.

Interpretace: Ztráta překročí hodnotu (Subscript[VaR, α])(L) s (malou) pravděpodobností 1-α. Typické hodnoty α=0.95, 0.99, v některých extrémních případech α=0.999. V praxi používané časové horizonty jsou jeden den, jeden týden, jeden měsíc. Někdy je však nutné použít i delší časové horizonty. Pro zjednodušení některých výpočtů se často předpokládá (1/2)<α<1.

Hullova interpretace hodnoty v riziku V na úrovni splehlivosti 1-X s časovým horizontem N dnů:

Hull, J. C.: Options, Futures, and Other Derivatives, 8th edition, p. 471: “Jsem si na X procent jist, že ztráta v příštích N dnech nebude vyšší než V dolarů.”

Zde X=α, V=(Subscript[VaR, α]).

Ještě jiná slovní interpretace: Hodnota v riziku je mez, nad kterou se ztráta dostane jen s velmi malou předepsanou pravděpodobností.

Poznámka

Pozor na nejednoznačnost! Někdy se v literatuře místo úrovně spolehlivost 1-α uvažuje úroveň spolehlivosti α. Smysl ale zůstává týž: Ztráta překročí VaR s malou pravděpodobností.

Podmíněná hodnota v riziku (Conditonal VaR, CVaR, Tail VaR, TVaR, expected shortfall)

Conditional value at risk (in case of continuous distributions also called Expected Shortfall), Očekávaná extrémní ztráta), Tail Conditional Expectation, Mean Excess Loss, Tail VaR) is the mean of the truncated distribution of L truncated at point (Subscript[VaR, α])(L) from the left:

(Subscript[CVaR, α])(L)=E{L|L≥(Subscript[VaR, α])(L)}.

Note that sometimes an analogue to the mean residual lifetime is considered:

(Subscript[CVaR, α])(L)-(Subscript[VaR, α])(L)=E{L-(Subscript[VaR, α])(L)|L≥(Subscript[VaR, α])(L)}.

For absolutely continuous distributions with the distribution function F(y) and density f(y) we have

(Subscript[CVaR, α])(L)=(1/(1-α))(∫_( (Subscript[Var, α])(L))^( ∞))y f(y)dy

or alternatively (after substitution F(y)=z, y=(F^(-1))(z), dz=f(y)dy)

(Subscript[CVaR, α])(L)=(1/(1-α))(∫_( α)^( 1))(F^(-1))(z)dz.

Poznámka (nejednoznačná terminologie)

V některých pramenech se v poslední době uvádí (podle mně nesprávně) poslední vzorec jako definice podmíněné hodnoty v riziku. Pro rozdělení, která nejsou absolutně spojitá (specielně diskrétní rozdělení) se mohou obě vyjádření lišit.

2019-12-12 Cvičení 11

Úlohy: Procedury pro výpočet VaR a CVaR, porovnání pro různá rozdělení

Ve všech následujících úlohách budeme pracovat s daty společnosti (Entity["Walt Disney", Entity["Company","Disney::59wfk"], "Entity[\"Company\", \"Disney::59wfk\"]", "company"]). Vypreparované ceny jsou v souboru  “pricesDisney.XML” (dále DIS) v pracovním adresáři.

Poznamenejme, že *.XML je pro Export ⟷ Import jeden ze spolehlivých universálních formátů. Je ale náročný na čas a na paměť. V rámci Wolfram Language je bezpečný Save ⟷ Get.

Úloha 13

Výsledky bodů 1 a 3 presentujte přehledně v tabulce.

1. Spočtěte výběrové charakteristiky DIS (průměr, SD, MAD, (Subscript[γ, 1]),(Subscript[γ, 2])). Ceny graficky znázorněte.

2. Z dat DIS spočtěte výnosy a ztráty (LossesDIS).

3. Z LossesDIS spočtěte výběrové charakteristiky DIS (průměr, SD, MAD, (Subscript[γ, 1]),(Subscript[γ, 2]))

4. Sestrojte vyhlazený histogram LossesDIS.

Řešení 13

13.1

SetDirectory[NotebookDirectory[]];
iDIS=Import["pricesDisney.XML"];

Střední absolutní odchylka je sice zabudovaná funkce (MeanDeviation), ale špatně se hledá. Navíc jejím argumentem může být jen seznam (náhodný výběr). Návody k hledání:

?*Mean*

?*Deviation*

Je však velmi jednoduché definovat vlastní funkci:

Clear[meanAbsoluteDeviation]
meanAbsoluteDeviation[randomSample_?VectorQ]:=Mean[Abs[randomSample-Mean[randomSample]]]

{meanAbsoluteDeviation[iDIS],MeanDeviation[iDIS]}

Grid[{{Mean,StandardDeviation,MeanDeviation,Skewness,Kurtosis},Through[charakteristiky={Mean,StandardDeviation,MeanDeviation,Skewness,Kurtosis}[iDIS]]},Frame->All]

ListLinePlot[iDIS,PlotLabel->"Prices Walt Disney"]

13.2

Clear[returns];
returns[prices_?VectorQ]:=((Rest[prices]-Most[prices])/(Most[prices]))

returnsDIS=returns[iDIS];
lossesDIS=-returnsDIS;
%//Length

13.3

Grid[{{Mean,StandardDeviation,MeanDeviation,Skewness,Kurtosis},Through[charakteristiky={Mean,StandardDeviation,MeanDeviation,Skewness,Kurtosis}[lossesDIS]]},Frame->All]

13.4

GraphicsRow@{SmoothHistogram[lossesDIS,PlotLabel->"Losses Walt Disney",PlotRange->All],SmoothHistogram[lossesDIS,PlotLabel->"Losses Walt Disney"]}

Úloha 14

1. Pro data LossesDIS vytvořte tabulku empiricky odhadnutých hodnot v riziku na úrovních 0.99, 0.95, 0.9 ve tvaru

Řešení 14

Grid[Prepend[Table[{1-α,Quantile[lossesDIS,α]},{α,{.99,.95,.9}}],{1-α,(Subscript[VaR, α])}],Frame->All]

nebo jednoduššeji (nevhodné pro větší tabulky)

Grid[{{1-α,(Subscript[VaR, α])},{0.01,Quantile[lossesDIS,0.99]},{0.05,Quantile[lossesDIS,0.95]},{0.1,Quantile[lossesDIS,0.9]}},Frame->All]

Úloha 15

1. Pro data LossesDIS vytvořte vyhlazené jádrové rozdělení, nazvěme jej iSmoothKernelDistributionLossesDIS.

2. Odhadněte parametry směsi rozdělení distribution pro data LossesDIS. Odpovídající rozdělení nechť je iMixtureDistributionLossesDIS.

3. Odhady střední hodnoty a směrodatné odchylky z 13.1 (špatný odkaz v zadání, správně má být 13.3) použijte pro definici normálního rozdělení s těmito parametry, nazvěme jej iNormalDistributionLossesDIS.

4. V jednom grafu znázorněte hustoty všech tří výše uvedených rozdělení.

Řešení 15

15.1

iSmoothKernelDistributionLossesDIS=SmoothKernelDistribution[lossesDIS]

15.2

Clear[𝒟mixtureOfNormal]
𝒟mixtureOfNormal[ϵ_,μ1_,μ2_,μ3_,σ_]:= MixtureDistribution[{ϵ,1-2ϵ,ϵ},{NormalDistribution[μ1,σ],NormalDistribution[μ2,σ],NormalDistribution[μ3,σ]}]

iMixtureDistributionLossesDIS=EstimatedDistribution[lossesDIS,𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]]

Ještě můžeme zkusit

iFindDistributionLossesDIS=FindDistribution[lossesDIS]

a graficky porovnat:

Plot[{PDF[iMixtureDistributionLossesDIS,x],PDF[iFindDistributionLossesDIS,x]},{x,-.1,.1},PlotRange->All,PlotTheme->"Detailed",PlotStyle->Thick]

Je patrné, že funkce FindDistribution dává v tomto případě nesmyslný výsledek.

15.3

iNormalDistributionLossesDIS=NormalDistribution[Mean[lossesDIS],StandardDeviation[lossesDIS]]

15.4

Pro představu o nosiči rozdělení pro účely grafického znázornění spočítáme

minmax={Min[lossesDIS],Max[lossesDIS]}
quantilesRange={Quantile[lossesDIS,0.01],Quantile[lossesDIS,0.99]}

Všimněte si, že maximální denní zisk byl téměř 16% a maximální denní ztráta menší než 10%.

Plot[{PDF[iSmoothKernelDistributionLossesDIS,x],PDF[iMixtureDistributionLossesDIS,x],PDF[iNormalDistributionLossesDIS,x]},{x,-.1,.1},PlotRange->All,PlotTheme->"Detailed",PlotStyle->Thick]

Test

Notes (tj. tento notebook s níže uvedeným názvem) s řešením zašlete e-mailem na mou adresu hurt@karlin.mff.cuni.cz.
Název notesu: 20191205_Jmeno_Prijmeni
Předmět zprávy: 20191205_Jmeno_Prijmeni
Nepoužívejte diakritiku.
Termín: 2019-12-05 12:30

Chyby

1.

ztráty =- výnosy. Místo se ztrátami jste počítali s výnosy. Další vyskytující se chyba byla v záměně cen a výnosů.

Několikrát se vyskytující chyba. V řadě firem se první taková chyba toleruje. Při druhé takové chybě čeká zaměstnance vyhazov. Totéž se týká případu, kdy místo 37.23445 uvedete 372.3445  a hájíte se tím, že je tam chyba jenom v desetinné čárce. Zaměstnavatele zvláště potěší, že aktiva firmy místo za 372.3445 mil. prodáte za 37.23445.

2.

Pozornost: MAD=Mean[Abs[disData-Mean[LossesDIS]. disData ⟶ LossesDIS

3.

vynosyDis[t_]:= (ceny[[t]]- ceny[[t+1]])/ ceny[[t]]

Pojmenování je matoucí. Jsou to ztráty. Reference na části výrazu (ceny[[t]]) a následné použití tabulky je zbytečné a neefektivní.

4. Moje chyba v zadání úlohy 15

3. Odhady střední hodnoty a směrodatné odchylky z 13.1 (špatný odkaz v zadání, správně má být 13.3) použijte pro definici normálního rozdělení s těmito parametry, nazvěme jej iNormalDistributionLossesDIS.

Někdo pochopil, že jde o chybu v zadání, někdo ne.

5.

V jednom případě jsem dostal nečitelný soubor byť s příponou *.nb. Příště před odesláním soubor uložte, zavřete, potom znovu otevřte, zkontrolujte, že je čitelný a teprve poté odešlete.

6.

Nevěřte bezmezně metodám strojového učení a umělé inteligence  (machine learning, Artificial Intelligence) - již výše uvedeno, že může vést k nesmyslům:

Clear[vynosy, dist];
SetDirectory[NotebookDirectory[]];
vynosy=Import["vwvynosy.xls"][[1,All,2]];
dist=FindDistribution[vynosy]
Show[Plot[PDF[dist,x],{x,-.15,.15},PlotStyle->{Thick,Red},PlotRange->All],SmoothHistogram[vynosy]]

Načtení dat

SetDirectory[NotebookDirectory[]];
iDIS=Import["pricesDisney.XML"];
Clear[returns];
returns[prices_?VectorQ]:=((Rest[prices]-Most[prices])/(Most[prices]))
returnsDIS=returns[iDIS];
lossesDIS=-returnsDIS;

Podmíněná hodnota v riziku (Conditonal VaR, CVaR, Tail VaR, TVaR, expected shortfall)

Conditional value at risk (in case of continuous distributions also called Expected Shortfall), Očekávaná extrémní ztráta), Tail Conditional Expectation, Mean Excess Loss, Tail VaR) is the mean of the truncated distribution of L truncated at point (Subscript[VaR, α])(L) from the left:

(Subscript[CVaR, α])(L)=E{L|L≥(Subscript[VaR, α])(L)}.

Note that sometimes an analogue to the mean residual lifetime is considered:

(Subscript[CVaR, α])(L)-(Subscript[VaR, α])(L)=E{L-(Subscript[VaR, α])(L)|L≥(Subscript[VaR, α])(L)}.

For absolutely continuous distributions with the distribution function F(y) and density f(y) we have

(Subscript[CVaR, α])(L)=(1/(1-α))(∫_( (Subscript[Var, α])(L))^( ∞))y f(y)dy

or alternatively (after substitution F(y)=z, y=(F^(-1))(z), dz=f(y)dy)

(Subscript[CVaR, α])(L)=(1/(1-α))(∫_( α)^( 1))(F^(-1))(z)dz.

Poznámka (nejednoznačná terminologie)

V některých pramenech se v poslední době uvádí (podle mně nesprávně) poslední vzorec jako definice podmíněné hodnoty v riziku. Pro rozdělení, která nejsou absolutně spojitá (specielně diskrétní rozdělení) se mohou obě vyjádření lišit.

Úloha 16

Pro data LossesDIS vytvořte tabulku empiricky odhadnutých hodnot v riziku, podmíněných hodnot v riziku na úrovních 0.99, 0.95, 0.9 ve tvaru

Řešení 16

SetDirectory[NotebookDirectory[]];
iDIS=Import["pricesDisney.XML"];
Clear[returns];
returns[prices_?VectorQ]:=((Rest[prices]-Most[prices])/(Most[prices]))
returnsDIS=returns[iDIS];
lossesDIS=-returnsDIS;

Clear[iVaR,iEDistribution]
(* Argumentem iVaR mohou být i data: *)
iVaR[distribution_,α_]:=Quantile[distribution,α]
iEDistribution[data_?VectorQ]:=EmpiricalDistribution[data]

iHodnotyvRiziku=Table[iVaR[lossesDIS,α],{α,{.99,.95,.9}}]

1. způsob

Z náhodného výběru se vyberou hodnoty ≥VaR a z těch se spočítá průměr:

Mean[Select[lossesDIS,#≥iVaR[lossesDIS,0.95]&]]

2. způsob

Spočítá se podmíněná střední hodnota:

Expectation[L\[Conditioned]L≥iVaR[lossesDIS,0.95],L\[Distributed]iEDistribution[lossesDIS]]

Clear[iCVaR]
iCVaR[data_?VectorQ,α_?(0<#<1&)]:=Mean[Select[data,#≥iVaR[data,α]&]]

Grid[Prepend[Table[{1-α,Quantile[lossesDIS,α],iCVaR[lossesDIS,α]},{α,{.99,.95,.9}}],{1-α,(Subscript[VaR, α]),(Subscript[CVaR, α])}],Frame->All]

Úloha 17

Z dat LossesDIS odhadněte parametry následujích rozdělení: 𝒟1=normální, 𝒟2=Studentovo, 𝒟3=Laplaceovo. Dále uvažujte rozdělení 𝒟4=vyhlazené jádrové a 𝒟5=empirické. Pro úrovně 0.99, 0.95, 0.9 vytvořte 3 tabulky typu

2019-12-19 Cvičení 12

FAQ

Úloha 17

Úloha 17

Z dat LossesDIS odhadněte parametry následujích rozdělení: 𝒟1=normální, 𝒟2=Studentovo, 𝒟3=Laplaceovo. Dále uvažujte rozdělení 𝒟4=vyhlazené jádrové a 𝒟5=empirické. Pro úrovně 0.99, 0.95, 0.9 vytvořte 3 tabulky typu

Řešení 17

SetDirectory[NotebookDirectory[]];
iDIS=Import["pricesDisney.XML"];
Clear[returns];
returns[prices_?VectorQ]:=((Rest[prices]-Most[prices])/(Most[prices]))
returnsDIS=returns[iDIS];
lossesDIS=-returnsDIS;

𝒟[i]

fiveDistributions={𝒟1=EstimatedDistribution[lossesDIS,NormalDistribution[μ,σ]],
𝒟2=EstimatedDistribution[lossesDIS,StudentTDistribution[μ,σ,ν]],
𝒟3=EstimatedDistribution[lossesDIS,LaplaceDistribution[μ,β]],
𝒟4=SmoothKernelDistribution[lossesDIS],
𝒟5=EmpiricalDistribution[lossesDIS]};

Alternativní definice CvaR: (Subscript[CVaR, α])(L)=(1/(1-α))(∫_((Subscript[Var, α])(L))^∞)y,f(y)dy. Nemůže fungovat pro 𝒟5 - nemá hustotu. Tam ale zafunguje definice CVaR: (Subscript[CVaR, α])(L)=E[L|L≥(Subscript[VaR, α])(L)].

Clear[CVaR,NCVaR]
CVaR[distribution_,α_]:=(1/(1-α))(∫_(Quantile[distribution,α])^∞)y PDF[distribution,y]\[DifferentialD]y
NCVaR[distribution_,α_]:=(1/(1-α))NIntegrate[y PDF[distribution,y],{y,Quantile[distribution,α],∞}]

Postup

α=0.95;
(* Hodnoty v riziku pro všechna rozdělení; ukázka, že iterátorem může být cokoli: *)var=Table[Quantile[i,α],{i,{𝒟1,𝒟2,𝒟3,𝒟4,𝒟5}}]
(* Funkcionální řešení *)
var=Map[Quantile[#,α]&,fiveDistributions]

cvar=Map[NCVaR[#,α]&,fiveDistributions]

Nefunguje pro empirické rozdělení, musíme ošetřit zvlášť. Pro první čtyři rozdělení (Most ze souboru fiveDistributions)  funkci NCVaR použít můžeme, pro empirické rozdělení vyjdeme z definice.

{varEmpirical=Quantile[𝒟5,.99],CVaREmpirical=Expectation[loss\[Conditioned]loss≥varEmpirical,loss\[Distributed]𝒟5]}

Pozor, musíme použít Flatten abychom dostali vektor.

(* Neregulární struktura: *) {Map[NCVaR[#,α]&,fiveDistributions//Most],Expectation[loss\[Conditioned]loss≥Quantile[𝒟5,α],loss\[Distributed]𝒟5]}
(* Vektor: *)
cvar={Map[NCVaR[#,α]&,fiveDistributions//Most],Expectation[loss\[Conditioned]loss≥Quantile[𝒟5,α],loss\[Distributed]𝒟5]}//Flatten

Pro všechny tři úrovně (v prvním sloupci VaR, ve druhém CVaR):

(tabulka=Table[{Map[Quantile[#,α]&,fiveDistributions],
{Map[NCVaR[#,α]&,fiveDistributions//Most],Expectation[loss\[Conditioned]loss≥Quantile[𝒟5,α],loss\[Distributed]𝒟5]}//Flatten},{α,αα={.99,.95,.90`}}])//TableForm

A všechno dohromady:

Column@Table[Grid[Prepend[{{"𝒟1","𝒟2","𝒟3","𝒟4","𝒟5"},tabulka[[i,1]],tabulka[[i,2]]}//Transpose,{"1-α = "<>ToString[1-αα[[i]]],"(Subscript[Var, α])","(Subscript[CVaR, α])"}],Frame->All],{i,1,3}]

Novoroční úloha

Termín odeslání na moji adresu 2020-01-07.

hurt@karlin.mff.cuni.cz

Předně se podívejte na Hodnoceni_YC

Uvažujte kupónovou obligaci s těmito parametry:
F=1000, c=0.027 (kupónová sazba), roční splatnost kupónu. Emitovaná obligace má v okamžiku emise dobu do splatnosti 30 let. Cílem úlohy je ale prozkoumat (tj. spočítat) její spravedlivou cenu pro doby do splatnosti  0≤t≤30 [rok] při hodnotících úrokových mírách daných výnosovou křivkou. 
Použijte aktuální výnosovou křivku US  (treasury yield curve) dostupné z Wolfram|Alpha: “yield curves”.  V protokolu uveďte datum stažení této výnosové křivky. Pro konstrukci můžete případně použít lineární interpolaci.   Použijte kalendářní konvenci “30/360”.

Vypracujte protokol (slohové cvičení) podle následující osnovy (užívejte odstavce, tj. Section, Subsection apod.):

1. Formulace úlohy
2. Vstupní data
3. Metoda řešení (použité funkce, algoritmy)
4. Výstup: graf zobrazující závislost spravedlivé ceny obligace (dirty price) na době do splatnosti 0≤t≤30.

?*Interpolat*

⟶ More ⟶” +“ ⟶ Computable data

Příklad (Hodnocení pri užití úrokových měr ze spotové výnosové křivky 2)

(RefLink[FinancialBond, paclet:ref/FinancialBond]) can be given a term structure of interest rates or a list of future spot rates:

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->{1->.02,2->.03,3->.04,5->.045,10->.06},"Settlement"->2/12,"DayCountBasis"->"ActualISDA"}]

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->{.03,.03,.04,.05,.05},"Settlement"->4/12,"DayCountBasis"->"Actual/365"}]

When dates are given along with a rate structure that varies with time, time 0 is assumed to occur at the nearest point before the settlement date that represents an integral number of years from the maturity date. In this case, time 0 occurs on December 31, 2010:

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->{2020,12,31},"CouponInterval"->1/2},{"InterestRate"->{1->.02,2->.03,3->.04,5->.045,10->.06},"Settlement"->{2011,2,4},"DayCountBasis"->"ActualISDA"}]

2020-01-09 Cvičení 13

FAQ

Novoroční úloha: komentáře

Novoroční úloha: možná řešení

Zadání

Termín odeslání na moji adresu 2020-01-07.

hurt@karlin.mff.cuni.cz

Předně se podívejte na Hodnoceni_YC

Uvažujte kupónovou obligaci s těmito parametry:
F=1000, c=0.027 (kupónová sazba), roční splatnost kupónu. Emitovaná obligace má v okamžiku emise dobu do splatnosti 30 let. Cílem úlohy je ale prozkoumat (tj. spočítat) její spravedlivou cenu pro doby do splatnosti  0≤t≤30 [rok] při hodnotících úrokových mírách daných výnosovou křivkou. 
Použijte aktuální výnosovou křivku US  (treasury yield curve) dostupné z Wolfram|Alpha: “yield curves”.  V protokolu uveďte datum stažení této výnosové křivky. Pro konstrukci můžete případně použít lineární interpolaci.   Použijte kalendářní konvenci “30/360”.

Vypracujte protokol (slohové cvičení) podle následující osnovy (užívejte odstavce, tj. Section, Subsection apod.):

1. Formulace úlohy
2. Vstupní data
3. Metoda řešení (použité funkce, algoritmy)
4. Výstup: graf zobrazující závislost spravedlivé ceny obligace (dirty price) na době do splatnosti 0≤t≤30.

?*Interpolat*

⟶ More ⟶” +“ ⟶ Computable data

Řešení

Nejsou uvedena kalendářní data, čas je tedy abstraktní, kalendářní konvence v tomto případě nehraje roli.

Konstrukce YC

yc20200106All={WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",1},"ComputableData"},PodStates->{"TreasuryYieldCurve:EconomicData__More"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",2},"ComputableData"},PodStates->{"TreasuryYieldCurve:EconomicData__More"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",3},"ComputableData"},PodStates->{"TreasuryYieldCurve:EconomicData__More"}]}

yc20200106=Map[First,yc20200106All[[2,All,2,All]]]
times={1/12,1/4,1/2,1,2,3,5,7,10,20,30};
timesIR={times,(yc20200106/100)}//Transpose
interpolatedIR=Interpolation[timesIR,InterpolationOrder->1]
Row@{"Discrete IR:\n",discreteIR={Range[31]-1,Map[interpolatedIR[#]&,Range[31]-1]}//Transpose}
discreteIR//Mean
Plot[interpolatedIR[t],{t,1,30}]

SetDirectory[NotebookDirectory[]]
Export["YieldCurve.XML",discreteIR]
Export["YieldCurve.XLSX",discreteIR]

Built-in function FinancialBond

Dirty values, two months interval:

fbTable2Months=Table[{t,FinancialBond[{"FaceValue"->1000,"Coupon"->.027,"Maturity"->30},{"InterestRate"->discreteIR,"Settlement"->t,"DayCountBasis"->"30/360"},{"FullValue"}]},{t,0,29.999999,1/6}]
%//ListLinePlot

Dirty values, continuous (Takes approx. 600 sec):

FinancialBond[{"FaceValue"->1000,"Coupon"->.027,"Maturity"->30},{"InterestRate"->discreteIR,"Settlement"->30,"DayCountBasis"->"30/360"},{"FullValue"}]

Plot[FinancialBond[{"FaceValue"->1000,"Coupon"->.027,"Maturity"->30},{"InterestRate"->discreteIR,"Settlement"->t,"DayCountBasis"->"30/360"},{"FullValue"}],{t,0,29.9999999}]

User-defined function - for discussion

Procedura pro spravedlivou (dirty) cenu. Zde  invertovaný čas n=T-t.

Clear[dirtyvalue];
(dirtyvalue[c_,F_,n_,i_]=Module[{v},
c F  (v^(n-⌊n⌋))(If[⌊n⌋-n==0,0,1]+ (∑_(j=1)^(⌊n⌋))(v^j))+F (v^n)/.v->(1/(1+i))]//FullSimplify[#,Assumptions->i>-1]&)

Plot[dirtyvalue[.027,1000,30-n,.02],{n,0,30},Exclusions->None]

Clear[dirtyvalue];
(dirtyvalue[c_,F_,n_,i_]=Module[{v},
c F  (v^(n-⌊n⌋))(If[⌊n⌋-n==0,0,1]+ (∑_(j=1)^(⌊n⌋))(v^j))+F (v^n)/.v->(1/(1+i))]//FullSimplify[#,Assumptions->i>-1]&)

(c F  (v^(n-⌊n⌋))(If[⌊n⌋-n==0,0,1]+ (∑_(j=1)^(⌊n⌋))(v^j))+F (v^n)/.v->(1/(1+i))/.n->T-t)//Simplify

***User-defined function - simple solution

iRR=Rest[discreteIR][[All,2]]

Clear[pvBondJH]
pvBondJH[coupon_,face_,interests_,maturity_]:=Module[{discounFactors,iT},
iT=
coupon Total[(1/(1+interests))]+(face/(1+interests[[maturity]]))]

pvBondJH[iC,iF,{r1,r2,r3,r4},2]

Příklad (Hodnocení pri užití úrokových měr ze spotové výnosové křivky 2)

Uses forward IR!!!

(RefLink[EffectiveInterest, paclet:ref/EffectiveInterest]) can be used to determine the actual discount rates that were used by (RefLink[FinancialBond, paclet:ref/FinancialBond]) when a term structure is given:

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->{1->.02,2->.03,3->.04,5->.045,10->.06},"Settlement"->3/12,"DayCountBasis"->"ActualISDA"}]
FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->EffectiveInterest[{1->.02,2->.03,3->.04,5->.045,10->.06}],"Settlement"->3/12,"DayCountBasis"->"ActualISDA"}]
EffectiveInterest[{1->.02,2->.03,3->.04,5->.045,10->.06}]

forwards=EffectiveInterest[{1->r1,2->r2,3->r3,5->r5}]

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->forwards,"Settlement"->3/12,"DayCountBasis"->"ActualISDA"}]

(RefLink[FinancialBond, paclet:ref/FinancialBond]) can be given a term structure of interest rates or a list of future spot rates:

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->{1->.02,2->.03,3->.04,5->.045,10->.06},"Settlement"->4/12,"DayCountBasis"->"ActualISDA"}]

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->{.03,.03,.04,.05,.05},"Settlement"->4/12,"DayCountBasis"->"Actual/365"}]

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->EffectiveInterest[{.02,.03,.04,.045,.05}],"Settlement"->4/12,"DayCountBasis"->"Actual/365"}]

When dates are given along with a rate structure that varies with time, time 0 is assumed to occur at the nearest point before the settlement date that represents an integral number of years from the maturity date. In this case, time 0 occurs on December 31, 2010:

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->{2020,12,31},"CouponInterval"->1/2},{"InterestRate"->{1->.02,2->.03,3->.04,5->.045,10->.06},"Settlement"->{2011,2,4},"DayCountBasis"->"ActualISDA"}]

Místo lokálních pravidel vytvořit dvojice :

{1->.02,2->.03,3->.04,5->.045,10->.06}//FullForm
spot={1->.02,2->.03,3->.04,5->.045,10->.06}/.Rule->List

iEff={1->.02,2->.03,3->.04,5->.045,10->.06}//EffectiveInterest
plotEff=%//ListLinePlot;
ListLinePlot[{iEff,spot}]

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->{2020,12,31},"CouponInterval"->1/2},{"InterestRate"->{{0,0.02`},{1,0.04009803921568622`},{2,0.0602922047318315`},{3,0.05003611551179721`},{4,0.05506009598071371`},{5,0.06312968216296877`},{6,0.06918120639622805`},{7,0.0752411453612658`},{8,0.08130945067919182`},{9,0.08738607431974055`}},"Settlement"->{2011,2,4},"DayCountBasis"->"ActualISDA"}]

Úloha 18 (CAPM Odhad faktoru β a bezrizikové ú. m.)

CAPM Odhad faktoru β pro následujících 12 firem. Tržní portfolio representováno DJIA.

Řešení 18 (Náznak)

https://en.wikipedia.org/wiki/Dow_Jones_Industrial _Average #Calculation

SetDirectory[NotebookDirectory[]];
returns12companies=Get["associationReturns12companies.save"]
djiAll=Import["^DJI (1).csv"]

{djiAll[[2]],djiAll//Last}

returns12companies//Keys

djiAll//Length
Map[Length,returns12companies]
{date1=returns12companies["SAP"]//First,date2=returns12companies["SAP"]//Last};
{date1=date1[[1,1,1;;3]],date2=date2[[1,1,1;;3]]}

Úloha 18 (CAPM Odhad faktoru β a bezrizikové ú. m.)

CAPM Odhad faktoru β

***Různé

InputForm

applied to the desired output form shows how to construct it.

Probability

law3σ[distribution_?DistributionParameterQ]:=With[{m=Mean[distribution],s=StandardDeviation[distribution]},
Probability[Abs[x-m]>3.s,x\[Distributed]distribution]]
{law3σ[NormalDistribution[]],law3σ[StudentTDistribution[0,1,4.1]]}

ProbabilityDistribution (Asymmetric Laplace Distribution)

Asymmetric Laplace Distribution (vhodné pro modelování měr výnosnosti)

Clear[AsymmetricLaplaceDistribution]
AsymmetricLaplaceDistribution[θ_,τ_,κ_]=ProbabilityDistribution[((κ(Sqrt[2]))/(τ(1+(κ^2))))Which[x≥θ,Exp[-((κ(Sqrt[2]))/τ)(x-θ)],
			   x<θ,Exp[-((Sqrt[2])/(κ τ))(θ-x)]],{x,-∞,∞}]

(* Distribution function *)CDF[AsymmetricLaplaceDistribution[θ,τ,κ],x]

Manipulate [Plot[PDF[AsymmetricLaplaceDistribution[θ,τ,κ],y]//Evaluate,{y,-6,6},PlotPoints->70,PlotRange->All],{{θ,0},-5,5,Appearance->"Labeled"},{{τ,4},0.1,5,Appearance->"Labeled"},{{κ,1/2},0.1,10,Appearance->"Labeled"},SaveDefinitions->True]

AsymmetricLaplaceDistribution[θ,τ,κ]//Mean

Through[{Mean,Variance,Skewness,Kurtosis}[AsymmetricLaplaceDistribution[θ,τ,κ]]]

***Efficient portfolios

n=6;SeedRandom[13];
xx={returns=RandomInteger[{10,30},n],risks=RandomInteger[{10,30},n]}//Transpose
%//Sort

Clear[efficientPortfolios]
efficientPortfolios:=k;xx/.{a___,b_,c___,d_,e___}:> {a,b,c}/;(b[[1]]≥d[[1]]∧b[[2]]<d[[2]])
%//Sort

membersDJI20191018=FinancialData["^DJI","Members"]

SP50020191108=FinancialData["SP500",{2017,1,1}]

GE20191108=FinancialData["GE",{2017,1,1}]
%//DateListPlot

membersDJI20191018[[15]]
DateListPlot[FinancialData[membersDJI20191018[[15]],{2017,1,1}]]

FinancialData["^DJI",{2017,1,1}]

SP500from20070101to20190826=FinancialData["SP500",{2007,1,1}]

Map[{DateList[#[[1]]][[1;;3]],#[[2]]}&,SP500from20070101to20190826["Path"][[3090;;]]]//DateListPlot