I want to compute the derivative of Black and Scholes option

Question

Good evening to everybody, my question is about computing the derivative of a Black and Scholes option and to draw it.

My code is

f\[ScriptCapitalD][μ_, σ_, t_] = LogNormalDistribution[μ t, σ Sqrt[t]];

PriceAtEpoch[S0_, μ_, σ_, t_] = S0 Mean[f\[ScriptCapitalD][μ, σ, t]]

PriceAtEpoch[Subscript[S, 0], μ, σ, t + 1] == 
 PriceAtEpoch[Subscript[S, 0], μ, σ, t] Exp[r]

{sol} = Solve[ForAll[{t, Subscript[S, 0]}, %], μ, Reals]

BlackScholesOptionPrice = 
 Assuming[Subscript[S, 0] > 0 && r > 0 && σ > 0 && t > 0 && \[ScriptCapitalK] > 0, 
  Exp[-r t] Expectation[optVal[\[ScriptF] Subscript[S, 0], 
      \[ScriptCapitalK]], \[ScriptF] \[Distributed] 
      f\[ScriptCapitalD][μ /. sol, σ, t]] // FullSimplify]

---

So how to compute the derivative (for example with respect to the volatility) for the call option?

---

Thank you for your feedback. I have tried to copy what is written in the book but it didn't work (see below the code)

{In[2]:= f = 1/(s_t b (2 pi)^.5) Exp [ - ((Log (s_t) - a^2)/(2 b^2))]


In[9]:= ((0.7071067811865475` E^(-((-a^2 + Log s_t)/(
  2 b^2))))/(b pi^0.5` s_t))

domain [f] = {s_t, 0, ∞} && {a ∈ Reals, b > 0}
 

Out[9]= (0.707107 E^(-((-a^2 + Log s_t)/(2 b^2))))/(b pi^0.5 s_t)

During evaluation of In[9]:= Rule::rhs: Pattern s_t appears on the right-hand side of rule domain[(0.707107 E^(-((-Power[<<2>>]+<<1>>)/(2 b^2))))/(b pi^0.5 s_t)]->{s_t,0,∞}&&{a∈Reals,b>0}. >>

Out[10]= {s_t, 0, ∞} && {a ∈ Reals, b > 0}

In[11]:= V_t = If[s_t > k, s_t - k, 0]

Out[11]= If[s_t > k, s_t - k, 0]

In[12]:= V = \[DoubleStruckE]^(-rT) Expect [V_t, f]

Out[12]= \[DoubleStruckE]^-rT Expect[V_t, (
  0.707107 E^(-((-a^2 + Log s_t)/(2 b^2))))/(b pi^0.5 s_t)]

In[15]:= Value = 
 V /. {a → Log [p] + (r - (σ^2)/2) T, 
   b → σ (T)^.5}

During evaluation of In[15]:= ReplaceAll::reps: {a→T (r-σ^2/2)+Log[p],b→T^0.5 σ} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

During evaluation of In[15]:= Set::wrsym: Symbol Value is Protected. >>

Out[15]= \[DoubleStruckE]^-rT Expect[V_t, (
   0.707107 E^(-((-a^2 + Log s_t)/(2 b^2))))/(
   b pi^0.5 s_t)] /. {a → T (r - σ^2/2) + Log[p], 
  b → T^0.5 σ}

In[16]:= Value /. {p → 104, k → 100, 
  r → .05, σ → .44, 
  T → 66/365}

During evaluation of In[16]:= ReplaceAll::reps: {p→104,k→100,r→0.05,σ→0.44,T→66/365} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

Out[16]= Value /. {p → 104, k → 100, 
  r → 0.05, σ → 0.44, 
  T → 66/365}

My question is still the same: I want to plot the first derivative of a call/put option.

Edit

Once again as, because of my inexperience, I have troubles with the code.

I want to compute the derivative of a Black and Scholes option and to plot it. After some iteration here is the code:

---

In[10]:= s = LogNormalDistribution[μ, σ]

Out[10]= LogNormalDistribution[μ, σ]

In[11]:= f = 
 1/(s * b (2 pi)^(1/2)) Expectation[ - ((Log (s) - a^2)/(2 b^2))]


During evaluation of In[11]:= Expectation::argmu: Expectation called with 1 argument; 2 or more arguments are expected. >>

Out[11]= Expectation[-((-a^2 + 
   Log LogNormalDistribution[μ, σ])/(
  2 b^2))]/(Sqrt[2] b Sqrt[pi] LogNormalDistribution[μ, σ])

In[12]:= ((0.7071067811865475` E^(-((-a^2 + Log s_t)/(
  2 b^2))))/(b pi^0.5` s_t))

domain [f] = {s, 0, ∞} && {a ∈ Reals, b > 0}
 

Out[12]= (0.707107 E^(-((-a^2 + Log s_t)/(2 b^2))))/(b pi^0.5 s_t)

Out[13]= {LogNormalDistribution[μ, σ], 
  0, ∞} && {a ∈ Reals, b > 0}

In[14]:= V_t = If[s_t > k, s_t - k, 0]

Out[14]= If[s_t > k, s_t - k, 0]

In[15]:= V = \[DoubleStruckE]^(-rT) Expect [V, f]

During evaluation of In[15]:= $RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of \[DoubleStruckE]^-rT. >>

Out[15]= Hold[V = \[DoubleStruckE]^-rT Expect[V, f]]

In[16]:= Value = 
 V /. {a → Log [p] + (r - (σ^2)/2) T, 
   b → σ (T)^.5}

During evaluation of In[16]:= $RecursionLimit::reclim2: Recursion depth of 
  1024 exceeded during evaluation of \[DoubleStruckE]^-rT. >>

Out[16]= Hold[
 Value = V /. {a → Log[p] + (r - σ^2/2) T, 
    b → σ T^0.5}]

In[17]:= Value /. {p → 104, k → 100, 
  r → .05, σ → .44, 
  T → 66/365}

During evaluation of In[17]:= ReplaceAll::reps: {p→104,k→100,r→0.05,σ→0.44,T→66/365} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

Out[17]= Value /. {p → 104, k → 100, 
  r → 0.05, σ → 0.44, 
  T → 66/365}

---

As you see it doesn't work. Is anybody able to tell me what I should do?

Answer 1

There are some obvious problems with your formulation .. the main one being that you define BlackScholesOptionPrice in terms of a function called optVal[] ... and you have not defined the latter. The second major issue is that you are just taking the expectation of a Lognormal random variable, whereas an option is the expectation of a censored Lognormal rv (censored below: leaving the upside gain).

There is an example in our Springer book, "Mathematical Statistics with Mathematica" at pp.70-71, setting up a Black-Scholes model in Mathematica from first principles, using the method you want. You can easily modify the example to use inbuilt functions such as Expectation that are now available in current versions, which should also work fine. You can download a free copy of the original printed book as a PDF file (or just Chapter 2) over here:

http://www.mathstatica.com/book/bookcontents.html

At the end of the example, it is also noted how to take the derivative that you seek.

---

Here is a worked example

Stock pries $S_T$ follow a Brownian motion. As per equation 2.38 above, the distribution model is:

$$\log S_T\sim N\left(a,b^2\right) \quad \quad \text{where} \quad \begin{align} a &= log{S(0)} + \left(r-\frac{\sigma ^2}{2}\right) T \\ b &=\sigma \sqrt{T} \\ \end{align} $$

To keep notation simple, let $S = S_T$.

The value of the option at expiry VT, with a strike of $k$, is:

VT = If[S > k, S-k, 0];

The value of a call option TODAY is $V0 = e^{-r t} E[\text{VT}]$:

V0 = Exp[-r T] Expectation[VT,Distributed[S,LogNormalDistribution[a,b]], Assumptions->k>0]

Finally, substitute in for $a$ and $b$, and we have:

 optionValue = V0 /. {a->Log[p] + (r-\[Sigma]^2/2) T, b->\[Sigma] Sqrt[T]};

All done.

The derivative wrt volatility is:

 D[optionValue, \[Sigma]] // Simplify