# R program to compute and also estimate the fair price
#     of an Asian Call Option

M = 1000
k = 8

# stock parameters:
a = 100
r = 0.03
T = 10
q = 105
sigma = 0.1

# Monte Carlo estimate:
h = 0.05
num = round(T/h)
vallist = rep(0,M)
for (j in 1:M) {
  X = rep(0,num+1)
  XiTk = rep(0,k)
  X[1] = a
  for (i in 1:num)
    X[i+1] = rnorm(1, mean = X[i] + r*X[i]*h, sd = sigma*X[i]*sqrt(h))
  for (i in 1:k)
    XiTk[i] = X[round(i*num/k)+1]
  Xbark = mean(XiTk)
  vallist[j] = exp(-r*T) * max(0, Xbark-q)
}
estmean = mean(vallist)
se = sd(vallist) / sqrt(M)

# Output results:

cat("Asian option, a=", a, " r=", r, " sigma=", sigma,
				" T=", T, " q=", q, " k=", k, "\n", sep="")
cat("Monte Carlo estimated price:", estmean, "+-", se, "\n")
cat("95% confidence interval: (", estmean-1.96*se, ",", estmean+1.96*se, ")\n")

# Plot the final simulation:
thetimes = (0 : (T/h))
plot(thetimes, X, type='l')