# Ornstein-Uhlenbeck process set.seed(123) d <- expression(-5 * x) s <- expression(3.5) sde.sim(X0=10,drift=d, sigma=s) -> X plot(X,main=”Ornstein-Uhlenbeck”)  

## (European) Option Surfaces OptionSurface <- function(EOres, label, fov=60) { axis.col <- “black” text.col <- axis.col ylab <- label xlab <- “Underlying” zlab <- “Volatility” y <- EOres ## clear scene: clear3d() clear3d(type=”bbox”) clear3d(type=”lights”) ## setup env: ## bg3d(color=”#887777″) bg3d(color=”#DDDDDD”) light3d() rgl.viewpoint(fov=fov) ##rgl.bg(col=”white”, fogtype=”exp2″) ##rgl.bg(col=”black”, fogtype=”exp2″) ##rgl.bg(col=”black”, fogtype=”exp”) ##rgl.bg(col=”white”, fogtype=”exp”) x <- (1:nrow(y)) z <- […]

## Option pricing Functions EuropeanOption <- function(type, underlying, strike, dividendYield, riskFreeRate, maturity, volatility) { UseMethod(“EuropeanOption”) } EuropeanOption.default <- function(type, underlying, strike, dividendYield, riskFreeRate, maturity, volatility) { type <- match.arg(type, c(“call”, “put”)) val <- europeanOptionEngine(type, underlying, strike, dividendYield, riskFreeRate, maturity, volatility) class(val) <- c(“EuropeanOption”, “Option”) val } AmericanOption <- function(type, underlying, strike, dividendYield, riskFreeRate, maturity, volatility, […]

## How to perform a Monte Carlo Simulation? ## First Step: # Write a function to generate the option’s innovations. # Use scrambled normal Sobol numbers: sobolInnovations = function(mcSteps, pathLength, init, …) { # Create Normal Sobol Innovations: innovations = rnorm.sobol(mcSteps, pathLength, init, …) # Return Value: innovations } ## Second Step: # Write a […]

# Monte Carlo Simulation / Simulation.r # basic parameters unit_cost <- 1.50 # per unit variable cost unit_price <- 4.00 # unit price per item unit_disp <- 0.20 # cost for the disposal of unsold inventories # simulate 3000 profits for 20K, 40K & 60K units each profit <- rep(NA, 3000) set.seed(1) for (i in […]

# R Monte Carlo Simulation Price Path Converging Volatility #number of simulation runs nSims = 1000 S = 100 r = 0 q = 0 volatility = 0.25 #drift term mu = r – q #every trading day increment dt = 1 / 365 vol_vec = vector() for (tDays in 1:365) { # standard normal […]

# MCS for single and multi assets asset.paths <- function(s0, mu, sigma, nsims = 10000, periods = c(0, 1) # time periods at which to simulate prices ) { s0 = as.vector(s0) nsteps = len(periods) dt = c(periods[1], diff(periods)) if( len(s0) == 1 ) { drift = mu – 0.5 * sigma^2 if( nsteps == […]

## model – # Define the Model Parameters for a Heston-Nandi Option: model = list(lambda = -0.5, omega = 2.3e-6, alpha = 2.9e-6, beta = 0.85, gamma = 184.25) S = X = 100 Time.inDays = 252 r.daily = 0.05/Time.inDays sigma.daily = sqrt((model$omega + model$alpha) / (1 – model$beta – model$alpha * model$gamma^2)) data.frame(S, X, […]

## model – # Define the Model Parameters for a Heston-Nandi Option: model = list(lambda = -0.5, omega = 2.3e-6, alpha = 2.9e-6, beta = 0.85, gamma = 184.25) S = X = 100 Time.inDays = 252 r.daily = 0.05/Time.inDays sigma.daily = sqrt((model$omega + model$alpha) / (1 – model$beta – model$alpha * model$gamma^2)) data.frame(S, X, […]

## hngarchSim – # Simulate a Heston Nandi Garch(1,1) Process: # Symmetric Model – Parameters: model = list(lambda = 4, omega = 8e-5, alpha = 6e-5, beta = 0.7, gamma = 0, rf = 0) ts = hngarchSim(model = model, n = 500, n.start = 100) par(mfrow = c(2, 1), cex = 0.75) ts.plot(ts, col […]