Fallstudie: Barrier Optionen

Nun ein Blick auf ein ernsteres Beispiel, das die Mächtigkeit der symbolischen Rechnungen auf dem Computer illustriert. Die Formeln erlauben Rabatte und basieren auf Rubenstein und Reiner (siehe "Breaking Down the Barriers", RISK, September 1991). Wir verwenden kontinuierlichen  Zinseszins in den Variablen.

Norm[x_] := (1/2) (Erf[(x/(Sqrt[2]))] + 1); 
    Norm[z_?NumberQ] := N[0.5 (Erf[(z/(Sqrt[2]))] + 1)];

Since the compilable subset does not contain the Erf function we will use an sigmoid-approximation instead.

MyErf[x_] := 2*(1/(E^(-((4*x)/Sqrt[Pi])) + 1) - 1/2)

Plot[{Erf[x], MyErf[x]}, {x, -5, 5}, 
      PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}];

Using this function instead of Erf we get the following Norm function

Norm[x_] := 1/2*(MyErf[x/Sqrt[2]] + 1); Norm[(z_)?NumberQ] =.

Definitionen

lambda[r_, q_, sd_] := 1/2 + (r - q)/sd^2

x[p_, k_, sd_, r_, q_, t_] := 
   Log[p/k]/(sd*Sqrt[t]) + lambda[r, q, sd]*sd*Sqrt[t];

xone[p_, h_, sd_, r_, q_, t_] := x[p, h, sd, r, q, t];

y[p_, k_, h_, sd_, r_, q_, t_] := 
   Log[h^2/(p*k)]/(sd*Sqrt[t]) + lambda[r, q, sd]*sd*Sqrt[t];

yone[p_, h_, sd_, r_, q_, t_] := x[h, p, sd, r, q, t];

a[r_, q_, sd_] := (r - q)/sd^2 - 1/2; 
b[r_, q_, sd_] := Sqrt[(r - q - sd^2/2)^2 + 2*r*sd^2]/sd^2; 
z[p_, h_, sd_, r_, q_, t_] := 
   Log[h/p]/(sd*Sqrt[t]) + b[r, q, sd]*sd*Sqrt[t];

intone[p_, k_, sd_, r_, q_, t_, phi_] := 
  phi*p*Exp[-q*t]*Norm[phi*x[p, k, sd, r, q, t]] - 
   phi*k*Exp[-r*t]*Norm[phi*x[p, k, sd, r, q, t] - phi*sd*Sqrt[t]]

inttwo[p_, k_, h_, sd_, r_, q_, t_, phi_] := 
  phi*p*Exp[-q*t]*Norm[phi*xone[p, h, sd, r, q, t]] - 
   phi*k*Exp[-r*t]*Norm[phi*xone[p, h, sd, r, q, t] - 
      phi*sd*Sqrt[t]]

intthree[p_, k_, h_, sd_, r_, q_, t_, phi_, eta_] := 
  phi*p*Exp[-q*t]*(h/p)^(2*lambda[r, q, sd])*
    Norm[eta*y[p, k, h, sd, r, q, t]] - 
   phi*k*Exp[-r*t]*(h/p)^(2*lambda[r, q, sd] - 2)*
    Norm[eta*y[p, k, h, sd, r, q, t] - eta*sd*Sqrt[t]]

intfour[p_, k_, h_, sd_, r_, q_, t_, phi_, eta_] := 
  phi*p*Exp[-q*t]*(h/p)^(2*lambda[r, q, sd])*
    Norm[eta*yone[p, h, sd, r, q, t]] - 
   phi*k*Exp[-r*t]*(h/p)^(2*lambda[r, q, sd] - 2)*
    Norm[eta*yone[p, h, sd, r, q, t] - eta*sd*Sqrt[t]]

intfive[reb_, p_, h_, sd_, r_, q_, t_, eta_] := 
  reb*Exp[-r*t]*(Norm[eta*xone[p, h, sd, r, q, t] - 
       eta*sd*Sqrt[t]] - (h/p)^(2*lambda[r, q, sd] - 2)*
      Norm[eta*yone[p, h, sd, r, q, t] - eta*sd*Sqrt[t]])

intsix[reb_, p_, h_, sd_, r_, q_, t_, eta_] := 
  reb*((h/p)^(a[r, q, sd] + b[r, q, sd])*
      Norm[eta*z[p, h, sd, r, q, t]] + 
     (h/p)^(a[r, q, sd] - b[r, q, sd])*
      Norm[eta*z[p, h, sd, r, q, t] - 2*eta*b[r, q, sd]*sd*Sqrt[t]])

Up/Down und Out Puts

Die Ergebnisse können numerisch und grafisch mit dem Original Rubenstein und Reiner Paper verglichen werden. Zunächst definieren wir die Funktionen für die Sensitivität des Wertes und Preises. Beachten Sie hier, daß alle Ableitungen, welche zur Berechnung von delta, gamma, vega etc. benötigt werden, automatisch vom Computer berechnet werden.

UpAndOutPut[reb_, p_, k_, h_, sd_, r_, q_, t_] :=
If[k > h,
inttwo[p, k, h, sd, r, q, t, -1] -
intfour[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1],
intone[p, k, sd, r, q, t, -1] -
intthree[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1]]

UpAndOutPutDelta[reb_, p_, k_, h_, sd_, r_, q_, t_] =
If[k > h,
Evaluate[D[inttwo[p, k, h, sd, r, q, t, -1] -
intfour[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1], p]],
Evaluate[D[intone[p, k, sd, r, q, t, -1] -
intthree[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1],p]]];

UpAndOutPutTheta[reb_, p_, k_, h_, sd_, r_, q_, t_] =
If[k > h,
Evaluate[-D[inttwo[p, k, h, sd, r, q, t, -1] -
intfour[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1], t]],
Evaluate[-D[intone[p, k, sd, r, q, t, -1] -
intthree[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1],t]]];

UpAndOutPutGamma[reb_, p_, k_, h_, sd_, r_, q_, t_] =
If[k > h,
Evaluate[D[inttwo[p, k, h, sd, r, q, t, -1] -
intfour[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1], {p,2}]],
Evaluate[D[intone[p, k, sd, r, q, t, -1] -
intthree[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1],{p,2}]]];

UpAndOutPutVega[reb_, p_, k_, h_, sd_, r_, q_, t_] =
If[k > h,
Evaluate[D[inttwo[p, k, h, sd, r, q, t, -1] -
intfour[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1], sd]],
Evaluate[D[intone[p, k, sd, r, q, t, -1] -
intthree[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1],sd]]];

UpAndOutPutRho[reb_, p_, k_, h_, sd_, r_, q_, t_] =
If[k > h,
Evaluate[D[inttwo[p, k, h, sd, r, q, t, -1] -
intfour[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1], r]],
Evaluate[D[intone[p, k, sd, r, q, t, -1] -
intthree[p, k, h, sd, r, q, t, -1, -1] +
intsix[reb, p, h, sd, r, q, t, -1],r]]];

Jetzt lassen wir den Computer arbeiten:

Z.B. - ein weniger triviales delta:

UpAndOutPutDelta[reb,  p,  k,  h,  sd,  r,  q,  t]

Numerische Verifikation gegen die Originalergebnisse

Zunächst können wir die Formeln vergleichen und stellen fest, daß unser exaktes Modell mit den Tabellen im Paper übereinstimmt. Die folgende Tabelle stimmt exakt mit den Ergebnissen eines Beispiels überein, wobei die Originalergebnisse eine limitierte Genauigkeit haben.

MatrixForm[
Table[N[UpAndOutPut[2,100,k,103,0.20,Log[1.1],Log[1.05],t]], 
{k, 90, 110, 10}, 
{t, 0.5, 1.5, 0.25}]]

Table[N[{2,100,k,103,0.20,Log[1.1],Log[1.05],t}], 
{k, 90, 110, 10}, 
{t, 0.5, 1.5, 0.25}] // InputForm

The numerical values in this matrix deviate slightly from the ones obtained if the exact Erf function is used in the definition of the Norm function.

Fall A
 K = H (barrier = strike) und q > r (div. yield greater than risk free rate.)

In den folgenden Graphiken untersuchen wir die Struktur der Formel und der dazugehörenden delta, gamma, und vega. 
Die erste Graphik zeigt die fast flache Struktur des Wertes als Oberfläche V(S,t).An der Barriere (hier 50) ist der Wert 0. Der Strike ist ebenso 50, die Volatilität ist 0.25 (25%), und der kontinuierliche Zinseszins ist r =10%. Das Ergebnis der Dividende ist 15%. Die Zeit läuft von 0 bis 1 Jahr. Der Aktienkurs variiert von 30  bis 50 (Barrier).

Plot3D[
UpAndOutPut[0, S, 50, 50, 0.25,0.1, 0.15, t],
{S, 30, 50}, {t, 0.001, 1}, PlotRange -> All];

Diese Graphik zeigt vega, symbolisch vom Computer berechnet. Es ist überall negativ.

Plot3D[
UpAndOutPutVega[0, S, 50, 50, 0.25,0.1, 0.15, t],
{S, 30, 50}, {t, 0.001, 1}, PlotRange -> All];

Da die Oberfläche des Wertes sehr flach ist, ist delta nahe -1, wie diese Graphik zeigt.

Plot3D[
UpAndOutPutDelta[0, S, 50, 50, 0.25,0.1, 0.15, t],
{S, 30, 50}, {t, 0.001, 1}, PlotRange -> All];

Gamma ist sehr klein und NEGATIV, wie der letzte Plot illustriert.

Plot3D[
UpAndOutPutGamma[0, S, 50, 50, 0.25,0.1, 0.15, t],
{S, 30, 50}, {t, 0.001, 1}, PlotRange -> All];

Fall B
 K = H (barrier = strike) und q = r (div. yield equal to risk free rate.)

Der Strike is ebenso 50, die Volatilität ist 0.25 (25%), und der kontinuierliche Zinseszins ist r =10%. Das Ergebnis der Dividende ist 10%. Die Zeit läuft von 0 bis 3 Jahren. Der Aktienkurs variiert von 30  bis 50 (Barrier).

Plot3D[
UpAndOutPut[0, S, 50, 50, 0.25,0.1, 0.1, t],
{S, 30, 50}, {t, 0.001, 3}, PlotRange -> All];

Diese Graphik zeigt vega, symbolisch vom Computer berechnet. Es ist überall negativ ( Machinengenauigkeit der gezeichneten Funktion ).

Plot3D[
UpAndOutPutVega[0, S, 50, 50, 0.25,0.1, 0.10, t],
{S, 30, 50}, {t, 0.001, 3}, PlotRange -> {-0.01, 0.01}];

Da die Oberfläche des Wertes sehr flach ist, ist delta nahe -1, wie diese Graphik zeigt (konstant im Preis, aber nicht über die Zeit)

Plot3D[
UpAndOutPutDelta[0, S, 50, 50, 0.25,0.1, 0.1, t],
{S, 30, 50}, {t, 0.001, 3}, PlotRange -> All];

Gamma ist Null:

Plot3D[
UpAndOutPutGamma[0, S, 50, 50, 0.25,0.1, 0.1, t],
{S, 30, 50}, {t, 0.001, 3}, PlotRange -> {-0.001, 0.001}];

MathCode C++

Load the package

Needs["MathCode`"]

Set a temporary directory

SetDirectory[".\\Tmp"]

Provide type information to MathCode C++

Declare[
      UpAndOutPutDelta[Real@ reb_, Real@ p_, Real@ k_, Real@ h_, 
          Real@ sd_, Real@ r_, Real@ q_, Real@ t_] -> Real]

We write a function that computes UpAndOutPutDelta over a 45×45 grid (the argument determine how many times to repeat the grid compuation for timing purposes)

\(computeGridValues[Integer@\ n_]\  -> \ Real[45, \ 45]\  := \ 
\ \ 
      Module[{
			Integer@\ {i, \ k, j}, 
			
          Real[45, 45]@\ xyMatrix}, 
		
        For[k = 1, k ≤ n, 
			
          xyMatrix = 
            Table[
					
              UpAndOutPutDelta[0. , 30 + \((50 - 30)\)/45\ i, \ 50, \ 50, \ 
                0.25, 0.1, \ 0.15, \ 0.001 + \((1 - 0.001)\)/45\ j], 

              					{i, 1, 45}, {j, 1, 45}], 
			\(k++\)]; 
		
        xyMatrix]\)

Generate C++ code and an executable file

BuildCode["Global`", EvaluateFunctions -> {UpAndOutPutDelta}]

Compute the time it takes to generate the following plot

Plot3D[
UpAndOutPutDelta[0., S, 50, 50, 0.25,0.1, 0.15, t],
{S, 30, 50}, {t, 0.001, 1}, PlotRange -> All,PlotPoints -> 45];//Timing

inttime = First[%]

Connect the MathCode C++ generated executable file to Mathematica

InstallCode[];

Compute the time for doing the above computation using external, compiled code (we repeat the grid-computation 100 times to get an accurate measure of the time)

\(start = AbsoluteTime[]; 
	computeGridValues[100]; 

    exttime = \((AbsoluteTime[] - start)\)/100  Second\)

The speed-up

inttime/exttime

We define a function that plots the result of the external computation (to check the result)

plotUpAndOutPutDelta[x_?MatrixQ] := 
      Show[SurfaceGraphics[
          Transpose[x], {PlotRange -> All, Axes -> True, 
            MeshRange -> {{30, 50}, {0.001, 1}}}]]

plotUpAndOutPutDelta[computeGridValues[1]]  //  Timing

Looks the same!

UninstallCode[]