Deutsches R-Forum

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)