Assmannhäuse Höllebesch!OK, kann auch Rotwein (ungekühlt) sein
Beim "Umschreiben" des Programms habe ich mich stark an Deine Vorlage gehalten, da müsste man meine Ergänzungen noch etwas glätten.
Auch müsste man - wenn die exakte Reproduzierbarkeit für Dich wichtig ist - noch etwas investieren.
Bitte beachten: das läuft so auf Linux (Ubuntu) - was bei Windows (wg. parallel) passiert, kann ich nicht sagen...
Code: Alles auswählen
library(Rcpp)
library(parallel)
#-------------------------------------------------------------------------------
cppFunction('double Simulation(double vStart, double SStart, double dt, double r,
double q, double theta, double kappa, double sigma,
NumericVector randv, NumericVector randS) {
int n = randv.size();
double v = vStart;
double S = SStart;
double vneu = 0;
double Sneu = 0;
double dt_vt = 0;
for(int i = 0; i < n; ++i) {
dt_vt = kappa*(theta - v)*dt + sigma*sqrt(v*dt)*randv[i];
vneu = v + dt_vt;
if (vneu <= 0)
vneu = vneu - 2*dt_vt;
Sneu = S*exp((r - q - v/2)*dt + sqrt(v*dt)*randS[i]);
v = vneu;
S = Sneu;
}
return S;
}')
#-------------------------------------------------------------------------------
einWert <- function(lnr, S, Tm, r, q, v0, theta, rho, kappa, sigma){
NSteps <- floor(Tm*252)
dt <- Tm/NSteps
rand.v <- rnorm(NSteps)
rand.h <- rnorm(NSteps)
sqrho <- sqrt(1 - rho^2)
rand.S <- rho * rand.v + sqrho * rand.h
Wert <- Simulation(v0, S, dt, r, q, theta, kappa, sigma, rand.v, rand.S)
return(Wert)
}
#-------------------------------------------------------------------------------
HestonSV_MC <- function(S, Tm, r, q, v0, theta, rho, kappa, sigma, NPaths, seed){
set.seed(seed)
# Hier mc.cores auf die verfügbaren Threads ändern!
Results <- mclapply(1:NPaths, FUN="einWert", S=S, Tm=Tm, r=r, q=q, v0=v0, theta=theta,
rho=rho, kappa=kappa, sigma=sigma, mc.cores = 40)
return(unlist(Results))
}
#-------------------------------------------------------------------------------
r <- 0.01
Tm <- 1
St <- HestonSV_MC(S=100, Tm=1, r=0.01, q=0.02, v0=0.04, theta=0.25, rho=-0.5,
kappa=4.0, sigma=1.0, NPaths=1000000, seed=606060)
Results <- data.frame(Strike=seq(80,120,10))
Results$PriceCall <- NA
Results$PricePut <- NA
for(i in 1:nrow(Results)){
Results$PriceCall[i] <- exp(-r*Tm)*mean(pmax(St - Results$Strike[i], 0))
Results$PricePut[i] <- exp(-r*Tm)*mean(pmax(Results$Strike[i] - St, 0))
}
print(Results)