別に金融工学が専門なわけではないのですが、ファイナンスは少しやったこともあるので先ほどふと書いてみようと思いついて書いてみました。とりあえずメモ程度に書きます。説明はおいおい書いていくつもりです（12月13日現在）。
d1部分の数式が少し違っていたので修正（12月14日）
コール・オプションのプレミアム価格（ヨーロピアンタイプ）
e ：ネイピア数
q ：原資産利回り(devidend yield)
t ：満期期間(time to maturity)
S ：原資産価格(spot price)
「ちょうちんマーク(x)」みたいなやつ：標準正規分布の累積確率密度関数
r ：安全利子率(interest rate)
K ：行使価格(exercise price)
σ ：ボラティリティ(volatility)
R言語のソースコード
bsmodel <- function(s, k, vol, t, r, q) {
if (t <= 0.0) {
data <-c(0.0, (s-k))
return(max(data))
} else {
return (exp(-q*t)*s*pnorm(d1(s, k, vol, t, r, q))-exp(-r*t)*k*pnorm(d2(s, k, vol, t, r, q)))
}
}
d1 <- function(s, k, vol, t, r, q) {
return ((log(s/k)+(r-q+vol^2/2)*t)/vol*sqrt(t))
}
d2 <- function(s, k, vol, t, r, q) {
return (d1(s, k, vol, t, r,q)-vol*sqrt(t))
}
実行結果（スポットレート ： S＝120 行使価格 ： K＝120 ボラティリティ ： vol＝0.1 期間 ： ｔ＝0.5 安全金利 ： r＝0.01 外国金利 ： ｑ＝0.03）
> bsmodel (120, 120, 0.1, 0.5, 0.01, 0.03) [1] 2.785521
計算結果が他の2サイト（下記参照）で算出したのと少し違っていまして…
とりあえず見返しもせずに作ったというのと、もしかしたら正規分布の累積確率密度関数の計算とはR言語の関数から呼び出したから、結果が少し違うのかもしれません…
修正（12月14日）：四捨五入すれば2.79なので、まあこんなもんでしょうか。
ソースコードはR言語ならそのままコピペすれば使えます。
参考文献
オプション取引入門講座 第9回 ブラック・ショールズ・モデル
ブラック・ショールズ・モデルって何だ？ その１ – 金融工学と経済学とプログラミング
Black Scholes Calculatorを作ってみた – 金融工学と経済学とプログラミング