Java tutorial
/** * Copyright (C) 2011 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.financial.model.finitedifference.applications; import org.apache.commons.lang.Validate; import com.opengamma.analytics.financial.model.finitedifference.BoundaryCondition; import com.opengamma.analytics.financial.model.finitedifference.ConvectionDiffusionPDE1DCoupledCoefficients; import com.opengamma.analytics.financial.model.finitedifference.CoupledFiniteDifference; import com.opengamma.analytics.financial.model.finitedifference.CoupledPDEDataBundle; import com.opengamma.analytics.financial.model.finitedifference.DirichletBoundaryCondition; import com.opengamma.analytics.financial.model.finitedifference.NeumannBoundaryCondition; import com.opengamma.analytics.financial.model.finitedifference.PDEFullResults1D; import com.opengamma.analytics.financial.model.finitedifference.PDEGrid1D; import com.opengamma.analytics.financial.model.finitedifference.PDEResults1D; import com.opengamma.analytics.financial.model.interestrate.curve.ForwardCurve; import com.opengamma.analytics.math.function.Function; import com.opengamma.analytics.math.function.Function1D; import com.opengamma.analytics.math.statistics.distribution.NormalDistribution; import com.opengamma.analytics.math.surface.FunctionalDoublesSurface; /** * Solves a coupled forward PDE (i.e. coupled Fokker-Plank) for the density of an asset when the process is CEV with vol levels determined by a * two state Markov chain. The densities, p(t,s,state1) & p(t,s,state2), are such that int_{0}^{\infty} p(t,s,stateX) ds gives the probability * of being in state X at time t, and (p(t,s,state1)+p(t,s,state2))*ds is the probability that the asset with be between s and s + ds at time t. */ public class TwoStateMarkovChainDensity { private static final double THETA = 1.0; private final ConvectionDiffusionPDE1DCoupledCoefficients _data1; private final ConvectionDiffusionPDE1DCoupledCoefficients _data2; private final Function1D<Double, Double> _initCon11; private final Function1D<Double, Double> _initCon12; public TwoStateMarkovChainDensity(final ForwardCurve forward, final double vol1, final double deltaVol, final double lambda12, final double lambda21, final double probS1, final double beta1, final double beta2) { this(forward, new TwoStateMarkovChainDataBundle(vol1, vol1 + deltaVol, lambda12, lambda21, probS1, beta1, beta2)); } public TwoStateMarkovChainDensity(final ForwardCurve forward, final TwoStateMarkovChainDataBundle data) { Validate.notNull(forward, "null forward"); Validate.notNull(data, "null data"); _data1 = getCoupledPDEDataBundle(forward, data.getVol1(), data.getLambda12(), data.getLambda21(), data.getBeta1()); _data2 = getCoupledPDEDataBundle(forward, data.getVol2(), data.getLambda21(), data.getLambda12(), data.getBeta2()); _initCon11 = getInitialCondition(forward, data.getP0()); _initCon12 = getInitialCondition(forward, 1.0 - data.getP0()); } PDEFullResults1D[] solve(final PDEGrid1D grid) { //BoundaryCondition lower = new FixedSecondDerivativeBoundaryCondition(0, grid.getSpaceNode(0), true); final BoundaryCondition lower = new NeumannBoundaryCondition(0.0, grid.getSpaceNode(0), true); //BoundaryCondition lower = new DirichletBoundaryCondition(0.0, grid.getSpaceNode(0));//TODO for beta < 0.5 zero is accessible and thus there will be non-zero //density there final BoundaryCondition upper = new DirichletBoundaryCondition(0.0, grid.getSpaceNode(grid.getNumSpaceNodes() - 1)); CoupledPDEDataBundle d1 = new CoupledPDEDataBundle(_data1, _initCon11, lower, upper, grid); CoupledPDEDataBundle d2 = new CoupledPDEDataBundle(_data2, _initCon12, lower, upper, grid); final CoupledFiniteDifference solver = new CoupledFiniteDifference(THETA, true); final PDEResults1D[] res = solver.solve(d1, d2); //handle this with generics final PDEFullResults1D res1 = (PDEFullResults1D) res[0]; final PDEFullResults1D res2 = (PDEFullResults1D) res[1]; return new PDEFullResults1D[] { res1, res2 }; } private Function1D<Double, Double> getInitialCondition(final ForwardCurve forward, final double initialProb) { //using a log-normal distribution with a very small Standard deviation as a proxy for a Dirac delta return new Function1D<Double, Double>() { private final double _volRootTOffset = 0.01; @Override public Double evaluate(final Double s) { if (s <= 0 || initialProb == 0) { return 0.0; } final double x = Math.log(s / forward.getSpot()); final NormalDistribution dist = new NormalDistribution(0, _volRootTOffset); return initialProb * dist.getPDF(x) / s; } }; } private ConvectionDiffusionPDE1DCoupledCoefficients getCoupledPDEDataBundle(final ForwardCurve forward, final double vol, final double lambda1, final double lambda2, final double beta) { final Function<Double, Double> a = new Function<Double, Double>() { @Override public Double evaluate(final Double... ts) { Validate.isTrue(ts.length == 2); double s = ts[1]; if (s <= 0.0) { //TODO review how to handle absorption s = -s; } return -Math.pow(s, 2 * beta) * vol * vol / 2; } }; final Function<Double, Double> b = new Function<Double, Double>() { @Override public Double evaluate(final Double... ts) { Validate.isTrue(ts.length == 2); final double t = ts[0]; double s = ts[1]; if (s < 0.0) { s = -s; } final double temp = (s < 0.0 ? 0.0 : 2 * vol * vol * beta * Math.pow(s, 2 * (beta - 1))); return s * (forward.getDrift(t) - temp); } }; final Function<Double, Double> c = new Function<Double, Double>() { @Override public Double evaluate(final Double... ts) { Validate.isTrue(ts.length == 2); final double t = ts[0]; double s = ts[1]; if (s < 0.) { s = -s; } double temp = (beta == 1.0 ? 1.0 : Math.pow(s, 2 * (beta - 1))); if (s < 0) { temp = 0.0; } return lambda1 + forward.getDrift(t) - vol * vol * beta * (2 * beta - 1) * temp; } }; return new ConvectionDiffusionPDE1DCoupledCoefficients(FunctionalDoublesSurface.from(a), FunctionalDoublesSurface.from(b), FunctionalDoublesSurface.from(c), -lambda2); } }