Advanced Topics of Financial Management

Interactive Demonstrations

Portfolio Theory

Returns of two perfectly negatively correlated assets (Výnosy dvou perfektně negativně korelovaných aktiv):

Example

Let us consider two assets A and B in the period of nine years with the corresponding returns (in per cent): 17, 13, 15, 20, 10, 16, 14, 12, 18 for the asset A and 13, 17, 15, 10, 20, 14, 16, 18, 12 for the asset B. The mean returns for both the assets based on these historical data are the same and equal to 15. The risks are also the same, 3.1225:

returnsA={17,13,15,20,10,16,14,12,18};
returnsB={13,17,15,10,20,14,16,18,12};
Text@Grid[{{"Mean A  ",Mean[returnsA],"   Mean B  ",Mean[returnsB]},{"Risk A  ",StandardDeviation[returnsA],"   Risk B  ",StandardDeviation[returnsB]}}//N,Frame->All]

If we invest all the unit wealth to any of the assets, we can expect a risky return 15 per cent. If we divide our unit wealth equally between the two assets, we have certain return 15 per cent over the given time interval with no risk. This is the case if the returns are perfectly negatively correlated, see the following Figure:

ListLinePlot[{returnsA,returnsB,(1/2)(returnsA+returnsB)},PlotTheme->"Detailed",PlotLabel->"Returns of the assets A, B, and the portfolio",PlotRange->{10,20}]

Efficient Frontier of the two assets portfolio (Eficientní hranice portfolia dvou aktiv):

In what follows, the doubled character denotes a vector or a matrix.

Example

Let us consider two assets 1, 2 with expected returns rr and the covariance matrix VV between returns:

rr={8,14};   VV=(({{9, 18ρ}, {18ρ, 36}}));

where ρ is the correlation between the returns. We will analyse the portfolio of the two assets with ρ as a parameter. Obviously, the risks of the assets 1 and 2 are 3 and 6, respectively. Let xx=((Subscript[x, 1]),(Subscript[x, 2])) denote a portfolio. Since (Subscript[x, 1])+(Subscript[x, 2])=1,  we can express the expected return on the portfolio

Clear[x1]
xx={x1,1-x1};
rp=rr.xx

The variance of the portfolio is

(σ2p=xx.VV.xx//Simplify)//TraditionalForm

For selected values of ρ, the dependence in risk-expected-return plane is illustrated in the following Figure for x1∈[0,1]:

A dynamic illustration (use slider to change the correlation coefficient):

Manipulate[ParametricPlot[{(Sqrt[36-72 x1+45 (x1^2)+36 x1 ρ-36 (x1^2) ρ]),14-6 x1},{x1,0,1},AxesOrigin->{0,8},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ρ]}],{{ρ,0},-1,1,Appearance->"Labeled"},SaveDefinitions->True]

Minimum Variance Portfolios with Risky Assets only (Portfolia pouze z rizikových aktiv s minimálním rizikem):

In risk-expected-return plane:

Minimum Variance Portfolios with Risky and Riskless Assets (Portfolia s minimálním rizikem při existenci bezrizikového aktiva):

Interest Rates and Yield Curves

Decomposition of the Interest Rates (Determinants of the Interest Rate (IR), 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

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

Here the only risk premium is for default.

We consider T-bonds and companies rated AAA, AA, A in years 1984 and 1991.

(* 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, migration 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. Chapman & Hall/CRC. Boca Raton 2003.

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

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

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

Another problem arising with a yield curves’ presentation is that the yields may be either declared or actually observed on the market. Here, by declared yields we mean the promised coupon rates (promised by the issuer) for usual fixed coupon bonds while the actually observed yields (i. e., the Yield To Maturity, YTM) are derived from the spot market price of the respective security.

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

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)

The following plots come from Wolfram|Alpha but at dates given in the headings :

Historical results (difficult to restore!):

December 15, 2014:

January 5, 2015:

Dec 3, 2015:

Dec 6, 2016:

Dec 5, 2017:

Nov 21, 2018:

Dec 3, 2018:

Actual:

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

EWGFM: European Working Group on Financial Modelling

Here we show that on an efficient market the yield curve cannot be flat (it means constant). Suppose on the contrary that the (annual) yields for all maturities are constant and equal to r(t)=r, 0≤t≤10 years. Let us consider three zero-coupon bonds (Subscript[B, 1]), (Subscript[B, 5]), and (Subscript[B, 10])  with maturities 1, 5, and 10 years and yields 1+r, (((1+r))^5), and (((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 weight w must fulfill 1*w+10*(1-w)=5 and thus 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,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,AxesLabel->{"i","PV"} ,PlotLabel->"PV(B)- PV(A)"]

Thus the investment B provides at least the same yield as the investment B. In an efficient market, supply and demand would drive the price of the five-year bond down and the prices of the one-year and ten-year bonds up.

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:

A possible realisation of the last table:

iVasicek7={3,4,5,6,7}/100.;
(tableVasicek7={{"i",Map[PaddedForm[#,{4,2}]&,iVasicek7]}//Flatten,{"(Subscript[PV, A])",Map[PaddedForm[#,{4,2}]&,Table[TimeValue[cfA[0.05],i1],{i1,iVasicek7}]//Map[Chop,#]&]}//Flatten,
{"(Subscript[PV, B])",Map[PaddedForm[#,{4,2}]&,Table[TimeValue[cfB[0.05,5/9],i1],{i1,iVasicek7}]//Map[Chop,#]&]}//Flatten,
{"Difference (Subscript[PV, B])- (Subscript[PV, A]) ",Map[PaddedForm[#,{4,2}]&,(Table[TimeValue[cfB[0.05,5/9],i1]/.i->i1,{i1,iVasicek7}]-Table[TimeValue[cfA[0.05],i1],{i1,iVasicek7}])//Map[Chop,#]&]}//Flatten})//TableForm;
tableVasicek7Grid=Text@Grid[tableVasicek7,Frame->All]

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.

Statistical and Numerical Methods  for Construction of the Yield Curves (Statistické a numerické metody konstrukce výnosových křivek)

Introduction

First we should point out that there is a substantial difference  between the statistical fitting and numerical approach interpolation. Approaching the problem from the statistical point of view, we try to construct a curve that ﬁts the observed data (quoted rates in our case) as close as possible with respect to a properly chosen criterion while the numerical methods are restricted to interpolation. This means that the interpolating function must take the same values at the observed quotations. From the point of financial analysis, the statistical approach sounds to be more reasonable. Market data are often influenced by wrong investors’ decisions and therefore some quoted prices or rates have no rational support. Nevertheless, some high posted managers of ﬁnancial institutions insist on the perfect match of quoted and estimated rates. So we will later present a method (fitting and interpolation by splines ) that fulﬁlls this requirement. From our point of view, however, a ﬁtted yield curve (not interpolated) is more useful for the ﬁnancial decision making because it eliminates  instantaneous mood or faults of the dealers and gives a more compact view on the term structure in question.

As we have already said, for the financial decision making and also for analysis we need information on the yields to maturities of the assets that are not available from the market.

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.

Warning interpolation 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];
(* poly3 is a polynomial of degree 21 *)
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 but 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. Resulting object is a piecewise polynomial function consisting of polynomials of the prescribed degree. 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]

Statistical approach (fitting, vyrovnání)

The principle

The principle of fitting: 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])) .

SetDirectory[NotebookDirectory[]];
Import["FM1e265_Fitting.png"]
Import["FM1e266_Fitting.png"]

Selected parametric models

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

Remind:

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

(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
"Fitted model:"
(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.

!Příklad (Nelson-Siegel)

(Subscript[g, N,S])(t;(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),c)= (Subscript[β, 0])+((((Subscript[β, 1])+(Subscript[β, 2]))(1-exp(-t/c)))/(t/c))-(Subscript[β, 2])(e^(-t/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]

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

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

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

Clear[f,R]
{R[t_]=(1/t)(∫_0^t)f[τ]\[DifferentialD]τ,D[t R[t],t]}

Motivation from the discrete case:

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

Nelson-Siegel forward curve:

D[t nelsonsiegelspot[(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),c,t],t]//Simplify//Expand//TraditionalForm

(* 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 (Nelson-Siegel)

(Subscript[g, N,S])(t;(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),c)= (Subscript[β, 0])+((((Subscript[β, 1])+(Subscript[β, 2]))(1-exp(-t/c)))/(t/c))-(Subscript[β, 2])(e^(-t/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
]

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

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

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

Clear[f,R]
{R[t_]=(1/t)(∫_0^t)f[τ]\[DifferentialD]τ,D[t R[t],t]}

Motivation comes from the discrete case:

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

Nelson-Siegel forward curve:

D[t nelsonsiegelspot[(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),c,t],t]//Simplify//Expand//TraditionalForm

(* 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}},PlotLabel->"Nelson-Siegel YC",PlotLegends->{"spot","forward"}],
{{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.

Note that the forward is not estimated directly. The estimated parameters of the spot are plugged  into formula for the forward instead.

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. Spot curve:

(Subscript[g, S,w])(t;(Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),(Subscript[β, 3]),(Subscript[c, 1]),(Subscript[c, 2]))= (Subscript[β, 0])+(1/t)[(Subscript[β, 1])(Subscript[c, 1])(1-(e^(-t/(Subscript[c, 1]))))+(Subscript[β, 2])((Subscript[c, 1])-(e^(-t/(Subscript[c, 1])))((Subscript[c, 1])+t))+(Subscript[β, 3])((Subscript[c, 2])-(e^(-t/(Subscript[c, 2])))((Subscript[c, 2])+t))]

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 is obtained as the integral mean of the forward:

Clear[swenssonspot];
(swenssonspot[b0_,b1_,b2_,b3_,c1_,c2_,t_]=(1/t)(∫_0^t)swenssonforward[b0,b1,b2,b3,c1,c2,τ]\[DifferentialD]τ)//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 extreme is determined by parameter c.

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

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}},PlotLabel->"Swensson YC",PlotLegends->{"spot","forward"}],
{{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

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

Note that again the forward is not estimated directly. The estimated parameters of the spot are plugged  into formula for the forward instead.

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]

Fitting by splines

The principle

Ilustrace je na datech   ir3

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

Procedure for fitting data by cubic splines 1

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}],PlotStyle->{Red,Thick,PointSize[0.02]}];
Show[tab1,plot2,PlotRange->Automatic]]

Examples

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]

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]

Example 3: Interest rates data

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 10 pairs. The following picture demonstrates that too small weight lead to numerical problems. But we usually do not fit the data by a straight line using cubic splines function.

x3=Take[Transpose[ir3][[1]],10];
u3=Take[Transpose[ir3][[2]],10];
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]],10];u3=Take[Transpose[ir3][[2]],10];pp3=Table[p3,{Length[x3]}];plotsplines[x3,u3,pp3] ,{{p3,0.000003,"1/λ"},0.00000000001,0.0001,0.000001,Appearance->"Labeled"},SaveDefinitions->True]

Procedure for fitting data by cubic splines 2

Source:

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

Test data

SeedRandom[249304];data249304=Table[{i,RiemannSiegelZ[i]+Sin[i]+RandomReal[NormalDistribution[0,.2]]},{i,0,25,.05}];
ListPlot[%]

The procedure

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

Observed yield curve 
(Completed with two additional observations since the procedure works only for the number of observations ≥24.)

ir33={{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`},{5840,5.56},{6205,5.57}};

In the following illustrations the smoothing constant is (10^λ) so that we manipulate with the power only!

Manipulate[smoothdata=CubicSplSmooth[data249304,10^lambda];
Show[ListPlot[data249304,PlotRange->{-5,3}],ListLinePlot[smoothdata,PlotStyle->Red]],{{lambda,2},-5,10,Appearance->"Labeled"}]

Manipulate[smoothdata=CubicSplSmooth[ir33,10^lambda];
Show[ListPlot[ir33,PlotRange->{3,6}],ListLinePlot[smoothdata,(* Mesh->All, *)PlotStyle->Red]],{{lambda,1},-5,10,Appearance->"Labeled"}]

Statistical approach neparametrický: jádrové odhady (kernel estimators)

Nonparameric approach: Kernel estimators (jádrové odhady)

Also known as Nonparametric Regression, Kernel Regression, Nonparametric Smoothing, Nonparametric fitting

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)

Source: Jeff Hamrick (2008):

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

Clear[analysis]
analysis[data1_,a_]:=Block[{input=data1,K,W,B,first,last,u,l,M,spline,designpoints,window},SeedRandom[1001];K[h_,arg_]:=If[-1≤(arg)/h≤1,(1/h)(3/4)(1-((arg/h))^2),0];W[data_,x_,h_]:=(K[h,#[[1]]-x])&/@data;B[data_,x_,h_,p_]:=Block[{d=Which[p==1,Table[{1,data[[i,1]]-x},{i,1,Length@data}],p==2,Table[{1,data[[i,1]]-x,(data[[i,1]]-x)^2},{i,1,Length@data}]],w=W[data,x,h]},LeastSquares[Sqrt[w] d,Sqrt[w] data[[All,-1]]]];first=First/@input;last=Last/@input;l=Quantile[first,0.01];u=Quantile[first,0.99];designpoints=Range[l,u,(u-l)/101];window[num_]:=(u-l)/(26-num);spline={{#},Sequence@@B[input,#,window[a],2]}&/@designpoints;M=Interpolation[spline];Show[ListPlot[input,PlotRange->{{l,u},Automatic}],Plot[M[x],{x,l,u},PlotStyle->{Red,Thick},PlotRange->All],ImageSize->{500,350}]]

Manipulate[
analysis[ir3,n],
{{n,16,"bandwidth"},0,100,0.5,Appearance->"Labeled"},SynchronousInitialization->False,TrackedSymbols->Manipulate,SaveDefinitions->True,ContinuousAction->False]

Using the yield curves for investment valuation
(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))}})

Stochastic Models of Interest Rates and Price Developments

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

Wiener process

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

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]

Geometric Brownian Motion Process

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

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

SeedRandom[13131];
Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[RandomFunction[GeometricBrownianMotionProcess[μ,σ,s0],{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],
{{s0,100},10,150,10,Appearance->"Labeled",ImageSize->Tiny},
{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},
SaveDefinitions->True,ControlPlacement->Left]

Fractional Brownian Motion

h ... Hurst coefficient

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

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

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)

ClearAll[generatorBasedOnTHeSumOf12Uniforms];
generatorBasedOnTHeSumOf12Uniforms[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,Min[pValues100]]

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

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

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

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

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

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

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

SDE: deterministic solutions via DSolve

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

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]//FullSimplify//First

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

Manipulate[Plot[(y[t,(Subscript[θ, 1]),(Subscript[θ, 2]),1])//Evaluate,{t,0,10},PlotStyle->Thick],{(Subscript[θ, 1]),0.1,10},{(Subscript[θ, 2]),0.1,10},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

Ito process 2

Example 2

μ,σ,θ,r0

example2ItoProcess[μ_,σ_,θ_,r0_]:=ItoProcess[\[DifferentialD]x[t]==θ(μ-x[t])\[DifferentialD]t+σ(Sqrt[x[t]])\[DifferentialD]w[t],x[t],{x,x0},t,w\[Distributed]WienerProcess[]]

Exercise. What is the name of this process?

example2SliceDistribution=SliceDistribution[example2ItoProcess[μ,σ,θ,r0],t]

Mean[example2SliceDistribution]

characteristicsofstochasticprocess[example2ItoProcess[μ,σ,θ,r0],s,t]

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

Financial Risk

The procedures as well as numbered references in

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

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

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

Coherent risk measure

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

Definitions of risk measures

Dispersion measure: Standard deviation (SD)

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

Dispersion measure: Mean absolute deviation (MAD)

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 technics and tricks in Mathematica coding.

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

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

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

The usual interpretation is that loss will exceed (Subscript[VaR, α])(X) 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.

Conditional Value at Risk (CVaR)

Conditional value at risk (also called Expected Shortfall, 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 interpretation.

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

ϕ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 Chandler):

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

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

Expectiles

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. 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 [***] 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[τ, ω])|)).

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)

Illustration

Data Apple

DateListPlot[
(* Prices *)Appleceny=FinancialData["AAPL",{{2011,1,1},{2018,4,16}}],Joined->True,PlotStyle->Thick,PlotLabel->"Apple Prices" ]
(* Returns *)Applevynosy=FinancialData["AAPL","Return",{{2011,1,1},{2018,4,16}},"Value"] ;
(* Losses *)
Appleztraty=-Applevynosy;
𝒟EmpiricaltApple=EmpiricalDistribution[Appleztraty]
(* SetDirectory[NotebookDirectory[]]; *)
(* Export["Appleztraty.txt",Appleztraty,"Table"]; *)

Characteristics of losses (Apple)

Row@{"Mean = ",Mean[Appleztraty]," SD = ",StandardDeviation[Appleztraty]," (Subscript[γ, 1]) = ",Skewness[Appleztraty]," (Subscript[γ, 2]) = ",Kurtosis[Appleztraty]}
Text@Row@{"Sample size = ",Length[Appleztraty],"   Min = ",Min[Appleztraty],"   Max =  ",Max[Appleztraty]}
Text@TableForm[Table[{1-α,Quantile[Appleztraty,α]},{α,{0.9,0.95,0.99}}],TableHeadings->{None,{"1 - α","VaR"}}]
SmoothHistogram[Appleztraty,PlotStyle->Thick,PlotLabel->"Apple Losses (Smooth Histogram)\n",PlotRange->All]

𝒟AppleLosses=EstimatedDistribution[Appleztraty,StudentTDistribution[μ,σ,ν]]
Show[Plot[PDF[𝒟AppleLosses,x],{x,-0.15,0.15},PlotStyle->{Thick,Red},PlotRange->All],Histogram[Appleztraty,Automatic,"PDF"]]

Row@{"p-value = ",   DistributionFitTest[Appleztraty,𝒟AppleLosses]}

Row@{"Test used above: ",DistributionFitTest[Appleztraty,Automatic,"AutomaticTest"]}

Nonparametric VaR (Neparametrická hodnota v riziku)

𝒟EmpiricaltApple=EmpiricalDistribution[Appleztraty]

Characteristics of losses (Apple, empirical distribution) and a comparison with parametric models

Row@{"Mean = ",Mean[𝒟EmpiricaltApple]," SD = ",StandardDeviation[𝒟EmpiricaltApple]," (Subscript[γ, 1]) = ",Skewness[𝒟EmpiricaltApple]," (Subscript[γ, 2]) = ",Kurtosis[𝒟EmpiricaltApple]}
Row@{Text@TableForm[Table[{1-α,Quantile[𝒟AppleLosses,α]},{α,{0.9,0.95,0.99}}],TableHeadings->{None,{"1 - α","VaR parametric"}}],"   ",Text@TableForm[Table[{1-α,Quantile[𝒟EmpiricaltApple,α]},{α,{0.9,0.95,0.99}}],TableHeadings->{None,{"1 - α","VaR nonparam."}}]}

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

Similar risk measures

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 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. 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, 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[τ, ω])|)).

Definitions of risk measures

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

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.

Conditional Value at Risk (CVaR)

This is valid only for continuous distributions!

Conditional value at risk (also called Expected Shortfall, 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 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

ϕ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 Chandler):

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

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

Expectiles

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. 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 [***] 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[τ, ω])|)).

Valuation of Securities

Fixed Income Securities

EffectiveInterest

Clear[p,t,r,i];
pp=Array[(Subscript[p, #])&,{5}];
tt=Array[(Subscript[t, #])&,{5}];
rr=Array[(Subscript[r, #])&,{5}];
ii=Array[(Subscript[i, #])&,{5}];

Compounding f-times yearly:

Clear[r,f];
{EffectiveInterest[r,1/f],
EffectiveInterest[r,1/4]}//TraditionalForm
(* continuous compounding *)
EffectiveInterest[r,0]//TraditionalForm

Set of the interest rates (IR):

EffectiveInterest[rr,1/f]

Average IR (geometric mean, Compound annual growth rate (CAGR) corresponding to a schedule of rates)

EffectiveInterest[rr]//TraditionalForm
{%/.Thread[rr-> {.05,.065,.04,.07,.085}],Mean[{.05,.065,.04,.07,.085}]}

Effective IR, compounding  f-times yearly:

{{0,3,5,6,10},ii}//Transpose

EffectiveInterest[{{0,3,5,6,10},ii}//Transpose,1/f]//TableForm//TraditionalForm

Convert a term structure of interest rates (yield curve or spot rates (Subscript[R, t])) to a list of implied forward rates (Subscript[f, t]) and the corresponding intervals over which they are valid:

(Subscript[R, 1])=(Subscript[f, 1])
(((1+(Subscript[R, t])))^t)=(((1+(Subscript[R, t-1])))^(t-1))(1+(Subscript[f, t]))

Thread[((1/Range[5])//Reverse)->rr]

Warning! For spot IR the number of the investment periods appear while for forward IR the starting time of the respective period is indicated (the beginning of the “bridge over” period).

EffectiveInterest[{1->0.01,2->0.015,3->0.02}]

EffectiveInterest[Thread[((1/Range[5])//Reverse)->rr]]//TraditionalForm

Thread[((Range[Length[rr]])^2)->rr]

EffectiveInterest[Thread[((Range[Length[rr]])^2)->rr]]//TraditionalForm

EffectiveInterest[{1/12->.01,3/12->.015,6/12->.02,1->.025,2->.03,3->.032,5->.035,7->.04,10->.045}]

Explanation:

(((1+(Subscript[R, (Subscript[t, 1])])))^(Subscript[t, 1]))(((1+(Subscript[f, (Subscript[t, 1]),(Subscript[t, 2])])))^((Subscript[t, 2])-(Subscript[t, 1])))=(((1+(Subscript[R, (Subscript[t, 2])])))^(Subscript[t, 2]))

{(((((((1+0.015))^(3/12))/(((1+0.01))^(1/12)))))^(1/(3/12-1/12)))-1,(((((((1+0.02))^(6/12))/(((1+0.015))^(3/12)))))^(1/(6/12-3/12)))-1}

Non-equidistant periods, the most frequent case in the practise:

(zz={1/12->.01,3/12->.015,6/12->.02,1->.025,2->.03,3->.032,5->.035,7->.04,10->.045})/.Rule->List

forwardrates=EffectiveInterest[spotrates={1/12->.01,3/12->.015,6/12->.02,1->.025,2->.03,3->.032,5->.035,7->.04,10->.045}]
ListPlot[{spotrates/.Rule->List,forwardrates},PlotStyle->PointSize[0.015],PlotTheme->"Detailed",PlotStyle->{{Blue,Thick},{Orange,Thick}},PlotLabel->"Spot and Forward interest rates",PlotLegends->{"spot","forward"}]

Transform Spot ⟶ Forward in general:

Thread[Range[Length[rr]]-> rr]

EffectiveInterest[Thread[Range[Length[rr]]-> rr]]

Najit nominalni urokovou miru odpovidajici efektivni urokove mire pri ctvrtletnim uroceni:

Find the nominal IR given the effective IR:

FindRoot[EffectiveInterest[r,1/4]==.05,{r,.2}]

TimeValue

Basic use

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

Details

Cashflow enables to enter calendar data:

Settlement date (valuation date) {2009,1,1}, valuation interest rate 0.08:

TimeValue[Cashflow[cf={{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],.08,{2009,1,1}]
GraphicsRow@{DateListPlot[cf,PlotStyle->Thick],
Plot[TimeValue[Cashflow[cf={{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2009,1,1}],{r,0,.2},PlotLabel->"Profil PV",PlotTheme->"Marketing"]}

IRR  for two cashflows with different settlement dates (Why the IRR’s are the same?

{FindRoot[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2000,1,1}],{r,0.1}],FindRoot[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2010,1,1}],{r,0.1}]}

TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2000,1,1}]
NSolve[TimeValue[Cashflow[{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2000,1,1}]==0,r]

But if the powers are not integers we see another picture. We change {2010,1,1} ⟶ {2010,1,15}. Why we cannot use NSolve?

timevalue2=TimeValue[Cashflow[{{{2010,1,15},-600},{{2011,1,1},300},{{2012,1,1},400}}],r,{2000,1,1}]
NSolve[timevalue2==0,r]

We must use FindRoot:

FindRoot[timevalue2==0,{r,0.2}]

Cashflow

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

Promenne urokove miry po obdobich

Variable IR in the respective periods:

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

Urokove miry ze spotove vynosove krivky

Spot IR (from the yield curve):

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

Annuity, AnnuityDue

Clear[PMT,i,T,v]

Present value (PV) of the annuity immediate (payments at the end of payments periods), PMT... regular payment, T... number of payments, i... valuation interest rate, v=1/(1+i) ... discount factor

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

PV of a perpetuity:

Limit[TimeValue[Annuity[PMT,T],i,0]//Simplify,T->∞]

Limit[((PMT-(((1+i))^(-T)) PMT)/i),T->∞,Assumptions->i>0]

Future value of the annuity immediate (payments at the end of payments periods):

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

Present value of the annuity due (payments at the beginning of payments periods):

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

(pvAnnuityDue/pvAnnuityImmediate)//Simplify

Future value of the annuity due (payments at the beginning of payments periods):

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

FinancialBond

FinancialBond

Notation: “FaceValue” F, “Coupon” coupon rate c (sometimes called interest), “Maturity”  T, “Settlement” valuation date t, “InterestRate” valuation interest rate (IR) r, market IR for a comparable investment, cost of capital, “CouponInterval” yearly frequency of coupon payments 1/f, (f=2 means twice a year, e.g.), “DayCountBasis” day count convention, “AccruedInterest” accrued interest (AI), “Value” pure, clean price (without AI) , “FullValue” fair, dirty price price (including AI)

Fair price = Clean price + Accrued interest
FullValue = Value + AccruedInterest

BEWARE OF SYNTAX!!! IF USING THIS AND OTHER SIMILAR FUNCTIONS GO TO HELP!!! (F1)

Properties:

({{"Value", price adjusted for accrued interest}, {"FullValue", unadjusted price}, {"AccruedInterest", accrued interest at settlement}, {"Duration", Macaulay duration}, {"ModifiedDuration", modified duration}, {"Convexity", convexity}, {"CouponPeriodDays", days in coupon period containing the settlement}, {"CouponToSettlementDays", days from previous coupon to settlement}, {"SettlementToCouponDays", days from settlement to next coupon}, {"NextCouponDate", next coupon date}, {"PreviousCouponDate", previous coupon date}, {"RemainingCoupons", remaining coupon payments}})

1 Example

Actual valuation interest rate is 6%. Calculate the clean price of the coupon bond with maturity 30 years, yearly coupon rate 5.75% and face value 1000 CZK?

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->30},{"InterestRate"->0.06,"Settlement"->0}]

Actual valuation interest rate is 6%. Calculate the fair price of the coupon bond with maturity 30 years, yearly coupon rate 5.75% and face value 1000 CZK? (Add property “FullValue” as the last parameter of the function FinancialBond.)

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->30},{"InterestRate"->0.06,"Settlement"->0},"FullValue"]

2 Example

Beware of terminology, in this example “yield”does not mean the yield of the bond (as YTM, Yield To Maturity) but the required gain of the investor.

(Clean) price of a 10-year semiannual coupon bond with a 5% yield (= Valuation IR), 9 months after the issue date:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->10,"CouponInterval"->(1/2)},{"InterestRate"->0.05,"Settlement"->(9/12)}]

Fair value of the same bond:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->10,"CouponInterval"->(1/2)},{"InterestRate"->0.05,"Settlement"->(9/12)},"FullValue"]

Accrued interest:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.0575,"Maturity"->10,"CouponInterval"->(1/2)},{"InterestRate"->0.05,"Settlement"->(9/12)},"AccruedInterest"]

3 Example

Price of a 4% quarterly coupon bond maturing December 31, 2030, and settling September 5, 2013 with  r=0.06:

FinancialBond[{"FaceValue"->1000,"Coupon"->0.04,"Maturity"->{2030,12,31},"CouponInterval"->(1/4)},{"InterestRate"->.06,"Settlement"->{2017,1,3}}(* ,"Value" *)]

Graph of the dependence on the valuation interest rate:
(Using Evaluate is almost necessary in this case.)

Plot[FinancialBond[{"FaceValue"->1000,"Coupon"->0.04,"Maturity"->{2030,12,31},"CouponInterval"->(1/4)},{"InterestRate"->r,"Settlement"->{2017,1,3}}]//Evaluate,{r,0.01,0.20},PlotPoints->200,AxesOrigin->{0,0}]

4a  Example (Accrued interest, different calendars)

Accrued interest under calendar conventions  Automatic, “30/360”, “Actual/365”, "Actual/360".

Map[{#,FinancialBond[{"FaceValue"->1000,"Coupon"->0.07,"Maturity"->{2030,12,31},"CouponInterval"->(1/2)},{"InterestRate"->0.06,"Settlement"->{2017,11,15},"DayCountBasis"->#},"AccruedInterest"]}&,{Automatic,"30/360","Actual/365","Actual/360"}]//TableForm

4b  Priklad (Accrued interest, different valuation dates)

Table[{settle,FinancialBond[{"FaceValue"->1000,"Coupon"->0.07,"Maturity"->{2030,12,31},"CouponInterval"->(1/2)},{"InterestRate"->0.06,"Settlement"->settle,"DayCountBasis"->"Actual/360"},"AccruedInterest"]},{settle,{{2010,11,12},{2010,12,12},{2010,12,30},{2010,12,31}}}];
Grid[Prepend[%,{"Settlement","Accrued Interest"}],Dividers->{{False,True},{False,True}}]

Seemingly looking paradox (AI greater than coupon) is due the used calendar, Actual/360 in this case.

ListPlot[Table[FinancialBond[{"FaceValue"->1000,"Coupon"->0.07,"Maturity"->{2030,12,31},"CouponInterval"->(1/2)},{"InterestRate"->0.06,"Settlement"->DayPlus[{2010,11,0},d],"DayCountBasis"->"Actual/360"},"AccruedInterest"],{d,1,100}]]

5 Example

Yield To Maturity (YTM, Yield To Maturity) is  called “Implied yield to maturity”) here.  Calculation of the YTM for the bond traded for 900 at 2013-5-3.

FindRoot[FinancialBond[{"FaceValue"->1000,"Coupon"->0.05,"Maturity"->{2018,6,31},"CouponInterval"->(1/4)},{"InterestRate"->ytm,"Settlement"->{2013,5,3}}]==900,{ytm,.1}]

6 Example

specifiValues1={F->1000,c->0.05,T->10,f->1,r->0.06,t->0};
specificValues2={F->1000,c->0.05,T->10,f->1,r->0.06,t->1/2};

It is possible to obtain general formulas:

FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->T,"CouponInterval"->(1/f)},{"InterestRate"->r,"Settlement"->t}]//Apart//Simplify//TraditionalForm
(* Pro konkretni hodnoty: *)
%/.specifiValues1
%%/.specificValues2

Value

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

FullValue

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

Here the expected result is c,F{f,t}:
({ft} is the fractional part of f,t)

fullvalue-value//FullSimplify
%/.specifiValues1
%%/.specificValues2

Accrued interest

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

7 Example (Duration, Convexity)

FinancialBond[{"FaceValue"->F,"Coupon"->c,"Maturity"->T,"CouponInterval"->(1/f)},{"InterestRate"->r,"Settlement"->0},"Duration"]//FullSimplify
(* pro roční frekvenci *)
(%/.f->1)//FullSimplify//TraditionalForm

Similarly for Convexity: Instead of  “Duration” use “Convexity”.

8 Example (Valuation using interest rates from the yield curve)

yieldCurve=Thread[Range[5]->Table[(Subscript[R, i]),{i,5}]]

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

Financial Derivatives

Introduction

Arbitráž a zajišťování (hedging)

All the models treated assume no-arbitrage principle, in other words, the absence of arbitrage opportunities. By an arbitrage opportunity we mean any of the two situations:
(1) At the same time, the same asset is sold at different prices at different places. Nowadays, this can hardly happen in the financial world since the information from one stock exchange is available on the stock exchange on the opposite side of the globe within a second.
(2) With zero investment at time 0 there is no probability of loss but there is a possibility of a riskless profit at time 1. More rigorously, an arbitrage opportunity in this case is a self-financing trading strategy with no initial investment, and a positive probability of positive profit and zero probability of negative profit later on.

Arbitrage opportunity is often characterized as a “money pump” and no-arbitrage
principle by the slogan: “There is no such thing like a free lunch.”.

Velmi často se setkáme s následujicí, poněkud nepřesnou terminologií.
Investor je v krátké pozici, short position, jestliže nějaké aktivum nevlastní, ale má závazky z něj plynoucí. Obvykle se jedná o emitenta dluhopisu nebo prodávajicího opce. V dlouhé pozici, long position je investor, ktery je držitelem (kupujícim) nějakého aktiva, ze kterého pro protistranu plynou nějaké závazky. V připadě opcí investorovi v krátké pozici nezbývá než trpně přihližet, jak se bude trh vyvíjet, a jak se protistrana zachová. V případě racionálně se chovajícího se investora v dlouhé pozici je jeho rozhodnutí determinováno trhem.

Zkráceně se někdy dlouhá pozice značí “+“, krátká “-“.

Příklad arbitráže

Nulová investice v čase 0 a nenulová pravděpodobnost zisku v čase 1 a nulová pravděpodobnost ztráty v čase 1.

Cena akcie v čase t je (Subscript[S, t]), t=0,1. V čase t=0 si půjčím částku B=-(Subscript[S, 0]), (úroková míra r) a za ní  si koupím akcii. Moje bilance (hodnota portfólia) v čase t=0 je (Subscript[S, 0])+B=0. V čase t=1 je hodnota portfólia (Subscript[S, 1])+(1+r)B. Jestliže P((Subscript[S, 1])+(1+r)B>0)>0 a zároveň P((Subscript[S, 1])+(1+r)B<0)=0, jedná se o arbitráž.

Bezriziková úroková míra (riskless, riskfree)

Abstraktní pojem, předpokládá se že existuje úroková míra (Subscript[r, 0]), která je stejná pro dlužníky, vypůjčovatele (borrowers) i věřitele (lenders). Investice při této ú. m. poskytují zaručený, jistý výnos. S touto úrokovou mírou se setkáváme v řadě modelů. V praktických výpočtech se obvykle nahrazuje (proxy) úrokovou mírou výnosu státních dluhopisů (poukázek) pro srovnatelnou dobu splatnosti. V současnosti při tzv. “kvantitavním uvolňování” (rozuměj tištění nekrytých bankovek) jak v USD tak EUR zóně je možné o bezrizikovosti s úspěchem pochybovat. To také platí při iracionálních zásazích centrální banky.

***Záleží na kontextu, stejná pro dlužníky, vypůjčovatele (borrowers) i věřitele (lenders) je jistě nereálné, navíc transakční náklady.

Forwardy a futures

Na rozdíl od opcí je smlouva mezi stranami závazná. Forwardy jsou individuálně uzavírané kontrakty o budoucím vypořádání, futures standardizované.

4.3.2 Termínové obchody: Forwardy a futures

Připomeňme nejprve mechanismus fungování opcí. U nich investor (držitel, kupující) v dlouhé pozici má právo, nikoli povinnost své právo zakotvené v opční smlouvě uplatnit. To znamená,že podkladové aktivum (underlying) může buď koupit (v případě CALL opce) nebo prodat (v případě PUT opce). Prodávající (pisatel) pak trpně čeká, jak se vyvíjí tržní situace, a jestli držitel své právo uplatní. Ten tak ovšem učiní v případě, že uplatněním opce dosáhne zisku nebo alespoň bude minimalizovat ztrátu vzniklou náklady na pořízení opce.

U kontraktů typu forward a futures není možná volba jako u opcí. Pro obě strany je smlouva závazná, tudíž tady v okamžiku uzavření kontraktu neexistuje cena kontraktu v pravém slova smyslu. Forwardový kontrakt je termínový obchod, smlouva mezi dvěma partnery (řekněme A a B), ve kterém se strana kupující A (v dlouhé pozici) zavazuje koupit od strany prodávající B (v krátké pozici) určitou komoditu v určeném množství za určitou cenu F fixovanou v okamžiku uzavření kontraktu v určitém budoucím termínu T (dodání, delivery). Vše je závazně domluveno dnes, dejme tomu v okamžiku  t = 0, a vlastní transakce proběhne ve sjednaném termínu v budoucnosti T, t<T. Obě strany jsou povinné své závazky splnit. Předpokládejme dále, že aktuální cena podkladového aktiva v okamžiku uzavření kontraktu je (Subscript[S, 0]).

Pozice obou stran v kontraktu je symetrická. V okamžiku vypořádání zisk jedné strany jest roven ztrátě druhé strany. Jestliže  smluvená cena kontraktu v čase vypořádání T je F, pak strana A obdrží v okamžiku T částku (Subscript[S, T])-F a strana B přirozeně obdrží  F-(Subscript[S, T]). To, co se obvykle v teorii neuvažuje, jsou transakční náklady. Pro individuálně uzavřené kontrakty mohou být nezanedbatelné.

Tyto kontrakty mají minimálně tyto dvě funkce. První: ochrana (uhájení) ceny komodity, která je předmětem kontraktu. Druhá: transfer vlastnictví; fixováním ceny dnes (v čase t) je kupující chráněn (jištěn) proti případné vyšší ceně v okamžiku vypořádání kontraktu (tj. okamžité (spotové) ceně na trhu v okamžiku dodání). Prodávající je chráněn (jištěn) proti případné nižší promtní ceně v okamžiku vypořádaní kontraktu (tj. okamžité ceně na trhu v okamžiku dodání).

Obchodování s forwardy není prakticky regulováno, jsou to obchody individuální, takže riziko z nesplacení (default) je nezanedbatelné. V tomto případě se rizikem selhání myslí to, že strana v dlouhé pozici cukne, nebude chtít koupit nebo strana v krátké pozici odmítne prodat. Toto riziko je z velké části eliminováno standardizací. Není možné postihnout všechny individuální obchody a převést je ke globálnímu obchodování, nicméně je mnoho komodit (obilí, ropa, strategické kovy, měny, . . . ), které je možné dovést ke standardizovanému obchodování. Standardizované obchody tohoto typu se provádějí na specializované burze a nazývají se futures.

Futures: Z matematického hlediska jsou forwards a futures prakticky totožné. Proto se na hodnocení futures vztahují stejná kritéria jako na forwardy. Dále budeme mluvit o hodnocení forwardů.

Uvedeme některé přístupy k hodnocení forwardových kontraktů. Nejprve uveďme hodnocení na základě principu neexistence arbitráže. Jako dříve předpokládáme, že cena podkladového aktiva v čase t=0 je (Subscript[S, 0]) a z držení aktiva neplynou žádné závazky např. ve formě dividend. Dále uvažujeme bezrizikovou hodnotící úrokovou míru r.

Strana B, krátká v kontraktu, jest povinnou doručit (deliver) částku (Subscript[S, T]) v okamžiku T. Strana B si vypůjčí částku (Subscript[S, 0]) v čase t=0 (čas uzavření kontraktu), koupí předmětné aktivum a v čase T obdrží F. Touto částkou (F) splatí půjčku, která v čase T je ve výši budoucí hodnoty (Subscript[S, 0])(e^(r,T)) , takže forwardová cena je

F=(Subscript[S, 0])(e^(r,T))

V okamžiku t=0 si strany nevyměňují žádné peníze, hodnota kontraktu je tudíž nulová.
V libovolném čase 0≤t≤T je však hodnota kontraktu

(Subscript[f, t])=(Subscript[S, t])-F,(e^(-r(,T-t)))

V takovém okamžiku může každá ze stran „uzamknout, lock” svou pozici tak, že uzavře pozici opačnou. Může tak fixovat zisk nebo ztrátu. Je-li (Subscript[S, t])<F,(e^(-r(T-t))), strana B si půjčí (Subscript[S, t]), koupí akcii, v okamžiku T splatí (Subscript[S, t])(e^(r(T-t))) a akcii prodá za F. Je-li naopak (Subscript[S, t])>F,(e^(-r(T-t))), strana A prodá akcii nakrátko, dostane za ni (Subscript[S, t]), tuto částku uloží a v okamžiku T vybere (Subscript[S, t])(e^(r(T-t)))>F. Koupí akcii za F a akcii vrátí.

4.3.9 Poznámka. (Ekvivalentní portfolio) Výplatní funkce forwardu v dlouhé pozici je (Subscript[S, T])-F. Dlouhá pozice ve forwardu je ekvivalentní v dlouhé evropské CALL a krátké pozici v evropské PUT, obě se stejnou dobou vypršení, stejnou realizační cenou K=F jako má forwardový kontrakt. Taková kombinace opcí má stejnou výplatní funkci jako forward, a tedy i v okamžiku uzavření kontraktu musí mít stejnou hodnotu. V čase t=0 (okamžik uzavření kontraktu) je však hodnota kontraktu nulová (nedochází k žádné výměně peněz), jinými slovy (Subscript[c, 0])-(Subscript[p, 0])=0. Ale z PUT-CALL parity (Subscript[c, 0])-(Subscript[p, 0])=(Subscript[S, 0])-F,(e^(-r,T)). Z toho tedy vyplývá vzorec pro F.

With[{K=50},GraphicsRow[{Plot[Max[0,S-K],{S,0,100},PlotStyle->{Thick,Red},PlotLabel->"CALL",AxesLabel->{"(Subscript[S, T])","(Subscript[c, T])"}],
Plot[Max[0,K-S],{S,0,100},PlotStyle->{Thick,Orange},PlotLabel->"PUT",AxesLabel->{"(Subscript[S, T])","(Subscript[p, T])"}],
Plot[Max[0,S-K]-Max[0,K-S],{S,0,100},PlotStyle->{Thick,Blue},PlotLabel->"Forward",AxesLabel->{"(Subscript[S, T])","(Subscript[f, T])"}]}]]

4.3.10 Poznámka. (Hodnocení pomocí Black–Scholesovy formule) Při rizikově neutrálním
oceňování a předpokladu, že cena akcie se řídí Black–Scholesovým vzorcem, je E((Subscript[S, T])|(Subscript[S, 0]))=(Subscript[S, 0])(e^(r,T)).  Výplatní funkce forwardu je (Subscript[S, T])-F. Střední hodnota této funkce podmíněná známou cenou akcie v okamžiku uzavření kontraktu t=0 a diskontovaná k okamžiku 0 je tedy (Subscript[S, 0])-F,(e^(-r,T)). Tento výraz je však roven nule, neboť v okamžiku uzavření kontraktu je cena nulová.
(WILMOTT Students, p. 99)

Options

Evropské a americké opce (vanilla options, obyčejné opce)

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

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

v okamžiku vypršení

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

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

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

v okamžiku realizace τ

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

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

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

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

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

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

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

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

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

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

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

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

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

Black-Scholesův vzorec

Risk Neutral Valuation (rizikově neutrální ptostředí, rizikově neutrální svět)

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

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

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

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

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

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

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

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

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

Then

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

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

Similarly for the variance:

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

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

Viz [1], str. 55-56.

Hodnota (cena) evropské opce

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

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

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

ct[[2,2]]//FullSimplify

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

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

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

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

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

Citlivosti (Comparative Statics, The Greeks)

Statics je sice vžitý, ale asi chybně použitý termín nevystihující podstatu. Greeks: charakteristiky (až na jednu) se označují řeckými písmeny.

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

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

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

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

Pro CALL resp. PUT dostaneme

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

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

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

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

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

Implied Volatility (Implikovaná volatilita)

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

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

?*ARCH*

?*AR*

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

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

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

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

Portfolia opcí

Příklad v [1], str. 62-63.

SetDirectory[NotebookDirectory[]]

Import["FM1aPortfolia_opci.jpg"]

Illustration of a cap:

With[{K=50,B=70},Plot[Max[0,S-K]-Max[0,S-K-B],{S,0,200},PlotTheme->"Business"]]