#------------------------------#
#                              #
#    Black-Scholes Formeln:    #
#    exakt vs. Monte Carlo     #
#                              #
#------------------------------#


# Input Parameter:

S0 = 100
K  = 100 
T  = 2
sigma = 20/100
r  = 3/100

N = 100000


# Exaktes Resultat, die Black-Scholes Formeln: 

dp = ( log(S0/K) + (r+sigma^2/2)*T ) / ( sigma*sqrt(T) )
dm = ( log(S0/K) + (r-sigma^2/2)*T ) / ( sigma*sqrt(T) )

bscall =  S0 * pnorm(dp)  - K*exp(-r*T) * pnorm(dm) 
bsput  = -S0 * pnorm(-dp) + K*exp(-r*T) * pnorm(-dm) 

bscall
bsput


# Mit Monte Carlo Simulation:

x = rnorm(N)
ST = S0 * exp( sigma*sqrt(T)*x + (r-sigma^2/2)*T ) 
Hcall = ifelse( ST > K , ST-K , 0 ) 
Hput  = ifelse( ST < K , K-ST , 0 )

Vcall = exp(-r*T) * sum(Hcall) / N
Vput  = exp(-r*T) * sum(Hput)  / N

Vcall
Vput


# Schauen wir uns an, wie schnell das Monte Carlo konvergiert: 

vcall = exp(-r*T) * cumsum(Hcall) / (1:N)
vput  = exp(-r*T) * cumsum(Hput)  / (1:N)

plot(vcall)

plot(vcall, ylim = c( 0.9*bscall, 1.1*bscall ) )
points( rep(bscall,N), col="red" )

plot(vput, ylim = c( 0.9*bsput, 1.1*bsput ) )
points( rep(bsput,N), col="red" )




#-------------------#
#    Beispiel 2     #
#-------------------#


# 2b) mit exponential-verteilten Zufallszahlen: 

N = 100000
x = rexp(N, rate = 1/2) 
F = 2*ifelse( x<log(4) , 1 , 0 )
res = sum(F)/N
res

# schauen wir uns die Konvergenz an:

res2 = cumsum(F)/(1:N)
plot(res2, ylim=c(0.9,1.1) ) 
points( rep(1,N), col="red" )


# 2c) mit gleichverteilten Zufallszahlen: 

N = 100000
x = runif(N, min=0, max=log(4) ) 
F = log(4) * exp(-x/2)
res = sum(F)/N
res

# schauen wir uns die Konvergenz an:

res2 = cumsum(F)/(1:N)
plot(res2, ylim=c(0.9,1.1) ) 
points( rep(1,N), col="red" )