ATFM Fall Term 2019/20

NMFM507 Advanced Topics of Financial Management

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

$LicenseExpirationDate

Introduction

Exams

The terms will be  in SIS.

2020-01-16, 2020-01-23, just one of 2020-02-06 and 2020-02-13: 11AM Lecture room K11.
Exam in written, duration about 1 hour, 2 hours booking recommended. Results available in SIS within 50 hours.

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.

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

[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.
http://www.karlin.mff.cuni.cz/~hurt/Financni_management_prednaska.nb

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

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

[15] Hurt, J.: Theory of Portfolio. Chapter In: Financial Managament. MATFYZPRESS (under preparation). http://www.karlin.mff.cuni.cz/~hurt/20190513_Teorie _portfolia _Template _mfpress _ 1.pdf

[16] Hurt, J.: Arbitrage Pricing Theory. Chapter In: Financial Managament. MATFYZPRESS (under preparation).
http://www.karlin.mff.cuni.cz/~hurt/20191109_APT _Template _mfpress _ 1.pdf

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

2019-10-03 Lesson 1

Theory of Portfolio

The whole chapter available at

Hyperlink[Button[" Theory of portfolio "], 
    "http://www.karlin.mff.cuni.cz/~hurt/20190513_Teorie_portfolia_Template_mfpress_1.pdf"]

Diversification, portfolio, efficient market, short sales, market and efficient portfolio

(* 
SetDirectory[NotebookDirectory[]];
Import["20161007_Template_mfpress_1_Stránka_1.png"]
Import["20161007_Template_mfpress_1_Stránka_2.png"]
Import["20161007_Template_mfpress_1_Stránka_3.png"]
Import["20161007_Template_mfpress_1_Stránka_4.png"]
Import["20161007_Template_mfpress_1_Stránka_5.png"]
*)

Covariance matrix

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 (ill-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.

h[n_]:=HilbertMatrix[n]
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]]]

Konec 2019-10-03

2019-10-10, 2019-10-17 Lesson 2, 3

Efficient frontier of two assets (Eficientní hranice portfolia dvou aktiv)

Hyperlink[Button[" Example of a portfolio of two assets "], 
    "http://www.karlin.mff.cuni.cz/~hurt/20171006_Example_2_assets.pdf"]

Formulation in Wolfram Language

Two assets, the correlation between the returns ρ. Doubled characters in formulas mean vectors or matrices. Thus xx, rr, or vv in a formula should be read as x, r, or V, respectively.

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

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

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

{rp=rr.xx//Simplify,σp2=xx.VV.xx//Expand,σp=(Sqrt[σp2])}

Minimize the variance (positive definite quadratic form ⟹ equal its derivative to zero):

D[σp2,x1]==0

(xxstar={(x1star=x1/.Solve[D[σp2,x1]==0,x1]//ExpandNumerator)//First,x2star=1-x1star//Simplify//First})//TraditionalForm

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

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

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]

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.

General solution (risky assets, short sales allowed, Mean-Variance)

See also the more recent hyperlink in Lesson 1.

Hyperlink[Button[" General solution (risky assets, short sales allowed, Mean-Variance) "], 
    "http://www.karlin.mff.cuni.cz/~hurt/20171010_Portfolio_General_Solution.pdf"]

General solution (Sharpe, riskless and risky assets, short sales allowed)

See also the more recent hyperlink in Lesson 1.

Hyperlink[Button[" Sharpe ratio\n General solution  (riskless and risky assets, short sales allowed) "],"http://www.karlin.mff.cuni.cz/~hurt/20171017_Portfolio_General_Solution2.pdf"]

Capital Asset Pricing Model

Security market line

Capital market line

Konec 2019-10-10

2019-10-24 Lesson 4

Beta Estimation from Monthly Returns

Here Market portfolio = 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]

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

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

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

linearModelsCompanies=Association[Thread[Rest[companies]->Map[LinearModelFit[Transpose[{returnsAssociation["SP500"],returnsAssociation[#]}],x,x]&,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 equation (Subscript[r, n])=(OverHat[(Subscript[α, n])])+(OverHat[(Subscript[β, n])])(Subscript[μ, M])  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] }

DJI 2018-10-29

membersDJI20181029=FinancialData["^DJI","Members"];
Text@Grid[Map[{#,FinancialData[#,"Name"]}&,membersDJI],Frame->All]

FinancialData[membersDJI20181029[[1]],{2015,

DJI 2017-10-31

membersDJI={"3M","American Express","Apple","Boeing","Caterpillar","Chevron","Cisco","Coca-Cola","Disney","DuPont","Exxon Mobil","General Electric","Goldman Sachs","Home Depot","Intel","IBM","Johnson & Johnson","JPMorgan Chase","McDonald's","Merck & Co.","Microsoft","Nike","Pfizer","Procter & Gamble","Travelers","UnitedHealth","United Technologies","Verizon Communications","Visa","Wal-Mart Stores"};
Clear[MMM,AXP,AAPL,BA,CAT,CVX,CSCO,KO,DIS,DD,XOM,GE,GS,HD,INTC,IBM,JNJ,JPM,MCD,MRK,MSFT,NKE,PFE,PG,TRV,UNH,UTX,VZ,V,WMT]
tickersDJI=Map[ToString,{MMM,AXP,AAPL,BA,CAT,CVX,CSCO,KO,DIS,DD,XOM,GE,GS,HD,INTC,IBM,JNJ,JPM,MCD,MRK,MSFT,NKE,PFE,PG,TRV,UNH,UTX,VZ,V,WMT}];
(* Add the returns of the market portfolio, i. e., SP500 *)tickersDJI=Join[{"SP500",tickersDJI}]//Flatten

Row@{"Starting {YY,MM}: ",startDate=DateList[][[1;;2]]-{5,0}," Ending {YY,MM}: ",endDate=DateList[][[1;;2]]}

DateObject[]
returnsMonthlyWithDates=Map[FinancialData[#,"Return",{startDate,endDate,"Month"}]&,tickersDJI];
returnsMonthlyValues=Map[FinancialData[#,"Return",{startDate,endDate,"Month"},"Value"]&,tickersDJI];

2017-10-24 The above operation timing 174 sec.

Technical note

returnsAssociation=Association[Thread[tickersDJI->returnsMonthlyValues]];
iKeys=Keys[returnsAssociation]
iLengths=Map[Length,Values[returnsAssociation]]

Sometimes it may happen that two or more companies have only one element. The above list looks like {59,59,59,1,59,...,59,1,59,59,...}. Obviously wrong. What are these companies? We use the following construction with a possible result:

Select[returnsAssociation,Length[#]==1&]
\[LeftAssociation]"HD"->Missing["NotAvailable"],"JNJ"->Missing["NotAvailable"]\[RightAssociation]

(* Select[returnsAssociation,Length[#]==1&] *)

It follows that at the moment of downloading the data of these two companies were not available.

End of the technical note

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

Next we perform the linear regression of the company return on the return of SP500 (proxy for the market portfolio) for all 30 companies. Explanation: The first element In tickersDJI are the returns SP500. We simply use Rest to work with dependent variables returns DJI only. The result is the association.

linearModelstickersDJI=Association[Thread[Rest[tickersDJI]->Map[LinearModelFit[Transpose[{returnsAssociation["SP500"],returnsAssociation[#]}],x,x]&,Rest[tickersDJI]]]]

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

parameterTabletickersDJI=Map[{#["ParameterTable"],#["RSquared"]}&,linearModelstickersDJI]

95 per cent confidence intervals for α, β:

(parameterConfidenceIntervalstickersDJI=Map[#["ParameterConfidenceIntervals"]&,linearModelstickersDJI[[1;;30]]])//Dataset

(* Estimated {α,β} for every company *) parameterstickersDJI=Map[#["BestFitParameters"]&,linearModelstickersDJI]

Recall CAPM  for asset n: (Subscript[r, n])-(Subscript[r, 0])=(Subscript[β, n])((Subscript[μ, M])-(Subscript[r, 0])). From the estimated regression equation (Subscript[r, n])=(OverHat[(Subscript[α, n])])+(OverHat[(Subscript[β, n])])(Subscript[μ, M])  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]]))&,parameterstickersDJI]),"\n\nAverage of the above:   " ,Mean[riskFreeRateEstimate] }

Histogram[riskFreeRateEstimate]

***Auxiliary calculations and actual prices

DateObject[]
Text@Grid[Prepend[ Map[{#,FinancialData[#,"Company"],FinancialData[#]}&,tickersDJI],{"Ticker","Name","Actual Price"}],Frame->All]

***Various functions of FinancialData

FinancialData["Coke","Company"]
cokePriceDaily=FinancialData["Coke","Price",{{2014,9},{2016,9}}];
cokeReturnDaily=FinancialData["Coke","Return",{{2014,9},{2016,9}}];
cokeReturnMonthly=FinancialData["Coke","Return",{{2014,9},{2016,9},"Month"}];
Text@"Other data available for the Coke company:"
FinancialData["Coke","Properties"]
GraphicsRow@{ DateListPlot[{cokeReturnDaily,cokeReturnMonthly},PlotRange->All,PlotLabel->"Coke Monthly and Daily Returns"],DateListPlot[cokePriceDaily,PlotLabel->"Coke Daily Prices"]}

***Relationship between SP500 and DJI based  on monthly returns

DateList[]

Row@{"Starting {YY,MM} = ",startDate=DateList[][[1;;2]]-{5,0},"   Ending {YY,MM} = ",endDate=DateList[][[1;;2]]}
indicesSP500DJI={"^SP500","^DJI"};
DateObject[]
Text@Grid[Prepend[ Map[{#,FinancialData[#,"Name"],FinancialData[#]}&,indicesSP500DJI],{"Ticker","Name","Actual Price"}],Frame->All]

DateObject[]
returnsMonthlyWithDates=Map[FinancialData[#,"Return",{startDate,endDate,"Month"}]&,indicesSP500DJI];
returnsMonthlyValues=Map[FinancialData[#,"Return",{startDate,endDate,"Month"},"Value"]&,indicesSP500DJI];

2015-11-06 The above operation timing 126 sec, 67 sec, ...

returnsAssociation=Association[Thread[indicesSP500DJI->returnsMonthlyValues]];
iKeys=Keys[returnsAssociation]
iLengths=Map[Length,Values[returnsAssociation]]

djiReturns=FinancialData["^DJI","Return",{startDate,endDate,"Month"},"Value"]

Sometimes it may happen that two or more companies have only one element. The above list looks like {59,59,59,1,59,...,59,1,59,59,...}. Obviously wrong. What are these companies? We use the following construction with a possible result:

Select[returnsAssociation,Length[#]==1&]
\[LeftAssociation]"HD"->Missing["NotAvailable"],"JNJ"->Missing["NotAvailable"]\[RightAssociation]

Select[returnsAssociation,Length[#]==1&]

(* Select[returnsAssociation,Length[#]==1&] *)

It follows that at the moment of downloading the data of these two companies were not available.

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

Next we perform the linear regression of the company return on the return of SP500 (proxy for the market portfolio) for all 30 companies. Explanation: The first element In indicesSP500DJI are the returns SP500. We simply use Rest to work with dependent variables returns DJI only. The result is the association.

linearModelsindicesSP500DJI=Association[Thread[Rest[indicesSP500DJI]->Map[LinearModelFit[Transpose[{returnsAssociation["SP500"],returnsAssociation[#]}],x,x]&,Rest[indicesSP500DJI]]]]

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

parameterTableindicesSP500DJI=Map[{#["ParameterTable"],#["RSquared"]}&,linearModelsindicesSP500DJI]

95 per cent confidence intervals for α, β:

(parameterConfidenceIntervalsindicesSP500DJI=Map[#["ParameterConfidenceIntervals"]&,linearModelsindicesSP500DJI[[1;;15]]])//Dataset
(parameterConfidenceIntervalsindicesSP500DJI=Map[#["ParameterConfidenceIntervals"]&,linearModelsindicesSP500DJI[[16;;30]]])//Dataset

(* Estimated {α,β} for every company *) parametersindicesSP500DJI=Map[#["BestFitParameters"]&,linearModelsindicesSP500DJI]

Recall CAPM  for asset n: (Subscript[r, n])-(Subscript[r, 0])=(Subscript[β, n])((Subscript[μ, M])-(Subscript[r, 0])). From the estimated regression equation (Subscript[r, n])=(OverHat[(Subscript[α, n])])+(OverHat[(Subscript[β, n])])(Subscript[μ, M])  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]]))&,parametersindicesSP500DJI]),"\n\nAverage of the above:   " ,Mean[riskFreeRateEstimate] }

Histogram[riskFreeRateEstimate]

Security Market Line

Security Risk Evaluation, Merill Lynch, Beta book

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

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

2019-11-07 Lesson 5

2019-11-14 Lesson 6

Reminder (Dedicated Portfolio)

A dedicated portfolio is a passively managed portfolio whose cash flows are designed to match the cash flows needed to fulfill a future obligation.

A Stochastic Model of Matching

Remark (Ruin Probability of an Insurance Company)

Yu Chen, Chun Su: Finite time ruin probability with heavy-tailed insurance and financial risks. Statistics & Probability Letters, Volume 76, Issue 16, 1 October 2006, Pages 1812-1820.

2019-11-21 Lesson 7

Interest Rates and Yield Curves

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

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:

(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

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

The transition matrix after 5 years

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

MatrixPower[transitionmatrix,5]//MatrixForm

2019-11-28 Lesson 8

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

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:

***Technical notes

DayCountConvention

?*onvention*

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

2019-12-05 Lesson 9

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

***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^degreeOfPolynomial),t]//Evaluate,{t,0,5475},PlotRange->All,PlotStyle->Thick],plot3],{{degreeOfPolynomial,1},1,Length[ir3]-1,Appearance->"Labeled"}]

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

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

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

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

***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 of 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

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 (vyrovnání, fitting, smoothing): splines 1: The principle

The principle and an example:

Hyperlink[Button[" Fitting by cubic splines "],"http://www.karlin.mff.cuni.cz/~hurt/20171219_Fitting_by_Cubic_Splines_FM1e2_1.pdf"]

***Statistical approach (vyrovnání, fitting, smoothing): splines 2: The procedure

SplinesJanHurt.nb
Reference 1: Späth, Helmuth: Algorithmen für elementare Ausgleichsmodelle. Oldenbourg Verlag, München Wien 1973 (pp. 59-61). Beware of a lot of errors in the formulas!!!

Reference 2 (added in proofs) : Fisher, M. and Zervos, D.: YieldCurve. In: Varian, H. R. (ed.) Computational Economics and Finance. Springer. New York 1996. (pp. 269-302).

Notice that with size n of the input data the system of 3n linear equations must be solved. Generally, this is simple task for Mathematica since the matrix of the system is five-diagonal, but rendering graphics of the spline functions usually takes a longer time. 

It is also useful to compare this method with the method of running medians (Tukey).

splinesJH::usage = "splinesJH[x, u, p, z] gives a list of cubic splines in the respective intervals given by the list x of independent observations, u being the list of dependent observations at points x, weights list p, and z as the argument of the spline functions. If p->∞, the method leads to usual numerical interpolation while for p->0 it leads to the least squares method."

plotsplines::usage = "plotsplines[x, u, p]  produces the plot of the spline function given the list x of independent variables, the list u of dependent variables with weights list p."

Clear[splinesJH];splinesJH[x_,u_,p_,z_]:=Module[{pp=p,dz,dyy,yy2,eq1,eq2,eq3,eq4,eq5,variables,sol,yyy,yyy2,dyyy,dyyy2,aa,bb,cc,dd,xk,n},
n=Length[x];
dx=Take[RotateLeft[x]-x,n-1];
yy=Array[y,n];
dyy=Take[RotateLeft[yy]-yy,n-1];
yy2=Array[y2,n];tt=Array[t,n];
eq1=Table[dx[[k-1]]*yy2[[k-1]]+dx[[k]]*yy2[[k+1]]+2*(dx[[k-1]]+dx[[k]])*yy2[[k]]==6*(dyy[[k]]/dx[[k]]-dyy[[k-1]]/dx[[k-1]]),{k,2,n-1}];
eq2={yy2[[1]]==0,yy2[[n]]==0};
eq3=Table[pp[[k]]*(u[[k]]-yy[[k]])==tt[[k]],{k,n}];
eq4={tt[[1]]==(yy2[[2]]-yy2[[1]])/dx[[1]],tt[[n]]==-((yy2[[n]]-yy2[[n-1]])/dx[[n-1]])};
eq5=Table[tt[[k]]==(yy2[[k+1]]-yy2[[k]])/dx[[k]]-(yy2[[k]]-yy2[[k-1]])/dx[[k-1]]==yy2[[k-1]]/dx[[k-1]]+yy2[[k+1]]/dx[[k]]-(1/dx[[k]]+1/dx[[k-1]])*yy2[[k]],{k,2,n-1}];
variables=Join[yy,yy2,tt];
sol=Solve[Join[eq1,eq2,eq3,eq4,eq5],variables];
yyy=Flatten[Partition[yy/.sol,1]];
yyy2=Flatten[Partition[yy2/.sol,1]];
dyyy=Take[RotateLeft[yyy]-yyy,n-1];
dyyy2=Take[RotateLeft[yyy2]-yyy2,n-1];
aa=dyyy2/(6*dx);
bb=1/2*Take[yyy2,n-1];
cc=dyyy/dx-1/6*dx*(4*bb+Take[RotateLeft[yyy2,1],n-1]);
dd=Take[yyy,n-1];xk=Take[x,n-1];
aa*(z-xk)^3+bb*(z-xk)^2+cc*(z-xk)+dd]

Clear[plotsplines];plotsplines[x_,u_,p_]:=Module[{s1,t,x1,tt,intervals,tab1,k,n,plot1,plot2},
(* Needs["PlotLegends`"];*)
s1=splinesJH[x,u,p,t];n=Length[x];x1=RotateLeft[x];tt=Table[t,{n}];intervals=Take[Transpose[{tt,x,x1}],n-1];tab1=Table[Plot[s1[[k]],Evaluate[intervals[[k]]],DisplayFunction->Identity,PlotStyle->Thick],{k,1,n-1}];(*plot1=Show[tab1,DisplayFunction->$DisplayFunction];*)plot2=ListPlot[Transpose[{x,u}],Joined->True,PlotStyle->{Red,Thick}];
Show[tab1,plot2,PlotRange->Automatic]]

Sorry, the above procedure has been created when the author was a junior beginner with Mathematica. The system of the equations may be created more effectively by direct symbolic calculation.

***Statistical approach (vyrovnání, fitting, smoothing): splines 3: Example 1

Since the interest rates data are too smooth we will illustrate the procedure on somewhat wilder data. Red are the data, blue the estimated spline function.

n=25;
SeedRandom[13];
uu=Range[n]+20 Sin[0.5 RandomReal[{0,1},n]-0.5];
xx=Range[n];
ppp=Table[10.500,{n}];
(* plotsplines[xx,uu,ppp]  *)
Manipulate[ppp=Table[p,{Length[xx]}];plotsplines[xx,uu,ppp],{{p,0.1,"1/λ"},0.0001,40.0001,2,Appearance->"Labeled"},SaveDefinitions->True]

***Statistical approach (vyrovnání, fitting, smoothing): splines 4: Example 2

xa=Range[40];ua=((Range[40])^3)/1000.-15 ((Range[40])^2)/1000+10 Sin[ ((Range[40])^2)];
pa=0.0001;
ppa=Table[pa,{40}];
(* plotsplines[xa,ua,ppa] *)
Manipulate[ppa=Table[pa,{Length[xa]}];plotsplines[xa,ua,ppa],{{pa,0.1,"1/λ"},0.0001,100.0001,5,Appearance->"Labeled"},SaveDefinitions->True]

Statistical approach (vyrovnání, fitting, smoothing): splines 5: Example 3
(Interest rates data)

Ilustration on the interest rates   ir3

Since there is almost impossible to see the differences between the fitted curve and the original data if we use all of them, we take just first n pairs. The following picture demonstrates that too small weight (i.e. large λ) leads to numerical problems. But we usually do not fit the data by a straight line using cubic splines function.

n=16;
x3=Take[Transpose[ir3][[1]],n];
u3=Take[Transpose[ir3][[2]],n];
(* 1/λ = *)p3=0.00000000010;
pp3=Table[p3,{Length[x3]}];plotsplines[x3,u3,pp3]

Next we see that sensitivity of the shape of the fitted curve on the weighting parameter.

Manipulate[x3=Take[Transpose[ir3][[1]],n];u3=Take[Transpose[ir3][[2]],n];pp3=Table[p3,{Length[x3]}];plotsplines[x3,u3,pp3] ,{{p3,0.00005,"1/λ"},0.0000000001,0.0001,0.000001,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.

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 (Subscript[K, h])(u)=(1/h)(3/4)(SuperPlus[(1-(u^2))])

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

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

2019-12-12 Lesson 10

Risk measures 1  (Míry rizika)

**Introduction

The procedures as well as numbered references are available in

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

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, ... .

Selected probability distributions of losses

?*Distribution

Data for Numerical Illustrations (Walt Disney)

Prices

iDISprices20070101=FinancialData["DIS",{2007,1,1},Method->"Legacy"]
%//Length
SetDirectory[NotebookDirectory[]];
Export["iDISprices20070101.XML",iDISprices20070101];

SetDirectory[NotebookDirectory[]];
iDISprices20070101=Import["iDISprices20070101.XML"];
%//DateListPlot

Returns

(iDISpricesWithoutDates20070101=iDISprices20070101[[All,2]])//Short

Extracted from the above data but here we just download them:

iDISReturns20070101=((Rest[iDISpricesWithoutDates20070101]-Most[iDISpricesWithoutDates20070101])/(Most[iDISpricesWithoutDates20070101]));
ListLinePlot[%,PlotRange->All]

Normal Distribution

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

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

Laplace distribution (double exponential)

Symmetric, higher, but constant kurtosis.

Density:

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]

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

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

Note that ν, the degrees of freedom need not be an integer!!! Symmetric distribution, arbitrary high kurtosis (or none)

PDF[StudentTDistribution[μ,σ,ν],x]//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,10},PlotLabel->"Kurtosis of Student t",AxesLabel->{"degrees of freedom",""},PlotStyle->Thick]

!!!Mixture of Distributions

Směs rozdělení (Mixture distribution)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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]

Smooth Kernel Distribution

SmoothKernelDistribution

SmoothHistogram[iDISpricesWithoutDates20070101,PlotLabel->"Prices Walt Disney"]

{SmoothHistogram[iDISReturns20070101,Automatic,"PDF",PlotLabel->"Returns Walt Disney"],SmoothHistogram[iDISReturns20070101,Automatic,"PDF",PlotLabel->"Returns Walt Disney",PlotRange->All]}

Manipulate[SmoothHistogram[iDISReturns20070101,{bw,"Epanechnikov"},"PDF",PlotLabel->"Returns Walt Disney",PlotRange->All],{bw,0.001,0.3,Appearance->"Labeled"}]

smoothDistr=SmoothKernelDistribution[iDISReturns20070101]

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

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.

Definitions of risk measures 1

Recall: 
ρ is the return (rate of return),
L=-ρ the loss.

Coherent risk measure:

SetDirectory[NotebookDirectory[]];
Import["20181218_Coherent_Risk_Measures.png"]

Dispersion measure: Standard deviation (volatility, směrodatná odchylka)

({{({{σ(ρ)=σ(L)=(Sqrt[var ρ ])=(Sqrt[var L])}})}})

in the financial context often called volatility.

Dispersion measure: Mean absolute deviation (MAD, střední absolutní odchylka)

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

Value at Risk at level 1-α (VaR, hodnota v riziku na úrovni 1-α)

The definitions of Value at Risk differ. Particularly, what concerns the significance level. Logic of VaR is to determine the value over which the loss will exceed with a small probability. Thus, in any case, our interest is in the right-hand tails of the the distributions of losses, or in the left-hand tails of profits (rates of return). In this course we suppose the significance level  1-α but in other sources often you can find α instead of 1-α.

Here is an example of the VaR in terms of  profits. The example comes from Nagy Gergely, the Graduate from Charles University, Financial & Insurance Mathematics. Licensed users may download the source code and take a benefit from the techniques and tricks in Mathematica coding.

http://demonstrations.wolfram.com/ValueAtRisk/

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

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

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

Often we can meet the objection that VaR is not a coherent risk measure, particularly, the axiom of subadditivity is not fulfilled.  This objection comes mostly from the academics but  there is a rational reason for this objection.

!!!Exercise. (Do not skip it! Important!) Find an instance when the axiom of subadditivity is not fulfilled for VaR.

Konec 2019-12-12

2019-12-19 Lesson 11

Risk measures 2  (Míry rizika)

Conditional Value at Risk (CVaR)

This is valid only for continuous distributions!

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

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

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

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

For absolutely continuous distributions we have

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

Average Value at Risk at level 1-α (AVaR)

Average value at risk is defined as

(Subscript[AVaR, α])(X)=(1/(1-α))(∫_α^1)VaR(Subscript[,, p])(X)dp=(1/(1-α))(∫_α^1)(F^(-1))(p)dp.

After substitution y=(F^(-1))(p) we get

(Subscript[AVaR, α])(X)=(1/(1-α))(∫_((Subscript[VaR, α])(X))^∞)ydF(y)=(Subscript[CVaR, α])(X).

Hence both characteristics coincide but they differ in the interpretation.

Weighted Average VaR (WVaR)

(Subscript[WAVaR, μ])(X)=(∫_0^1)AVaR(Subscript[,, λ])(X)μ(dλ),

where μ is a probability measure on [0,1].

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

Median value at risk is defined as the median of the conditional distribution given that the loss exceeded value at risk.

(Subscript[MVaR, α])(X)=Median{X|X≥(Subscript[VaR, α])(X)}.

Therefore, (Subscript[MVaR, α])(X) is a value satisfying

P(X≥ (Subscript[MVaR, α])(X)|X≥(Subscript[VaR, α])(X))=((P(X≥ (Subscript[MVaR, α])(X),X≥(Subscript[VaR, α])(X)))/(P(X≥(Subscript[VaR, α])(X))))=(1/2)

or

P(X≥ (Subscript[MVaR, α])(X),X≥(Subscript[VaR, α])(X))=(1/2)(1-α).

Since P(X≥(Subscript[VaR, α])(X))=1-α>(1/2)(1-α), (Subscript[MVaR, α])(X)≥(Subscript[VaR, α])(X), and therefore

P(X≥ (Subscript[MVaR, α])(X),X≥(Subscript[VaR, α])(X))=P(X≥ (Subscript[MVaR, α])(X))=(1/2)(1-α)

so that P(X<(Subscript[MVaR, α])(X))=((1+α)/2) and

(Subscript[MVaR, α])(X)=(Subscript[VaR, (1+α)/2])(X).

Sometimes an analogy to the median residual lifetime is considered:

(Subscript[MVaR, α])(X)-(Subscript[VaR, α])(X)=Median{X-(Subscript[VaR, α])(X)|X≥(Subscript[VaR, α])(X)}.

Quantile Value at Risk

Beware of the symbol β-quantile. It means 100β percent  quantile, "-" is not meant as minus sign!

Quantile Value at Risk (Quantile Shortfall, Tail Conditional Quantile) is defined as

(Subscript[QVaR, α,β])(X)=β-quantile of the conditional distribution{X|X≥(Subscript[VaR, α])(X)}.

It means that the probability of loss greater than QVaR given that it exceeded VaR will be 1-β.

Thus (Subscript[QVaR, α,β])(X) satisfies the relation

P(X≥ (Subscript[QVaR, α,β])(X)|X≥(Subscript[VaR, α])(X))=((P(X≥ (Subscript[QVaR, α,β])(X),X≥(Subscript[VaR, α])(X)))/(P(X≥(Subscript[VaR, α])(X))))=1-β,

and, since it has sense to consider (Subscript[QVaR, α,β])(X)≥(Subscript[VaR, α])(X) only, we have

P(X≥ (Subscript[QVaR, α,β])(X))=(1-β)(1-α),

so that (due to 1-(1-β)(1-α)=α+β-α β)

P(X< (Subscript[QVaR, α,β])(X))=α+β-α β.

It follows

(Subscript[QVaR, α,β])(X)=(F^(-1))(α+β-α β)=(α+β-α β)-quantile of X.

Spectral measures

Suppose that ϕ is a non-decreasing non-negative real function defined on [0,1], (∫_0^1)ϕ(p)dp=1, sometimes called spectrum. (ϕ is simply non-decreasing PDF with support [0,1].The quantity

(Subscript[M, ϕ])(X)=(∫_0^1)(F^(-1))(p)ϕ(p)dp

is called spectral risk measure. If F is absolutely continuous with the probability density function f and strictly increasing on its support then after substitution (F^(-1))(p)=y we obtain

(Subscript[M, ϕ])(X)=(∫_(-∞)^∞)y ϕ(F(y))f(y)dy.

In fact, ϕ is a weighting function giving higher weights to higher losses.

Examples of spectra:

Dual-power risk aversion function

ν//FullForm

ϕdual[u_,ν_]:=ν (u^(ν-1)) (* ν≥1 *)
Manipulate[Plot[ϕdual[u,ν],{u,0,1},PlotLabel->"ϕdual"],{{ν,1.5},1,10,Appearance->"Labeled"},SaveDefinitions->True]

Proportional Hazard risk aversion function

ϕhazard[u_,γ_]:=(1/γ)(((1-u))^(1/γ-1)) (* γ≥1 *)
Manipulate[Plot[ϕhazard[u,γ],{u,0,1},PlotLabel->"ϕhazard"],{{γ,1.5},1,10,Appearance->"Labeled"},SaveDefinitions->True]

Wang's risk aversion function

(Subscript[ϕ, Wang])(u,δ)=exp(δ (Φ^(-1))(u)-(δ^2)/2).

Beware of the fact that in literature this is often presented with the mistake as

(Subscript[ϕ, Wang])(u,δ)=exp(-δ (Φ^(-1))(u)-(δ^2)/2)

with wrong sign at the first term of the exponential. Check:

Integrate[Exp[δ Quantile[NormalDistribution[0,1],u]-(δ^2)/2],{u,0,1}]

ϕWang[u_,δ_]:=Exp[δ Quantile[NormalDistribution[0,1],u]-(δ^2)/2] (* δ>0 *)
Manipulate[Plot[ϕWang[u,δ],{u,0,1},PlotLabel->"ϕWang"],{{δ,1.5},0,10,Appearance->"Labeled"},SaveDefinitions->True]

Exponential risk aversion function

ϕArrowPratt[u_,k_]:=((k Exp[-k(1-u)])/(1-Exp[-k]))
Manipulate[Plot[ϕArrowPratt[u,k],{u,0,1},PlotLabel->"ϕArrowPratt"],{{k,1},0.01,10,Appearance->"Labeled"},SaveDefinitions->True]

where k is the Arrow-Pratt degree of absolute risk aversion (ARA).

Let U be a utility function. Coefficient of absolute risk aversion or Arrow–Pratt measure of absolute risk-aversion (ARA) as a function of the wealth u is defined as

(Subscript[ARA, U])(u)=-((U''(u))/(U'(u))).

Here we have an example of the constant absolute risk aversion:

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

***See also Wolfram Demonstration by Chandler):

http://demonstrations.wolfram.com/ConstantRiskAversionUtilityFunctions/

Expectiles

For more details see

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

We start with motivation given in [***] Yao (1996) despite the idea originates in [***] Newey and Powell (1987). First we give an alternative approach for defining quantiles through an optimization problem. For 0<α<1 define

(Subscript[R, α])(x)=(1-α)|x| ,   x≤0,
(Subscript[R, α])(x)=α|x| ,                 x>0.

The α-quantile (Subscript[q, α]) of X may be defined as

(Subscript[q, α])=arg (min_a) E,(Subscript[R, α])(X-a).

Obviously, for α=(1/2) the corresponding quantile is the median.

Obviously α=P(X≤(Subscript[q, α])). In order to have an analogy with expectiles, write

α=((E,I{X≤(Subscript[q, α])})/(E,1)),

where I means the indicator function.

To define expectiles, for 0<ω<1 let us define

(Subscript[Q, ω])(x)=(1-ω)(x^2) ,   x≤0,
(Subscript[Q, ω])(x)=ω (x^2) ,                 x>0.

The ω-expectile (Subscript[τ, ω]) of X is defined as the minimizer

(Subscript[τ, ω])=arg (min_a) E,(Subscript[Q, ω])(X-a),

see [***] Yao (1996), e.g. For ω=(1/2)  we have (Subscript[τ, 1/2])=E,X. Since (Subscript[Q, ω])(.) has a continuous first derivative, (Subscript[τ, ω]) satisfies the equation

E,(Subscript[L, ω])(X-(Subscript[τ, ω]))=0

where

(Subscript[L, ω])(x)=(1-ω)x ,   x≤0,
(Subscript[L, ω])(x)=ω x ,                 x>0.

Hence (Subscript[τ, ω]) satisfies

ω=((E,[|X-(Subscript[τ, ω])|I{X≤(Subscript[τ, ω])}] )/(E|X-(Subscript[τ, ω])|)).

* 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

Expectiles (expektily)

For more details see

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

We start with motivation given in [22] Yao (1996) despite the idea originates in [18] Newey and Powell (1987). First we give an alternative approach for defining quantiles through an optimization problem. Define

(Subscript[R, α])(x)=(1-α)|x| ,   x≤0,
(Subscript[R, α])(x)=α|x| ,                 x>0.

The α-quantile (Subscript[q, α]) of X may be defined as

(Subscript[q, α])=arg (min_a) E,(Subscript[R, α])(X-a).

Obviously α=P(X≤(Subscript[q, α])). In order to have an analogy with expectiles, write

α=((E,I{X≤(Subscript[q, α])})/(E,1)),

where I means the indicator function.

To define expectiles, let us define

(Subscript[Q, ω])(x)=(1-ω)(x^2) ,   x≤0,
(Subscript[Q, ω])(x)=ω (x^2) ,                 x>0.

The ω-expectile of X is defined as the minimizer

(Subscript[τ, ω])=arg (min_a) E,(Subscript[Q, ω])(X-a),

see [22] Yao (1996), e.g. For ω=(1/2)  we have (Subscript[τ, 1/2])=E,X. Since (Subscript[Q, ω])(.) has a continuous first derivative, (Subscript[τ, ω]) satisfies the equation

E,(Subscript[L, ω])(X-(Subscript[τ, ω]))=0

where

(Subscript[L, ω])(x)=(1-ω)x ,   x≤0,
(Subscript[L, ω])(x)=ω x ,                 x>0.

Hence (Subscript[τ, ω]) satifies

ω=((E,[|X-(Subscript[τ, ω])|I{X≤(Subscript[τ, ω])}] )/(E|X-(Subscript[τ, ω])|)).

Expectiles (Calculation)

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 the above formula for  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]

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

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.

Stochastic Models of Interest Rates and Price Developments 1

Introductory Notes

Calibration

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

**General Procedure for Characteristics

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 version folows.

!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];RandomFunction[WienerProcess[μ,σ]/.{μ->1,σ->0.05},{0,1,0.01}] (* generates the object TemporalData *)
%//Normal//First (* visible observations *)
%//Length (* number of observations *)

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

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]

2020-01-09 Lesson 13

Stochastic Models of Interest Rates and Price Developments 2

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,θ,(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]

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

***Ito process 1

ItoProcess

dX(t)==a(X(t),t)dt+b(X(t),t).dW(t)

is represented as ItoProcess[{a,b},x,t]]. It may be multidimensional.

Alternatively, it may be defined by a stochastic differential equation (SDE): ItoProcess[SDEquations, expr, x,t,w\[Distributed]dproc]  represents an Ito process specified by a stochastic differential equation(s) sdeqns, output expression expr, with state x and time t, driven by w following the stochastic process dproc.

Row@{Hyperlink[Button["Ito Process"],"paclet:ref/ItoProcess"],Spacer[20],Hyperlink[Button["Stochastic Differential Equation Processes"],"paclet:guide/StochasticDifferentialEquationProcesses"],Spacer[20],Hyperlink[Button["Heston Model"],"paclet:example/HestonModel"]}

Built-in functions ending with “Process”:

?*Process

?*PoissonProcess

Example 1

example1ItoProcess=ItoProcess[\[DifferentialD]x[t]==-x[t]\[DifferentialD]t+(Sqrt[1+x[t]^2])\[DifferentialD]w[t],x[t],{x,1},t,w\[Distributed]WienerProcess[]]

characteristicsofstochasticprocess[example1ItoProcess,s,t]

SeedRandom[131311];
ListLinePlot[Table[RandomFunction[example1ItoProcess,{0.,5.,0.01}],{5}],PlotStyle->Thick,PlotTheme->"Detailed",PlotRange->All]

example1SliceDistribution=SliceDistribution[example1ItoProcess,t]

Mean[example1SliceDistribution]

CovarianceFunction[example1ItoProcess,s,t]

Plot3D[CovarianceFunction[example1ItoProcess,s,t]//Evaluate,{s,0,5},{t,0,5},ColorFunction->"Rainbow"]

***SDE: deterministic solutions via DSolve

Here we deal with the deterministic part of the above SDE’s. We just skip the stochastic term.

See also: S. M. Iacus (University of Milan), Chicago, R/Finance 2011, April 29th

Statistical data analysis of financial time series and option pricing in R

Geometric Brownian motion

gBm

DSolve[{x'[t]==μ x[t],x[0]==x0},x[t],t]

Cox-Ingersoll-Ross (CIR)

Clear[y];
y[t_,θ1_,θ2_,x0_]=(x[t]/.DSolve[{x'[t]==(θ1-θ2 x[t]),x[0]==x0},x[t],t])//First//Apart//FullSimplify

y[t,(Subscript[θ, 1]),(Subscript[θ, 2]),1]

Manipulate[Plot[(y[t,(Subscript[θ, 1]),(Subscript[θ, 2]),1])//Evaluate,{t,0,10},PlotStyle->Thick],{{(Subscript[θ, 1]),1},0.1,10,Appearance->"Labeled"},{{(Subscript[θ, 2]),2},0.1,10,Appearance->"Labeled"},SaveDefinitions->True]

Chan-Karolyi-Longstaff-Sanders (CKLS)

The same as CIR, the difference is only in the stochastic term.

Nonlinear mean reversion (Aït-Sahalia)

DSolve[{x'[t]==((Subscript[α, -1])/(x[t]))+(Subscript[α, 0])+(Subscript[α, 1])x[t]+(Subscript[α, 2])((x[t])^2),x[0]==x0},x[t],t]//ToRadicals

Double Well potential (bimodal behaviour, highly nonlinear)

DSolve[{x'[t]==(x[t]-((x[t])^3)),x[0]==x0},x[t],t]

Test

D[-((E^t)/(Sqrt[-1+(E^(2 t))+(1/(x0^2))])),t]==-((E^t)/(Sqrt[-1+(E^(2 t))+(1/(x0^2))]))-(((-((E^t)/(Sqrt[-1+(E^(2 t))+(1/(x0^2))]))))^3)

Jacobi diffusion (political polarization)

DSolve[{x'[t]==-θ(x[t]-(1/2)),x[0]==x0},x[t],t]//FullSimplify

Ohrenstein-Uhlenbeck (OU)

It differs from gBm only in the stochastic term.

Radical Ohrenstein-Uhlenbeck

DSolve[{x'[t]==(θ/(x[t]))-x[t],x[0]==x0},x[t],t]//FullSimplify

Hyperbolic diffusion (dynamics of sand)

DSolve[{x'[t]==(θ/(x[t]))-x[t],x[0]==x0},x[t],t]//FullSimplify

!!! 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.

**Black-Karasinski model

({{({{d log,X(t)=(θ(t)-α(t) log,X(t))dt+σ(t)dW(t)}})}})

In practise often α(t) and σ(t) are supposed to be constant so that the model reads:

({{({{d log,X(t)=(θ(t)-α)dt+σdW(t)}})}})

**Black-Derman-Toy model

({{({{d log,X(t)=(θ(t)-α(t) log,X(t))dt+σ(t)dW(t)}})}})

is a special case of Black-Karasinski where α(t)=σ'(t)/σ(t).

Hodnota firmy, Value of the Firm, The Cost of the Capital

[11] Brigham, E. F.: Fundamentals of Financial Management. Dryden Press. 6th edition Forth 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"]

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.