You're welcome.

Since you're a Mathematica novice, here's a simple version to get you started.

Examples:

cblack[100, 105, 0.05, 0.40, 0.25]

5.83223

Timing[cblackavg[100, 105, 0.05, 0.40, 0.25]]

{0., 0.839571}

As advertised, the timing is essentially instantaneous.

With slightly more complex code, you can find cblackavg with an arbitrary number of accurate digits, although the timing will increase under high precision. You can also display more digits in the answer under Machine precision, which is the precision of my code here.

=====================================================

Finally, you can ask Mathematica to

try to do it symbolically (analytically). Trying

Integrate[cblack[F, F, 0, sig, t], {t, 0, T}]

will return the analytic ATM case you've already done in terms of Erf[..].

But with Integrate[cblack[F, K, 0, sig, t], {t, 0, T}] or even

Integrate[Erf[Sqrt[t] + a/Sqrt[t]], {t, 0, T}], Mathematica trys, but gives up after a few minutes, returning essentially the unevaluated integral. The fail on that last one suggests a general analytic expression (in terms of special functions with no infinite sums or integrations) may be very difficult to find, if at all.

```
Cumnormal[x_] := (1 + Erf[x/Sqrt[2]])/2
CN = Cumnormal;
d1[S_, K_, r_, sig_, T_] :=
(Log[S/K] + (r + sig^2/2) T)/(sig Sqrt[T])
d2[S_, K_, r_, sig_, T_] :=
(Log[S/K] + (r - sig^2/2) T)/(sig Sqrt[T])
cblack[F_, K_, r_, sig_, T_] := E^(-r T)
(F CN[d1[F, K, 0, sig, T]] - K CN[d2[F, K, 0, sig, T]])
cblackavg[F_, K_, r_, sig_, T_] :=
NIntegrate[cblack[F, K, r, sig, t], {t, 0, T}]
```