Option pricing in R

# see more: http://www.reakkt.com/2012/03/diverging-on-option-prices.html

library(fOptions)
library(RQuantLib)

# S - spot
# K - strike price
# rf - risk free rate
# TTM - time to maturity
# sigma - volatility
BSCall <- function(S,K,rf,TTM,sigma) {

  d1 <- (log(S/K)+(rf+sigma^2/2)*TTM) / (sigma*sqrt(TTM))

  d2 <- d1-sigma*sqrt(TTM)

  pnorm(d1)*S-pnorm(d2)*K*exp(-rf*TTM)

}

BSPut <- function(S,K,rf,TTM,sigma) {

  K*exp(-rf*TTM)-S+BSCall(S,K,rf,TTM,sigma)

}

# q - dividend yield
BSCall.dividend <- function(S,K,rf,q,TTM,sigma) {

  F = S * exp((rf-q)*TTM)

  d1 <- (log(F/K)+(sigma^2/2)*TTM) / (sigma*sqrt(TTM))

  d2 <- d1-sigma*sqrt(TTM)

  exp(-rf*TTM)*( (F*pnorm(d1))-K*pnorm(d2) )

}

BSPut.dividend <-  function(S,K,rf,q,TTM,sigma) {

  F = S * exp((rf-q)*TTM)

  d1 <- (log(F/K)+(sigma^2/2)*TTM) / (sigma*sqrt(TTM))

  d2 <- d1-sigma*sqrt(TTM)

  exp(-rf*TTM)*( (K*pnorm(-d2))-F*pnorm(-d1) )

}

BSCall(S=5,K=7,rf=0.045,TTM=0.5,sigma=0.27)

BSPut(S=10,K=7,rf=0.045,TTM=0.5,sigma=0.27)

BSCall.dividend(S=5,K=7,rf=0.045,q=0.02,TTM=0.5,sigma=0.27)

BSPut.dividend(S=10,K=7,rf=0.045,q=0.02,TTM=0.5,sigma=0.27)

###

xSpot   = 601042.9378
xStrike = 254417.3841
xRf     = 0.14395813
xDiv    = 0.11951086
xTTM    = 6.22885299
xSigma  = 1.04225573

BSPut.dividend(S=xSpot,K=xStrike,rf=xRf,q=xDiv,TTM=xTTM,sigma=xSigma)

GBSOption(TypeFlag="p",S=xSpot,X=xStrike,Time=xTTM,r=xDiv,b=xRf,sigma=xSigma)@price

EuropeanOption(type="put",underlying=xSpot,strike=xStrike,dividendYield=xDiv,riskFreeRate=xRf,maturity=xTTM,volatility=xSigma)$value