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.

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