ATFM Fall Term 2022/23

NMFM507 Advanced Topics of Financial Management

NMFM507 Pokročilé partie finančního managementu

Technical paragraph

chttp://www.karlin.mff.cuni.cz/~hurt/
chttp://www2.karlin.mff.cuni.cz/~hurt/

$LicenseExpirationDate

E-mails

2023-01-08

1. Chapter Risk measures (PDF): Typos on p. 21 corrected: symbol q at the beginning should be read as u. Other minor changes had been made.

2. Important illustrations, amendments, and minor changes in Mathematica notebook PPFM2022ZS.nb in 2022-12-21 Lesson 12 and 2023-01-04 Lesson 13.

3. Pay attention to the construction of portfolios briefly mentioned during the last lecture Lesson 13. Particularly, think about a construction of  a portfolio if you have a given payoff function (expressed in a graphical form).

4. From a plot of sample data of prices and/or returns you may be asked for to guess an appropriate type of the corresponding process and/or distribution.

Please, do not hesitate ask me in case of ambiguities or if you suffer from any other problem (if it is not a secrete, of course).

Introduction

Introductory and technical notes

1. The ideas and methods presented in this course are useful in a world with prevailing rational economic behavior. This is unfortunately not the case of the European Union and many other countries worldwide nowadays.

2a. Theoretical parts will be available at my web page in PDF files. Numerical illustrations will be presented in Wolfram Language (WL, often known as Mathematica), *.nb files. The use of WL chosen for this course enables to make interactive presentations, provide and analyze up to date information from financial markets and, in the main, efficiently to simply model varieties of of financial, insurance, and economic problems.

2b. Some knowledge of WL is an advantage but not the necessity. The installation instructions as well as simple hints for use of WL in this course without any previous knowledge of the system follow.

3. The students who did not pass the courses Introduction to Finance and Financial Management in Bc syllabus are recommended to have a look at the Basics of Financial Mathematics available here:

http://www2.karlin.mff.cuni.cz/~hurt/BasicsOfFinMat.pdf

It deals with basic notions, notation and terminology used in this course.  A small part of the material is written partly in Czech but the understanding of the essence of the text is language independent.

4. Installation instructions for Mathematica:

http://www2.karlin.mff.cuni.cz/~hurt/Mathematica_Installation_Students.pdf

5. Majority of students here understand Czech. Please do not hesitate to ask in Czech or Slovak but not in Swahili (I only understand “tutaonana” which means “See you later.”

How to read and understand Wolfram Language files

2. Designed for students who did not pass the course Mathematica for the beginners:
2.1. To study and analyze examples in the Mathematica notebook you need only the knowledge how to:
2.1.1. Open the cells (double click on the blue bracket of the cell right).Mathematica notebook is divided in Sections, Subsections, Subsubsections, Input cells (by default grouped with the results, output cells), ... .  These cells are automatically grouped in higher levels. See also the following Subsubsection Cell Groups and Outlining.
2.1.2. Initialization of the cells or starting calculations is realized by   ⇧+⌤ or ⌤ on the numerical keyboard. (On Mac ⇧+⏎) Symbols, variables, that are not defined in the current notebook or are unknown in the system shine blue.

Further very useful pieces of information are on the Slavík pages (only Czech):
https://www2.karlin.mff.cuni.cz/~slavik/info.html

and also: An Elementary Introduction to the Wolfram Language

Cell Groups and Outlining

double-click a group bracket — open or close a cell group
double-click cell brackets — close all other cells in a group
triple-click — select a function and its arguments
Subscript (Ctrl+ _): (Subscript[x, n])  ▪  Special Character  ▪  Convert to StandardForm  ▪  Why the Coloring?  ▪  
Superscript (Ctrl+ ^): (x^y)
Alt+7 — text cell (with default stylesheet) (Format ▶ Style ▶ ...)
Alt+4 — section-heading cell (with default stylesheet)

Exams

The terms will be  in SIS.  Preliminary terms:

2023-01-10, 2023-01-31, 2023-02-07: **AM Lecture room K*.
Exam in written, duration about 1 hour, 2 hours booking recommended.

References

[1] Dupačová, J., Hurt, J., Štěpán, J.: Stochastic Modeling in Economics and Finance. Kluwer Academic Publishers. Dordrecht 2002.

[2] Cipra, T.: Finanční ekonometrie. Ekopress, Praha 2008. (Czech)

[3] Cipra, T.: Matematika cenných papírů. Kamil Mařík Professional Publishing. Praha 2013. (Czech)

[4] Hurt, J.: Yield curves with Mathematica 6.0. In: Wolfram Technology Conference 2007. http://library.wolfram.com/infocenter/Conferences/6956/ . Champaign (IL) 2007.

[5] Hurt, J.: Risk measures in finance. In: 2008 International Mathematica User Conference. http://library.wolfram.com/infocenter/Conferences/7230/.  Champaign (IL) 2008.

[6] Hurt, J.: Risk measures in finance revisited. In: Wolfram Technology Conference 2010. http://library.wolfram.com/infocenter/Conferences/7861/  . Champaign (IL) 2010.

[7] Sciencedirect:  http://www.sciencedirect.com/

[8] Mathematica Cookbook

[9] Luenberger, D.: Investment Science. Oxford University Press. New York 1998.

[10] Hurt, J.: Finanční management. Přednáška MFF UK.

[11] Brigham, E. F.: Fundamentals of Financial Management. Dryden Press. 6th edition Fort Worth 1992. Chapter 7, The Cost of Capital.

[12] Cipra, T.: Financial and Insurance Formulas. Springer-Verlag. Berlin Heidelberg 2010.

[13] Cipra, T.: Riziko ve financích a pojišťovnictví: Basel III a Solvency II. Ekopress. Praha 2015. (Czech)

[14] Hurt, J.: NMFM507_Advanced_Topics_of_Financial_Management_Interactive_Demonstrations. To appear on my web page soon.

[17] McNeil, A. J., Frey, R., Embrechts, P.: Quantitative Risk Management. Princeton University Press. Princeton 2005.

[18] Witzany, J.: Credit Risk Management. Springer. Cham, 2017.

Paragraphs marked by * should be interpreted this way:

*   Exercise, understanding necessary.

**  Important information but its misunderstanding may not be fatal.

*** Inspiratory subjects or work in progress. May be useful in future problem solving.

***Final Exam Requirements  Financial & Insurance Mathematics

http://www.karlin.mff.cuni.cz/~hurt/20140521_Final _Exam _Requirements _MSc _Financial _ 2.docx

or in Czech http://www.karlin.mff.cuni.cz/~garant/stud/nmgr_ob_fpm _szz.shtml

20221-10-05 Lesson 1

FAQ

Q1

A1

Reminders

Prices and returns (definitions)

(Subscript[P, 0]),(Subscript[P, 1]),... (prices of an asset typically stocks, shares , ceny aktiva) in times 0,1,...  are nonnegative random variables (limited liability, princip omezeného ručení ) 
(Subscript[p, t])=ln (Subscript[P, t]) ... logarithmic price, arbitrary sign

(Subscript[R, t])=((Subscript[P, t+1])-(Subscript[P, t]))/(Subscript[P, t])=(Subscript[P, t+1])/(Subscript[P, t])-1 ... rate of return (ROR, míra výnosnosti, míra zisku), return (výnos), arbitrary sign

(Subscript[r, t]^✶)=(Subscript[p, t+1])-(Subscript[p, t]) ... difference of the logarithmic prices

1+(Subscript[R, t])=(Subscript[P, t+1])/(Subscript[P, t]) ... These random variables are nonnegative

(R^((T)))=(∏_(t=0)^(T-1))(1+(Subscript[R, t]))-1=(Subscript[P, T-1])/(Subscript[P, 0]) ...  return in period [0,T]

(((1+(R^((T)))))^(1/T))-1 ...  annualized return in period [0,T]  (geometric mean minus 1)

Theory of Portfolio 1

The whole chapter on Portfolio Theory available at

http://www2.karlin.mff.cuni.cz/~hurt/20221002_Theory_Of_Portfolio.pdf

Actual part end: p. 6 up to 1.3.1 Example

End 2022-10-05

2022-10-12 Lesson 2

FAQ

Q1

A1

Q2 Short sales (Prodeje nakrátko) - Illustration?

A2 Short sales (Prodeje nakrátko)

English

A simple example:

In the basics of portfolio theory we consider the time horizon one period.  The investor has a disposable fond W and realizes his trades on a stock exchange through his broker. At time 0  the investor borrows stocks from him supposing to give them back at time 1. The investor sells these assets for the price (Subscript[P, 0]) and thus his/her disposable funds increase to W+(Subscript[P, 0]). He/she invests this amount. At time 1 the investor has to buy stocks for the actual price (Subscript[P, 1])  and gives them back to the broker.

Short sale strategy is a speculative one.  The investor bets on the price decrease. He/she hopes that (Subscript[P, 1])<(Subscript[P, 0]) will take the place.  In such a case he/she gains  (Subscript[P, 0])-(Subscript[P, 1]) (we do not consider the time value of money for simplicity). The case  (Subscript[P, 1])>(Subscript[P, 0]) results in loss  (Subscript[P, 1])-(Subscript[P, 0]). Short sales are very risky, the possible loss may be unlimited.

Poznamenejme, že objem W, který investor hodlá investovat nebývá ve většině případů celý jeho majetek. Může tedy ztrátu hradit majetkem, který mu ještě zůstal. Makléři ještě navíc vyžadují nějakou pojistku (vinkulovaný vklad), podobně jako je tomu při pobchodování s futures aj., kde se skládají “marže” (margins, neplést s bankovními přirážkami) .

Czech

Zde podrobněji ilustruji na na jednoduchém příkladu.

Připomeňme, že v základech teorie portfolia uvažujeme časový horizont jedno období. Investor má bohatství (disponibilní fondy) W. Investor obchoduje na burze prostřednictvím svého makléře (broker). V čase 0 si od něj investor půjčí akcie s tím, že mu je vrátí v čase 1. Tyto akcie prodá za aktuální (promtní) cenu  (Subscript[P, 0]), čímž se mu zvýší disponibilní zdroje na W+(Subscript[P, 0]). Ty investuje. V čase 1 musí zakoupit akcie za cenu (Subscript[P, 1]) platnou v čase 1 aby je vrátil svému makléřovi.

Strategie prodejů nakrátko je spekulativní. Investor sází na pokles ceny. Doufá tedy, že bude (Subscript[P, 1])<(Subscript[P, 0]). V takovém případě má zisk (Subscript[P, 0])-(Subscript[P, 1]) (pro jednoduchost neuvažujeme časovou hodnotu peněz). V případě, že (Subscript[P, 1])>(Subscript[P, 0]), utrpí ztrátu (Subscript[P, 1])-(Subscript[P, 0]). Obchody nakrátko jsou velmi riskantní - ztráta může být neomezená, v řadě zemí jsou prozřetelně legislativně zakázány.

Poznamenejme, že objem W, který investor hodlá investovat nebývá ve většině případů celý jeho majetek. Může tedy ztrátu hradit majetkem, který mu ještě zůstal. Makléři ještě navíc vyžadují nějakou pojistku (vinkulovaný vklad), podobně jako je tomu při pobchodování s futures aj., kde se skládají “marže” (margins, neplést s bankovními přirážkami) .

Theory of Portfolio 2

http://www2.karlin.mff.cuni.cz/~hurt/20221002_Theory_Of_Portfolio_Notation.pdf

1.3.1 Example (Formulation in Wolfram Language)

Two assets, (Subscript[r, 1])=8[%], (Subscript[r, 2])=14[%], (Subscript[σ, 11])=9[((%)^2)], (Subscript[σ, 22])=36[((%)^2)], the correlation between the returns ρ so that (Subscript[σ, 12])=(Subscript[σ, 21])=18ρ[((%)^2)]. Doubled characters in commands mean vectors or matrices. Thus xx, rr, or VV in an input command should be read as x, r, or V, respectively.

r=(({{8}, {14}}))   [%]                   V=(({{9, 18ρ}, {18ρ, 36}}))  [((%)^2)]

Since x=((((Subscript[x, 1]),(Subscript[x, 2])))^T) is a portfolio, (Subscript[x, 1])+(Subscript[x, 2])=1 so we may simplify x=((((Subscript[x, 1]),1-(Subscript[x, 1])))^T). In WL notation:

Clear[x1,xx,ρ];
xx={x1,1-x1};rr={8,14}; VV={{9,18ρ},{18ρ,36}};

The expected return, variance, and the risk of the portfolio, respectively:

Row@{"(Subscript[r, p]) = ",rp=rr.xx//Simplify,"   (Subscript[σ, p]^2) = ",σp2=xx.VV.xx//Expand,"   (Subscript[σ, p]) = ",σp=(Sqrt[σp2])}

Here we seek for a minimum-variance portfolio(MVP) without any prescribed expected return (called global minimum-variance portfolio, GMVP)  so that we will minimize the variance (Subscript[σ, p]^2) (for |ρ|<1 positive definite quadratic form ⟹ equal its derivative to zero):

D[σp2,x1]==0

Row@{"GPMV portfolio: (Subscript[x, GMVP]^⋆) = ", (xxstar={(x1star=x1/.Solve[D[σp2,x1]==0,x1]//ExpandNumerator)//First,x2star=1-x1star//Simplify//First})//TraditionalForm}

The corresponding expected return, variance, and the risk of the optimal portfolio:

{xxstar.rr,xxstar.VV.xxstar,(Sqrt[xxstar.VV.xxstar])}//Simplify//TraditionalForm

Next we illustrate the risks and expected returns in the risk- expected-return plane for (Subscript[x, 1])∈[0,2].  (If you see the black window you must enable dynamic updating. (Blue button in the upper part of the screen or from Menu ⟶ Evaluation).)

Manipulate[ParametricPlot[{(Sqrt[36-72 x1+45 (x1^2)+36 x1 ρ-36 (x1^2) ρ]),14-6 x1},{x1,0,2},AxesOrigin->{0,(*8*),0},AspectRatio->3/4,AxesLabel->{Style["(Subscript[σ, p])",Italic,16],Style["(Subscript[r, p])",Italic,16]},AxesStyle->Directive[14],PlotPoints->50,PlotStyle->{Black,Thick,Dashing[0.01+0.01ρ]},PlotRange->{{0,7},{2,14}}, PlotLabel->Text@Row@{"(Subscript[x, 1]) = ",((2(ρ-2))/(4ρ-5)),"     (Subscript[σ, p]) = ",6(Sqrt[(((ρ^2)-1)/(4ρ-5))]),"\nRed point represents the global minimum variance portfolio."}, Epilog->{Red,PointSize[0.05],Point[{6(Sqrt[(((ρ^2)-1)/(4ρ-5))]),((44ρ-46)/(4ρ-5))}]}],{{ρ,0},-1,1,Appearance->"Labeled"},SaveDefinitions->True]

1.3.3 Remark

A more general case:

Clear[σ1,σ2,ρ,x1]
xx={x1,1-x1};(vv={{((σ1)^2),ρ σ1 σ2},{ρ σ1 σ2,((σ2)^2)}})//MatrixForm
derivative=D[xx.vv.xx,x1]//Simplify
Row@{"GMVP = ",{(x1opt=SolveValues[derivative==0,x1][[1]])//Simplify,x2opt=1-x1opt//Simplify},""}
Row@{"Particular cases:\nρ =  1 GMVP = ",({x1opt,x2opt}/.ρ->1)//Simplify,"\nρ = -1 GMVP = ",({x1opt,x2opt}/.ρ->-1)//Simplify}

End 2022-10-12

2022-10-19 Lesson 3

FAQ

Q1

A1

Theory of Portfolio 3

Data SP500, Coke, Coca-Cola, VW, Apple

Data of the 4 companies from

Actual Data of the 4 companies

For the next purposes we take data from last 5 years (monthly returns):

DateList[]

Row@{"Starting {YY,MM}: ",startDate=DateList[][[1;;2]]-{5,0},"   Ending {YY,MM}: ",endDate=DateList[][[1;;2]]}
companies={"SP500","Coke","KO","F:VOW","AAPL"};
DateObject[]
Text@Grid[Prepend[ Map[{#,FinancialData[#,"Name"],FinancialData[#]}&,companies],{"Ticker","Name","Actual Price"}],Frame->All]

Monthly returns and prices

DateObject[]
returnsMonthlyWithDates=Map[FinancialData[#,"Return",{startDate,endDate,"Month"},Method->"Legacy"]&,companies];
returnsMonthlyValues=Map[FinancialData[#,"Return",{startDate,endDate,"Month"},"Value",Method->"Legacy"]&,companies];
pricesMonthlyWithDates=Map[FinancialData[#,{startDate,endDate,"Month"},Method->"Legacy"]&,companies];

SetDirectory[NotebookDirectory[]];
Save["returnsMonthlyWithDates",returnsMonthlyWithDates]
Save["returnsMonthlyValues",returnsMonthlyValues]
Save["pricesMonthlyWithDates",pricesMonthlyWithDates]

(* This is only for the case if the direct download of the data from Wolfram fails. *)
(* Reading data from the files: *) 

SetDirectory[NotebookDirectory[]]
returnsMonthlyWithDates=Get["returnsMonthlyWithDates"];
returnsMonthlyValues=Get["returnsMonthlyValues"];
pricesMonthlyWithDates=Get["pricesMonthlyWithDates"];

Grid[{{"Expected returns",SpanFromLeft,SpanFromLeft},companies,rr=Map[Mean,returnsMonthlyValues]},Frame->All]
Grid[{{"Standard deviations",SpanFromLeft,SpanFromLeft},companies,Map[StandardDeviation,returnsMonthlyValues]},Frame->All]

We see that the selected companies did not fulfil the assumptions of the efficient market: Coke versus AAPL, e.g.

DateListPlot[returnsMonthlyWithDates,PlotLegends->companies,PlotRange->All,PlotStyle->Thick]

Recommendation: If possible, avoid a graphical representation of longer time series of returns. Usually  we are only able to identify peaks or irregularities on a plot. Graphs of prices are more instructive.

(* Plot of monthly prices *)
DateListPlot[pricesMonthlyWithDates//Rest,PlotLegends->(companies//Rest),PlotLabel->"Prices", PlotRange->All,PlotStyle->Thick]

Before we go further,  we look on the covariance structure of the returns. Since it is rather difficult to analyze the covariances, it is a good idea to look also on correlations and variances

covariancematrix=Covariance[Transpose[returnsMonthlyValues]];

So far rarely used is to figure out the covariance structure in the form of a covariance- correlation matrix. In such a matrix, the elements above the diagonal are covariances, the variances are on the diagonal, and the correlations are under the diagonal. (Note that the diagonal of a correlation matrix is redundant, it contains ones.) Correlations usually give better insight into the dependency structure of the observations.

(* Procedure for Colored Covariance-Correlation Matrix *)
Remove[coloredMatrix]
coloredMatrix[dataMatrix_(* rows observations, columns variables *)]:=
Module[{m=dataMatrix,mix,covupper,covlower,aux},
covupper=UpperTriangularize[Covariance[m]];
covlower=LowerTriangularize[Correlation[m],-1];
mix=covupper+covlower;
aux=ReplacePart[mix,{i_,i_}:>Style[mix[[i,i]],Red]];
aux=ReplacePart[aux,{i_,j_}:>Style[mix[[i,j]],Blue]/;i<j];
aux=ReplacePart[aux,{i_,j_}:>Style[mix[[i,j]],Green]/;i>j];
Column@{"Covariance\!\(TraditionalForm\`\(\(\\ \)\(-\)\)\) Correlation matrix:\n(Variances Red, Covariances Blue, Correlations Green)","\n",aux//MatrixForm}]

companies
coloredMatrix[Transpose[returnsMonthlyValues]]

Prepend[Prepend[(corrWithSP500=Correlation[returnsMonthlyValues//Transpose]),companies]//Transpose,{" ",companies}//Flatten]//Grid[#,Frame->All]&

1.3.1 General Solution (Risky Assets, Short sales Allowed)

companies

(i) Minimize risk with prescribed expected return

Define Lagrange function:

Clear[xx,rr,vv,μ,λ1,λ2,lagrange]
lagrange[xx_,rr_,vv_,μ_,λ1_,λ2_]:=(1/2)xx.vv.xx+λ1(1-Total[xx])+λ2(μ-rr.xx)

expectedreturns=Mean[Transpose[returnsMonthlyValues]];
(covariancematrix=Covariance[Transpose[returnsMonthlyValues]])//MatrixForm;
xx=Table[x[i],{i,1,Length[expectedreturns]}]
Grid[{companies,expectedreturns,Diagonal[covariancematrix]}//Transpose,Frame->All]

We see that this set of assets is not efficient (it does not fulfil some assumptions of an efficient market).

lag=lagrange[xx,expectedreturns,covariancematrix,μ,λ1,λ2]//Simplify
(* The derivative *)
derlag=Map[D[lag,#]&,xx]
xxλ=SolveValues[derlag==0,xx][[1]]

λ1 and λ2 are calculated from the initial conditions for arbitrary prescribed expected return:

{λ1opt,λ2opt}=Solve[{Total[xxλ]==1,xxλ.expectedreturns==μ},{λ1,λ2}][[1]]
xxoptimal=xxλ/.{λ1opt,λ2opt}//Simplify

Framed@ Row@{"MVP(μ) = ",xxoptimal,"\nExpected return = ",er=xxoptimal.expectedreturns//Simplify,"\nVariance(μ) = ",var=xxoptimal.covariancematrix.xxoptimal//Simplify,"\nRisk(μ) = ",(Sqrt[var])}

Plot[var,{μ,0,0.018},AxesLabel->{(Subscript[r, p]),(Subscript[σ, p]^2)},PlotLabel->"Representation of MVP",PlotStyle->Thick,AspectRatio->1]

⚠⚠⚠  Numerical complexity of linear problems (from the simplest to the most difficult): Solve a system of linear equations ⟶Find an inverse matrix ⟶ Calculate the determinant ⟶ Find the largest eigenvalue(s) ⟶Find all the eigenvalues. It means that if you can reach your goal by solving a system of linear equations then do it, ..., if it is sufficient to know the largest eigenvalue, use an algorithm  for finding it etc.

In theory we usually present the results int an explicit form using the inverse matrix in our cases, e. g. Solving the systems of linear equations  is more simple to perform, particularly if solving by hand.

$ImportFormats

(ii) Maximize Sharpe ratio

From Schwarz inequality it follows that the optimal portfolio is a multiple of the solution to the equation V,x=r.  (See 1.3.2 in PDF file.) Afterwards we just normalize to get (1^T)x=1:

xxSharpe=xx/.Solve[covariancematrix.xx==expectedreturns,xx][[1]]
xxSharpe=(xxSharpe/(Total[xxSharpe]))
%//Total

Framed@ Row@{"Sharpe = ",xxSharpe,"\nExpected return Sharpe = ",erSharpe=xxSharpe.expectedreturns//Simplify,"\nVariance Sharpe = ",varSharpe=xxSharpe.covariancematrix.xxSharpe//Simplify,"\nRisk Sharpe = ",(Sqrt[varSharpe])}

Compare with MVP with the same expected return, μ=0.052277668167643156`:

Framed@ Row@{"MVP(μ) = ",xxoptimal,"\nExpected return = ",er=xxoptimal.expectedreturns//Simplify,"\nVariance(μ) = ",var=xxoptimal.covariancematrix.xxoptimal//Simplify,"\nRisk(μ) = ",(Sqrt[var])}/.μ->0.052277668167643156`

(RowDefault["MVP(μ) = ", {-3.0233784726925395`,0.4068835894457534`,1.3355185719707465`,0.2827721966651403`,1.998204114610899`}, "\nExpected return = ", 0.052277668167643156`, "\nVariance(μ) = ", 0.018919075549744473`, "\nRisk(μ) = ", 0.13754663045580023`])

By the way, from the economical point of view it is nonsense to create such an portfolio. Nobody will lend you SP500 in such an enormous amount.

(iii) Maximize risk-adjusted expected return

Follow 1.1.3, p. 14 from PDF file. First define the Lagrange function, calculate its derivative, and equate it with 0.
⚠ Remind that here the Lagrange multiplier isγ while λ is an investor’s level of aversion!!!

Clear[λ,γ]
(* 1.20 *)
lagrangeRiskAdjusted=expectedreturns.xx-(1/2)λ xx.covariancematrix.xx+γ(1-Total[xx])//Simplify
(* 1.21 *)
lagRiskAdjusted=Map[D[lagrangeRiskAdjusted,#]&,xx]//Simplify
(* 1.22 *)
xxRiskAdjusted=NSolveValues[lagRiskAdjusted==0,xx ]//Simplify//First

(* Normalize i.e. find γ:*)
γnormed=NSolve[Total[xxRiskAdjusted]==1,γ]//Simplify//First

(* with risk aversion λ *)
xxoptimalRiskAdjusted=xxRiskAdjusted/.γnormed//Simplify

Prepend [Table[{λ,xxoptimalRiskAdjusted,xxoptimalRiskAdjusted.expectedreturns}//Flatten,{λ,0.01,300,10}],{"λ","SP500","Coke","KO","F:VOW","AAPL","(Subscript[r, p])"}]//Grid

Re

Optimal portfolio: risky assets only and short sales allowed (use of the built-in functions)

(* This is only for the case if the direct download of the data from Wolfram fails. *)
(* Reading data from the files: *) 

SetDirectory[NotebookDirectory[]]
returnsMonthlyWithDates=Get["returnsMonthlyWithDates"];
returnsMonthlyValues=Get["returnsMonthlyValues"];
pricesMonthlyWithDates=Get["pricesMonthlyWithDates"];

Grid[{{"Expected returns",SpanFromLeft,SpanFromLeft},companies,rr=Map[Mean,returnsMonthlyValues]},Frame->All]
Grid[{{"Standard deviations",SpanFromLeft,SpanFromLeft},companies,Map[StandardDeviation,returnsMonthlyValues]},Frame->All]
expectedreturns=rr;
vvWithSP500=covariancematrix=Covariance[Transpose[returnsMonthlyValues]];
Clear[x,xx,r];iN=Length[rr];xx=Array[x,iN]
portfolioreturn=rr.xx;
portfoliovariance=xx.vvWithSP500.xx//Simplify;

(i) minimize risk

Target expected return μ. In recent versions the problem may be solved explicitly.

Clear[optimalportfolio,optimalportfolioSD]
solution=Minimize[{portfoliovariance,portfolioreturn==μ,Total[xx]==1},xx]//Simplify
optimalportfolioSD[μ_]=(Sqrt[solution[[1]]])
optimalportfolio[μ_]=xx/.solution[[2]]

Manipulate[ParametricPlot[{optimalportfolioSD[μ],rr.optimalportfolio[μ]},{μ,Min[rr],Max[rr]+0.05},PlotStyle->Thick,AspectRatio->5/8,AxesLabel->{Style["(Subscript[σ, p])",Italic,14],Style["(Subscript[r, p])",Italic,14]},AxesStyle->Directive[12],PlotPoints->50,PlotRange-> {{0.039,0.08},{0.007,0.035}},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],optimalportfolio[μ]},Frame->All],"\nExpected return = ",rr.optimalportfolio[μ],"   Risk = ",optimalportfolioSD[μ],"\nRed point represents the above portfolio.\nGreen point represents the Sharpe portfolio."}, Epilog->{Red,PointSize[0.05],Point[{optimalportfolioSD[μ],rr.optimalportfolio[μ]}],Green,PointSize[0.05],Point[{0.0782,0.0327}]}],{{μ,Median[rr]},Min[rr],Max[rr],Appearance->"Labeled"},SaveDefinitions->True]

(ii) maximize Sharpe ratio

(* Sharpe ratio *)
solutionSharpe=Maximize[{(portfolioreturn/(Sqrt[portfoliovariance])),Total[xx]==1},xx]
Row@{"Expected return = ",erSharpe=(xxS=xx/.solutionSharpe[[2]]).rr,"   Risk = ",riskSharpe= (Sqrt[xxS.vvWithSP500.xxS])}

(* Square of Sharpe ratio *)
Maximize[{((((portfolioreturn))^2)/portfoliovariance),Total[xx]==1},xx]

(iii) maximize risk adjusted return

Clear[adjustedReturn]
adjustedReturn[rr_(*expected returns*),vv_(*covariance matrix*),λ_:0.0001(*risk aversion*)]:=
Module[{iN=Length[rr],xx=Array[x,iN],aux,solution},
aux=rr.xx-λ xx.vv.xx//Simplify;
aux=Maximize[{aux,Total[xx]==1},xx];
solution=xx/.aux[[2]];
Row@{"λ = ",λ," Optimal portfolio: ",solution,"\nExpected return: ",solution.rr,"   Risk: ",(Sqrt[solution.vv.solution])}
]
adjustedReturn[rr,vvWithSP500]
adjustedReturn[rr,vvWithSP500,10]

Scenarios:

Manipulate[ParametricPlot[{adjustedReturn[rr,vvWithSP500,λ1][[6]],adjustedReturn[rr,vvWithSP500,λ1][[4]]},{λ1,0.0001,40},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],adjustedReturn[rr,vvWithSP500,λ][[2]]},Frame->All],"\nExpected return = ",adjustedReturn[rr,vvWithSP500,λ][[4]],"   Risk = ",adjustedReturn[rr,vvWithSP500,λ][[6]],"\nRed point represents the above portfolio."}, Epilog->{Red,PointSize[0.05],Point[{adjustedReturn[rr,vvWithSP500,λ][[6]],adjustedReturn[rr,vvWithSP500,λ][[4]]}]}],{{λ,10},0.0001,40,Appearance->"Labeled"},SaveDefinitions->True]

End 2022-10-19

2022-10-26 Lesson 4

FAQ

Q1

A1

Theory of Portfolio 4

Optimal portfolio: risky assets only and short sales not allowed  (use of the built-in functions)

All the problems may be solved by optimization methods (mathematical programming, quadratic programming in this case: a quadratic objective function and linear constraints in the form of equalities or inequalities assuming nonnegative variables).

Preliminary calculations

companies={"SP500","Coke","KO","F:VOW","AAPL"};
(* The following download is only in the case if the direct download of the data from Wolfram fails. *)
(* Reading data from the files: *) 
SetDirectory[NotebookDirectory[]];
returnsMonthlyWithDates=Get["returnsMonthlyWithDates"];
returnsMonthlyValues=Get["returnsMonthlyValues"];
pricesMonthlyWithDates=Get["pricesMonthlyWithDates"];
expectedreturns=rr=Mean[Transpose[returnsMonthlyValues]];
(covariancematrix=vv=vvWithSP500=Covariance[Transpose[returnsMonthlyValues]]);
Row@{"Actual Date: ",returnsMonthlyWithDates//Last//Last//First//DateObject}
xx=Table[x[i],{i,1,Length[expectedreturns]}]
iN=Length[xx];
portfolioreturn=rr.xx//Simplify
portfoliovariance=xx.vv.xx//Simplify

(i) minimize risk

rr
rr//Sort

Target expected return μ. In the previous procedures we just add condition x>=0 to the constraints.

The output of the following calculations should be read this way:
If 0.007183567833648296`≤μ≤0.007496239848598976` then x[1]=11.977603352131746` -1528.1547562914243` μ,
if 0.007496239848598976`<μ≤0.007708965150447653` then x[1]=1.3517985117753657` -110.67011666342634` μ,
...,
if 0.0142622308653579`≤μ≤0.026461659008839532` then x[1]=0,
in other cases (True) x[1] is not defined,
if 0.02520985677872382`≤μ≤0.026461659008839532` then x[2]=3.536534908291886` -133.64751269414003` μ,
... .

Clear[optimalportfolio,optimalportfolioSD,μ]
solution=Minimize[{portfoliovariance,portfolioreturn==μ,Total[xx]==1,Thread[xx>=0]},xx]//Simplify;
optimalportfolioSD[μ_]=(Sqrt[solution[[1]]]);
optimalportfolio[μ_]=xx/.solution[[2]]

Manipulate[ParametricPlot[{optimalportfolioSD[μ],rr.optimalportfolio[μ]},{μ,Min[rr],Max[rr]+0.05},PlotStyle->Thick,AspectRatio->5/8,AxesLabel->{Style["(Subscript[σ, p])",Italic,14],Style["(Subscript[r, p])",Italic,14]},AxesStyle->Directive[12],PlotPoints->50,PlotRange-> {{0.039,0.091},{0.007,0.035}},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],optimalportfolio[μ]},Frame->All],"\nExpected return = ",rr.optimalportfolio[μ],"   Risk = ",optimalportfolioSD[μ],"\nRed point represents the above portfolio.\nGreen point represents the Sharpe portfolio."}, Epilog->{Red,PointSize[0.05],Point[{optimalportfolioSD[μ],rr.optimalportfolio[μ]}],Green,PointSize[0.05],Point[{0.0671,0.0279}]}],{{μ,Median[rr]},Min[rr],Max[rr],Appearance->"Labeled"},SaveDefinitions->True]

(ii) maximize Sharpe ratio

(* Sharpe ratio *)
solutionSharpe=Maximize[{(portfolioreturn/(Sqrt[portfoliovariance])),Total[xx]==1,Thread[xx>=0]},xx]
Row@{"Expected return = ",erSharpe=(xxS=xx/.solutionSharpe[[2]]).rr,"   Risk = ",riskSharpe= (Sqrt[xxS.vvWithSP500.xxS])}

(* Square of Sharpe ratio *)
Maximize[{((((portfolioreturn))^2)/portfoliovariance),Total[xx]==1,Thread[xx>=0]},xx]

(iii) maximize risk adjusted return

Remove[adjustedReturn]
adjustedReturn[rr_(*expected returns*),vv_(*covariance matrix*),λ_:0.0001(*risk aversion*)]:=
Module[{iN=Length[rr],xx=Array[x,iN],aux,solution},
aux=rr.xx-λ xx.vv.xx//Simplify;
aux=Maximize[{aux,Total[xx]==1,Thread[xx>=0]},xx];
solution=xx/.aux[[2]];
{"Optimal portfolio: ",solution,"   Expected return: ",solution.rr,"   Risk: ",(Sqrt[solution.vv.solution])}
]

Row@adjustedReturn[rr,vv,0]

Row@adjustedReturn[rr,vv,30]

ListLinePlot[Table[{adjustedReturn[rr,vv,λ1][[6]],adjustedReturn[rr,vv,λ1][[4]]},{λ1,0,20,1}],AxesLabel->{"(Subscript[σ, p])","(Subscript[r, p])"},PlotStyle->Thick,PlotRange->{Automatic,{0.00718357,0.0264617}}]

Scenarios:

1.3.4 General Solution (Riskless andRisky Assets, Short sales Allowed)

http://www2.karlin.mff.cuni.cz/~hurt/20221002_Theory_Of_Portfolio.pdf

pages 14 to19 in the above material

Efficient frontier

(* SetDirectory[NotebookDirectory[]];
Import["concepts7_optimal_portfolio.gif"] *)

Note that on the next picture the point representing the riskless return is above the expected return of the global minimum variance portfolio which is in the contradiction with an assumption of the efficient market.

Clear[Δ,A,B,iC]
σ[r_]=(Sqrt[(1/Δ)(A (r^2)-2 B r+iC)])

Simplify[D[σ[r],r],Assumptions->Δ>0]//TraditionalForm

End 2022-10-26

2022-11-02 Lesson 5

FAQ

Q1

A1

CAPM

The whole chapter on Capital Asset Pricing Model available at

http://www2.karlin.mff.cuni.cz/~hurt/20221026_CAPM.pdf

Example

Initialization

All cells should be initialized:

companies={"SP500","Coke","KO","F:VOW","AAPL"};
SetDirectory[NotebookDirectory[]];
returnsMonthlyWithDates=Get["returnsMonthlyWithDates"];
returnsMonthlyValues=Get["returnsMonthlyValues"];
pricesMonthlyWithDates=Get["pricesMonthlyWithDates"];
Row@{"Monthly returns from the period ",returnsMonthlyWithDates[[1,1,1]]//DateObject," to ",returnsMonthlyWithDates[[Length[companies]]]//Last//First//DateObject}

returnsAssociation=Association[Thread[companies->returnsMonthlyValues]];
Keys[returnsAssociation]
Map[Length,Values[returnsAssociation]]

DateListPlot[returnsMonthlyWithDates,PlotLegends->companies,PlotRange->All,PlotStyle->Thick]

Linear Regression Model

Market return = Return on S&P 500

linearModelsCompanies=Association[Thread[Rest[companies]->Map[LinearModelFit[Transpose[{returnsAssociation["SP500"],returnsAssociation[#]}],(Subscript[μ, M]),(Subscript[μ, M])]&,Rest[companies]]]]

Estimates of α, β, their standard errors, t-statistic, p-value, (R^2):

parameterTableCompanies=Map[{#["ParameterTable"],#["RSquared"]}&,linearModelsCompanies]

95 per cent confidence intervals for α, β (first row in the respective company for α, the second for β):

(parameterConfidenceIntervalsCompanies=Map[#["ParameterConfidenceIntervals"]&,linearModelsCompanies])//Dataset

parametersCompanies=Map[#["BestFitParameters"]&,linearModelsCompanies]

Recall CAPM  for asset n: (Subscript[r, n])-(Subscript[r, 0])=(Subscript[β, n])((Subscript[μ, M])-(Subscript[r, 0])). From the estimated regression equations (Subscript[r, n])=(OverHat[(Subscript[α, n])])+(OverHat[(Subscript[β, n])])(Subscript[μ, M]) , n="Coke","KO","F:VOW","AAPL" we may estimate the risk free rate as (OverHat[(Subscript[r, 0])])=(OverHat[(Subscript[α, n])])/(1-(OverHat[(Subscript[β, n])])):

Row@{"Risk free rate estimates (OverHat[(Subscript[r, 0])]):  ",(riskFreeRateEstimate=Map[((#[[1]])/(1-#[[2]]))&,parametersCompanies]),"\nAverage of the above:   " ,Mean[riskFreeRateEstimate] }

Histogram[riskFreeRateEstimate]

Security Market Line (SML)

Security Risk Evaluation, Merill Lynch, Beta book

Estimation of the respective α’s and β’s is based on 60 monthly observations (5 years history).

The following link is likely not valid nowadays:

http://newarkwww.rutgers.edu/RBS/ALPHA_Jan05.pdf

Indices 2004

SetDirectory[NotebookDirectory[]];
Import["20151103_Merrill_Lynch_Security_Risk_Evaluation_3.png"]

Indices 2006

SetDirectory[NotebookDirectory[]];
Import["20161102_Merrill_Lynch_Security_Risk_Evaluation_2007_3.png"]

Coke & other companies 2004

SetDirectory[NotebookDirectory[]];
Import["20151103_Merrill_Lynch_Security_Risk_Evaluation_63.png"]

Coke & other companies 2006

SetDirectory[NotebookDirectory[]];
Import["20161102_Merrill_Lynch_Security_Risk_Evaluation_2007_66.png"]

References concerning betas

http://pages.stern.nyu.edu/~adamodar/New_Home_Page/datafile/Betas.html

Excel file:

http://www.stern.nyu.edu/~adamodar/pc/datasets/betas.xls

Since the above server is not reliable, you can use the following link

http://www.karlin.mff.cuni.cz/~hurt/20211108_adamodar.betas.xls

SetDirectory[NotebookDirectory[]];
i1=Import["20211108_adamodar.betas.xls","XLS"];
i1[[2]][[10;;]]//TableForm

20211108_adamodar.betas.xls

Capital Market Line (CML)

End 2022-11-02

2022-11-09 Lesson 6 The Dean’s day

2022-11-16 Lesson 7

FAQ

Q1

A1

Arbitrage Pricing Theory (APT), Arbitrážní cenový model

The whole chapter on Arbitrage Pricing Theory available at

http://www2.karlin.mff.cuni.cz/~hurt/20221031_APT.pdf

*McNeil Factor Analysis (FA)

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

(sectorsDJI=Map[{#,FinancialData[#,"Name"],FinancialData[#,"Sector"]}&,membersDJI])//TableForm

FinancialData["PM",{"Name","Sector"}]

companiesMcNeil={"PM","KO","KODK","HPQ","INTC","MSFT","IBM","MCD","WMT","DIS"};
(compMcN=FinancialData[companiesMcNeil,{"Symbol","Name","Sector"}])//Grid[#,Frame->All]&

intersection=Intersection[membersDJI,compMcN[[All,1]]]
complement=Complement[compMcN[[All,1]],intersection]

sectorsMissinginDJI=FinancialData[complement,{"Symbol","Name","Sector"}]//Grid[#,Frame->All]&

(*
startDate={2014,1,1};endDate={2019,12,31};
returnsMonthlyValues=Map[FinancialData[#,"Return",{startDate,endDate,"Month"},"Value",Method->"Legacy"]&,companiesMcNeil];
pricesMonthlyWithDates=Map[FinancialData[#,{startDate,endDate,"Month"},Method->"Legacy"]&,companiesMcNeil];
pricesSP500=FinancialData["SP500",{startDate,endDate,"Month"},Method->"Legacy"];
*)

(*
SetDirectory[NotebookDirectory[]];
Save["returnsMonthlyValues20140101to20191231",returnsMonthlyValues]
Save["pricesSP50020140101to20191231",pricesSP500]
*)

SetDirectory[NotebookDirectory[]];
returnsMonthlyValues=Get["returnsMonthlyValues20140101to20191231"];
pricesSP500=Get["pricesSP50020140101to20191231"];

Append[returnsMonthlyValues,pricesSP500]//Dimensions
DateListPlot[pricesSP500,PlotLabel->"Prices ^SP500",PlotStyle->Thick,PlotTheme->"Detailed"]

startDate={2014,1,1};endDate={2019,12,31};

DateListPlot[FinancialData["NYSE:KODK",startDate],PlotLabel->"Prices of Eastman Kodak"]

Principal Components Analysis (PCA)

Original see [17] McNeil ..., p. 113. HWP = Hewlett-Packard, EK = Eastman Kodak, MO = Philil Morris.

Binarize[]

correlationMcNeil=rr=Correlation[returnsMonthlyValues//Transpose];
Row@{"Eigenvalues of the correlation matrix (variances of the PC):\n", eigenvaluesMcNeil=λλ=Eigenvalues[correlationMcNeil]}
eigenvectorsMcNeil=Eigenvectors[correlationMcNeil];
Row@{"Percentage of the total variability explained by the first m principal components:\n", eigenvaluesMcNeilAccumulated=10 Accumulate[eigenvaluesMcNeil]}

Test of the spectral decomposition (decomposition minus the correlation matrix should be zero matrix)

Transpose[eigenvectorsMcNeil].DiagonalMatrix[eigenvaluesMcNeil].eigenvectorsMcNeil-correlationMcNeil//Chop

Factor analysis (FA)

factorLoadings=(Sqrt[eigenvaluesMcNeil])*eigenvectorsMcNeil (* in rows *)

(* 3 common factors *)
(bb1=factorLoadings)[[1;;3]]

(* Orthogonal rotation to make the loadings of the first factor positive *)bb1=(factorLoadings.(-IdentityMatrix[Length[factorLoadings]]))[[1;;3]]

companiesMcNeil={"PM","KO","KODK","HPQ","INTC","MSFT","IBM","MCD","WMT","DIS"};GraphicsRow@{BarChart[bb1[[1]],BarOrigin->Left,PlotLabel ->"(1^st) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}],BarChart[bb1[[2]],BarOrigin->Left,PlotLabel ->"(2^nd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]}

BarChart[bb1[[3]],BarOrigin->Left,PlotLabel ->"(3^rd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]

***Iteration

(* Reduced corr matrix *)

rr1=rr-DiagonalMatrix[Diagonal[Transpose[bb1].bb1]]
λλ1=Eigenvalues[rr1]
xx1=Eigenvectors[rr1];

bb10=(Sqrt[λλ1])*xx1 (* in rows *)

bb2=bb10[[1;;3]].(-IdentityMatrix[Length[factorLoadings]])

GraphicsRow@{BarChart[bb2[[1]],BarOrigin->Left,PlotLabel ->"(1^st) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}],BarChart[bb2[[2]],BarOrigin->Left,PlotLabel ->"(2^nd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]}

BarChart[bb2[[3]],BarOrigin->Left,PlotLabel ->"(3^rd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]

2022-11-23 Lesson 8

FAQ

Q1

A1

Matching Assets and Liabilities (Sladění aktiv a pasiv)

The whole chapter on Matching Assets and Liabilities available at

http://www2.karlin.mff.cuni.cz/~hurt/20221031_ALM.pdf

End 2022-11-23

2022-11-30 Lesson 9

FAQ

Q1

A1

Interest Rates and Yield Curves 1

The whole chapter on Interest Rates and Yield Curves available at

http://www2.karlin.mff.cuni.cz/~hurt/20221123InterestRatesYC.pdf

Additional reference
[18] Witzany, J.: Credit Risk Management. Springer. Cham, 2017.

Decomposition of the Interest Rates (IR) (Determinants of the Interest Rate, Dekompozice úrokové míry)

Determinants of the IR

In a free economy, interest rates, as a price of money, are mainly determined by market supply and demand, and partly mastered by the government or central bank via money supply policy. Interest rates vary with economic environment, market position, used financial instrument, and time. The economic units which are willing to pay higher interest rates for the  funds
(= borrowed money in this case) expect higher returns on their investments. The returns are usually measured by the rate of return defined by:

Rate of return=((Ending price + Cash income - Beginning price)/(Beginning price)),

sometimes quoted in per cent.

Every investment should be valued from the point of view of return, risk, inflation, and liquidity. The firm with higher return will pay higher interest for funds (money). With the rate of return 25 per cent the firm pays 20 per cent interest with pleasure. Another firm, with the rate of return 20 per cent, would
not pay 20 per cent interest since then it would not have reason to develop any activity. More risky investment should be more expensive than an investment with (almost) certain return, in terms of interest rates. Inflation also makes funds more expensive. If the inflation is high, the funds may be not accessible. Short term funds (money borrowed for short time) are usually cheaper than long term funds (money borrowed for long time). Short term interest rates more or less reflect the actual state of the economy while long term interest rates reflect expectations, rational or less rational. Situation is more complicated, however, see the concept of yield curves next in this part. Denote r the rate of interest comprising all the factors mentioned above. In this context, r is also called cost of capital.

Taxation. Almost all incomes coming from investment are subject to taxes. The few exemptions are returns on some government or municipal bonds, e.g. Thus the taxation reduces the returns. Moreover, the taxes are often different for various types of investment and sometimes are progressive, i.e., the higher the return, the higher the taxes. Thus any investment should be carefully valued with respect to tax consideration.

In 2014 the United States the marginal corporate income tax rate was 39.1 percent (consisting of the 35% federal rate plus a combined state rate). See

https://en.wikipedia.org/wiki/Corporate_tax_in _the _United _States

Decomposition

1+r=(1+(Subscript[r, 0]))(1+(Subscript[r, infl])) (1+(Subscript[r, default]))  (1+(Subscript[r, liquid])) (1+(Subscript[r, mat]))

r ... quoted or nominal interest rate (IR) in a general sense (nominální, kótovaná úroková míra v zobecněném slova smyslu)

Here (Subscript[r, 0]) denotes the risk-free, riskless IR if we do not consider inflation, (dependent on domestic and international economical conditions, time preference of consumption),  (Subscript[r, inﬂ]) inflation premium (prémie za inﬂaci) expected rate of inflation, (Subscript[r, default]) default risk premium, credit risk premium (prémie za riziko nesplacení, kreditní riziko) is the premium charged for the default risk, that is the risk that the debtor will not pay either principal or interest or both, completely or partly, deliberately or unwittingly. The term  (Subscript[r, liquid]), liquidity premium (prémie za likviditu) stands for the risk that an asset in question is not readily convertible into cash without considerable cost. In other words, the respective asset cannot be easily converted into cash for its market value. Finally (Subscript[r, mat]) maturity risk premium (prémie za riziko v okamžiku splatnosti, riziko z vypršení) is the premium for the risk produced by possible changes of interest rates during the life of an asset. There are two types of the maturity risk. Consider bonds, e.g. For long-term bonds, it is the interest rate risk; if the market interest rate rises, the prices of bonds go down. This kind of premium appears when the interest rates are more volatile. For short-term bonds, it is the reinvestment rate risk; if these bills become due and the actual (market) interest rates are low, the reinvestment will result in interest income loss. It expresses the investor’s impotency what to do with the gained cash. This is also the reason why banks charge extra fees for the early repayment of a loan. At maturity the investment opportunities may substantially differ from those existing at the beginning.

Additive form

Sometimes the decomposition is given in the additive form

which is a good approximation if the components of r are sufficiently small since the cross-factors of type (Subscript[r, 0])(Subscript[r, infl]) are small of twice higher order than the original components.

Danger of the additive form

We have already seen that the only correct average of the interest rates over several periods is based on  the geometric mean. Here we give a brutal example showing that the arithmetic mean may lead to a nonsense.

Suppose the price development in times 0,1,2 is (Subscript[P, 0])=100,(Subscript[P, 1])=50,(Subscript[P, 2])=100. Then, according to the definition of  returns (Subscript[R, t]), (Subscript[R, 1])=((Subscript[P, 1])-(Subscript[P, 0]))/(Subscript[P, 0])=-1/2, (Subscript[R, 2])=((Subscript[P, 2])-(Subscript[P, 1]))/(Subscript[P, 1])=1. The arithmetic mean (Subscript[(OverBar[R]), A])=((Subscript[R, 1])+(Subscript[R, 2]))/2=1/4 or 25 per cent and the geometric mean calculated from (((1+(Subscript[(OverBar[R]), G])))^2)=(1+(Subscript[R, 1]))(1+(Subscript[R, 2])) is (Subscript[(OverBar[R]), G])=0 which sounds more realistic.

Riskless IR with expected inflation

Riskless IR  (Subscript[r, 0infl]) with the expected inflation:

(1+(Subscript[r, 0infl]))=(1+(Subscript[r, 0]))(1+(Subscript[r, infl]))

Example (calculation (Subscript[r, default]))

Companies rated AAA, AA, A in years 1984 and 1991. We consider default risk only.

(* in per cent *)
Clear[r,r0infl,rdefault];
Solve[1+(r/100)==(1+(r0infl/100))(1+(rdefault/100)),rdefault]

Grid[{{" ","(Subscript[r, 1991])","\!\(\*SubscriptBox[\(r\), \(default\\ 1991\)]\)","(Subscript[r, 1984])","\!\(\*SubscriptBox[\(r\), \(default\\ 1984\)]\)"},
{"T-bonds",8,0,12.34,0},
{"AAA",8.9,((100 (r-r0infl))/(100+r0infl))/.{r->8.9,r0infl->8},12.50,((100 (r-r0infl))/(100+r0infl))/.{r->12.50,r0infl->12.34}},
{"AA",9.1,((100 (r-r0infl))/(100+r0infl))/.{r->9.1,r0infl->8},12.70,((100 (r-r0infl))/(100+r0infl))/.{r->12.70,r0infl->12.34}},{"A",9.4,((100 (r-r0infl))/(100+r0infl))/.{r->9.4,r0infl->8},13.41,((100 (r-r0infl))/(100+r0infl))/.{r->13.41,r0infl->12.34}}},Dividers->{{True,True},{True,True}},Frame->True]

Rating

Example (Matrix of transition probabilities, transition matrix, matice pravděpodobností přechodu)

Credit Rating Transition Probabilities

Rating categories, credit classes (in Markov chains’ terminology states of the chain) Standard & Poor’s: AAA, AA, A, BBB, BB, B, C, D. (In practice, we meet more detailed subdivision: AA+, AA-, &c.). AAA means the best credit quality (extremely reliable with regard to financial obligations), AA very good credit quality (very reliable), ..., C close to or already bankrupt, D payment default on some financial obligation has actually occurred (bankrupt firms).

Source: Bluhm & al.: Credit Risk Modeling

spRatingCategories=Binarize[]

Let (Subscript[X, t]) be a random variable taking values AAA, AA, A, BBB, BB, B, CCC, D at time t. Transition matrix (matice pravděpodobností přechodu) among credit classes is a square matrix P=((Subscript[p, i,j])), i,j= AAA, AA, A, BBB, BB, B, CCC, D, where (Subscript[p, i,j])=P((Subscript[X, t+1])=j|(Subscript[X, t])=i).  Time is usually in years. In literature the last row is often missing because it is always (0,0,...,0,1). State (class) D is so called absorbing state (pohlcující stav). If a firm falls in D, it cannot quit.

If the transition probabilities remain constant in time (a homogenous Markov chain, e. g.) then the transition matrix from time t to time t+n is simply the n-th power (P^((n)))=(P^n).

Dimitrios Kavvathas 
Goldman Sachs Group, Inc.
Estimating Credit Rating Transition Probabilities for Corporate Bonds,
AFA 2001 New Orleans Meetings 
http://web.sakarya.edu.tr/~adurmus/basel_II/SSRN-id252517 %20 rating.pdf

The corresponding transition matrix P:

(transitionmatrix={{0.9446`,0.0477`,0.006`,0.001`,0.0007,0.0001`,0.`,0.`},{0.0049`,0.9179`,0.0694`,0.0059`,0.0007,0.0009`,0.0002`,0.0001`},{0.001`,0.0187`,0.9198999999999999`,0.052199999999999996`,0.0055000000000000005`,0.0025`,0.0001`,0.0001`},{0.0004`,0.0032`,0.0547`,0.8819`,0.0496`,0.0083`,0.0007,0.0011`},{0.0004`,0.0014,0.0084`,0.0675`,0.8388`,0.0713`,0.0072`,0.005`},{0.0001`,0.0007,0.0027`,0.006`,0.0562`,0.8479,0.0481`,0.0384`},{0.0012`,0.0001`,0.0051`,0.0046`,0.017`,0.0891`,0.5401`,0.3428`},{0.`,0.`,0.`,0.`,0.`,0.`,0.`,1.`}})//MatrixForm

The transition matrix after 3 years (P^3)

Map[PaddedForm[#,{4,4}]&,MatrixPower[transitionmatrix,3],{2}]//MatrixForm

The transition matrix after 5 years (P^5)

Map[PaddedForm[#,{4,4}]&,MatrixPower[transitionmatrix,5],{2}]//MatrixForm

MatrixPower[transitionmatrix,5]//MatrixForm

Additional reference
Witzany, J.: Credit Risk Management. Springer. Cham, 2017.

Pomocný materiál od odstavce 1.1.3 Časová struktura ... až do konce kapitoly:

https://www.karlin.mff.cuni.cz/~hurt/20201116InterestRatesYCT1.pdf

Term Structure of Interest Rates (Časová struktura úrokových měr)

Term structure of interest rates explains the relationships between the yields of comparable types of investment for different maturities. Usually it is applied to bonds but generally it may be used for an arbitrary investment.  By comparability we mean a set of corporate bonds rated AAA, say. The term structure is changing in time so we must clearly state the respective date τ to which the term structure relates. Sometimes we speak on interest rates implied in the term structure at the particular time τ. Here, for the sake of simplified notation we set τ=0.

(Subscript[R, t]) spot (okamžité),  (Subscript[f, t]) forward (forwardové) interest rates

(((1+(Subscript[R, t])))^t)=(∏_(j=1)^t)(1+(Subscript[f, j])),  t=1,...,T

Yield Curves (YC, Výnosové křivky)

YC ...  dependence of the interest rates on maturity (závislost úrokové míry na době splatnosti): r=r(t), or R=(Subscript[R, t]), plot of this function (graf této závislosti)

Generally, a yield curve (výnosová křivka, YC) plots interest rates paid on interest bearing securities against the time to maturity. Sometimes we also speak on the term structure of the interest rates. Usually it is applied to zero-coupon bonds but similarly it may be used for an arbitrary investment. We must emphasize that we have to consider only comparable investments. Yield curves differ both in time and with the type of investment. Thus at the same time we may plot yield curves for government zero coupon bonds for the maturities 1, 3, 6, 9, 12 months getting a completely different picture for AA rated firm’s bonds for the same maturities, symbolically Y,(Subscript[C, T-bills,2016])≠Y,(Subscript[C, AA,2016]). We should also take into account the risk factors (cf. decomposition of interest rate) and also comparable taxation conditions. Even for the same type of securities (like T-bills), the shape of the yield curve differs in time, i.e., the shape is different in years 2015 and 2016, say, ceteris paribus. Symbolically, for an AA rated company Y,(Subscript[C, AA,2015])≠Y,(Subscript[C, AA,2016]). This feature may be explained by many factors, like the change in spot riskless rate, inflation, and other exogenous factors. Another important feature is the internal need of the issuer for short, medium, or long financial funds.

(Subscript[R, t])-(Subscript[R, 1]) ...  yield spread (výnosové rozpětí)

***Definitions of NPV, D, C, IRR

Note that in last versions instead of the users’ defined functions we can make use of the built-in Timevalue &c.

(* npv, duration, convexity,  IRR *)
Clear[npv,dur,con,altIRR];
npv[cf_,i_]:=cf.(1/(1+i))^(Range[Length[cf]]-1)
dur[cf_,i_]:=(cf*(Range[Length[cf]]-1)).(1/(1+i))^(Range[Length[cf]]-1)/npv[cf,i]
con[cf_,i_]:=Module[{times=Range[Length[cf]]-1},(cf*times*(times+1)).(1/(1+i)^times)]/npv[cf,i]
altIRR[cf_]:=Module[{a},a=N[Solve[npv[cf,x]==0,x]];Cases[x/.a,z_Real/;z≥0]]
(* Profitability index *)
Clear[pi];
pi[cf_,i_]:=-npv[Join[{0},cf//Rest]//Flatten,i]/First[cf]//N
(* Test whether the project cf is a type I project, i. e., the outflows precede the inflows *)
Clear[typeI];
typeI[cf_]:=Module[{au},
	au=DeleteCases[Sign[cf],0];
	Sort[au]==au
				]
(* Discounted Payback method *)
Clear[payback];
payback[cf_,i_]:=Module[{dcf,au,ak,k},
If[Not[typeI[cf]],Print["    The project\n ",cf,  "\n    is not of type I (outflows do not precede inflows).   "],
	dcf=cf*(1/(1+i))^(Range[Length[cf]]-1);
	au=FoldList[Plus,0,dcf]//Rest;
	ak=Select[au,#>0&]//First;
	k=Position[au,ak]-1//Flatten;
	If[k=={},∞,k-1-au[[k]]/dcf[[k+1]]//N]]]
Clear[ieff];
ieff[inom_,n_]:=(1+inom/n)^n

Example (Comparison of the yield curves at different reference dates, Porovnání výnosových křivek v různých referenčních časech, Treasury YC USA)

December 15, 2014:

January 5, 2015:

Dec 10, 2015:

Dec 6, 2016:

Dec 5, 2017:

Notes:

Nov 21, 2018:

Dec 3, 2018:

Sep 18, 2019:

Feb 3, 2020:

Nov 17, 2020:

Mar 25, 2022:

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

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

Nov 28, 2022:

End 2022-11-30

2022-12-07 Lesson 10

FAQ

Q1

A1

Interest Rates and Yield Curves 2

***Technical notes

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

{WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",1},"FormattedData"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",2},"FormattedData"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",3},"FormattedData"}]}

{WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",1},"ComputableData"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",2},"ComputableData"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",3},"ComputableData"}]}

https://www.imf.org/external/pubs/ft/wp/2010/wp10258.pdf

Example (Oldřich Alfons Vasicek, EWGFM Oct 2010 Prague), in Czech

Na eficientním trhu nemůže být výnosová křivka konstantní. Předpokládejme, že je tomu naopak a r(t)=0.05=:r, 1≤t≤10, čas v rocích. Předpokládejme počáteční investici W=100. Můžeme si vybrat ze dvou možností:

Portfolio (investice) A: pětiletý bezkuponový dluhopis s výše uvedeným výnosem r=0.05. Odpovídající peněžní tok cfA je

cfA=Cashflow[{-100,0,0,0,0,100(((1+r))^5)}/.r->0.05]

Portfolio (investice) B: kombinace jednoročního bezkuponového dluhopisu ((Subscript[B, 1])) a desetiletého bezkuponového dluhopisu ((Subscript[B, 10])), přičemž taková kombinace má pětiletou průměrnou dobou splatnosti (průměrná doba splatnosti ne ve smyslu durace!). Nechť w je váha (Subscript[B, 1]) a tedy 1-w váha  (Subscript[B, 10]). w tedy musí splňovat rovnici

({{({{1,w+10(1-w)=5}})}})

w0=w/.Solve[1 w+10(1-w)==5,w]//First

Odpovídající peněžní tok cfB je tedy

Clear[cfB];cfB[w_,r_]:=Cashflow[{-100,w 100(1+r),0,0,0,0,0,0,0,0,(1-w)100(((1+r))^10)}]
cfBw0=cfB[w0,0.05]

Současné hodnoty při hodnotící úrokové míře i:

{npvA=TimeValue[cfA,i],npvBw0=TimeValue[cfBw0,i]}//Simplify

Rozdíl mezi současnou hodnotou portfolia B a současnou hodnotou portfolia A v závisloti na úrokové míře:

Plot[npvBw0-npvA,{i,0,0.1},PlotStyle->Thick,PlotTheme->"Detailed",PlotLabel->"PV(B)- PV(A)"]

Investice (portfolio) B poskytuje tedy vždy minimálně stejný výnos jako investice A. Na eficientním trhu nabídka a poptávka budou stlačovat cenu dluhopisu A dolů a cenu kombinace B nahoru.

Table[{i,npvBw0-npvA}/.i->i1,{i1,{0,0.03,0.04,0.05,0.06,0.07,0.1}}]//TableForm

!!! Example (Oldřich Alfons Vasicek, EWGFM Oct 29, 2010 Prague), in English

EWGFM: European Working Group on Financial Modelling

Here we show that on an efficient market the yield curve cannot be constant (flat). Suppose on the contrary that the (annual) yields for all maturities are constant and equal r(t)=r, 0≤t≤10 years. Let us consider three zero coupon bonds (Subscript[B, 1]), (Subscript[B, 5]), (Subscript[B, 10]) with maturities 1,5,10 years and yields 1+r, (((1+r))^5), (((1+r))^10), respectively. Let W be an initial investment. We create two portfolios. Portfolio A consists of (Subscript[B, 5]) only, portfolio B consists of a combination of (Subscript[B, 1]) and (Subscript[B, 10]) with weights w and 1-w such that the average time to maturity (not the duration!) of such a portfolio is 5 years, the same as of A=(Subscript[B, 5]). Hence the weights in portfolio B fulfil

({{({{1,w+10(1-w)=5}})}})

so that w=5/9.

We choose the following numerical values:

Current yields of all maturities: r=5% (flat yield curve)
Initial investment: W=$100
Portfolio A: 5-year bond
Portfolio B: Combination of 1-year bond and 10-year bond with 5 year average time to maturity ⟹ weights w and 1-w fulfil

({{({{1,w+10(1-w)=5}}).}})

so that w=5/9.

Cash flows for the two portfolios and general r and w:

Clear[cfA,cfB]
cfA[r_]:=Cashflow[{-100,0,0,0,0,100 (((1+r))^5) }]
cfB[r_,w_]:=Cashflow[{-100,100(1+r)w,0,0,0,0,0,0,0,0,100(((1+r))^10)(1-w)}]

The corresponding present values with general r and i (the valuation (market) interest rate) and w=5/9:

{TimeValue[cfA[r],i],TimeValue[cfB[r,5/9],i]}

The plots of the  present values dependent on the valuation (market) interest rate i:

Plot[({TimeValue[cfA[r],i],TimeValue[cfB[r,5/9],i]}/.r->0.05)//Evaluate,{i,0,0.1},PlotStyle->Thick,(*PlotTheme->"Detailed",*)AxesLabel->{"i","PV"} ,PlotLabel->"Present values of the portfolios A and B",PlotLegends->{"PV(A)","PV(B)"}]

The plot of the difference of the present values dependent on the valuation (market) interest rate i:

Plot[(TimeValue[cfB[r,5/9],i]-TimeValue[cfA[r],i]/.r->0.05)//Evaluate,{i,0,0.1},PlotStyle->Thick,(*PlotTheme->"Detailed",*)AxesLabel->{"i","PV"} ,PlotLabel->"PV(B)- PV(A)"]

We see that for any valuation interest rate i portfolio B outperforms portfolio A in the sense of higher present values. Numerically for selected values of i:

Consequently, in an efficient market supply and demand would push the price of (Subscript[B, 5]) down and the prices of (Subscript[B, 1]) and (Subscript[B, 10]) up. The yield curve would not stay flat.

***Exercise

Místo “durace” ve výše uvedeném příkladu použijte skutečnou duraci.

Instead of the average time to maturity in the above Example use the duration (the characteristic of a cash flow).

***Wolfram|Alpha Yield Curves (Technical Notes)

{WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",1},"FormattedData"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",2},"FormattedData"}],WolframAlpha["yield curves",{{"TreasuryYieldCurve:EconomicData",3},"FormattedData"}]}

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

Today

p20151202=ListLinePlot[{{1,0.14`},{3,0.19`},{6,0.39`},{12,0.52`},{24,0.93`},{36,1.25`},{60,1.66`},{84,2.01`},{120,2.23`},{240,2.64`},{360,3.`}},PlotStyle->Red,AxesLabel->{"Month","Yield"},PlotLabel->"US Yield Curve"]

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

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

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

p2014xxxx=ListLinePlot[{{1,0.02`},{3,0.04`},{6,0.11`},{12,0.22`},{24,0.6`},{36,1.03`},{60,1.58`},{84,1.9`},{120,2.12`},{240,2.45`},{360,2.74`}},PlotStyle->Thick,AxesLabel->{"Month","Yield"},PlotLabel->"US Yield Curve"]

Show[p20151202,p2014xxxx]

Statistical and Numerical Methods for Construction of the Yield Curves

N comparable securities 1,...,N,
observed, quoted yields (interest rates) (Subscript[y, 1]),...,(Subscript[y, N]) with maturities (Subscript[T, 1]),...,(Subscript[T, N])

Let T≠ (Subscript[T, n]) , n=1,...,N.

For the financial decision making it is necessary to know also yields for non-observable maturities (different of quoted ones) at least in the interval  [(Subscript[T, 1]),(Subscript[T, N])].

Numerical approach (interpolation)

The principle of the interpolation: Find a function y=f(t) (called interpolation function) which must fulfill (Subscript[y, n])=f((Subscript[T, n])), n=1,...,N.

Example (Beware of the axes origins in the plots!)

Data ir3 below are the real data in the form of pairs {maturity (in days), spot interest rate (in per cent)}, N=22.

plot3=ListPlot[(ir3={{7,4.74},{14,4.728},{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.49},{5110,5.51},{5475,5.54}}),PlotStyle->{Red,PointSize[0.015]}]

The simplest (and perhaps the best in most cases) is the linear interpolation, i.e., jointing the neighbouring points by line segments.

int31=Interpolation[ir3,InterpolationOrder->1];
plotint31=Plot[int31[t],{t,7,10000},PlotRange->All,PlotStyle->Thick];
Show[plot3,plotint31]

In detail:

The interpolation by polynomials results in a catastrophe:

poly3=InterpolatingPolynomial[ir3,t];
poly3plot=Plot[poly3,{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All,PlotStyle->Thick];
Show[plot3,poly3plot,PlotRange->All]

The interpolation polynomial is of order N-1=21 in our case. The coefficients of the polynomial are not so small as they look; take into account that t∈[7,290] so that (200^11)= 2.048*(10^25), e.g.

poly3//Expand//TraditionalForm

Interpolation by the polynomials of lower order, given by the option InterpolationOrder are illustrated in the following plot. We see that the smaller interpolation order the better results we get.

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,Appearance->"Labeled"},SaveDefinitions->True]

Check it by fitting by the least squares method. If the degree of the polynomial is number of observed pairs minus one, the resulting function is the interpolating polynomial.

(t^(Range[Length[ir3]]-1))

Fit[ir3,(t^(Range[Length[ir3]]-1)),t]

Show[Plot[Fit[ir3,(t^(Range[Length[ir3]]-1)),t]//Evaluate,{t,0,5475},PlotRange->All,PlotStyle->Thick],plot3]

***BEGIN
Skip it, not yet finished.

Table[Show[Plot[Fit[ir3,(t^k),t]//Evaluate,{t,0,5475},PlotRange->All,PlotStyle->Thick],plot3],{k,18,Length[ir3]}]

Manipulate[Show[Plot[Fit[ir3,(t^(Range[degreeOfPolynomial+1]-1)),t]//Evaluate,{t,0,5475},PlotRange->All,PlotStyle->Thick],plot3],{{degreeOfPolynomial,1},1,Length[ir3]-1,Appearance->"Labeled"},SaveDefinitions->True]

*** END

Statistical approach (fitting, smoothing, vyrovnání): parametric case

In what follows we use this structure of observed interest rates, time in days:

plot3=plotir3=ListPlot[(ir3={{7,4.74},{14,4.728},{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.49},{5110,5.51},{5475,5.54}}),PlotStyle->{Red,PointSize[0.015]}]

The principle of fiting: Find a function y=g(t) such that g((Subscript[T, n])) does not differ from (Subscript[y, n]) too much. The phrase “not to differ too much” is meant in sense of the sum of squares of the deviations (Subscript[y, n])-g((Subscript[T, n])) , e.g.

Example (fitting by a cubic polynomial by least squares method)

g(t; (Subscript[β, 0]), (Subscript[β, 1]),(Subscript[β, 2]),(Subscript[β, 3]))=(Subscript[β, 0])+(Subscript[β, 1])t+(Subscript[β, 2])(t^2)+(Subscript[β, 3])(t^3)

(Subscript[β, i]) unknown parameters

(cubic3fit=LinearModelFit[ir3,{t,(t^2),(t^3)},t])//Normal
Show[plot3,Plot[cubic3fit[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]}]]

***Example (Bradley-Crane model, BC model)

(Subscript[g, B,C])(t;α,β,γ)=α,(t^β)(e^(γ,t))

α, β, γ unknown parameters

This function is nonlinear in parameters so that we cannot us a linear regression model directly. We may, however, transform it to the function which is linear in parameters by a logarithmic transformation and use the methods for the linear regression:

(y^★)= log g(t;α,β,γ)=log α +β log ,t+γ,t, after reparametrization  (α^★):=log α it becomes

(y^★)=(α^★)+β log ,t+γ,t.

We estimate the parameters (α^★),β,γ and then we plug them into the original model:

y=exp((y^★))=exp((α^★)+β log ,t+γ,t)=α,(t^β)(e^(γ,t)).

(* Logarithms of the quoted yields (y^★) *)
(ir3log=Map[{#//First,Log[#//Last]}&,ir3])//Transpose//TableForm
(bradley3logfit=NonlinearModelFit[ir3log,αstar+β Log[t]+γ t,{αstar,β,γ},t]//Normal//Exp//Simplify)
Show[plot3,Plot[bradley3logfit,{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotStyle->Thick]]

Next we estimate the model directly using build-in function NonlinearModelFit. We have to use some initial estimates, without them the method does not converge. Here we use the initial estimates from the above linearized model.

(bradley3fit=NonlinearModelFit[ir3,α (t^β)(E^(γ t)),{{α,4},{β,-.01},{γ,.00005}} (* {α,β,γ} *),t])//Normal
Show[plot3,Plot[bradley3fit[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotStyle->Thick,PlotRange->All]]

Remark. We have seen that the BC model is unsatisfactory for the ir3 data. From the methodological point of view, it is good to see how to attack more complicated models than BC is, however.

!Example (Nelson-Siegel)

(Subscript[g, N,S])(t;(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),c)= ... (see PDF material).

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 b,0+b,1 at time 0 and has asymptote at level b,0; a local extreme is determined by parameter c.

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},PlotStyle->{Thick}],
{{b0,5.8,"b0"},5,6,0.1,Appearance->"Labeled"},
{{b1,-1,"b1"},-2,-0.5,0.1,Appearance->"Labeled"},
{{b2,-3,"b2"},-5,-1,0.2,Appearance->"Labeled"},
{{c,400,"c"},100,1000,10,Appearance->"Labeled"},SaveDefinitions->True
]

Remember: r ... spot rate, f ... forward rate, then

R(t)=(1/t)(∫_(,, 0)^t)f(τ)dτ

(the integral mean, integrální průměr)

Motivation from the discrete case:

({{({{log (1+(Subscript[R, t]))}})=(1/t)(∑_(j=1)^t)log(1+(Subscript[f, j]))}})

(* Forward curve *)
Clear[nelsonsiegelforward];
nelsonsiegelforward[b0_,b1_,b2_,c_,t_]:=b0+b1 (E^(-t/c))+b2 (t/c) (E^(-t/c))

Both spot and forward curves:

Manipulate[Plot[{nelsonsiegelspot[b0,b1,b2,c,t],nelsonsiegelforward[b0,b1,b2,c,t]},{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3},PlotStyle->{{Blue,Thick},{Orange,Thick}}],
{{b0,5.8,"b0"},5,6,0.1,Appearance->"Labeled"},
{{b1,-1,"b1"},-2,-0.5,0.1,Appearance->"Labeled"},
{{b2,-3,"b2"},-5,-1,0.2,Appearance->"Labeled"},
{{c,400,"c"},100,1000,10,Appearance->"Labeled"},SaveDefinitions->True
]

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,Thick}],Plot[curvenel3,{t,0,10000},PlotRange->All,PlotStyle->Thick],PlotRange->All]

!Example (Swensson)

This is a generalization of the Nelson-Siegel model. It has two more parameters.

(Subscript[g, S,w])(t;(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),(Subscript[β, 3]),(Subscript[c, 1]),(Subscript[c, 2]))= see the derivation below.

We will start from the simpler form of the forward rate:

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

***BEGIN Editing notes

b0+b1 (E^(-t/c1))+b2 (t/c1) (E^(-t/c1))+b3 (t/c2)(E^(-t/c2))/.{b0->(Subscript[β, 0]),b1->(Subscript[β, 1]),b2->(Subscript[β, 2]),c1->(Subscript[γ, 1]),b3->(Subscript[β, 3]),c2->(Subscript[γ, 2])}
%//TeXForm

***END Editing notes

Swensson spot rate obtained as the integral mean:

Clear[swenssonspot];
(swenssonspot[b0_,b1_,b2_,b3_,c1_,c2_,t_]=(1/t)(∫_0^t)swenssonforward[b0,b1,b2,b3,c1,c2,τ]\[DifferentialD]τ//Simplify)//TraditionalForm

swenssonspot[b0,b1,b2,b3,c1,c2,t]//Expand

The parameters have the following meaning: the function starts at level b,0+b,1 at time 0 and has asymptote at level b,0; a local extremes are determined by parameters c1, c2.

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

Manipulating spot:

Manipulate[Plot[swenssonspot[b0,b1,b2,b3,c1,c2,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3},PlotStyle->{Thick}],
{{b0,5.8,"b0"},5,6,0.1,Appearance->"Labeled"},
{{b1,-1,"b1"},-2,-0.5,0.1,Appearance->"Labeled"},
{{b2,-3,"b2"},-5,-1,0.2,Appearance->"Labeled"},
{{b3,-3,"b3"},-5,-1,0.2,Appearance->"Labeled"},
{{c1,200,"c1"},100,1000,10,Appearance->"Labeled"},
{{c2,5000,"c2"},100,10000,10,Appearance->"Labeled"},SaveDefinitions->True
]

Remeber: r ... spot rate, f ...forward rate, then

R(t)=(1/t)(∫_0^t)f(τ)dτ

Both spot and forward curves

Manipulate[Plot[{swenssonspot[b0,b1,b2,b3,c1,c2,t],swenssonforward[b0,b1,b2,b3,c1,c2,t]},{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3},PlotStyle->{{Blue,Thick},{Orange,Thick}}],
{{b0,5.8,"b0"},5,6,0.1,Appearance->"Labeled"},
{{b1,-1,"b1"},-2,-0.5,0.1,Appearance->"Labeled"},
{{b2,-3,"b2"},-5,-1,0.2,Appearance->"Labeled"},
{{b3,-3,"b3"},-5,-1,0.2,Appearance->"Labeled"},
{{c1,200,"c1"},100,1000,10,Appearance->"Labeled"},
{{c2,5000,"c2"},100,10000,10,Appearance->"Labeled"},SaveDefinitions->True
]

Fitting Swensson: on the graph, estimated spot is blue and estimated forward is orange.

plot3=ListPlot[(ir3={{7,4.74},{14,4.728},{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.49},{5110,5.51},{5475,5.54}}),PlotStyle->{Red,PointSize[0.015]}]

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,PlotStyle->Thick],Plot[curveswe3forward,{t,0,10000},PlotRange->All,PlotStyle->{Orange,Thick}],PlotRange->All]

!!Statistical approach: Splines - the principle

See PDF file.

***Splines: Procedure 1

Source: Adapted from

https://mathematica.stackexchange.com/questions/33206/implementation-of-smoothing-splines-function/33262

CubicSplSmooth[data_,lambda_]:=Module[{M,Knots,X,Dsq,a},M=Length@data;Knots=Flatten@{1,1,1,Range@M,M,M,M};X=Table[Evaluate@N@BSplineBasis[{3,Knots},n,t],{t,1,M},{n,0,M+1}];Dsq=Differences[X,2];a=LinearSolve[Transpose[X].X+lambda*Transpose[Dsq].Dsq,Transpose[X].data,Method->"Multifrontal"];Return[X.a]];

*** BEGIN

SeedRandom[249304];data=Table[{i,RiemannSiegelZ[i]+Sin[i]+RandomVariate[NormalDistribution[0,.2]]},{i,0,25,0.05}];
%//ListPlot

*** END

Manipulate[smoothdata=CubicSplSmooth[ir3,10^lambda];Show[ListLinePlot[smoothdata,Mesh->All,PlotStyle->Red],ListPlot[ir3,PlotRange->{4,6}]],{{lambda,1},-5,4,Appearance->"Labeled"},SaveDefinitions->True]

Splines: Procedure 2

Source: Adapted from

https://mathematica.stackexchange.com/questions/33206/implementation-of-smoothing-splines-function/33262

SmoothingSplineFunction[dat_?MatrixQ,p:(_?NumericQ|Automatic):Automatic]:=Module[{n=Length[dat],pv=p,cc,dc,del,h,qg,qm,rh,tm,uv,xa,ya},{xa,ya}=Transpose[dat];h=Differences[xa];rh=1/h;del=Differences[ya] rh;qm=SparseArray[{Band[{1,1}]->Most[rh],Band[{1,2}]->-ListCorrelate[{1,1},rh],Band[{1,3}]->Rest[rh]},{n-2,n}];tm=SparseArray[{Band[{2,1}]->Most[Rest[h]],Band[{1,1}]->ListCorrelate[{2,2},h],Band[{1,2}]->Drop[h,-2]},{n-2,n-2}];qg=qm.Transpose[qm];If[pv===Automatic,pv=1/(1+Tr[tm]/(6 Tr[qg]))];uv=LinearSolve[6 (1-pv) qg+pv tm,Differences[del]];dc=ya-6 (1-pv) Differences[ArrayPad[Differences[ArrayPad[uv,1]]/h,1]];Interpolation[Transpose[{List/@xa,dc,Append[Differences[dc]/h-h ListCorrelate[{2,1},ArrayPad[pv uv,1]],pv Last[uv] Last[h]-(Subtract@@Take[dc,-2])/Last[h]]}],InterpolationOrder->3,Method->"Hermite"]]

Manipulate[smoothdata=SmoothingSplineFunction[ir3,1/10^lambda];Show[ListLinePlot[smoothdata,Mesh->All,PlotStyle->Red],ListPlot[ir3,PlotRange->{4,6}]],{{lambda,1},0,14,Appearance->"Labeled"},SaveDefinitions->True]

***Statistical approach (vyrovnání, fitting, smoothing): Nonparametric Regression, Kernel Regression

See also Cambell, J. Y. et al.: Econometrics of Financial Markets, e.g.

We wish to estimate the relation between two variables X and Y based on the observations ((Subscript[X, 1]),(Subscript[Y, 1])),...,((Subscript[X, N]),(Subscript[Y, N])) satisfying the relationship

(Subscript[Y, n])=m((Subscript[X, n]))+(Subscript[ε, n]),   n=1,...,N.

where m is an unknown function, (Subscript[ε, n])’s are iid with zero mean. Obviously E(Y|X)=m(X).  A crucial role in kernel regression plays a kernel, here denoted K which is simply a symmetric probability density function. In practise we use a rescaling kernel

(Subscript[K, h])(x)=(1/h)K((x/h))

for some h>0 called bandwidth or window width (šířka okénka). Note that (Subscript[K, h]) is also a symmetric density. The Nadaraya-Watson kernel estimator of m is

(OverHat[m])(x)=(((∑_(n=1)^N)(Subscript[K, h])(x-(Subscript[X, n]))(Subscript[Y, n]))/((∑_(n=1)^N)(Subscript[K, h])(x-(Subscript[X, n])))).

The motivation for the above formula:

There are plenty of kernels and a vast number of recommendations in various situations. We will restrict our illustrations to two frequently used kernels: Gauss and Epanechnikov. These are also often used in kernel density estimation.

Recommended choice of the bandwidth h from the practical point of view :  look on plots and make experiments.

Recommended choice of the bandwidth h from the theoretical point of view by Silverman, often only of a theoretical interest. See also Help to

SmoothHistogram

Silverman, B. W. 1992. Density Estimation for Statistics and Data Analysis. London: Chapman & Hall. ISBN 9780412246203

The optimal width is the width that would minimize the mean integrated squared error if the data were Gaussian and a Gaussian kernel were used, so it is not optimal in any global sense. In fact, for multimodal and highly skewed densities, this width is usually too wide and oversmooths the density (Silverman 1992).
The Silverman optimal bandwidth h is:

h=((0.9m)/(n^(1/5)))    with    m=min((Sqrt[var X]),((IQR(X))/1.349)),

where n is the number of observations on X and IQR(X) its interquartile range.

Implementation:
Input: data in the form {{(Subscript[X, 1]),(Subscript[Y, 1])},...,{(Subscript[X, N]),(Subscript[Y, N])}}, kernel (a predefined function), h (bandwidth)
Output: formula (5)

Example.

(Now you can use built-in function SmoothKernelDistribution.)

Clear[kernelEstimator]
kernelEstimator[data_,h_,kernel_,x_]:=Module[{xx=data[[All,1]],yy=data[[All,2]]},
((Map[kernel,(x-xx)/h].yy)/(Total[Map[kernel,(x-xx)/h]]))]

Epanechnikov kernel K(u)=(3/4)(SuperPlus[(1-(u^2))])

Clear[kernelEpanechnikov]
kernelEpanechnikov[x_]:=Which[-1<x<1,(3/4)(1-(x^2)),
							   True,0]

Properties of Epanechnikov kernel:

(* Is it a density? *)
Integrate[kernelEpanechnikov[x],{x,-∞,∞}]

(* Moments *)
distrEpanechnikov=ProbabilityDistribution[kernelEpanechnikov[x],{x,-∞,∞}];
{Mean[distrEpanechnikov],Variance[distrEpanechnikov]}

kernelGauss[x_]=PDF[NormalDistribution[],x]

Plot[{kernelEpanechnikov[x],kernelGauss[x]},{x,-3,3},PlotTheme->"Detailed",PlotStyle->Thick]

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`}};
listPlotir3=ListPlot[ir3,PlotStyle->Red];

Manipulate[Show[Plot[{kernelEstimator[ir3,bwE,kernelEpanechnikov,x],kernelEstimator[ir3,bwG,kernelGauss,x]},{x,7,5475},PlotRange->{4.4,5.6},PlotStyle->{{Blue,Thick},{Orange,Thick}}],listPlotir3],
{{bwE,600,"bw Epanechnikov"},100,2000,10,Appearance->"Labeled"},
{{bwG,200,"bw Gauss"},100,1000,10,Appearance->"Labeled"},
SaveDefinitions->True]

A fruitful demonstration (despite not dealing with the interest rates) is

http://demonstrations.wolfram.com/EstimatingTheLocalMeanFunction/

Using the yield curves for investment appraisal
(Užití výnosových křivek pro hodnocení investic)

Based on the present value of the cash flow C,F using spot rates R=((Subscript[R, 1]),...,(Subscript[R, T]))  taken from the yield curve

({{({{PV(C,F,R)}})=(∑_(t=0)^T)((C,(Subscript[F, t]))/(((1+(Subscript[R, t])))^t))}})

End 2022-12-07

2022-12-14 Lesson 11

The whole chapter on Interest Return and Risk available at

http://www2.karlin.mff.cuni.cz/~hurt/20221214_ReturnRisk.pdf

FAQ

Q1

A1

Return and Risk (Výnos a riziko) 1

**Introduction

The procedures as well as numbered references are partly available in

Hyperlink[Button[" Risk Measures Revisited "],"http://www.karlin.mff.cuni.cz/~hurt/WTC2010JanHurtRiskMeasuresRevisited2withReference.nb"]

but the information in this notebook is self-contained.

Notation

(Subscript[S, 0]),(Subscript[S, 1]),... prices of a financial asset (random variables)

Rate Of Return (ROR, míra výnosnosti) or shortly return (výnos):

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

loss (ztráta):

(Subscript[L, t])=-(Subscript[ρ, t]).

In literature we can often find very misleading assumptions that returns or losses are always nonnegative, hidden in the assumptions about lognormal distribution of returns (losses) &c.  It is obviously nonsense (Explain why). There are very rare exceptions in insurance industry, however, that we seek for risk measures even for positive random variables. Nevertheless, the methodology is analogous.

In what follows we alternatively denote loss as L, (Subscript[L, t]), X, ... .

1.1 Selected probability distributions of losses

Data for Numerical Illustrations (Walt Disney, updated 2022-12-06)

Prices and returns

iDISprices20080101=FinancialData["DIS",{{2008,1,1},{2022,12,6}},Method->"Legacy"]
iDISreturns20080101=FinancialData["DIS","Return",{{2008,1,1},{2022,12,6}},Method->"Legacy"]
Row@{"Number of observations = ",%%//Length}
SetDirectory[NotebookDirectory[]]
Save["iDISprices20080101",iDISprices20080101];Save["iDISreturns20080101",iDISreturns20080101];

SetDirectory[NotebookDirectory[]]
iDISprices20080101=Get["iDISprices20080101"];
iDISreturns20080101=Get["iDISreturns20080101"];
iDISreturnsWithoutDates20080101=iDISreturns20080101[[All,2]];
GraphicsRow@{DateListPlot[iDISprices20080101,PlotLabel->"Prices Walt Disney",PlotTheme->"Marketing"],DateListPlot[iDISreturns20080101,PlotLabel->"Returns Walt Disney",PlotRange->All,PlotTheme->"Marketing"]}

SmoothHistogram[iDISreturnsWithoutDates20080101,PlotRange->All,PlotStyle->Thick,PlotLabel->"Returns of Walt Disney (kernel density estimator)"]

A general procedure for distribution characteristics of the returns

Input of this procedure is the list of desired distributions (even if it consists of just one distribution). The output are the estimates of the unknown parameters (making use of the built-in function FindDistributionParameters that by default provides the estimates by default maximum likelihood), the moments, and the plot of the density of  estimated distribution together with a normal density.

Clear[distributionCharacteristics]
distributionCharacteristics[listOfDistrubutions_,data_]:=Module[{distr, returnsParameters,mean,sd,sk,ku},
Table[Grid[{{ToString[distr],SpanFromLeft},
{"Estimated parameters",returnsParameters=FindDistributionParameters[data,distr]},
{"Mean, Standard Deviation,\nSkewness, Kurtosis",{mean,sd,sk,ku}={Mean[distr],StandardDeviation[distr],Skewness[distr],Kurtosis[distr]}/.returnsParameters},
{NormalDistribution[{mean,sd}/.List->Sequence],SpanFromLeft},
{Plot[{(PDF[distr/.returnsParameters,x])//Evaluate,PDF[NormalDistribution[{mean,sd}/.List->Sequence],x]},{x,-0.15,0.15},PlotStyle->Thick,PlotLegends->{ToString[distr[[0]]],"Normal"},PlotRange->All],SpanFromLeft}},Alignment->Center ,Frame->All],{distr,listOfDistrubutions}]//Column]
(* %//Row*)

1.1.1 Normal Distribution

Often it does not fit the financial data well. Symmetry and low kurtosis may serve as obstacles.

Clear[μ,σ]
Text@Grid[{{"Mean ",Mean[NormalDistribution[μ,σ]]},{"SD  ",StandardDeviation[NormalDistribution[μ,σ]]},{"Skewness     ",Skewness[NormalDistribution[μ,σ]]},{"Kurtosis",Kurtosis[NormalDistribution[μ,σ]]}},Frame->All,Alignment->Left]

Walt Disney data:

Needs["HypothesisTesting`"]
Row@{"95 per cent confidence interval for the mean of returns: ",MeanCI[iDISreturnsWithoutDates20080101]}
distributionCharacteristics[{NormalDistribution[μ,σ]},iDISreturnsWithoutDates20080101]

1.1.2 Lognormal Distribution

Clear[μ,σ]
Text@Grid[{{"Mean ",Mean[LogNormalDistribution[μ,σ]]},{"SD  ",StandardDeviation[LogNormalDistribution[μ,σ]]},{"Skewness     ",Skewness[LogNormalDistribution[μ,σ]]//TraditionalForm},{"Kurtosis",Kurtosis[LogNormalDistribution[μ,σ]]//TraditionalForm}},Frame->All,Alignment->Left]

1.1.3 Laplace distribution (double exponential)

Symmetric: higher, but constant kurtosis equal to 6.

Clear[μ,σ]
Row@{"Density of LaplaceDistribution[μ,σ]:\n", PDF[LaplaceDistribution[μ,σ],x]}
Plot[PDF[LaplaceDistribution[0,1],x],{x,-6,6}]

Text@Grid[{{"Mean ",Mean[LaplaceDistribution[μ,σ]]},{"SD  ",StandardDeviation[LaplaceDistribution[μ,σ]]},{"Skewness     ",Skewness[LaplaceDistribution[μ,σ]]},{"Kurtosis",Kurtosis[LaplaceDistribution[μ,σ]]}},Frame->All,Alignment->Left]

Walt Disney data:

distributionCharacteristics[{LaplaceDistribution[μ,σ]},iDISreturnsWithoutDates20080101]

***1.1.4 Asymmetric Laplace distribution (ALD)

Note: It is not necessary to know the formulas for the density &c., but it is necessary to know the approximate ranges of skewness, kurtosis, ... .

The Laplace Distribution and Generalizations
A Revisit with New Applications
Samuel Kotz
Department of Engineering Management & System Engineering
The George Washington University, Washington, DC 20052
Tomasz J. Kozubowski
Department of Mathematics
University of Nevada, Reno, NV 89557
Krzysztof Podg´orski
Department of Mathematical Sciences
Indiana University–Purdue University, Indianapolis, IN 46202.
January 24, 2001

See [21] Trindade for an introduction and statistical inference to ALD. Also the standardized ALD is discussed there.

See also http://www.mathestate.com/tools/Financial/s07b/Laplace.html

Probability density function (definition)

The density of ALD is

f(x;θ,τ,κ)=((κ (Sqrt[2]))/(τ(1+(κ^2))))exp{-((κ(Sqrt[2]))/τ)(x-θ)},    x≥θ,
f(x;θ,τ,κ)=((κ (Sqrt[2]))/(τ(1+(κ^2))))exp{-((Sqrt[2])/(κ τ))(θ-x)},        x<θ,

where real θ is a location parameter (parametr polohy), τ>0 is a scale parameter (parametr měřítka), and κ>0 is a shape parameter (parametr tvaru).

Definition of ALD in Mathematica and illustration

Clear[AsymmetricLaplaceDistribution];Unprotect[Mean, Variance, StandardDeviation, Skewness, Kurtosis,
         PDF, CDF, ExpectedValue, DistributionParameterQ,
         DistributionDomainQ,DistributionDomain,
         CharacteristicFunction,Quantile,RandomReal];
AsymmetricLaplaceDistribution/: PDF[AsymmetricLaplaceDistribution[θ_,τ_,κ_],x_]:=((κ(Sqrt[2]))/(τ(1+(κ^2))))Which[x≥θ,Exp[-((κ(Sqrt[2]))/τ)(x-θ)],
			   x<θ,Exp[-((Sqrt[2])/(κ τ))(θ-x)]]
AsymmetricLaplaceDistribution/: CDF[AsymmetricLaplaceDistribution[θ_,τ_,κ_],x_]:=Which[x≤θ,(((E^(((Sqrt[2]) (x-θ))/(κ τ))) (κ^2))/(1+(κ^2))),
						 x>θ,((1-(E^(((Sqrt[2]) (-x+θ) κ)/τ))+(κ^2))/(1+(κ^2)))];
AsymmetricLaplaceDistribution/: Mean[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=θ-(((-1+(κ^2)) τ)/((Sqrt[2]) κ))
AsymmetricLaplaceDistribution/: 
Variance[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=(((1+(κ^4)) (τ^2))/(2 (κ^2)))
AsymmetricLaplaceDistribution/: 
StandardDeviation[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=Sqrt[(((1+(κ^4)) (τ^2))/(2 (κ^2)))]
AsymmetricLaplaceDistribution/: 
Skewness[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=((2 (1-(κ^6)))/((((κ^4)+1))^(3/2)))
AsymmetricLaplaceDistribution/: 
Kurtosis[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=((9+6 (κ^4)+9 (κ^8))/(((1+(κ^4)))^2))
AsymmetricLaplaceDistribution/: Quantile[AsymmetricLaplaceDistribution[θ_,τ_,κ_],q_]:=
	Which[0<q≤((κ^2)/(1+(κ^2))),θ+((κ τ Log[q (1+(1/(κ^2)))])/(Sqrt[2])),
		  ((κ^2)/(1+(κ^2)))<q<1,θ-((τ Log[-(-1+q) (1+(κ^2))])/((Sqrt[2]) κ))]
AsymmetricLaplaceDistribution/: RandomReal[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=Quantile[AsymmetricLaplaceDistribution[θ,τ,κ],Random[]]
Protect[Mean, Variance, StandardDeviation, Skewness, Kurtosis,
         PDF, CDF, ExpectedValue, DistributionParameterQ,
         DistributionDomainQ,DistributionDomain,
         CharacteristicFunction,Quantile,RandomReal];

Illustrations

Plot[(Quantile[AsymmetricLaplaceDistribution[θ,τ,κ],q]/.{θ->0,τ->4,κ->1/2}),{q,0.01,0.99},AxesLabel->{"α","Quantile"},PlotLabel-> "Quantiles of ALD",PlotStyle->Thick]

Due to the three parameters the distribution is apparently flexible. Nevertheless its kurtosis is always high (≥6), which might be unsatisfactory in some practical situations.

{Skewness[AsymmetricLaplaceDistribution[θ,τ,κ]],Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]]}//Simplify//TraditionalForm

Text@Row[{"Minimum kurtosis is ", FindMinimum[Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]],{κ,2}]//First//Round, " and is achieved for κ = ",κ/.(FindMinimum[Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]],{κ,2}]//Last)//Round,"."}]

Plot[{Skewness[AsymmetricLaplaceDistribution[θ,τ,κ]],Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]]},{κ,0,3},AxesLabel->{"κ"," "},PlotStyle->Thick,PlotLegends->{"Skewness","Kurtosis"},AxesOrigin->{0,0},PlotRange->{-2,9},PlotStyle->Thick]

Probability density function (illustration)

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

Study specific parameters for ALD and simulation

In specific numerical illustrations we will consider the distributions with the same mean and standard deviation for comparisons. We have chosen AsymmetricLaplaceDistribution[0, 4, 1/2] as a standard (sometimes called reference distribution in the sequel). The corresponding moments are:

Text@Grid[{{"Mean ",Mean[AsymmetricLaplaceDistribution[0.,4,0.5]]},{"SD  ",StandardDeviation[AsymmetricLaplaceDistribution[0.,4,0.5]]},{"Skewness     ",Skewness[AsymmetricLaplaceDistribution[0.,4,0.5]]},{"Kurtosis",Kurtosis[AsymmetricLaplaceDistribution[0.,4,0.5]]}},Frame->All,Alignment->Left]

Simulation:

t1=Table[RandomReal[AsymmetricLaplaceDistribution[0.,4,0.5]],{10000}];
Text@Grid[{{"Mean ",Mean[t1]},{"SD  ",StandardDeviation[t1]},{"Skewness     ",Skewness[t1]},{"Kurtosis",Kurtosis[t1]}},Frame->All,Alignment->Left]
Histogram[t1]

1.1.5 Student t-distribution with location μ and scale σ parameters, ν degrees of freedom

Note that ν, the degrees of freedom need not be an integer as is usual in standard statistical methods!!! Symmetric distribution, arbitrary high kurtosis (or none)

Clear[μ,σ,ν]
PDF[StudentTDistribution[μ,σ,ν],x]//Simplify//TraditionalForm

Text@Grid[{{"Mean ",Mean[StudentTDistribution[μ,σ,ν]]},{"SD  ",StandardDeviation[StudentTDistribution[μ,σ,ν]]},{"Skewness     ",Skewness[StudentTDistribution[μ,σ,ν]]},{"Kurtosis",Kurtosis[StudentTDistribution[μ,σ,ν]]}},Frame->All,Alignment->Left]

Kurtosis[StudentTDistribution[μ,σ,4.1]]

Plot[Kurtosis[StudentTDistribution[μ,σ,ν]],{ν,4.1,20},PlotLabel->"Kurtosis of Student t",AxesLabel->{"degrees of freedom",""},PlotStyle->Thick]

Manipulate[Column[{
Text@Style[Row[{"Mean:          ",Mean[StudentTDistribution[0,1,ν]],"       Variance: ",Variance[StudentTDistribution[0,1,ν]]}],Bold,14],
Text@Style[Row[{"Skewness:  ",Skewness[StudentTDistribution[0,1,ν]],"       Kurtosis:   ",Kurtosis[StudentTDistribution[0,1,ν]]}],Bold,14],
Show[Plot[PDF[StudentTDistribution[0,1,ν],x],{x,-4,4},AxesOrigin->{0,0},PlotStyle->{Black,Thick},PlotRange->All,ImageSize->Large,PlotLabel->"The density of Student t-distribution (black)\nThe density of N(μ,σ) with the same mean and variance (red)"],Plot[PDF[NormalDistribution[μ,σ]/.{μ->Mean[StudentTDistribution[0,1,ν]],σ->StandardDeviation[StudentTDistribution[0,1,ν]]},x]//Evaluate,{x,-4,4},PlotStyle->{Red,Thick}]]}],
{{ν,5},2.1,30,Appearance->"Labeled"},SaveDefinitions->True]

Walt Disney data:

distributionCharacteristics[{StudentTDistribution[μ,σ,ν]},iDISreturnsWithoutDates20080101]

*1.1.7 Meixner distribution

Suppose m∈ℝ, a>0, d>0, –π<b<π.

Clear[a,b,m,d]
Row@{"Density:  ",PDF[MeixnerDistribution[a, b, m, d],x]}//TraditionalForm
Row@{"Support:   ",
DistributionDomain[MeixnerDistribution[a, b, m, d]]}//TraditionalForm

Text@Grid[{{"Mean ",Mean[MeixnerDistribution[a, b, m, d]]},{"SD  ",StandardDeviation[MeixnerDistribution[a, b, m, d]]},{"Skewness     ",Skewness[MeixnerDistribution[a, b, m, d]]},{"Kurtosis",Kurtosis[MeixnerDistribution[a, b, m, d]]}},Frame->All,Alignment->Left]

distributionCharacteristics[{MeixnerDistribution[a, b, m, d]},iDISreturnsWithoutDates20080101]

Kurtosis[MeixnerDistribution[a, b, m, d]]

Manipulation: Initial values of the parameters are those as estimated for iDISReturns20080101 data

Manipulate[Column[{
Text@Style[Row[{"Mean:          ",Mean[MeixnerDistribution[a, b, m, d]],"        Standard Deviation: ",StandardDeviation[MeixnerDistribution[a, b, m, d]]}],Bold,14],
Text@Style[Row[{"Skewness:  ",Skewness[MeixnerDistribution[a, b, m, d]],"        Kurtosis: ",Kurtosis[MeixnerDistribution[a, b, m, d]]}],Bold,14],
Show[Plot[PDF[MeixnerDistribution[a, b, m, d],x],{x,-0.15,0.15},AxesOrigin->{0,0},PlotStyle->{Black,Thick},PlotRange->All,ImageSize->Large,PlotLabel->(*"Červeně je vyznačena hustota N(μ,σ) se stejnými parametry."*)"The density of Meixner distribution (black)\nThe density of N(μ,σ) with the same mean and variance (red)"],Plot[PDF[NormalDistribution[μ,σ]/.{μ->Mean[MeixnerDistribution[a, b, m, d]],σ->StandardDeviation[MeixnerDistribution[a, b, m, d]]},x]//Evaluate,{x,-0.15,0.15},PlotStyle->{Red,Thick}]]}],
(* Estimated parameters for iDISReturns20080101 data: *) (* {a->0.05992549278174418`,b->0.02867076820572458`,m->0.0004443642663303386`,d->0.16738388235624066`} *)
{{a,0.0599},0.03,0.2,Appearance->"Labeled"},
{{b,0.02867},-2,2,Appearance->"Labeled"},
{{m,0.0004443},-1,1,Appearance->"Labeled"},
{{d,0.1673},0.1,3,Appearance->"Labeled"},SaveDefinitions->True]

1.1.8 Smooth Kernel Distribution

Let start with a usual histogram. If (Subscript[X, 1]),...(Subscript[X, n]) is a random sample from a continuous distribution with density f, then the histogram is simply the plot of the function

(Subscript[f, n])(x)=(1/(n,h))(∑_(i=1)^n)I(x-(h/2)<(Subscript[X, i])<x+(h/2))=(1/(n,,h))(∑_(i=1)^n)I(-,(1/2)<(((Subscript[X, i])-x)/h)<(1/2))

where I stands for an indicator function and h is the width of bin. This is a crude density estimator (check that the integral is equal to 1). Illustration:

Histogram[iDISreturnsWithoutDates20080101,{0.0025},"PDF"]

Instead of the indicator function we may take an arbitrary symmetric probability density function K called kernel. In practise we use a rescaling kernel

(Subscript[K, h])(x)=(1/h)K((x/h))

for some h>0 called bandwidth or window width (šířka okénka.  Note that (Subscript[K, h]) is also a symmetric density. The kernel density estimator of f is

(OverHat[f])(x)=(1/(n,h))(∑_(i=1)^n)K((((Subscript[X, i])-x)/h)).

With h approaching 0 the individual observations are stressed while for higher values the shape of the kernel is prevailing on the image.  On the basis of such a density we are able to define a probability distribution. This is done by the built-in function SmoothKernelDistribution which may be used as any other *Distribution. Possible choices of the bandwidth and the kernel (default is “Gaussian”; another frequently used is “Epanechnikov”) are in Help:

paclet:ref/SmoothKernelDistribution

Walt Disney data:

𝒟KernelDISreturnsWithoutDates20080101=SmoothKernelDistribution[iDISreturnsWithoutDates20080101]
Row@{"Support: ",DistributionDomain[𝒟KernelDISreturnsWithoutDates20080101]}
Row@{"Density f(x) = ",PDF[𝒟KernelDISreturnsWithoutDates20080101,x]}

Plot[PDF[𝒟KernelDISreturnsWithoutDates20080101,x],{x,-0.1,0.1},PlotRange->All]

Show[Plot[PDF[𝒟KernelDISreturnsWithoutDates20080101,x],{x,-0.1,0.1},PlotRange->All],Histogram[iDISreturnsWithoutDates20080101,Automatic,"PDF"]]

Use with defaults:

distributionCharacteristics[{𝒟KernelDISreturnsWithoutDates20080101},iDISreturnsWithoutDates20080101]

Here we use SmoothHistogram because it is faster than SmoothKernelDistribution and sufficient for our purposes.

GraphicsRow@{SmoothHistogram[iDISreturnsWithoutDates20080101,Automatic,"PDF",PlotLabel->"Returns Walt Disney"],SmoothHistogram[iDISreturnsWithoutDates20080101,Automatic,"PDF",PlotLabel->"Returns Walt Disney, PlotRange ⟶ All",PlotRange->All]}

Manipulate with bandwidth

Manipulate[SmoothHistogram[iDISreturnsWithoutDates20080101,{bw,"Epanechnikov"},"PDF",PlotLabel->"Returns Walt Disney",PlotRange->All,PlotStyle->Thick],{{bw,0.035},0.0001,0.3,Appearance->"Labeled"},  SaveDefinitions->True]

smoothDistr=SmoothKernelDistribution[iDISreturnsWithoutDates20080101]

Grid[{{"Mean ",Mean[smoothDistr]},{"SD  ",StandardDeviation[smoothDistr]},{"Skewness     ",Skewness[smoothDistr]},{"Kurtosis",Kurtosis[smoothDistr]}},Frame->All,Alignment->Left]

!!! 1.1.9 Mixture of Distributions (Směs rozdělení)

Let (Subscript[p, 1]),...,(Subscript[p, K]) be weights such that  (Subscript[p, k])≥0, (∑_(k=1)^K)(Subscript[p, k])=1 and let (Subscript[F, 1]),...,(Subscript[F, K]) be distribution functions. Then

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

is a distribution function and is called the distribution function of the mixture  (Subscript[F, 1]),...,(Subscript[F, K]) with the weights (Subscript[p, 1]),...,(Subscript[p, K]). The corresponding distributions is called a mixture distribution. Similarly for the densities:

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

Such a distribution is highly flexible.  Wide classes of distributions appearing in finance can be modelled with a mixture, particularly those with heavy tails (fat tails, těžké chvosty) or multimodal distributions. Problems with parameter’s estimation may arise, however.

A simple mixture consists of three normals with different means and the same variance and symmetric weigths ϵ,1-2ϵ,ϵ, 0<ϵ<1/2:

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

A (Mathematica^®) (user defined) function for this case:

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

Example of a trimodal distribution (Subscript[μ, 1])=-3, (Subscript[μ, 2])=0, (Subscript[μ, 3])=4 , σ=(2/3), ϵ=(1/20)

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

Use of the builtin function MixtureDistribution

Manipulate[Plot[PDF[MixtureDistribution[{ϵ,1-2ϵ,ϵ},{NormalDistribution[μ1,σ],NormalDistribution[μ2,σ],NormalDistribution[μ3,σ]}],x],{x,-5,5},PlotLabel->"The density of mixture",PlotRange->All,PlotStyle->Thick],
{{ϵ,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,PlotStyle->Thick,AxesOrigin->{0,0}, 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]

Very flexible, catches such phenomena like leptokurtosis. Highly recommendable!!!

𝒟mixture=MixtureDistribution[{1/20,17/20,2/20},{NormalDistribution[-3,1/2],NormalDistribution[],NormalDistribution[3.5,1/2]}];
Plot[PDF[𝒟mixture,x],{x,-5,5},PlotLabel->"The density of mixture\n",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]

Text@Grid[{{"Mean ",Mean[𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]]},{"SD  ",StandardDeviation[𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]]},{"Skewness     ",Skewness[𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]]},{"Kurtosis",Kurtosis[𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]]}},Frame->All,Alignment->Left]

Walt Disney data:

distributionCharacteristics[{𝒟mixtureOfNormal[ϵ,μ1,μ2,μ3,σ]},iDISreturnsWithoutDates20080101]

Example: Wolfram Demonstration

http://demonstrations.wolfram.com/ModelingReturnDistributions/

Exercise

Analyze carefully and comment whether the symmetric distributions used in the above demonstration fit properly the  DJI companies returns data.

End 2022-12-14

2022-12-21 Lesson 12

FAQ

Q1

A1

Return and Risk (Výnos a riziko) 2

1.2 Modelling of Prices and Returns

Introductory Notes

Calibration

Calibration (calibrating instruments) in financial literature is in fact estimation in statistical terminology.

A General Procedure for Characteristics of Stochastic Processes

This procedure should be initialized before further calculations:

In Mathematica, the Normal and related distributions have the second parameter standard deviation σ and not the variance (σ^2)!!!

Clear[characteristicsofstochasticprocess];characteristicsofstochasticprocess[process_,s_,t_]:=Module[{a},Text@Grid[{{"Process:   ",   Style["  "<>( process//HoldForm//ToString)<>"  ",Bold]},{"Mean   ",Mean[process[t]]//FullSimplify//TraditionalForm},
{"Variance   ",Variance[process[t]]//FullSimplify//TraditionalForm},
{"Skewness   ",Skewness[process[t]]//FullSimplify//TraditionalForm},
{"Kurtosis   ",Kurtosis[process[t]]//FullSimplify//TraditionalForm},
{"Covariance function   ",CovarianceFunction[process,s,t]//FullSimplify//TraditionalForm},{"Slice distribution   ",SliceDistribution[process,t]//FullSimplify//TraditionalForm}},Frame->All]]

1.2.1 Random Walk

Prices (Subscript[P, 0]),(Subscript[P, 1]),... (alternatively (Subscript[S, 0]),(Subscript[S, 1]),...)
Random walk hypothesis for logarithmic prices (Subscript[p, t])=log (Subscript[P, t]):

(Subscript[p, t+1])-(Subscript[p, t])=μ+(Subscript[ϵ, t+1]), t=0,1,...

ϵ's independent  identically distributed (iid) , E,(Subscript[ϵ, t])=0

Random walk with μ=0:

Clear[p,μ,σ,s,t];
characteristicsofstochasticprocess[RandomWalkProcess[p],s,t]

Clear[simulationRandomWalkProcess]
simulationRandomWalkProcess[p_,length_,numberOfTrajectories_,randomSeed_:20221221]:=
Module[{},SeedRandom[randomSeed];ListLinePlot[Table[RandomFunction[RandomWalkProcess[p],{0,length}],{numberOfTrajectories}],PlotStyle->Thick,
PlotLabel->Text@Row@{"Random Walk: Parameter p = ",p,"\nLength = ",length,"   No. of trajectories = ",numberOfTrajectories}]]

GraphicsGrid[{{simulationRandomWalkProcess[0.5,15,1],
simulationRandomWalkProcess[0.5,250,1]},
{simulationRandomWalkProcess[0.5,250,10],
simulationRandomWalkProcess[0.5,1000,10]}},ImageSize->1000]

The continuous version follows.

!Wiener process

?*Process

Clear[μ,σ,s,t];
characteristicsofstochasticprocess[WienerProcess[μ,σ],s,t]

RandomFunction is a discrete version of a continuous stochastic process. In case of Wiener process it is simply the random walk. Here is an example: time interval [0,1], step 0.01, results in 101 data points.

Simulation of one trajectory. For an ω given by the random number generator we plot a discretized version of W(t,ω) as a function of t.

SeedRandom[13131];ListLinePlot[RandomFunction[WienerProcess[μ,σ]/.{μ->1,σ->0.1},{0,1,0.01}],ImageSize->{350,300},PlotLabel->"Trajectory of Wiener process, drift = 1, volatility = 0.1"]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[WienerProcess[μ,σ],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,5},PlotTheme->"Detailed",PlotLabel->"Wiener process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0},-1,1,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.001,2,0.1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

Fifty trajectories, i. e.,  (Subscript[ω, 1]),...,(Subscript[ω, 50]):

SeedRandom[13131];ListLinePlot[RandomFunction[WienerProcess[μ,σ]/.{μ->1,σ->0.5},{0,1,0.01},50],ImageSize->{350,300},PlotStyle->Thick]

!1.2.3 Geometric Brownian Motion Process

First see the text for motivation, a discrete version,

The state X(t) of a geometric Brownian motion satisfies an Ito differential equation

dX(t)==μ X(t)dt +σ X(t) dW(t).

Here  σ>0, μ∈ℝ.

Clear[μ,σ,S0,s,t];
characteristicsofstochasticprocess[GeometricBrownianMotionProcess[μ,σ,S0],s,t]

SeedRandom[13131];ListLinePlot[RandomFunction[GeometricBrownianMotionProcess[μ,σ,(Subscript[S, 0])]/.{μ->1,σ->0.05,(Subscript[S, 0])->100},{0,1,0.01}],ImageSize->{350,300}]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[GeometricBrownianMotionProcess[μ,σ,(Subscript[S, 0])],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,300},PlotTheme->"Detailed",PlotLabel->"Geometric Brownian Motion\n"]],
Style["Time horizon",Bold],
{{T,5},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0},-1,1,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.001,2,0.1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
Style["Initial value (Subscript[S, 0])",Bold],
{{(Subscript[S, 0]),100},10,150,10,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

!1.2.4 Fractional Brownian Motion Process

Basics

Clear[μ,σ,S0,s,t,H];
characteristicsofstochasticprocess[FractionalBrownianMotionProcess[μ,σ,H],s,t]

ListLinePlot[Evaluate[Table[SeedRandom[131313];RandomFunction[FractionalBrownianMotionProcess[0.5,0.25,H],{0,1,0.01}],{H,{1/5,1/2,3/4}}]],ImageSize->Large,PlotTheme->"Detailed",PlotStyle->Thick,PlotLegends->{"H=1/5","H=1/2","H=3/4"},PlotLabel->"Fractional Brownian Motion"]

Details

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[FractionalBrownianMotionProcess[μ,σ,h],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,5},PlotTheme->"Detailed",PlotLabel->"Fractional Brownian Motion Process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0},-1,1,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.001,2,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Hurst coefficient",Bold],
{{h,0.5},0.05,.95,0.1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

!1.2.6 Vasicek model (Ornstein–Uhlenbeck process)

The state x(t) of an Ornstein–Uhlenbeck process satisfies an Ito differential equation dx(t)==θ (μ-x(t))dt +σ dw(t), where w follows a standard (RefLink[WienerProcess, paclet:ref/WienerProcess])[].

Basics

Clear[μ,σ,S0,s,t,H];
characteristicsofstochasticprocess[OrnsteinUhlenbeckProcess[μ,σ,θ,r0],s,t]

Example (interest rates in per cents per month):

characteristicsofstochasticprocess[OrnsteinUhlenbeckProcess[1.5,2,2,1.5],s,t]

ListLinePlot[Evaluate[Table[SeedRandom[131313];RandomFunction[OrnsteinUhlenbeckProcess[1.5,σ,2,1.5],{0,24,1}],{σ,{1,2,4}}]],ImageSize->Large,PlotTheme->"Detailed",PlotStyle->Thick,PlotLegends->{"σ=1","σ=2","σ=4"},PlotLabel->"Vasicek model",PlotRange->All]

Details

The state X(t) of a Ornstein–Uhlenbeck process satisfies an Ito differential equation

dX(t)==θ (μ-X(t))dt +σ dW(t),

with long-term mean μ, volatility σ, and mean reversion speed θ. Part μ-X(t) is the mean reverting drift pulled to a level μ at rate θ. Here θ>0, σ>0, μ∈ℝ. A drawback of the solution is that X(t) may attain negative values.

In more details (beware of the different notation):

Oldrich Alfons Vasicek (1977): spot rate r(t), α>0

dr=α(γ-r)dt+ρdz

In contrast to the random walk (the Wiener process), which is an unstable process
and after a long time will diverge to infinite values, the Ornstein-Uhlenbeck
process possesses a stationary distribution. The instantaneous drift α(γ-r)
represents a force that keeps pulling the process towards its long-term mean γ
with magnitude proportional to the deviation of the process from the mean.
The stochastic element, which has a constant instantaneous variance (ρ^2), causes
the process to fluctuate around the level γ in an erratic, but continuous, fashion.
The conditional expectation and variance level are ...

Clear[μ,σ,s,t,θ,(Subscript[r, 0])(* initial value of arbitrary sign *)];
characteristicsofstochasticprocess[OrnsteinUhlenbeckProcess[μ,σ,θ,(Subscript[r, 0])],s,t]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[OrnsteinUhlenbeckProcess[μ,σ,θ,r0],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,5},PlotTheme->"Detailed",PlotLabel->"Ornstein Uhlenbeck Process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0.04},0.01,0.11,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.01,2,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Rate θ",Bold],
{{θ,0.2},0.05,0.5,0.05,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
Style["Initial value (Subscript[r, 0])",Bold],
{{r0,0.02},0.01,0.3,0.001,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

! 1.2.7 Cox-Ingersoll-Ross Process (CIR)

CoxIngersollRossProcess

The state X(t) of a Cox-Ingersoll-Ross process (CIR) satisfies an Ito differential equation

dX(t)==θ (μ-X(t))dt+σ (Sqrt[X(t)]) dW(t),

with long-term mean μ, volatility σ, speed of adjustment θ, and (possible) initial condition (Subscript[x, 0])>0. Here σ>0, θ μ>0, usually μ>0. The last condition means that both θ and μ must be nonzero and of the same sign. As compared with the Ornstein-Uhlenbeck, CIR can never attain negative values. It can be shown that X(t)>0 a. s. if (σ^2)≤2θ,μ.

Clear[μ,σ,s,t,θ,r0(* initial value positive *)];
characteristicsofstochasticprocess[CoxIngersollRossProcess[μ,σ,θ,r0],s,t]

Limit of the variance:

Limit[(((σ^2) (E^(-2 θ t)) ((E^(θ t))-1) (2 r0+μ ((E^(θ t))-1)))/(2 θ)),t->∞,Assumptions->θ>0&&r0>0&&μ>0]

Skewness and kurtosis:

Text@Row@{"The limits of skewness:  ",Limit[(((Sqrt[2]) (3 r0+μ ((E^(θ t))-1)))/((2 r0+μ ((E^(θ t))-1)) (Sqrt[((θ (μ+((2 r0)/((E^(θ t))-1))))/(σ^2))]))),t->∞,Assumptions->θ>0&&r0>0&&μ>0]//TraditionalForm,Spacer[20],"and kurtosis:   ",Limit[((3 (σ^2) (1-((4 (r0^2))/(((2 r0+μ ((E^(θ t))-1)))^2))))/(θ μ))+3,t->∞,Assumptions->θ>0&&r0>0&&μ>0]}

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[CoxIngersollRossProcess[μ,σ,θ,r0],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-0.05,0.05},PlotTheme->"Detailed",PlotLabel->"Cox Ingersoll Ross Process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0.06},0.01,0.11,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,0.05},0.01,1,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Rate θ",Bold],
{{θ,0.15},0.05,0.5,0.05,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
Style["Initial value (Subscript[r, 0])",Bold],
{{r0,0.02},0.01,0.1,0.001,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

** Calculation of Risk Measures

**General Formulas for Calculations

Asymmetric Laplace distribution (ALD)

See [21] Trindade for an introduction and statistical inference to ALD. Also the standardized ALD is discussed there.

See also http://www.mathestate.com/tools/Financial/s07b/Laplace.html

Probability density function (definition)

The density of ALD is

f(x;θ,τ,κ)=((κ (Sqrt[2]))/(τ(1+(κ^2))))exp{-((κ(Sqrt[2]))/τ)(x-θ)},    x≥θ,
f(x;θ,τ,κ)=((κ (Sqrt[2]))/(τ(1+(κ^2))))exp{-((Sqrt[2])/(κ τ))(θ-x)},        x<θ,

where real θ is a location parameter, τ>0 is a scale parameter, and κ>0 is a shape parameter.

Definition of ALD in Mathematica

Clear[AsymmetricLaplaceDistribution];Unprotect[Mean, Variance, StandardDeviation, Skewness, Kurtosis,
         PDF, CDF, ExpectedValue, DistributionParameterQ,
         DistributionDomainQ,DistributionDomain,
         CharacteristicFunction,Quantile,RandomReal];
AsymmetricLaplaceDistribution/: PDF[AsymmetricLaplaceDistribution[θ_,τ_,κ_],x_]:=((κ(Sqrt[2]))/(τ(1+(κ^2))))Which[x≥θ,Exp[-((κ(Sqrt[2]))/τ)(x-θ)],
			   x<θ,Exp[-((Sqrt[2])/(κ τ))(θ-x)]]
AsymmetricLaplaceDistribution/: CDF[AsymmetricLaplaceDistribution[θ_,τ_,κ_],x_]:=Which[x≤θ,(((E^(((Sqrt[2]) (x-θ))/(κ τ))) (κ^2))/(1+(κ^2))),
						 x>θ,((1-(E^(((Sqrt[2]) (-x+θ) κ)/τ))+(κ^2))/(1+(κ^2)))];
AsymmetricLaplaceDistribution/: Mean[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=θ-(((-1+(κ^2)) τ)/((Sqrt[2]) κ))
AsymmetricLaplaceDistribution/: 
Variance[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=(((1+(κ^4)) (τ^2))/(2 (κ^2)))
AsymmetricLaplaceDistribution/: 
StandardDeviation[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=Sqrt[(((1+(κ^4)) (τ^2))/(2 (κ^2)))]
AsymmetricLaplaceDistribution/: 
Skewness[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=((2 (1-(κ^6)))/((((κ^4)+1))^(3/2)))
AsymmetricLaplaceDistribution/: 
Kurtosis[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=((9+6 (κ^4)+9 (κ^8))/(((1+(κ^4)))^2))
AsymmetricLaplaceDistribution/: Quantile[AsymmetricLaplaceDistribution[θ_,τ_,κ_],q_]:=
	Which[0<q≤((κ^2)/(1+(κ^2))),θ+((κ τ Log[q (1+(1/(κ^2)))])/(Sqrt[2])),
		  ((κ^2)/(1+(κ^2)))<q<1,θ-((τ Log[-(-1+q) (1+(κ^2))])/((Sqrt[2]) κ))]
AsymmetricLaplaceDistribution/: RandomReal[AsymmetricLaplaceDistribution[θ_,τ_,κ_]]:=Quantile[AsymmetricLaplaceDistribution[θ,τ,κ],Random[]]
Protect[Mean, Variance, StandardDeviation, Skewness, Kurtosis,
         PDF, CDF, ExpectedValue, DistributionParameterQ,
         DistributionDomainQ,DistributionDomain,
         CharacteristicFunction,Quantile,RandomReal];

Illustrations

Plot[(Quantile[AsymmetricLaplaceDistribution[θ,τ,κ],q]/.{θ->0,τ->4,κ->1/2}),{q,0.01,0.99},AxesLabel->{"α","Quantile"}]

Due to the three parameters the distribution is apparently flexible. Nevertheless its kurtosis is always high (≥6), which might be unsatisfactory in some practical situations.

Print["Minimal kurtosis is ", FindMinimum[Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]],{κ,2}]//First//Round, " and is achieved for κ = ",κ/.(FindMinimum[Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]],{κ,2}]//Last)//Round,"."]

Plot[{Skewness[AsymmetricLaplaceDistribution[θ,τ,κ]],Kurtosis[AsymmetricLaplaceDistribution[θ,τ,κ]]},{κ,0,3},AxesLabel->{"κ"," "}]

Probability density function (illustration)

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

Study specific parameters for ALD and simulation

In specific numerical illustrations we will consider the distributions with the same mean and standard deviation for comparisons. We have chosen AsymmetricLaplaceDistribution[0, 4, 1/2] as a standard (sometimes called reference distribution in the sequel). The corresponding moments are:

Text@Grid[{{"Mean ",Mean[AsymmetricLaplaceDistribution[0.,4,0.5]]},{"SD  ",StandardDeviation[AsymmetricLaplaceDistribution[0.,4,0.5]]},{"Skewness     ",Skewness[AsymmetricLaplaceDistribution[0.,4,0.5]]},{"Kurtosis",Kurtosis[AsymmetricLaplaceDistribution[0.,4,0.5]]}},Frame->All,Alignment->Left]

Simulation

t1=Table[RandomReal[AsymmetricLaplaceDistribution[0.,4,0.5]],{10000}];
Text@Grid[{{"Mean ",Mean[t1]},{"SD  ",StandardDeviation[t1]},{"Skewness     ",Skewness[t1]},{"Kurtosis",Kurtosis[t1]}},Frame->All,Alignment->Left]
Histogram[t1]

VaR

Information about admissible values of the parameters can be transferred through the option assuming.

Clear[VaR];
VaR[distribution_,α_,assuming___]:=Assuming[{assuming},Quantile[distribution,α]//Simplify]

CVar

Clear[CVaR];
CVaR[distribution_,α_,assuming___]:=(1/(1-α))Assuming[{assuming},(∫_(VaR[distribution,α])^∞)y PDF[distribution,y]\[DifferentialD]y]//PowerExpand//Simplify

Clear[NCVaR];
NCVaR[distribution_,α_]:=(1/(1-α))NIntegrate[y PDF[distribution,y],{y,VaR[distribution,α],∞}]

***Average Value at Risk

Clear[AVaR];
AVaR[distribution_,α_,assuming___]:=(1/(1-α))Assuming[{assuming},(∫_α^1)VaR[distribution,p]\[DifferentialD]p]//PowerExpand//Simplify
Clear[NAVaR];
NAVaR[distribution_,α_]:=(1/(1-α))NIntegrate[VaR[distribution,p],{p,α,1}]

***Median Value at Risk (Median Shortfall, Tail Conditional Median)

Clear[MVaR];
MVaR[distribution_,α_,assuming___]:=Assuming[{assuming},Quantile[distribution,(1+α)/2]//PowerExpand//Simplify]

Spectral measures

First we start with the definitions of some spectra:

Dual-power risk aversion function

ϕdual[u_,ν_]:=ν (u^(ν-1)) (* ν≥1 *)

Proportional Hazard risk aversion function

ϕhazard[u_,γ_]:=(1/γ)(((1-u))^(1/γ-1)) (* γ≥1 *)

Wang's risk aversion function

ϕWang[u_,δ_]:=Exp[-δ Quantile[NormalDistribution[0,1],u]-(δ^2)/2] (* δ>0 *)

Exponential risk aversion function

ϕArrowPratt[u_,k_]:=((k Exp[-k(1-u)])/(1-Exp[-k]))

where k is the Arrow-Pratt degree of absolute risk aversion (ARA, see also Wolfram Demonstration by [3] Chandler):

-((D[ϕArrowPratt[u,k],{u,2}])/(D[ϕArrowPratt[u,k],u]))

Spectral measure symbolic

Clear[SpectralRisk];
SpectralRisk[distribution_,ϕ_,assuming___]:=Assuming[{assuming},(∫_(-∞)^∞)y ϕ[CDF[distribution,y]]PDF[distribution,y]\[DifferentialD]y]

Spectral measure numeric

Clear[NSpectralRisk];
NSpectralRisk[distribution_,ϕ_]:=NIntegrate[y ϕ[CDF[distribution,y]]PDF[distribution,y],{y,-∞,∞}]

Expectiles

Since this characteristic is not widely used we go more in details. We present two methods for calculating expectiles. First method is based on direct minimization of (1). It works reliably but it is very slow in some cases.

Clear[Expectile];
Expectile[distribution_,ω_]:=a/.NMinimize[(1-ω)(∫_(-∞)^a)(((x-a))^2)PDF[distribution,x]\[DifferentialD]x+ω(∫_a^∞)(((x-a))^2)PDF[distribution,x]\[DifferentialD]x,a][[2]]

Basic illustrations N(17/4,4):

Text@Grid[Prepend[{{0.75,Expectile[NormalDistribution[17/4,4],0.75]//Timing}//Flatten,{0.75,Expectile[NormalDistribution[17/4,4.],0.75]//Timing}//Flatten,{0.95,Expectile[NormalDistribution[17/4,4.],0.95]//Timing}//Flatten},{"ω","Time","Expectile"}],Frame->All]

The second method is based on formula (2), the first order condition. We take the first derivative of E,(Subscript[Q, ω])(x-a) and equate it to zero. For that, we use FindRoot with the 0.5 and 1.5 multiples of the mean of the distribution as starting points.

(* The derivative *)
D[(1-ω)(∫_(-∞)^a)(((x-a))^2)PDF[distribution,x]\[DifferentialD]x+ω(∫_a^∞)(((x-a))^2)PDF[distribution,x]\[DifferentialD]x,a]//Simplify

Clear[Expectile2];
Expectile2[distribution_,ω_]:=a/.FindRoot[ω (∫_a^∞)2 (a-x) PDF[distribution,x]\[DifferentialD]x-(-1+ω) (∫_(-∞)^a)2 (a-x) PDF[distribution,x]\[DifferentialD]x,{a,0.5 Mean[distribution],1.5Mean[distribution]},MaxIterations->200]

Comparison of timings of the Expectile and Expectile2

We start with Asymmetric Laplace Distribution:

Expectile[AsymmetricLaplaceDistribution[0.,4.,0.5],0.75]//Timing

Procedure based on NMinimize takes enormously long time. Expectile2 procedure is substantially less time consuming for both rational and real input parameters.

Text@"ω = 0.75"
Text@Grid[Prepend[{Expectile2[AsymmetricLaplaceDistribution[0,4,1/2],0.75]//Timing,Expectile2[AsymmetricLaplaceDistribution[0.,4.,0.5],0.75]//Timing},{"Time","Expectile2"}],Frame->All]

Comparison of NMinimize (procedure Expectile) and FindRoot (procedure Expectile2) for N(4.24264,5.83095) and FindRoot (procedure Expectile2) for (ALD(0,4,1/2):

Off[Power::"infy"];Off[∞::"indet"];Off[FindRoot::"nlnum"];Text@Grid[Prepend[timecomparisontable=Table[{ω,Expectile[NormalDistribution[3. (Sqrt[2]),(Sqrt[34.])],ω]//Timing,Expectile2[NormalDistribution[3. (Sqrt[2]),(Sqrt[34.])],ω]//Timing,Expectile2[AsymmetricLaplaceDistribution[0.,4.,0.5],ω]//Timing}//Flatten,{ω,omegas={0.5,0.6,0.7,0.75,0.8,0.85,0.9,0.95,0.99}}],{"ω","Timing","Expectile\n Normal","Timing","Expectile2\n  Normal","Timing","Expectile2\n    ALD"}],Frame->All]

We see that for the normal distribution, the times for both procedures are almost the same independently of ω. In case of asymmetric Laplace distribution and Expectile2 procedure we observe that with increasing ω the timing increases with increasing ω. Above example shows that procedure Expectile (direct minimization) is extremely time consuming for the latter distribution.

***Distribution Specific Calculations

Normal distribution

VaR

VaRNormal[α_,μ_,σ_]=VaR[NormalDistribution[μ,σ],α]

CVaR

CVaRNormal[α_,μ_,σ_]=CVaR[NormalDistribution[μ,σ],α,Re[(σ^2)]>0]

AVaR

AVaRNormal[α_,μ_,σ_]=AVaR[NormalDistribution[μ,σ],α]//PowerExpand//Simplify

Spectral measure

Recall that throughout the paper if we speak on  specific distributions, we consider distributions with mean and standard deviation equal to

{3 (Sqrt[2]),(Sqrt[34])}//N

We start with the spectral risk for a normal distribution.

Manipulate[Plot[NSpectralRisk[NormalDistribution[3. (Sqrt[2]),σ],ϕdual[#,ν]&],{ν,1,10},AxesOrigin->{0,0},PlotRange->{0,20},AxesLabel->{"ν","Spectral risk"}],{{σ,5.83},4,10},SaveDefinitions->True]

Manipulate[Plot[NSpectralRisk[NormalDistribution[3. (Sqrt[2]),σ],ϕArrowPratt[#,k]&],{k,0.1,10},AxesOrigin->{0,0},PlotRange->{0,20},AxesLabel->{k,"σ"}],{{σ,5.83},4,10},SaveDefinitions->True]

The reader will simply discover that such plots are rather boring. More interesting and more important are the ratios of the risk measures.

Expectiles of the normal distribution

Text@Grid[Prepend[tableExpectile2Normal=Table[{ω,Expectile2[NormalDistribution[3. (Sqrt[2]),5.83],ω]},{ω,0.85,0.99,0.01}],{"ω","Expectile2"}],Frame->All]

ListPlot[tableExpectile2Normal,Joined->True,AxesOrigin->{0.85,8},AxesLabel->{"ω","Expectile"}]

Laplace distribution

VaR

VaRLaplace[α_,μ_,β_]=VaR[LaplaceDistribution[μ,β],α,β>0&&(1/2)<α<1]

CVaR

CVaRLaplace[α_,μ_,β_]=CVaR[LaplaceDistribution[μ,β],α,β>0&&(1/2<α<1)]

For (1/2)<α<1 we have 2-2α<1 so that μ-β Log[2-2α]>μ and thus

CVaRLaplace[α_,μ_,β_]=Assuming[β>0,(((∫_(μ-β Log[2-2 α])^∞)(((E^(-((y-μ)/β))) y)/(2 β))\[DifferentialD]y)/(1-α))]//FullSimplify

Observe that direct calculation of the AVaR (which coincides with CVaR) is straightforward and does need any manual adjustment.

AVaR

AVaRLaplace[α_,μ_,β_]=AVaR[LaplaceDistribution[μ,β],α,β>0&&(1/2<α<1)]

Spectral measure

We choose the parameters to have the same mean and standard deviation as in the normal case. We have

{Mean[LaplaceDistribution[μ,β]],StandardDeviation[LaplaceDistribution[μ,β]]}

{μ,(Sqrt[2]) β}

so that a solution to these equations returns the equivalent parameters.

Solve[{μ==3. (Sqrt[2]),(Sqrt[2]) β==5.83},{μ,β}]

For ϕdual the computation mostly takes longer time than for ϕArrowPratt, e.g. However, the absolute exponential risk aversion is, fortunately, most popular in practice. Thus it would be convenient to give some examples with the Arrow-Pratt spectrum only.

Nevertheless, we start with ϕdual:

Text@Grid[Prepend[Table[{ν,NSpectralRisk[LaplaceDistribution[4.24264,4.12243],ϕdual[#,ν]&]},{ν,1,10}],{"ν","Spectral Risk   "}]//Transpose,Frame->True]

And now for the absolute exponential risk aversion.

Manipulate[Plot[NSpectralRisk[LaplaceDistribution[4.24264,β],ϕArrowPratt[#,k]&],{k,0.1,10},AxesOrigin->{0,0},PlotRange->{0,20}],{{β,4.12243},1,10},SaveDefinitions->True]

Expectiles

Expectile2[LaplaceDistribution[4.24264,4.12243],0.75]//Timing

On an  old HP notebook the following calculation appeared to be too ambitious.

tableExpectile2Laplace=Table[{ω,Expectile2[LaplaceDistribution[4.24264,4.12243],ω]},{ω,0.85,0.99,0.01}]

Aborted after 5349 seconds.

Asymmetric Laplace distribution (ALD)

VaR

VaRALD[α_,θ_,τ_,κ_]=VaR[AsymmetricLaplaceDistribution[θ,τ,κ],α]

Manipulate[Plot[VaRALD[α,θ,τ,κ],{α,1/2,1},AxesOrigin->{1/2,0},PlotRange->{0,20},AxesLabel->{"α","VaR"}],{{θ,0},-3,3},{{τ,4},1/2,10},{{κ,1/2},0.1,4},SaveDefinitions->True]

CVaR

See Section Calculations for explanation.

CVaRALD[α_,θ_,τ_,κ_]=Which[0<α≤((κ^2)/(1+(κ^2))),(1/(1-α))(1/2) (1+(κ^2)) τ (-2 (-1+α) θ+(((Sqrt[2]) τ (1+(κ^2) (-1+α-α (Log[α]+Log[1+(1/(κ^2))]))))/κ)),
((κ^2)/(1+(κ^2)))<α<1,-(1/(1-α))(((-1+α) (1+(κ^2)) τ ((Sqrt[2]) θ κ+τ-τ Log[-(-1+α) (1+(κ^2))]))/((Sqrt[2]) κ))];

Manipulate[Plot[CVaRALD[α,θ,τ,κ],{α,1/2,1},AxesOrigin->{1/2,0},AxesLabel->{"α","CVaR"}],{{θ,0},-3,3},{{τ,4},1/2,10},{{κ,1/2},0.1,4},SaveDefinitions->True]

AVaR

See Section Calculations for details.

AVaRALD[α_,θ_,τ_,κ_]=Which[0<α≤((κ^2)/(1+(κ^2))),-(1/(2 (-1+α) κ))(2 θ κ-2 α θ κ+(Sqrt[2]) τ-(Sqrt[2]) (κ^2) τ+(Sqrt[2]) α (κ^2) τ-(Sqrt[2]) α (κ^2) τ (Log[α]+Log[1+(1/(κ^2))])),
((κ^2)/(1+(κ^2)))<α<1,((2 θ κ+(Sqrt[2]) τ-(Sqrt[2]) τ Log[1-α]-(Sqrt[2]) τ Log[1+(κ^2)])/(2 κ))];

Plot[AVaRALD[α,0,4,1/2],{α,0,1},AxesLabel->{"α","AVaR"}]

Spectral measure

NSpectralRisk[AsymmetricLaplaceDistribution[0,4,1/2],ϕdual[#,1]&]//Timing

Table[NSpectralRisk[AsymmetricLaplaceDistribution[0,4,1/2],ϕArrowPratt[#,k]&],{k,0.1,10}]

Manipulate[Plot[NSpectralRisk[AsymmetricLaplaceDistribution[θ,τ,κ],ϕArrowPratt[#,k]&],{k,0.1,10},AxesOrigin->{0,0},AxesLabel->{"k","Spectral risk"}],{{θ,0},-3,3},{{τ,4},1/2,10},{{κ,1/2},0.1,4},SaveDefinitions->True]

Expectiles

First illustration:

Expectile2[AsymmetricLaplaceDistribution[0,4,1/2],0.75]//Timing

(tableExpectile2ALD=Table[{ω,Expectile2[AsymmetricLaplaceDistribution[0.,4.,1./2],ω]}//Flatten,{ω,0.85,0.99,0.01}])

Text@Grid[Prepend[{{0.87`,9.257299332722239`},{0.88`,9.54517970636458`},{0.89`,9.859036271596631`},{0.9`,10.204110341986302`},{0.91`,10.587353105732982`},{0.9199999999999999`,11.018251652751681`},{0.9299999999999999`,11.51021881284805`},{0.94`,12.08307975792722`},{0.95`,12.767884865996166`},{0.96`,13.617254460402796`},{0.97`,14.731159451739403`},{0.98`,16.337714446735852`},{0.99`,19.181023514893155`}},{"ω","Expectile2"}],Frame->True
]

The plot of the expectiles for the reference distribution

ListPlot[tableExpectile2ALD,AxesLabel->{"ω","Expectiles"},AxesOrigin->{0.85,8},Joined->True]

**Value at Risk (VaR, hodnota v riziku)

at confidence level 1-α is the α-quantile of the loss distribution of L, shortly (Subscript[VaR, α])(L):

({{({{P(L≤(Subscript[VaR, α])(L)}}))=α}})

so that

({{({{P(L>(Subscript[VaR, α])(L))=1-α}})}})

Typically α=0.95, 0.99 &c., so that the probability in the above formula is small.

**Parametric VaR (Parametrická hodnota v riziku)

See the text

**Nonparametric VaR (Neparametrická hodnota v riziku)

See the text

**Conditional VaR (CVaR, podmíněná hodnota v riziku)

Expected shortfall (očekávaná extrémní ztráta), Tailed VaR

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

Nonparametric estimate of CVaR

See the text

**Similar risk measures (příbuzné míry rizika)

Average VaR

Weighted Average VaR

Median VaR

Quantile VaR

***Valuation of Securities 1 (Hodnocení investic s pevným výnosem)

https://www.karlin.mff.cuni.cz/~hurt/20201213SecuritiesValuationT1.pdf

Skip:
1.1.2 Exercise, 1.1.5 Remark
Important note:
For Callable bonds (1.1.6 Remark) see also the text in this notebook.

Coupon bonds

See PDF materials

Hodnocení obligací (dluhopisů)

Obligace s maturitou T roků (všude předpokládáme T∈ℕ) a nominální hodnotou F  emitovaná v čase 0 vyplácí kupóny ve výši C v časech 1,...,T. Nominální hodnota je splatná taktéž v čase T. Pro držitele obligace po celou dobu její životnosti je odpovídající peněžní tok:

(* Import["Cashflowofabond.pdf"]//First *)

Obligace je však obchodovatelný cenný papír, nutno znát jeho hodnotu v libovolném čase do splatnosti  t<T. Situace je poněkud komplikována administrativními důvody (v tomto případě pochopitelnými, nutno evidovat řádově tisíce držitelů): Emitent stanoví datum ex-kupónu (Subscript[t, ex]), tj. datum, od kterého je dluhopis obchodován bez nejbližšího následujícího kupónu. (The ex coupon date is the first day the bond starts trading without the coupon attached to it.) Datum ex-kupónu je obvykle několik dnů před výplatou kupónu, typicky 14 dnů. Na obrázku je hodnotící čas t∉ℕ, t<(Subscript[t, ex]). Držiteli obligace tedy náleží kupóny v časech ⌊t⌋+1,⌊t⌋+2,...,T-1,T, celkem tedy T-⌊t⌋ kupónů. Od času (Subscript[t, ex])  jsou vypláceny kupóny v časech ⌊t⌋+2,...,T-1,T. V dalším uvedeme spravedlivé ceny obligace ve všech možných případech.

(* Import["CashflowofabondValuation.eps"] *)

V dalším užíváme elementární vztahy mezi celou ⌊⌋, zlomkovou částí {} a stropem ⌈⌉, zdůrazněme T celé.

t=⌊t⌋+{t} ,  ⌊T-t⌋=T-⌈t⌉,   ⌊t⌋+1-t=1-{t}.

Spravedlivá cena obligace: t∉ℕ, t<(Subscript[t, ex])

P,V=C,(v^(⌊t⌋+1-t))+C,(v^(⌊t⌋+2-t))+...+(Cv^(T-1-t))+(C+F)(v^(T-t))=
C,(v^(1-{t}))((1-(v^(T-⌊t⌋)))/(1-v))+F,,(v^(T-t))

Spravedlivá cena cena obligace: t∉ℕ, (Subscript[t, ex])≤t≤⌊t⌋+1

P,V=C,(v^(⌊t⌋+2-t))+C,(v^(⌊t⌋+3-t))+...+(Cv^(T-1-t))+(C+F)(v^(T-t))=
C,(v^(2-{t}))((1-(v^(T-⌊t⌋-1)))/(1-v))+F,,(v^(T-t))

Spravedlivá cena cena obligace: t∈ℕ, v čase t se pochopitelně žádný kupón nevyplácí

P,V=C,v+C,(v^2)+...+(Cv^(T-1-t))+(C+F)(v^(T-t))=
C,v((1-(v^(T-t)))/(1-v))+F,,(v^(T-t))

Čistá cena (t=0, kupóny se vyplácejí v okamžicích t=1,2,...,T):

P,(Subscript[V, p])=C,v+C,(v^2)...+(C+F)(v^T)=C,(∑_(j=1)^T)(v^j)+F,(v^T)=((C,,v((v^T)-1))/(v-1))+F,(v^T)=(C/i)+(1/(((1+i))^T))(F-(C/i))

Example (Explicit formulas)

Označení: “FaceValue” nominální hodnota F, ”Coupon” kuponova sazba c, “Maturity” doba do splatnosti od okamziku vydani T, “Settlement” hodnotici okamzik t, “InterestRate” hodnotici urokova mira r, trzni urokova mira pro srovnatelnou investici, cenu kapitalu, “CouponInterval” rocni frekvence vyplaty kuponu 1/f, (f=2 znamena dvakrat rocne), “DayCountBasis” kalendářová konvence, “AccruedInterest” alikvotni urok, “Value” clean value, “FullValue” dirty value.

Poznámka: 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}})

Clear[F,c,T,f,r,t]
FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->T,"CouponInterval"->(1/f)},{"InterestRate"->r,"Settlement"->t}]//Apart//TraditionalForm

Value (Clean price)

(value=FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->T,"CouponInterval"->(1/f)},{"InterestRate"->r,"Settlement"->t},"Value"])//Apart//FullSimplify//TraditionalForm

FullValue (Dirty price)

(fullvalue=FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->T,"CouponInterval"->(1/f)},{"InterestRate"->r,"Settlement"->t},"FullValue"])//Apart//TraditionalForm

Alikvótní úrok

(accrued=FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->T,"CouponInterval"->(1/f)},{"InterestRate"->r,"Settlement"->t},"AccruedInterest"])//Apart//TraditionalForm

což je vlastně c,F{f,t}

fullvalue-value//FullSimplify//TraditionalForm

fullvalue-value==accrued//FullSimplify

FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->5},{"InterestRate"->r,"Settlement"->0}]
%//Apart//TraditionalForm

Manipulate[Plot[({value,fullvalue}/.{F->1,c->ic,f->1,T->iT,r->ir})//Evaluate ,{t,0,iT},Exclusions->"Singularities",PlotStyle->Thick,AxesLabel->{"Time to Maturity","Value"},PlotLabel->"Value of coupon bond\nPar Value F = 1, yearly coupons",PlotLegends->{"Clean price","Fair price"}],
Style["Maturity",Bold],
{{iT,10},1,20,1,Appearance->"Labeled",ImageSize->Tiny},
Style["Market IR",Bold],
{{ir,0.05},0.001,0.3,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Coupon rate",Bold],
{{ic,0.1},0.01,0.2,0.01,Appearance->"Labeled",ImageSize->Tiny},
(* Delimiter,
Style["Initial value (Subscript[S, 0])",Bold],
{{(Subscript[S, 0]),100},10,150,10,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&}, *)
SaveDefinitions->True,ControlPlacement->Left]

Callable bonds

See PDF materials

Svolatelný dluhopis (vypověditelný, callable bond)

Emitent může dluhopis splatit (odkoupit, svolat) před uplynutím normální doby splatnosti. Zpravidla tak učiní v případě, že má přístup k “levnějším” penězům než musí platit v úrocích emitovaného dluhopisu. Pro držitele dluhopisu (investora) to bývá nevýhodné (emitent obvykle svolává v případě, že tržní ú. m. je nižší než kupónová sazba, investor tedy obvykle na trhu nenajde stejně výhodnou investiční příležitost) , takže emitent musí nabídnout určitou kompenzaci. Ta je obvykle zakotvena ve svolávacím plánu (call schedule).

Je-li T doba do splatnosti, nominální hodnota F, pak emitent může obvykle svolat v okamžicích (Subscript[T, -k])=T-k, k=1,...,K    se splatnými částkami (Subscript[F, -1]),...,(Subscript[F, -K]). Z podstaty věci vyplývá, že by mělo platit (Subscript[F, -k])>F. K exaktnímu hodnocení musíme též znát pravděpodobnosti svolání (Subscript[p, -1]),...(Subscript[p, K]) v okamžicích (Subscript[T, -1]),...,(Subscript[T, K]). Je přirozené, že tyto pravděpodobnosti se mění v závislosti na tržní situaci a jsou závislé na hodnotícím okamžiku t: (Subscript[p, -1])(t),...(Subscript[p, K])(t). Struktura CALL opce je tedy poměrně komplikovaná.

V praxi se uvádějí svolávací ceny v procentech: (Subscript[F, 0])=100%,(Subscript[F, -1])=102%,(Subscript[F, -2])=104%,....

Mluví se o též dluhopisu s vnořenou opcí (embedded option). Opci může uplatnit emitent (je v dlouhé pozici), takže

(cena svolatelného dluhopisu) = (cena ekvivalentního dluhopisu bez možnosti svolání) - (cena CALL opce).

***Price of a callable bond is always lower than the price of a straight bond because the call option adds value to an issuer. Yield on a callable bond is higher than the yield on a straight bond.

Viz též [3], str. 11, 90 a dále.

Poznámka (puttable bond)

V tomto případě má držitel právo prodat dluhopis emitentovi, opět za předem dohodnutých podmínek (prodejní opce). Exisují také dluhopisy s oběma opcemi, CALL opce a PUT opce se vzájemně nevylučují.

Amortized bonds

See PDF materials

Amortizované obligace

Emitent nesplatí v okamžiku vypršení nominální hodnotu, ale “amortizuje” ji v konstatních kupónech. Místo standardního peněžního toku pro běžnou kupónovou obligaci {C,C,...,C,C+F} tak investor dostává anuitu {(SuperStar[C]),(SuperStar[C]),...,(SuperStar[C])}. Amortizovaná obligace je, ceteris paribus, méně riskantní (proč?), takže “ekvivalentní” kupón bude nižší než vypočtený z přesné ekvivalence peněžních toků (samozřejmě závislých na hodnotící ú. m.).

Bonds under uncertainty

Simple Bonds under Uncertainty I

Suppose that a zero coupon bond pays F with probability p and pays η,F with probability 1-p at maturity one year , where η ∈[0, 1) is called the recovery rate (míra pokrytí, míra zotavení), and  λ:=1-η  is called the loss rate (ztrátový poměr). The case η=0  is equivalent to the default (nesplacení). Suppose that the valuation interest rate is r. The fair value of the bond is the expected present value:

(OverBar[PV]) =(1/(1+r))[p,F+(1-p)η,F]

Simple Bonds under Uncertainty 2

See PDF material.

Uvažujme pro změnu obligaci s nulovým kuponem splatnou za jeden rok s nominální hodnotou F=1, ale ve skutečnosti vyplatí X, kde X je náhodná veličina, 0≤X≤1. Nechť r je hodnotici urokova mira. Spravedlivá hodnota obligace je pak

(OverBar[PV]) =(1/(1+r))E,X

Risk premiums

See PDF material.

Podobnou myšlenku můžeme užít i pro hodnoceni půjček, resp. při stanovení prémie z rizika nesplacení (default risk premium). Předpokládejme, že půjčka ve výši PV = 1 je splatná za 1 rok, ale ve skutečnosti dlužník splatí náhodnou částku X, 0≤X≤1, úrok v to nepočítaje. Otázkou je, jaká by měla být prémie za riziko z nesplacení (Subscript[r, default]) účtovaná věřitelem, aby dostal v průměru bezrizikový úrok (Subscript[r, 0]). Jistá budoucí hodnota půjčky je 1+(Subscript[r, 0]), takže věřitel by měl požadovat rizikovou přirážku vyplývajici ze vztahu

1+(Subscript[r, 0])=(1+(Subscript[r, 0]))(1+(Subscript[r, default]))E,X

z čehož

(Subscript[r, default])=(1/(E,X))-1.

2023-01-04 Lesson 13

1. Numerical illustrations will be completed in this notebook and you will be informed by e-mail in a due time. Pay attention, it will contain valuable information for exam purposes.
2. State exams requirements: Ask for a meeting with a responsible supervisor.

Risk measures: Expectiles: The last Section of
http://www2.karlin.mff.cuni.cz/~hurt/20221214_ReturnRisk.pdf

Valuation of Securities 1 (Fixed Income)

http://www2.karlin.mff.cuni.cz/~hurt/20230102SecuritiesValuation.pdf

Valuation of Securities 2 (Derivatives)

http://www2.karlin.mff.cuni.cz/~hurt/20230102FinancialDerivatives.pdf

FAQ

Due the time pressure the questions will be exceptionally answered in written by request send by e-mail.

Q1

A1

***Valuation of Securities 2 (Finanční deriváty)

http://www.karlin.mff.cuni.cz/~hurt/20201220FinancialDerivativesT1.pdf

Numbering of formulas and the references needn’t be correct in this material.

Skip:
p. 3: Natural Boundaries for option prices
p. 4: 1.1.3 Exercise
pp. 6-10: The Binomial Option Pricing Model
p. 18: skip all

Without derivations, just the meaning, interpretation, and natural properties:
pp. 11-14: Comparative Statics - The Greeks

***Options (Opce)

Opakování

Evropské a americké opce (vanilla options)

r ... bezriziková úr. míra (nekdy značena  (Subscript[r, 0]))
(Subscript[S, t]) ... cena podkladového aktiva (akcie) v čase t
(Subscript[c, t]) resp. (Subscript[p, t]) ... cena evropské CALL resp. PUT opce
K ... realizační cena (strike)
T ... čas vypršení evropské opce (expiry)

Výplatní funkce (payoff function) evropské opce

v okamžiku vypršení

CALL: (SuperPlus[((Subscript[S, T])-K)])

PUT: (SuperPlus[(K-(Subscript[S, T]))])

PUT - CALL parita (platí pouze pro evropské opce!)

({{({{(Subscript[S, t])+(Subscript[p, t])=(Subscript[c, t])+K,(e^(-r(T-t)))}})}})

Mnemotechnická pomůcka pro zapamatování: státní podnik = císařsko-královské

Důkaz: viz též [1], str. 54. Uvažujme portfolio +PUT +akcie -CALL, vše na stejnou akcii, opce evropské se stejnou dobou vypršení. Hodnota portfolia v čase t≤T je

({{(Subscript[Π, t])=(Subscript[S, t])+(Subscript[p, t])-(Subscript[c, t]).}})

V čase vypršení T je hodnota tohoto portfolia

({{({{(Subscript[Π, T])=(Subscript[S, T])+(SuperPlus[(K-(Subscript[S, T]))])-(SuperPlus[((Subscript[S, T])-K)]).}})}})

V případě (Subscript[S, T])≤K je (Subscript[Π, T])=(Subscript[S, T])+K-(Subscript[S, T])-0=K, v případě (Subscript[S, T])>K je (Subscript[Π, T])=(Subscript[S, T])+0-((Subscript[S, T])-K)=K. To znamená, že portfolio je bezrizikové vedoucí k jistému zisku K. Budoucí hodnota je  (Subscript[Π, T])=K  a současná hodnota je tedy (Subscript[Π, t])=K,(e^(-r(T-t))).

Současná hodnota výplatní funkce (payoff function) evropské opce v čase t≤T

Diskontujeme k času t:

CALL:   (SuperPlus[(e^(-r(T-t)))((Subscript[S, T])-K)])

PUT:   (SuperPlus[(e^(-r(T-t)))(K-(Subscript[S, T]))])

Black-Scholesův vzorec

Risk Neutral Valuation

We assume that the stochastic process of the prices of an asset  (Subscript[{(Subscript[S, t])}, t≥0]) is governed by a geometrical Brownian motion, hence

({{({{ℒ((Subscript[S, T])|(Subscript[S, t]))=L,N(ln (Subscript[S, t])+(μ-(1/2)(σ^2))(T-t),(σ^2)(T-t))}})}})

where μ is the rate of return on the underlying asset and σ its volatility.

Rizikově neutrální hodnocení: Investoři se chovají racionálně, čím větší riziko (volatilita), tím vyšší očekávaný výnos požadují. Čím vyšší averze k riziku, tím vyšší očekávaný výnos je požadován. Vyloučením rizikových preferencí dojdou investoři ke shodě na jediném rizikově neutrálním výnosu μ=r. Pozor, v žádném případě to neznamená, že by výnos podkladového aktiva byl roven r !!!

Alternativní princip vedoucí k témuž závěru: Očekávaný výnos z držení akcie v hodnotě (Subscript[S, t])  v období (t,T) by měl být stejný jako bezrizikový výnos z investice  ve výši (Subscript[S, t])  při úrokové míře r:

({{({{(SuperStar[E])((Subscript[S, T])|(Subscript[S, t]))=FV((Subscript[S, t]),r)=(Subscript[S, t])(e^(r(T-t)))}}).}})

Hvězdičkou zde i v dalším značíme skutečnost, že za výnos bereme rizikově-neutrální výnos r.

Neuvažujeme-li žádné rizikové preference investorů, je

({{({{(SuperStar[ℒ])((Subscript[S, T])|(Subscript[S, t]))=L,N(ln (Subscript[S, t])+(r-(1/2)(σ^2))(T-t),(σ^2)(T-t))}}).}})

Then

({{({{(SuperStar[E])((Subscript[S, T])|(Subscript[S, t]))=}})}})

Mean[LogNormalDistribution[Log[(Subscript[S, t])]+(r-(1/2)(σ^2))(T-t),σ(Sqrt[T-t])]]//Simplify//TraditionalForm

Similarly for the variance:

({{({{(SuperStar[var])((Subscript[S, T])|(Subscript[S, t]))=}})}})

Variance[LogNormalDistribution[Log[(Subscript[S, t])]+(r-(1/2)(σ^2))(T-t),σ(Sqrt[T-t])]]//Simplify//TraditionalForm

See [1], pp. 55-56.

Hodnota (cena) evropské opce

Jako přirozené se jeví vzít za hodnotu evropské opce v čase t podmíněnou střední hodnotu současné hodnoty výplatní funkce při dané hodnotě podkladového aktiva (Subscript[S, t]) při čemž střední hodnota je uvažovaná při bezrizikové (tj. rizikově neutrální) úrokové míře r:

({{({{(Subscript[c, t])=(e^(-r(T-t)))(SuperStar[E])((SuperPlus[((Subscript[S, T])-K)])|(Subscript[S, t]))}, {(Subscript[p, t])=(e^(-r(T-t)))(SuperPlus[(SuperStar[E])(K-(Subscript[S, T]))])|(Subscript[S, t]))}})}})

(* 65 sec!!! *)ct=(E^(-r(T-t)))Expectation[Max[0,((Subscript[S, T])-K)],(Subscript[S, T])\[Distributed]LogNormalDistribution[Log[(Subscript[S, t])]+(r-(1/2)(σ^2))(T-t),σ(Sqrt[T-t])]]

ct[[2,2]]//FullSimplify

Tuto formuli lze převést do běžně uváděného tvaru (Black-Scholes formula)

({{({{(Subscript[c, t])=(Subscript[S, t])Φ((Subscript[d, 1]))-K (e^(-r(T-t)))Φ((Subscript[d, 2])),
(Subscript[d, 1,2])=((log((Subscript[S, t])/K)+(r±(σ^2)/2)(T-t))/(σ(Sqrt[T-t]))).}})}})

Hodnota PUT opce se získá z PUT - CALL parity: ({{({{(Subscript[p, t])=(Subscript[c, t])+K,(e^(-r(T-t)))}})-(Subscript[S, t]).}})

**Poznámka (opce na akcie vyplácející dividendu)

Vyplácí-li akcie dividendu s dividendovou intenzitou d, pak se v předchozím vzorci nahradí r→r-d. Vyplácením dividendy se snižuje hodnota firmy a tudíž i kapitálový výnos.

**Citlivosti (Comparative Statics, The Greeks)

Statics je sice vžitý, ale asi chybně použitý termín nevystihující podstatu.

Nechť (Subscript[V, t]) označuje hodnotu CALL resp. PUT opce, (Subscript[V, t])=(Subscript[c, t]) resp. (Subscript[V, t])=(Subscript[p, t]). Z Black-Scholesova vzorce je patrné, že (Subscript[V, t]) je funkcí pěti proměnných (parametrů):

({{({{(Subscript[V, t])=(Subscript[V, t])((Subscript[S, t]),K,σ,T-t,r)}}).}})

Derivace podle jednotlivých proměnných se nazývají citlivosti a protože se označují řeckými písmeny, jsou též známé jako “The Greeks”. Jejich výpočet je přímočarý, ale v některých případech technicky dosti komplikovaný. Uvedeme jen jednu míru citlivosti delta, derivaci podle aktuální ceny akcie:

({{({{Δ=((∂(Subscript[V, t]))/(∂(Subscript[S, t])))}}).}})

Pro CALL resp. PUT dostaneme

({{({{(Subscript[Δ, c])=Φ((Subscript[d, 1]))   resp.   (Subscript[Δ, p])=(Subscript[Δ, c])-1=-Φ(}})-(Subscript[d, 1])).}})

Charakteristiky Δ se používá k tzv. delta-zajišťování, delta hedging. Uvažujme portfolio skládající se z dlouhé pozice jedné akcie a krátké pozice v A kusech CALL opcí na tutéž akcii. Hodnota portfolia je tedy (Subscript[Π, t])=(Subscript[S, t])-A,(Subscript[c, t]). Chceme najít A takové, aby hodnota portfolia byla invariantní vůči (malým) změnám ceny akcie, tj.

({{({{((∂,(Subscript[Π, t]))/(∂(Subscript[S, t])))}})=0.}})

Hledané A=1/(Subscript[Δ, c]), nazývané zajišťovací poměr, hedge ratio. Vzhledem k tomu, že se delta mění v čase, musí se portfolio často měnit → dynamické zajišťování, dynamic hedge. To může být nákladné kvůli transakčním nákladům. Tudíž vhodné pro makléře a tvůrce trhu, kteří mají relativně nízké transakční náklady. Všimněme si, že čím nižší (Subscript[S, t]) ve srovnání s K, tím vyšší je zajišťovací poměr. Podobnou mírou citlivosti je elasticita ceny CALL opce vzhledem k ceně akcie definovaná jako (Subscript[E, c])=(Subscript[Δ, c])(Subscript[S, t])/(Subscript[c, t]). Vždy je (Subscript[E, c])>1. (Dokažte za cvičení.) Proto je CALL opce vždy riskantnější než její podkladové aktivum. Analogie k Δ je durace.

**Další míry citlivosti viz [1], str. 60-61.

*Implied Volatility (Implikovaná volatilita)

Ve výše uvedeném vzorci pro (Subscript[V, t]) v okamžiku t vždy známe hodnoty (Subscript[S, t]),K,T-t,r. Při aplikaci Black-Scholesova vzorce pro neznámou volatilitu obvykle použijeme výběrovou směrodatnou odchylku z historických (známých) dat logaritmických výnosů ln((Subscript[S, t+1])/(Subscript[S, t])) pro nějaké rozumné historické hodnoty. Sofistikovanější metody jsou založené na modelech volatility typu ARCH, GARCH, jejich zobecněních a modifikacích. Viz zabudované funkce (nezahrnují zdaleka všechny dnes používané modely):

?*ARCH*

?*AR*

Můžeme ale též přistoupit k odhadu volatility jiným způsobem. V daném časovém okamžiku t, (teď, hodnotící časový okamžik), známe nejenom  hodnoty (Subscript[S, t]),K,T-t,r, ale i tržní hodnotu opce (Subscript[V, t]^M). Po substituci známých hodnot do Black-Scholesovy formule dostaneme rovnici pro neznámou σ:

({{(Subscript[V, t]^M)=({{(Subscript[V, t])((Subscript[S, t]),K,σ,T-t,r)}}).}})

Tato rovnice nemusí mít žádné řešení (řešení je nutné hledat numericky), může mít právě jedno (nejpříjemnější situace, i když i tam můžeme býti překvapeni) nebo více řešení - z těch je pak zapotřebí citlivě vybírat.

Důležité! Volatilita a implikovaná volatilita viz [1], str. 61-62.

Illustrations

Technické záležitosti týkající se funkce FinancialDerivative:

Označení: “StrikePrice” realizacni cena K, “Expiration” doba do splatnosti, doba do vyprseni T, pokud se uvede “ReferenceTime” t, pak je to T-t,  “InterestRate” bezrizikova urokova mira (Subscript[r, 0]), “Volatility” volatilita ceny podkladoveho aktiva σ (underlying asset), “CurrentPrice” soucasna cena podkladoveho aktiva (Subscript[S, 0]) resp. (Subscript[S, t]) . “Dividend” je dividendova intenzita q. Vsechny udaje se vztahuji ke stejnemu casovemu obdobi, typicky 1 rok.

(* vsechny zabudovane derivaty, zde zvane "Contracts" *)
FinancialDerivative[](* //Length *)

POZOR NA SYNTAXI!!! PRI RESENI ULOH ZKOPIRUJTE Z HELPU!!! Typicky ma FinancialDerivative 3 az 4 argumenty, kazdy z nich je seznam (list).

1 Priklad (hodnoty opci)

FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50, "Dividend"->0.05}]

FinancialDerivative[{"American","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50, "Dividend"->0.05}]

**2 Priklad (hodnoty opci a jejich citlivosti, The Greeks)

V ... price of the option, Δ=(∂/(∂(Subscript[S, t])))V,  Γ=((∂^2)/(∂(Subscript[S, t]^2)))V,  Ρ=(∂/(∂r))V,  Θ=(∂/(∂t))V,  𝒱=(∂/(∂σ))V

Opce z predchoziho prikladu, vystupem hodnota opce a citlivosti

FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50, "Dividend"->0.05}, {"Value", "Greeks"}]

3 Priklad (Grafické znázornění závislosti ceny opce na ceně podkladového aktiva)

Manipulate[Plot[FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->T},  {"InterestRate"-> 0.1, "Volatility" ->σ, "CurrentPrice"-> s, "Dividend"->0.0}],{s,40,60},AxesOrigin->{Automatic,0},PlotRange->{0,20},PlotStyle->{Thick,Red},AxesLabel->{"(Subscript[S, t])",Automatic}],{{T,0.05},0.0001,1,Appearance->"Labeled"},{{σ,0.5},0.01,1,Appearance->"Labeled"},SaveDefinitions->True]

Manipulate[Plot[{FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->T},  {"InterestRate"-> 0.1, "Volatility" ->σ, "CurrentPrice"-> s, "Dividend"->0.0}],FinancialDerivative[{"European","Put"}, {"StrikePrice"-> 50.00, "Expiration"->T},  {"InterestRate"-> 0.1, "Volatility" ->σ, "CurrentPrice"-> s, "Dividend"->0.0}]},{s,40,60},AxesOrigin->{Automatic,0},PlotRange->{0,20},PlotStyle->Thick,AxesLabel->{Style["(Subscript[S, t])",18],Automatic},PlotLegends->{"CALL","PUT"}],{{T,0.05},0.0001,1,Appearance->"Labeled"},{{σ,0.5},0.01,1,Appearance->"Labeled"},SaveDefinitions->True]

Manipulate[Plot[FinancialDerivative[{"American","Call"}, {"StrikePrice"-> 50.00, "Expiration"->T},  {"InterestRate"-> 0.1, "Volatility" ->σ, "CurrentPrice"-> s, "Dividend"->0.0}],{s,40,60},AxesOrigin->{Automatic,0},PlotRange->{0,25},PlotStyle->Thick,AxesLabel->{"(Subscript[S, t])",Automatic}],{{T,0.05},0.0001,1,Appearance->"Labeled"},{{σ,0.5},0.01,1,Appearance->"Labeled"},SaveDefinitions->True]

4 Příklad (Grafické znázornění závislosti na ceně podkladového aktiva, různé doby do splatnosti)

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}}]}//Flatten,{st,30,60},PlotLabel->"Values of CALL options for various expiry times",PlotStyle->Thick,PlotStyle->Thick,AxesLabel->{Style["(Subscript[S, t])",18],Automatic}]

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},PlotStyle->Thick,PlotStyle->Thick,AxesLabel->{Style["(Subscript[S, t])",16],Automatic}],
{{t,0,"Doba do vyprseni"},0,1,Appearance->"Labeled"},SaveDefinitions->True]

***BEGIN*****************************************************
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"}]

***END***********************************************

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]

Cena tohoto instrumentu podle Black-Scholesova vzorce je tedy (musíme přidat bezrizikovou úrokovou míru r,0, volatilitu σ 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,  v následující ilustraci K=1, B=2, (Subscript[r, 0])=0.05 jsou fixované, nezávisle proměnná je cena podkladového aktiva (Subscript[S, t]):

Manipulate[Plot[cenaCAPnaCALL[st,1,2,0.05,σ,Tt],{st,0,5},PlotRange->{0,2},PlotStyle->Thick,AxesLabel->{Style["(Subscript[S, t])",16],Automatic}],{{σ,0.3},0.1,0.5,Appearance->"Labeled"},{{Tt,1},0.01,2,Appearance->"Labeled"},SaveDefinitions->True]

FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50}, {"Value", "Greeks"}]

FinancialDerivative[{"American","Call"}, {"StrikePrice"-> 50.00, "Expiration"->1},  {"InterestRate"-> 0.1, "Volatility" -> 0.5, "CurrentPrice"-> 50}, {"Value", "Greeks"}]

***Stochastic Models of Interest Rates and Price Developments

Introductory Notes

Calibration

Calibration (calibrating instruments) in financial literature is in fact estimation in statistical terminology.

A General Procedure for Characteristics of Stochastic Processes

This procedure should be initialized before further calculations:

In Mathematica, the Normal and related distributions have the second parameter standard deviation σ and not the variance (σ^2)!!!

Clear[characteristicsofstochasticprocess];characteristicsofstochasticprocess[process_,s_,t_]:=Module[{a},Text@Grid[{{"Process:   ",   Style["  "<>( process//HoldForm//ToString)<>"  ",Bold]},{"Mean   ",Mean[process[t]]//FullSimplify//TraditionalForm},
{"Variance   ",Variance[process[t]]//FullSimplify//TraditionalForm},
{"Skewness   ",Skewness[process[t]]//FullSimplify//TraditionalForm},
{"Kurtosis   ",Kurtosis[process[t]]//FullSimplify//TraditionalForm},
{"Covariance function   ",CovarianceFunction[process,s,t]//FullSimplify//TraditionalForm},{"Slice distribution   ",SliceDistribution[process,t]//FullSimplify//TraditionalForm}},Frame->All]]

Random walk

Prices (Subscript[P, 0]),(Subscript[P, 1]),... 
Random walk hypothesis for logarithmic prices (Subscript[p, t])=log (Subscript[P, t]):

(Subscript[p, t+1])-(Subscript[p, t])=μ+(Subscript[ϵ, t+1]), t=0,1,...

ϵ's independent  identically distributed (iid) , E,(Subscript[ϵ, t])=0

The continuous and discretized versions folow.

!Wiener process

?*Process

Clear[μ,σ,s,t];
characteristicsofstochasticprocess[WienerProcess[μ,σ],s,t]

RandomFunction is a discrete version of a continuous stochastic process. In case of Wiener process it is simply the random walk. Here is an example: time interval [0,1], step 0.01, results in 101 data points.

Simulation of one trajectory. For an ω given by the random number generator we plot a discretized version of W(t,ω) as a function of t.

SeedRandom[13131];ListLinePlot[RandomFunction[WienerProcess[μ,σ]/.{μ->1,σ->0.1},{0,1,0.01}],ImageSize->{350,300},PlotLabel->"Trajectory of Wiener process, drift = 1, volatility = 0.1"]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[WienerProcess[μ,σ],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,5},PlotTheme->"Detailed",PlotLabel->"Wiener process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0},-1,1,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.001,2,0.1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

Fifty trajectories, i. e.,  (Subscript[ω, 1]),...,(Subscript[ω, 50]):

SeedRandom[13131];ListLinePlot[RandomFunction[WienerProcess[μ,σ]/.{μ->1,σ->0.5},{0,1,0.01},50],ImageSize->{350,300},PlotStyle->Thick]

!Geometric Brownian Motion Process

First see the text for motivation, a discrete version,

The state X(t) of a geometric Brownian motion satisfies an Ito differential equation

dX(t)==μ X(t)dt +σ X(t) dW(t).

Here  σ>0, μ∈ℝ.

Clear[μ,σ,(Subscript[S, 0]),s,t];
characteristicsofstochasticprocess[GeometricBrownianMotionProcess[μ,σ,(Subscript[S, 0])],s,t]

SeedRandom[13131];ListLinePlot[RandomFunction[GeometricBrownianMotionProcess[μ,σ,(Subscript[S, 0])]/.{μ->1,σ->0.05,(Subscript[S, 0])->100},{0,1,0.01}],ImageSize->{350,300}]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[GeometricBrownianMotionProcess[μ,σ,(Subscript[S, 0])],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,300},PlotTheme->"Detailed",PlotLabel->"Geometric Brownian Motion\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0},-1,1,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.001,2,0.1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
Style["Initial value (Subscript[S, 0])",Bold],
{{(Subscript[S, 0]),100},10,150,10,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

Fractional Brownian Motion

h ... Hurst coefficient

Clear[μ,σ,s,t,h];
characteristicsofstochasticprocess[FractionalBrownianMotionProcess[μ,σ,h],s,t]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[FractionalBrownianMotionProcess[μ,σ,h],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,5},PlotTheme->"Detailed",PlotLabel->"Fractional Brownian Motion Process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0},-1,1,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.001,2,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["Hurst coefficient",Bold],
{{h,0.5},0.05,.95,0.1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

Distribution of Z(t)-Z(s), {Z(t),t≥0} the standard fBm

difffBm=TransformedProcess[z[t]-z[s],z\[Distributed]FractionalBrownianMotionProcess[h],t];

Variance[difffBm[t]]

Jump Diffusion Processes, Illustration of CLT, 3σ Rule

Poisson Process

Let {N(t),t≥0} be a homogeneous Poisson process with intensity λ, i. e., a stochastic process with independent increments, N(0)=0 w. p. 1, and N(t)-N(s)\[Distributed]Poisson(λ|t-s|). The trajectories are nondecreasing with jumps of magnitude 1.

Clear[λ,s,t];
characteristicsofstochasticprocess[PoissonProcess[λ],s,t]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[PoissonProcess[λ],{0,T,0.01}],ImageSize->{350,300},PlotRange->{0,2 λ T},PlotTheme->"Detailed",PlotStyle->Thick,PlotLabel->"Poisson Process\n",InterpolationOrder->0]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Intensity λ",Bold],
{{λ,20},1,100,10,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

Compound Poisson Process

Let {N(t),t≥0} be a homogeneous Poisson process with intensity λ, independent of iid (Subscript[Z, 1]),(Subscript[Z, 2]),... . Then

({{{({{X(t):=(∑_(i=1)^(N(t)))(Subscript[Z, i])}}), t≥0}}})

is a compound Poisson process. Note that (Subscript[Z, i])'s may be of arbitrary sign.

Built-in function:CompoundPoissonProcess[λ, distribution of jumps]

Example

Clear[λ,μ,ν,s,t];
characteristicsofstochasticprocess[CompoundPoissonProcess[λ,ExponentialDistribution[ν]],s,t]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[CompoundPoissonProcess[λ,ExponentialDistribution[ν]],{0,T,0.01}],ImageSize->{350,300},PlotRange->{0,100 λ T},PlotTheme->"Detailed",PlotStyle->Thick,PlotLabel->"Compound Poisson Process\n",InterpolationOrder->0]],
Style["time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["jumps",Bold],
{{ν,0.01},0.01,0.2,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["rate parameter λ",Bold],
{{λ,20},10,100,10,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

***Insurance Example (see Mathematica Help)

Jumps follow iid Pareto type II distribution with the density function

PDF[ParetoDistribution[k,α,μ],x]//TraditionalForm

For α>2

Text@Row@{"Mean:   ", Mean[ParetoDistribution[k,α,μ]][[1,1,1]]//TraditionalForm,Spacer[20],"Variance:    ",Variance[ParetoDistribution[k,α,μ]][[1,1,1]]//TraditionalForm}

Aggregate claims from a risk follow a compound Poisson process with Poisson parameter 200. The claim amount distribution is a Pareto distribution with minimum value parameter 300,shape parameter 3,and location parameter 0. The insurer has effected excess of loss reinsurance with retention level 300. Simulate the claims process for four years:

The mean of the corresponding Pareto distribution is

Mean[TransformedDistribution[Max[0,x-300], x\[Distributed]ParetoDistribution[300,3,0.]]]

and the standard deviation

StandardDeviation[TransformedDistribution[Max[0,x-300], x\[Distributed]ParetoDistribution[300,3,0.]]]

claimsProcess=CompoundPoissonProcess[200, TransformedDistribution[Max[0,x-300], x\[Distributed]ParetoDistribution[300,3,0]]];

characteristicsofstochasticprocess[claimsProcess,s,t]

SeedRandom[13131];
Manipulate[ListLinePlot[RandomFunction[claimsProcess,{0,4},n],InterpolationOrder->0,PlotStyle->Thick,PlotTheme->"Detailed"],
Style["Number of trajectories",Bold],
{{n,8},2,100,2,Appearance->"Labeled",ImageSize->Tiny},
SaveDefinitions->True,ControlPlacement->Left]

Slice distributions for the first four years:

Histogram[Table[RandomVariate[TruncatedDistribution[{0,50000},claimsProcess[t]],5 10^3],{t,4}],150,"PDF",ChartLegends->(StringJoin["t = ",ToString[#]]&/@Range[4])]

Mean and standard deviation of the reinsurer’s aggregate claims for the first four years:

Table[Mean[claimsProcess[t]], {t,4}]

Table[StandardDeviation[claimsProcess[t]], {t,4}]//N

**Illustration of the Central Limit Theorem

Define a random number generator (n=12 in practise)

Clear[generator];
generator[n_]:=Module[{i},((∑_(i=1)^n)RandomVariate[UniformDistribution[{0,1}]]-n/2)/(Sqrt[n/12])]

Functional programming approach, results are the same but this generator is much more faster:

Clear[generator2];
generator2[n_?(IntegerQ[#]&&Positive[#]&)]:=(Total[RandomReal[{0,1},n]]-n/2)/(Sqrt[n/12])

Illustration

SeedRandom[13131];
{generator[12],generator[12],generator[12]}
SeedRandom[13131];
{generator2[12],generator2[12],generator2[12]}

SeedRandom[13131];
tab10000=Table[generator2[12],{10000}];//Timing
Show[{Plot[PDF[NormalDistribution[],x],{x,-4,4},PlotStyle->Thick],Histogram[tab10000,100,"PDF"]}]
Text@Row@{"p-value:   ",DistributionFitTest[tab10000]}

{Mean[tab10000],StandardDeviation[tab10000],Skewness[tab10000],Kurtosis[tab10000]}

More generally for more samples (numberOfSamples) of the same size (sizeOfSample) . Implicitely, random seed is set to 13131.

tableOfpValues[sizeOfSample_,numberOfSamples_,seedrandom_:13131]:=(SeedRandom[seedrandom];Map[DistributionFitTest,Table[Table[generator2[12],{sizeOfSample}],{numberOfSamples}]])

(* 10 samples of size 10000 *)
tableOfpValues[10000,10]

(* 100 samples of size 10000 *)
pValues100=tableOfpValues[10000,100]

Histogram[pValues100,10]
Text@Row@{"Minimum of p-values:  ",Min[pValues100],Spacer[20],"Maximum of p-values:   ",Max[pValues100]}
(* Number of p-values ≤ 0.05 and ≥0.95 *)
Text@Row@{"Out of ",Length[pValues100],"  p-values ",Count[pValues100,_?(#≤0.05&)]," are less or equal 0.05 and ",Count[pValues100,_?(#≥0.95&)],"  p-values are greater or equal 0.95."}

pos=Position[pValues100,0.00008228157715383714`]

*3σ Rule and 2σ Rule

Probability[Abs[X-μ]<3σ,X\[Distributed]NormalDistribution[μ,σ]]//N

Probability[Abs[X-μ]<2σ,X\[Distributed]NormalDistribution[μ,σ]]//N

Compound Poisson Process with drift

Let {N(t),t≥0} be a homogeneous Poisson process with intensity λ, independent of iid (Subscript[Z, 1]),(Subscript[Z, 2]),... . Then

({{{({{X(t):=(∑_(i=1)^(N(t)))(Subscript[Z, i])}}), t≥0}}})

is a compound Poisson process and

({{Y(t)=(Subscript[y, 0])+c t+X(t)}})

is a compound Poisson process with drift c, (Subscript[y, 0]) an initial value.

***Example

Clear[c,λ,x,μ,s,t];
CompoundPoissonProcessWithDrift[c_,λ_,μ_,y0_]:=TransformedProcess[y0+c t+x[t],x\[Distributed]CompoundPoissonProcess[λ,ExponentialDistribution[μ]],t]

characteristicsofstochasticprocess[CompoundPoissonProcessWithDrift[c,λ,μ,y0],s,t]

Skewness and kurtosis:

{Skewness[CompoundPoissonProcessWithDrift[c,λ,μ,y0][t]],Kurtosis[CompoundPoissonProcessWithDrift[c,λ,μ,y0][t]]}//FullSimplify[#,Assumptions->μ>0∧λ>0∧t>0]&

Gauss-Poisson Process

({{Y(t)=W(t)+N(t)}})

where {W(t),t>0} is a standard Wiener process and {N(t),t≥0} is a homogeneous Poisson process, W and N independent.

Slice distribution:

({{({{P(Y(t)≤x)=(∑_(k=1)^∞)P(Y(t)≤x,N(t)=k)=}})(∑_(k=1)^∞)P(W(t)≤x-k,N(t)=k)=}, {(∑_(k=1)^∞)P(W(t)≤x-k)P(N(t)=k)=(∑_(k=1)^∞)(e^(-λt))((((λt))^k)/(k!))Φ(x-k)}})

***Gauss-Compound-Poisson Jump Process with Drift

({{Y(t)=(Subscript[y, 0])+μ t+σ W(t)+X(t)}})

where {W(t),t≥0} is a standard Wiener process, {X(t),t≥0} a compound Poisson process, (Subscript[Z, 1]),(Subscript[Z, 2]),... iid independent of {N(t),t≥0}, W and X independent. Recall that (Subscript[Z, i])'s may be of arbitrary sign.

Example

Clear[λ,y0,μ,σ,s,t];
GaussCompoundPoissonProcessWithDrift[λ_,μ_,σ_,ν_,y0_]:=TransformedProcess[y0+μ t+σ w[t]+x[t],{w\[Distributed]WienerProcess[],x\[Distributed]CompoundPoissonProcess[λ,ExponentialDistribution[ν]]},t]

characteristicsofstochasticprocess[GaussCompoundPoissonProcessWithDrift[λ,μ,σ,ν,y0],s,t]

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[GaussCompoundPoissonProcessWithDrift[λ,μ,σ,ν,y0],{0,T,0.01}],ImageSize->{350,300},PlotRange->All(*{-10 λ T,100 λ T}*),PlotTheme->"Detailed",PlotStyle->Thick,PlotLabel->"Gauss Compound Poisson Process\n",InterpolationOrder->0]],
Style["time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["drift",Bold],
{{μ,2},0.1,100,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["volatility",Bold],
{{σ,50},0.1,300,0.1,Appearance->"Labeled",ImageSize->Tiny},
Style["jumps",Bold],
{{ν,0.01},0.01,0.2,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["rate parameter λ",Bold],
{{λ,20},1,100,10,Appearance->"Labeled",ImageSize->Tiny},
Style["initial value",Bold],
{{y0,10},1,200,1,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

!!! Ornstein-Uhlenbeck Process, Vasicek Model

The state X(t) of a Ornstein–Uhlenbeck process satisfies an Ito differential equation

dX(t)==θ (μ-X(t))dt +σ dW(t),

with long-term mean μ, volatility σ, and mean reversion speed θ. Part μ-X(t) is the mean reverting drift pulled to a level μ at rate θ. Here θ>0, σ>0, μ∈ℝ. A drawback of the solution is that X(t) may attain negative values.

In more details (beware of the different notation):

Oldrich Alfons Vasicek (1977): spot rate r(t), α>0

dr=α(γ-r)dt+ρdz

In contrast to the random walk (the Wiener process), which is an unstable process
and after a long time will diverge to infinite values, the Ornstein-Uhlenbeck
process possesses a stationary distribution. The instantaneous drift α(γ-r)
represents a force that keeps pulling the process towards its long-term mean γ
with magnitude proportional to the deviation of the process from the mean.
The stochastic element, which has a constant instantaneous variance (ρ^2), causes
the process to fluctuate around the level γ in an erratic, but continuous, fashion.
The conditional expectation and variance level are ...

Clear[μ,σ,s,t,θ,r0(* initial value of arbitrary sign *)];
characteristicsofstochasticprocess[OrnsteinUhlenbeckProcess[μ,σ,θ,r0],s,t]//Simplify

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[OrnsteinUhlenbeckProcess[μ,σ,θ,r0],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-5,5},PlotTheme->"Detailed",PlotLabel->"Ornstein Uhlenbeck Process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0.04},0.01,0.11,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,1},0.01,2,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Rate θ",Bold],
{{θ,0.2},0.05,0.5,0.05,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
Style["Initial value (Subscript[r, 0])",Bold],
{{r0,0.02},0.01,0.3,0.001,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

!!! Cox-Ingersoll-Ross Process (CIR)

CoxIngersollRossProcess

The state X(t) of a Cox-Ingersoll-Ross process (CIR) satisfies an Ito differential equation

dX(t)==θ (μ-X(t))dt+σ (Sqrt[X(t)]) dW(t),

with long-term mean μ, volatility σ, speed of adjustment θ, and (possible) initial condition (Subscript[x, 0])>0. Here σ>0, θ μ>0, usually μ>0. The last condition means that both θ and μ must be nonzero and of the same sign. As compared with the Ornstein-Uhlenbeck, CIR can never attain negative values. It can be shown that X(t)>0 a. s. if (σ^2)≤2θ,μ.

Clear[μ,σ,s,t,θ,r0(* initial value positive *)];
characteristicsofstochasticprocess[CoxIngersollRossProcess[μ,σ,θ,r0],s,t]

Limit of the variance:

Limit[(((σ^2) (E^(-2 θ t)) ((E^(θ t))-1) (2 r0+μ ((E^(θ t))-1)))/(2 θ)),t->∞,Assumptions->θ>0&&r0>0&&μ>0]

Skewness and kurtosis:

Text@Row@{"The limits of skewness:  ",Limit[(((Sqrt[2]) (3 r0+μ ((E^(θ t))-1)))/((2 r0+μ ((E^(θ t))-1)) (Sqrt[((θ (μ+((2 r0)/((E^(θ t))-1))))/(σ^2))]))),t->∞,Assumptions->θ>0&&r0>0&&μ>0]//TraditionalForm,Spacer[20],"and kurtosis:   ",Limit[((3 (σ^2) (1-((4 (r0^2))/(((2 r0+μ ((E^(θ t))-1)))^2))))/(θ μ))+3,t->∞,Assumptions->θ>0&&r0>0&&μ>0]}

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[CoxIngersollRossProcess[μ,σ,θ,r0],{0,T,0.01}],ImageSize->{350,300},PlotRange->{-0.05,0.05},PlotTheme->"Detailed",PlotLabel->"Cox Ingersoll Ross Process\n"]],
Style["Time horizon",Bold],
{{T,1},1,10,1/250.,Appearance->"Labeled",ImageSize->Tiny},
Style["Drift",Bold],
{{μ,0.06},0.01,0.11,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Diffusion",Bold],
{{σ,0.05},0.01,1,0.01,Appearance->"Labeled",ImageSize->Tiny},
Style["Rate θ",Bold],
{{θ,0.15},0.05,0.5,0.05,Appearance->"Labeled",ImageSize->Tiny},
Delimiter,
Style["Initial value (Subscript[r, 0])",Bold],
{{r0,0.02},0.01,0.1,0.001,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

!!! Ho-Lee model

({{({{dX(t)=θ(t)dt+σdW(t)}})}})

!Hull-White model

({{({{dX(t)=(θ(t)-α(t)X(t))dt+σ(t)(X^β)(t)dW(t)}})}})

The mean reversion level is given by θ(t)/α(t). For Vasicek model β=0, other parameters constant. For CIR, β=1/2.

***Hodnota firmy, Value of the Firm, The Cost of the Capital

[11] Brigham, E. F.: Fundamentals of Financial Management. Dryden Press. 6th edition Fort Worth 1992. Chapter 7, The Cost of Capital.

Hyperlink[Button[" Hodnota firmy, cena kapitálu "],"http://www.karlin.mff.cuni.cz/~hurt/Value_of_the_firm_WACC.pdf"]

There are some misprints in the above material but not severe for the understanding the principal idea.

Despite the explanation in the above reference is rather short, about 15 rows, it takes some time to understand the essence of the topic thoroughly.

***Misc

Skip:
1.1.5 Asymmetric Laplace Distribution, 1.1.13 Mnohorozměrné ..., 1.1.14 Mnohorozměrné ..., 1.1.15 Kopuly, 1.2 Selected Statistical Metdos ..., 1.3.2 Itô process, 1.3.6 Variance Gamma, 1.4.4 Example, 1.4.7 Example, 1.4.8 Remark, 1.4.11 Remark, the proof of Lemma 1.4.14.

Just for overview: 1.4.7 Příbuzné míry rizika

Returns & Covariance matrix {SP500, Coke, Coca-Cola, VW, Apple}

companies={"SP500","Coke","KO","F:VOW","AAPL"};
SetDirectory[NotebookDirectory[]]
returnsMonthlyWithDates=Get["returnsMonthlyWithDates"];
returnsMonthlyValues=Get["returnsMonthlyValues"];

Map[Length,returnsMonthlyValues]

Grid[{{"Expected returns",SpanFromLeft,SpanFromLeft},companies,rr=Map[Mean,returnsMonthlyValues]},Frame->All]
Grid[{{"Standard deviations",SpanFromLeft,SpanFromLeft},companies,Map[StandardDeviation,returnsMonthlyValues]},Frame->All]

We see that the selected companies did not fulfil the assumptions of the efficient market: Coke versus AAPL, e.g.

DateListPlot[returnsMonthlyWithDates,PlotRange->All,PlotLegends->companies,PlotLabel->"Monthly Returns of 5 companies (61 months)"]

Kovarianční matice

Column@{"Covariance matrix with SP500:",(vvWithSP500=Covariance[returnsMonthlyValues//Transpose])//MatrixForm}
Column@{"Covariance matrix without SP500:",(vvWithoutSP500=Covariance[returnsMonthlyValues//Rest//Transpose])//MatrixForm}
{Det[vvWithSP500],Det[vvWithoutSP500]}

Korelační matice

Grid[Map[{#,FinancialData[#,"MarketCap"]}&,{"SP500","Coke","KO","F:VOW","AAPL"}]//Transpose,Frame->All]

companies={"SP500","Coke","KO","F:VOW","AAPL"};

Prepend[Prepend[(corrWithSP500=Correlation[returnsMonthlyValues//Transpose]),companies]//Transpose,{" ",companies}//Flatten]//Grid[#,Frame->All]&
(corrWithoutSP500=Correlation[returnsMonthlyValues//Rest//Transpose])//MatrixForm
{Det[corrWithSP500],Det[corrWithoutSP500]}

Optimal portfolio: risky assets only and short sales allowed

rr//Sort

Clear[x,xx,r];iN=Length[rr];xx=Array[x,iN]
portfolioreturn=rr.xx;
portfoliovariance=xx.vvWithSP500.xx//Simplify;

(i) minimize risk

Target expected return μ. In recent versions the problem may be solved explicitly.

Clear[optimalportfolio,optimalportfolioSD]
solution=Minimize[{portfoliovariance,portfolioreturn==μ,Total[xx]==1},xx]//Simplify
optimalportfolioSD[μ_]=(Sqrt[solution[[1]]])
optimalportfolio[μ_]=xx/.solution[[2]]

Manipulate[ParametricPlot[{optimalportfolioSD[μ],rr.optimalportfolio[μ]},{μ,Min[rr],Max[rr]+0.05},PlotStyle->Thick,AspectRatio->5/8,AxesLabel->{Style["(Subscript[σ, p])",Italic,14],Style["(Subscript[r, p])",Italic,14]},AxesStyle->Directive[12],PlotPoints->50,PlotRange-> {{0.039,0.08},{0.007,0.035}},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],optimalportfolio[μ]},Frame->All],"\nExpected return = ",rr.optimalportfolio[μ],"   Risk = ",optimalportfolioSD[μ],"\nRed point represents the above portfolio.\nGreen point represents the Sharpe portfolio."}, Epilog->{Red,PointSize[0.05],Point[{optimalportfolioSD[μ],rr.optimalportfolio[μ]}],Green,PointSize[0.05],Point[{0.0782,0.0327}]}],{{μ,Median[rr]},Min[rr],Max[rr],Appearance->"Labeled"},SaveDefinitions->True]

(ii) maximize Sharpe ratio

(* Sharpe ratio *)
solutionSharpe=Maximize[{(portfolioreturn/(Sqrt[portfoliovariance])),Total[xx]==1},xx]
Row@{"Expected return = ",erSharpe=(xxS=xx/.solutionSharpe[[2]]).rr,"   Risk = ",riskSharpe= (Sqrt[xxS.vvWithSP500.xxS])}

(* Square of Sharpe ratio *)
Maximize[{((((portfolioreturn))^2)/portfoliovariance),Total[xx]==1},xx]

(iii) maximize risk adjusted return

Clear[adjustedReturn]
adjustedReturn[rr_(*expected returns*),vv_(*covariance matrix*),λ_:0.0001(*risk aversion*)]:=
Module[{iN=Length[rr],xx=Array[x,iN],aux,solution},
aux=rr.xx-λ xx.vv.xx//Simplify;
aux=Maximize[{aux,Total[xx]==1},xx];
solution=xx/.aux[[2]];
{"Optimal portfolio: ",solution,"   Expected return: ",solution.rr,"   Risk: ",(Sqrt[solution.vv.solution])}
]
adjustedReturn[rr,vvWithSP500,10]

Scenarios:

Manipulate[ParametricPlot[{adjustedReturn[rr,vvWithSP500,λ1][[6]],adjustedReturn[rr,vvWithSP500,λ1][[4]]},{λ1,0.0001,40},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],adjustedReturn[rr,vvWithSP500,λ][[2]]},Frame->All],"\nExpected return = ",adjustedReturn[rr,vvWithSP500,λ][[4]],"   Risk = ",adjustedReturn[rr,vvWithSP500,λ][[6]],"\nRed point represents the above portfolio."}, Epilog->{Red,PointSize[0.05],Point[{adjustedReturn[rr,vvWithSP500,λ][[6]],adjustedReturn[rr,vvWithSP500,λ][[4]]}]}],{{λ,10},0.0001,40,Appearance->"Labeled"},SaveDefinitions->True]

Optimal portfolio: risky assets only and short sales not allowed

All these problems may be solved by optimization methods (mathematical programming)

(i) minimize risk

Target expected return μ. In the previous procedures we just add condition x>=0 to the constraints.

Clear[optimalportfolio,optimalportfolioSD]
solution=Minimize[{portfoliovariance,portfolioreturn==μ,Total[xx]==1,Thread[xx>=0]},xx]//Simplify;
optimalportfolioSD[μ_]=(Sqrt[solution[[1]]]);
optimalportfolio[μ_]=xx/.solution[[2]]

Manipulate[ParametricPlot[{optimalportfolioSD[μ],rr.optimalportfolio[μ]},{μ,Min[rr],Max[rr]+0.05},PlotStyle->Thick,AspectRatio->5/8,AxesLabel->{Style["(Subscript[σ, p])",Italic,14],Style["(Subscript[r, p])",Italic,14]},AxesStyle->Directive[12],PlotPoints->50,PlotRange-> {{0.039,0.08},{0.007,0.035}},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],optimalportfolio[μ]},Frame->All],"\nExpected return = ",rr.optimalportfolio[μ],"   Risk = ",optimalportfolioSD[μ],"\nRed point represents the above portfolio.\nGreen point represents the Sharpe portfolio."}, Epilog->{Red,PointSize[0.05],Point[{optimalportfolioSD[μ],rr.optimalportfolio[μ]}],Green,PointSize[0.05],Point[{0.0671,0.0279}]}],{{μ,Median[rr]},Min[rr],Max[rr],Appearance->"Labeled"},SaveDefinitions->True]

(ii) maximize Sharpe ratio

(* Sharpe ratio *)
solutionSharpe=Maximize[{(portfolioreturn/(Sqrt[portfoliovariance])),Total[xx]==1,Thread[xx>=0]},xx]
Row@{"Expected return = ",erSharpe=(xxS=xx/.solutionSharpe[[2]]).rr,"   Risk = ",riskSharpe= (Sqrt[xxS.vvWithSP500.xxS])}

(* Square of Sharpe ratio *)
Maximize[{((((portfolioreturn))^2)/portfoliovariance),Total[xx]==1,Thread[xx>=0]},xx]

(iii) maximize risk adjusted return

Clear[adjustedReturn]
adjustedReturn[rr_(*expected returns*),vv_(*covariance matrix*),λ_:0.0001(*risk aversion*)]:=
Module[{iN=Length[rr],xx=Array[x,iN],aux,solution},
aux=rr.xx-λ xx.vv.xx//Simplify;
aux=Maximize[{aux,Total[xx]==1,Thread[xx>=0]},xx];
solution=xx/.aux[[2]];
{"Optimal portfolio: ",solution,"   Expected return: ",solution.rr,"   Risk: ",(Sqrt[solution.vv.solution])}
]
adjustedReturn[rr,vvWithSP500,10]

Scenarios:

Manipulate[ParametricPlot[{adjustedReturn[rr,vvWithSP500,λ1][[6]],adjustedReturn[rr,vvWithSP500,λ1][[4]]},{λ1,0.0001,40},PlotLabel->Text@Row@{Grid[{{"Optimal portfolio",SpanFromLeft,SpanFromLeft},Map[ToString,xx],adjustedReturn[rr,vvWithSP500,λ][[2]]},Frame->All],"\nExpected return = ",adjustedReturn[rr,vvWithSP500,λ][[4]],"   Risk = ",adjustedReturn[rr,vvWithSP500,λ][[6]],"\nRed point represents the above portfolio."}, Epilog->{Red,PointSize[0.05],Point[{adjustedReturn[rr,vvWithSP500,λ][[6]],adjustedReturn[rr,vvWithSP500,λ][[4]]}]}],{{λ,10},0.0001,40,Appearance->"Labeled"},SaveDefinitions->True]

***Slovní formulace vzorců

Zákon o dani z příjmu §32, zrychlené odepisování

https://www.zakonyprolidi.cz/cs/1992-586#p32

(3) Při zrychleném odpisování majetku zvýšeného o jeho technické zhodnocení se odpisy stanoví

a) v roce zvýšení zůstatkové ceny jako podíl dvojnásobku zvýšené zůstatkové ceny majetku a přiřazeného koeficientu zrychleného odpisování platného pro zvýšenou zůstatkovou cenu,

b) v dalších zdaňovacích obdobích jako podíl dvojnásobku zůstatkové ceny majetku a rozdílu mezi přiřazeným koeficientem zrychleného odpisování platným pro zvýšenou zůstatkovou cenu a počtem let, po které byl odpisován ze zvýšené zůstatkové ceny.

(4) Při zkrácení doby odpisování podle § 30 odst. 1 u již odpisovaného hmotného majetku, kdy rozdíl mezi přiřazeným koeficientem pro zrychlené odpisování a počtem let, po které byl již hmotný majetek odpisován, je menší než 2, se pro stanovení ročního odpisu použije vždy rozdíl mezi přiřazeným koeficientem pro zrychlené odpisování a počtem let, po které byl již hmotný majetek odpisován, ve výši 2; přitom do počtu let se zahrnují pouze celá zdaňovací období podle § 21a.

***1.3.1 Example (Continued)

Here we make use of the built-in function Minimize. Without success, aborted by force, from Menu -> Evaluation -> Abort Evaluation after more then 1 hour.

Clear[x1,xx,ρ,μ];
xx={x1,x2};rr={8,14}; VV={{9,18ρ},{18ρ,36}};

(* endless *)
example131=Minimize[{xx.VV.xx,Total[xx]==1,rr.xx==μ},xx]

*** END

Using explicit results:

b1={1,1};VVi=Inverse[VV]//Simplify;iA=b1.VVi.b1;iB=b1.VVi.rr;iC=rr.VVi.rr;Δ=iA iC-(iB^2)//Simplify;
(* Basis portfolios*)
xx1=(1/iA)VVi.b1//Simplify;xx2=(1/iB)VVi.rr//Simplify;
δ=(1/Δ)iA(iC-μ iB)//FullSimplify;
xxoptimal=δ xx1+(1-δ)xx2//Simplify
{meanoptimal,varianceoptimal}={rr.xxoptimal,xxoptimal.VV.xxoptimal}//FullSimplify

***

***

returnsAssociation=Association[Thread[companies->returnsMonthlyValues]];
Keys[returnsAssociation]
Map[Length,Values[returnsAssociation]]

DateListPlot[returnsMonthlyWithDates,PlotLegends->companies,PlotRange->All,PlotStyle->Thick]

***Covariance matrix: properties

Covariance matrix of the returns is usually estimated from the historical stock exchange data. An unpleasant feature of the estimated matrix is that it is

Ill-conditioned or badly conditioned  matrix (špatně podmíněná matice)

Vaguely speaking: Small changes in the original matrix cause large changes in its inverse. It is well known that the Hilbert matrix is very badly conditioned.

Not finished yet:

h[n_]:=HilbertMatrix[n]
hDisturbed[n_,ϵ_]:=Module[{h=HilbertMatrix
h[5]//MatrixForm

Here are the determinants of the Hilbert matrix of order 20, approximate and exact:

{detNum=N[Det[h[20]],20],detExact=Det[h[20]]}

The same for the inverse matrix:

{detNumInv=N[Det[Inverse[h[20]]],20],detExactInv=Det[Inverse[h[20]]]}
detExact detExactInv

{detNum detNumInv,detExact detExactInv}

{Det[N[h[20]]] ,Det[Inverse[N[h[20]]]]}
Det[N[h[20]]] Det[Inverse[N[h[20]]]]

Det[h[20]] Det[Inverse[h[20]]]

Efficient frontier

(* SetDirectory[NotebookDirectory[]];
Import["concepts7_optimal_portfolio.gif"] *)

Note that on the next picture the point representing the riskless return is above the expected return of the global minimum variance portfolio which is in the contradiction with an assumption of the efficient market.

Clear[Δ,A,B,iC]
σ[r_]=(Sqrt[(1/Δ)(A (r^2)-2 B r+iC)])

Simplify[D[σ[r],r],Assumptions->Δ>0]//TraditionalForm

***

***

returnsAssociation=Association[Thread[companies->returnsMonthlyValues]];
Keys[returnsAssociation]
Map[Length,Values[returnsAssociation]]

DateListPlot[returnsMonthlyWithDates,PlotLegends->companies,PlotRange->All,PlotStyle->Thick]

***Covariance matrix: properties

Covariance matrix of the returns is usually estimated from the historical stock exchange data. An unpleasant feature of the estimated matrix is that it is

Ill-conditioned or badly conditioned  matrix (špatně podmíněná matice)

Vaguely speaking: Small changes in the original matrix cause large changes in its inverse. It is well known that the Hilbert matrix is very badly conditioned.

Not finished yet:

h[n_]:=HilbertMatrix[n]
hDisturbed[n_,ϵ_]:=Module[{h=HilbertMatrix
h[5]//MatrixForm

Here are the determinants of the Hilbert matrix of order 20, approximate and exact:

{detNum=N[Det[h[20]],20],detExact=Det[h[20]]}

The same for the inverse matrix:

{detNumInv=N[Det[Inverse[h[20]]],20],detExactInv=Det[Inverse[h[20]]]}
detExact detExactInv

{detNum detNumInv,detExact detExactInv}

{Det[N[h[20]]] ,Det[Inverse[N[h[20]]]]}
Det[N[h[20]]] Det[Inverse[N[h[20]]]]

Det[h[20]] Det[Inverse[h[20]]]

Efficient frontier

(* SetDirectory[NotebookDirectory[]];
Import["concepts7_optimal_portfolio.gif"] *)

Note that on the next picture the point representing the riskless return is above the expected return of the global minimum variance portfolio which is in the contradiction with an assumption of the efficient market.

Clear[Δ,A,B,iC]
σ[r_]=(Sqrt[(1/Δ)(A (r^2)-2 B r+iC)])

Simplify[D[σ[r],r],Assumptions->Δ>0]//TraditionalForm

20220105 E-mail

Mail studentům, kteří jsou přihlášeni na předmět Pokročilé partie finančního managementu (NMFM507)

Na mé webové stránce je aktualisovaný notes PPFM2021ZS.nb.

V posledních dvou přednáškách si zvláště všímejte chování cen finančních instrumentů pomocí funkce Manipulate (při studiu ji užijte aktivně). Pro posluchače, kteří se vyhýbali přednáškám sděluji, že i bez znalosti systému Mathematica mohou bez problémů aktivně notes prostudovat. Stačí buňky, které na levé straně svítí modře aktivovat buď  Shift+Enter nebo Enter na numerické klávesnici.

V případě jakýchkoli nejasností mně bezprostředně kontaktujte.

jh


Tento email v SIS / This message in SIS:
https://is.cuni.cz/studium/omne/index.php?do=email&doe=email_detail&mail_id=72187245

McNeil Factor Analysis (FA)

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

(sectorsDJI=Map[{#,FinancialData[#,"Name"],FinancialData[#,"Sector"]}&,membersDJI])//TableForm

FinancialData["PM",{"Name","Sector"}]

companiesMcNeil={"PM","KO","KODK","HPQ","INTC","MSFT","IBM","MCD","WMT","DIS"};
(compMcN=FinancialData[companiesMcNeil,{"Symbol","Name","Sector"}])//Grid[#,Frame->All]&

intersection=Intersection[membersDJI,compMcN[[All,1]]]
complement=Complement[compMcN[[All,1]],intersection]

sectorsMissinginDJI=FinancialData[complement,{"Symbol","Name","Sector"}]//Grid[#,Frame->All]&

(*
startDate={2014,1,1};endDate={2019,12,31};
returnsMonthlyValues=Map[FinancialData[#,"Return",{startDate,endDate,"Month"},"Value",Method->"Legacy"]&,companiesMcNeil];
pricesMonthlyWithDates=Map[FinancialData[#,{startDate,endDate,"Month"},Method->"Legacy"]&,companiesMcNeil];
pricesSP500=FinancialData["SP500",{startDate,endDate,"Month"},Method->"Legacy"];
*)

(*
SetDirectory[NotebookDirectory[]];
Save["returnsMonthlyValues20140101to20191231",returnsMonthlyValues]
Save["pricesSP50020140101to20191231",pricesSP500]
*)

SetDirectory[NotebookDirectory[]];
returnsMonthlyValues=Get["returnsMonthlyValues20140101to20191231"];
pricesSP500=Get["pricesSP50020140101to20191231"];

Append[returnsMonthlyValues,pricesSP500]//Dimensions
DateListPlot[pricesSP500,PlotLabel->"Prices ^SP500",PlotStyle->Thick,PlotTheme->"Detailed"]

DateListPlot[FinancialData["NYSE:KODK",startDate],PlotLabel->"Prices of Eastman Kodak"]

Principal Components Analysis (PCA)

correlationMcNeil=rr=Correlation[returnsMonthlyValues//Transpose];
Row@{"Eigenvalues of the correlation matrix (variances of the PC):\n", eigenvaluesMcNeil=λλ=Eigenvalues[correlationMcNeil]}
eigenvectorsMcNeil=Eigenvectors[correlationMcNeil];
Row@{"Percentage of the total variability explained by the first m principal components:\n", eigenvaluesMcNeilAccumulated=10 Accumulate[eigenvaluesMcNeil]}

Test of the spectral decomposition (decomposition minus the correlation matrix should be zero matrix)

Transpose[eigenvectorsMcNeil].DiagonalMatrix[eigenvaluesMcNeil].eigenvectorsMcNeil-correlationMcNeil//Chop

Factor analysis (FA)

factorLoadings=(Sqrt[eigenvaluesMcNeil])*eigenvectorsMcNeil (* in rows *)

(* 3 common factors *)
(bb1=factorLoadings)[[1;;3]]

(* Orthogonal rotation to make the loadings of the first factor positive *)bb1=(factorLoadings.(-IdentityMatrix[Length[factorLoadings]]))[[1;;3]]

GraphicsRow@{BarChart[bb1[[1]],BarOrigin->Left,PlotLabel ->"(1^st) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}],BarChart[bb1[[2]],BarOrigin->Left,PlotLabel ->"(2^nd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]}

BarChart[bb1[[3]],BarOrigin->Left,PlotLabel ->"(3^rd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]

***Iteration

(* Reduced corr matrix *)

rr1=rr-DiagonalMatrix[Diagonal[Transpose[bb1].bb1]]
λλ1=Eigenvalues[rr1]
xx1=Eigenvectors[rr1];

bb10=(Sqrt[λλ1])*xx1 (* in rows *)

bb2=bb10[[1;;3]].(-IdentityMatrix[Length[factorLoadings]])

GraphicsRow@{BarChart[bb2[[1]],BarOrigin->Left,PlotLabel ->"(1^st) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}],BarChart[bb2[[2]],BarOrigin->Left,PlotLabel ->"(2^nd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]}

BarChart[bb2[[3]],BarOrigin->Left,PlotLabel ->"(3^rd) common factor",ChartLabels->{Placed[companiesMcNeil,Center]}]