St<-20
K<-21
t<-0.25
r<-0.01
sigma<-0.2

## First we use the Black-Scholes Formula

d1<-(log(St/K) + (r+(sigma^2)/2)*t)/(sigma*sqrt(t))

d2<-d1-sigma*sqrt(t)

V<-pnorm(d1)*St - K*exp(-r*t)*pnorm(d2)

V

n<-2000
## Now we use a Monte Carlo approach
set.seed(27)

norms<-rnorm(n,0,1)
s = St*exp((r-(sigma^2)/2)*t+sigma*sqrt(t)*norms)

V2<-mean(exp(-r*t)*pmax(s-K,0))

V2

all.equal(V, V2, tolerance=0.01)