Asian option numerical pricing method generates a negative time value

I use R to write a function which simulates price path and calculates the value of an arithmetic Asian option. I found sometimes the value of the option can be lower than its intrinsic value, i.e., negative time value. I am wondering if this is a feature of arithmetic Asian option, or is there anything wrong with my code.

###############################################################################
OptionMC = function (style, type, S, X, r, q, tDays, volatility, nSims) {
#   drift term
mu = r - q
#   underlying price, S can be a vector with previous prices
S0 = tail(S, n=1)
dt = 1 / 252
#   time to maturity express in fraction of years
tMaturity = tDays / 252
#   standard normal distribution random number
z = rnorm(tDays*nSims, mean=0, sd=1)

#   generate log-normal return matrix
return_matrix = matrix(exp((mu - 0.5 * volatility ^ 2) * dt + volatility * sqrt(dt) * z), ncol=nSims)
#   return value: price path matrix
path = rbind(matrix(rep(S,nSims),ncol=nSims), S0*apply(return_matrix,2,cumprod))
#   calculate option delta, generate a new price path, with the same brownian motion term "z"
path_delta = rbind(matrix(rep(S,nSims),ncol=nSims), S0*1.0001*apply(return_matrix,2,cumprod))
path_delta[length(S),] = path_delta[length(S),] * 1.0001

#   different option style
if (style == "vanilla") {
#       plain vanilla option, use the final price as settlement price
ref = tail(path, n=1)
ref_delta = tail(path_delta, n=1)
} else if (style == "arithmetic") {
#       arithmetic mean of the entire price path, fixed strike price
ref = apply(path, 2, mean)
ref_delta = apply(path_delta, 2, mean)
} else if (style == "geometric") {
#       geometric mean of the entire price path, fixed strike price
ref = apply(path, 2, function(x){exp(mean(log(x)))})
ref_delta = apply(path_delta, 2, function(x){exp(mean(log(x)))})
} else {
}

#   option's value, average of difference between terminal value and strike price, then discount back to present value
value = exp(- r * tMaturity) * sum(pmax(ifelse(type=="call",1,-1) * (ref - X), 0)) / nSims

#   option's value with "path_delta", i.e., values calculated for calculating delta
value_delta = exp(- r * tMaturity) * sum(pmax(ifelse(type=="call",1,-1) * (ref_delta - X), 0)) / nSims
#   delta is defined as the change of option value w.r.t. change of underlying price (S0)
delta = (value_delta - value) / (S0 * 0.0001)

#   return value:
list(path=path, value=value, delta=delta)
}

• Does it get negative even for a large nSims? – James Sep 15 '14 at 15:01
• @James Yes, it does. I have tried 100,000 runs, but it still gets a negative time value. Thanks. – 2607 Sep 18 '14 at 5:21