net.finmath.montecarlo.process.EulerSchemeFromProcessModel Maven / Gradle / Ivy
/*
* (c) Copyright Christian P. Fries, Germany. Contact: [email protected].
*
* Created on 18.11.2012
*/
package net.finmath.montecarlo.process;
import java.util.ArrayList;
import java.util.Map;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import net.finmath.concurrency.FutureWrapper;
import net.finmath.montecarlo.IndependentIncrements;
import net.finmath.montecarlo.model.ProcessModel;
import net.finmath.stochastic.RandomVariable;
/**
* This class implements some numerical schemes for multi-dimensional multi-factor Ito process.
*
*
* It features the standard Euler scheme and the standard predictor-corrector Euler scheme
* for Y, then applies the state space transform \( X = f(Y) \).
*
*
*
* For the standard Euler scheme the process \( Y = (Y_{1},\ldots,Y_{d}) \) is discretized as
* \[
* Y_{j}(t_{i+1}) = Y_{j}(t_{i}) + \mu_{j}(t_{i}) \Delta t_{i} + \lambda_{1,j}(t_{i}) \Delta W_{1}(t_{i}) + \ldots + \lambda_{m,j} \Delta W_{m} \text{.}
* \]
*
*
* The parameters have to be provided by a {@link net.finmath.montecarlo.model.ProcessModel}:
*
* - \( f \) -
applyStateSpaceTransform
* - \( Y(t_{0}) \) -
getInitialState
* - \( \mu \) -
getDrift
* - \( \lambda_{j} \) -
getFactorLoading
*
*
*
* Using the state space transform \( (t,x) \mapsto \exp(x) \), it is possible to create a log-Euler scheme, i.e.,
* \[
* X(t_{i+1}) = X(t_{i}) \cdot \exp\left( (r(t_{i}) - \frac{1}{2} \sigma(t_{i})^2) \Delta t_{i} + \sigma(t_{i}) \Delta W(t_{i}) \right) \text{.}
* \]
* for a process \( \mathrm{d} X = r X \mathrm{d}t + \sigma X \mathrm{d}W \).
*
*
* The dimension \( d \) is called numberOfComponents here.
* The value \( m \) is called numberOfFactors here.
* The default for numberOfFactors is 1.
*
* @author Christian Fries
* @see net.finmath.montecarlo.model.ProcessModel The interface definition contains more details.
* @version 1.4
*/
public class EulerSchemeFromProcessModel extends MonteCarloProcessFromProcessModel {
private static boolean isUseMultiThreadding;
static {
// Default value is true
isUseMultiThreadding = Boolean.parseBoolean(System.getProperty("net.finmath.montecarlo.process.EulerSchemeFromProcessModel.isUseMultiThreadding","true"));
}
public enum Scheme {
EULER,
PREDICTOR_CORRECTOR,
EULER_FUNCTIONAL,
PREDICTOR_CORRECTOR_FUNCTIONAL
}
private final IndependentIncrements stochasticDriver;
private final Scheme scheme;
// Used locally for multi-threadded calculation.
private ExecutorService executor;
/*
* The storage of the simulated stochastic process.
*/
private transient RandomVariable[][] discreteProcess = null;
private transient RandomVariable[] discreteProcessWeights;
/**
* Create an Euler discretization scheme.
*
* @param model The model (the SDE specifcation) used to generate the (sampling of the) stochastic process.
* @param stochasticDriver The stochastic driver of the process (e.g. a Brownian motion).
* @param scheme The scheme to use. See {@link Scheme}.
*/
public EulerSchemeFromProcessModel(final ProcessModel model, final IndependentIncrements stochasticDriver, final Scheme scheme) {
super(stochasticDriver.getTimeDiscretization(), model);
this.stochasticDriver = stochasticDriver;
this.scheme = scheme;
}
/**
* Create an Euler discretization scheme.
*
* @param model The model (the SDE specification) used to generate the (sampling of the) stochastic process.
* @param stochasticDriver The stochastic driver of the process (e.g. a Brownian motion).
*/
public EulerSchemeFromProcessModel(final ProcessModel model, final IndependentIncrements stochasticDriver) {
super(stochasticDriver.getTimeDiscretization(), model);
this.stochasticDriver = stochasticDriver;
Scheme scheme = Scheme.EULER_FUNCTIONAL; // Default, unless applyStateSpaceTransformInverse is not provided
try {
model.applyStateSpaceTransformInverse(null, 0, 0, null);
}
catch(UnsupportedOperationException e) {
// If the inverse is not specified we cannot perform the function EULER
scheme = Scheme.EULER;
}
catch(Exception e) {};
this.scheme = scheme;
}
/**
* This method returns the realization of the process at a certain time index.
*
* @param timeIndex Time index at which the process should be observed
* @return A vector of process realizations (on path)
*/
@Override
public RandomVariable getProcessValue(final int timeIndex, final int componentIndex) {
// Thread safe lazy initialization
synchronized(this) {
if (discreteProcess == null || discreteProcess.length == 0) {
doPrecalculateProcess();
}
}
if(discreteProcess[timeIndex][componentIndex] == null) {
throw new NullPointerException("Generation of process component " + componentIndex + " at time index " + timeIndex + " failed. Likely due to out of memory");
}
// Return value of process
return discreteProcess[timeIndex][componentIndex];
}
/**
* This method returns the weights of a weighted Monte Carlo method (the probability density).
*
* @param timeIndex Time index at which the process should be observed
* @return A vector of positive weights
*/
@Override
public RandomVariable getMonteCarloWeights(final int timeIndex) {
// Thread safe lazy initialization
synchronized(this) {
if (discreteProcessWeights == null || discreteProcessWeights.length == 0) {
doPrecalculateProcess();
}
}
// Return value of process
return discreteProcessWeights[timeIndex];
}
/**
* Calculates the whole (discrete) process.
*/
private void doPrecalculateProcess() {
if (discreteProcess != null && discreteProcess.length != 0) {
return;
}
final int numberOfPaths = this.getNumberOfPaths();
final int numberOfFactors = this.getNumberOfFactors();
final int numberOfComponents = this.getNumberOfComponents();
// Allocate Memory
discreteProcess = new RandomVariable[getTimeDiscretization().getNumberOfTimeSteps() + 1][getNumberOfComponents()];
discreteProcessWeights = new RandomVariable[getTimeDiscretization().getNumberOfTimeSteps() + 1];
// Set initial Monte-Carlo weights
discreteProcessWeights[0] = stochasticDriver.getRandomVariableForConstant(1.0 / numberOfPaths);
// Set initial value
final RandomVariable[] initialState = getInitialState();
final RandomVariable[] currentState = new RandomVariable[numberOfComponents];
for (int componentIndex = 0; componentIndex < numberOfComponents; componentIndex++) {
currentState[componentIndex] = initialState[componentIndex];
discreteProcess[0][componentIndex] = applyStateSpaceTransform(0, componentIndex, currentState[componentIndex]);
}
/*
* Evolve the process using an Euler scheme.
* The evolution is performed multi-threadded.
* Each component of the vector runs in its own thread.
*/
executor = Executors.newCachedThreadPool();
// Evolve process
for (int timeIndex2 = 1; timeIndex2 < getTimeDiscretization().getNumberOfTimeSteps()+1; timeIndex2++) {
final int timeIndex = timeIndex2;
// Generate process from timeIndex-1 to timeIndex
final double deltaT = getTime(timeIndex) - getTime(timeIndex - 1);
// Fetch drift vector
final RandomVariable[] drift;
try {
drift = getDrift(timeIndex - 1, discreteProcess[timeIndex - 1], null);
}
catch(final Exception e) {
throw new RuntimeException(e + " - drift calculaton failed at time index " + timeIndex + " (time=" + getTime(timeIndex - 1) + ") . See cause of this exception for details.", e);
}
// Fetch brownianIncrement vector
final RandomVariable[] brownianIncrement = stochasticDriver.getIncrement(timeIndex - 1);
// Calculate new realization
final ArrayList> discreteProcessAtCurrentTimeIndex = new ArrayList<>(numberOfComponents);
for (int componentIndex2 = 0; componentIndex2 < numberOfComponents; componentIndex2++) {
final int componentIndex = componentIndex2;
final RandomVariable driftOfComponent = drift[componentIndex];
// Check if the component process has stopped to evolve
if (driftOfComponent == null) {
discreteProcessAtCurrentTimeIndex.add(componentIndex, null);
continue;
}
final Callable worker = new Callable() {
@Override
public RandomVariable call() {
if(scheme == Scheme.EULER_FUNCTIONAL || scheme == Scheme.PREDICTOR_CORRECTOR_FUNCTIONAL) {
currentState[componentIndex] = applyStateSpaceTransformInverse(timeIndex - 1, componentIndex, discreteProcess[timeIndex - 1][componentIndex]);
}
final RandomVariable[] factorLoadings = getFactorLoading(timeIndex - 1, componentIndex, discreteProcess[timeIndex - 1]);
// Check if the component process has stopped to evolve
if (factorLoadings == null) {
return null;
}
// Apply drift
if(driftOfComponent != null) {
currentState[componentIndex] = currentState[componentIndex].addProduct(driftOfComponent, deltaT); // mu DeltaT
}
// Apply diffusion
currentState[componentIndex] = currentState[componentIndex].addSumProduct(factorLoadings, brownianIncrement); // sigma DeltaW
// Transform the state space to the value space and return it.
return applyStateSpaceTransform(timeIndex, componentIndex, currentState[componentIndex]);
}
};
/*
* Optional multi-threadding (asyncronous calculation of the components)
*/
Future result = null;
try {
if(isUseMultiThreadding) {
result = executor.submit(worker);
} else {
result = new FutureWrapper<>(worker.call());
}
} catch (final Exception e) {
throw new RuntimeException("Euler step failed at time index " + timeIndex + " (time=" + getTime(timeIndex) + "). See cause of this exception for details.", e);
}
// The following line will add the result of the calculation to the vector discreteProcessAtCurrentTimeIndex
discreteProcessAtCurrentTimeIndex.add(componentIndex, result);
}
// Fetch results and move to discreteProcess[timeIndex]
for (int componentIndex = 0; componentIndex < numberOfComponents; componentIndex++) {
try {
final Future discreteProcessAtCurrentTimeIndexAndComponent = discreteProcessAtCurrentTimeIndex.get(componentIndex);
if(discreteProcessAtCurrentTimeIndexAndComponent != null) {
discreteProcess[timeIndex][componentIndex] = discreteProcessAtCurrentTimeIndexAndComponent.get().cache();
} else {
discreteProcess[timeIndex][componentIndex] = discreteProcess[timeIndex-1][componentIndex];
}
} catch (final InterruptedException | ExecutionException e) {
throw new RuntimeException("Euler step failed at time index " + timeIndex + " (time=" + getTime(timeIndex) + "). See cause of this exception for details.", e.getCause());
}
}
if (scheme == Scheme.PREDICTOR_CORRECTOR || scheme == Scheme.PREDICTOR_CORRECTOR_FUNCTIONAL) {
// Apply corrector step to realizations at next time step
final RandomVariable[] driftWithPredictor = getDrift(timeIndex - 1, discreteProcess[timeIndex], null);
for (int componentIndex = 0; componentIndex < getNumberOfComponents(); componentIndex++) {
final RandomVariable driftWithPredictorOfComponent = driftWithPredictor[componentIndex];
final RandomVariable driftWithoutPredictorOfComponent = drift[componentIndex];
if (driftWithPredictorOfComponent == null || driftWithoutPredictorOfComponent == null) {
continue;
}
// Calculated the predictor corrector drift adjustment
final RandomVariable driftAdjustment = driftWithPredictorOfComponent.sub(driftWithoutPredictorOfComponent).div(2.0).mult(deltaT);
// Add drift adjustment
currentState[componentIndex] = currentState[componentIndex].add(driftAdjustment);
// Re-apply state space transform
discreteProcess[timeIndex][componentIndex] = applyStateSpaceTransform(timeIndex, componentIndex, currentState[componentIndex]);
} // End for(componentIndex)
} // End if(scheme == Scheme.PREDICTOR_CORRECTOR)
// Set Monte-Carlo weights
discreteProcessWeights[timeIndex] = discreteProcessWeights[timeIndex - 1];
} // End for(timeIndex)
try {
executor.shutdown();
}
catch(final SecurityException e) {
// @TODO Improve exception handling here
}
}
/**
* @return Returns the numberOfPaths.
*/
@Override
public int getNumberOfPaths() {
return stochasticDriver.getNumberOfPaths();
}
/**
* @return Returns the numberOfFactors.
*/
@Override
public int getNumberOfFactors() {
return stochasticDriver.getNumberOfFactors();
}
/**
* @return Returns the independent increments interface used in the generation of the process
*/
@Override
public IndependentIncrements getStochasticDriver() {
return stochasticDriver;
}
/**
* @return Returns the scheme.
*/
public Scheme getScheme() {
return scheme;
}
@Override
public EulerSchemeFromProcessModel clone() {
return new EulerSchemeFromProcessModel(getModel(), getStochasticDriver(), scheme);
}
@Override
public MonteCarloProcess getCloneWithModifiedModel(ProcessModel model) {
return new EulerSchemeFromProcessModel(model, getStochasticDriver(), scheme);
}
@Override
public MonteCarloProcess getCloneWithModifiedData(final Map dataModified) {
final ProcessModel newModel = (ProcessModel) dataModified.getOrDefault("model", getModel());
if(dataModified.containsKey("seed") && dataModified.containsKey("stochasticDriver")) {
throw new IllegalArgumentException("Simultaneous specification of stochasticDriver and seed.");
}
final IndependentIncrements newStochasticDriver;
if(dataModified.containsKey("seed")) {
newStochasticDriver = getStochasticDriver().getCloneWithModifiedSeed((int)dataModified.get("seed"));
}
else if(dataModified.containsKey("stochasticDriver")) {
newStochasticDriver = (IndependentIncrements) dataModified.getOrDefault("stochasticDriver", stochasticDriver);
}
else {
newStochasticDriver = stochasticDriver;
}
final Scheme newScheme = (Scheme) dataModified.getOrDefault("scheme", scheme);
return new EulerSchemeFromProcessModel(newModel, newStochasticDriver, newScheme);
}
@Override
public Object getCloneWithModifiedSeed(final int seed) {
return new EulerSchemeFromProcessModel(getModel(), getStochasticDriver().getCloneWithModifiedSeed(seed));
}
@Override
public String toString() {
return "EulerSchemeFromProcessModel [stochasticDriver=" + stochasticDriver + ", scheme=" + scheme + ", executor="
+ executor + "]";
}
}