Java tutorial
/** * Copyright (C) 2012 - 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.ConvectionDiffusionPDE1DCoefficients; import com.opengamma.analytics.financial.model.finitedifference.ConvectionDiffusionPDE1DStandardCoefficients; import com.opengamma.analytics.financial.model.finitedifference.DirichletBoundaryCondition; import com.opengamma.analytics.financial.model.finitedifference.ExponentialMeshing; import com.opengamma.analytics.financial.model.finitedifference.HyperbolicMeshing; import com.opengamma.analytics.financial.model.finitedifference.MeshingFunction; import com.opengamma.analytics.financial.model.finitedifference.PDE1DDataBundle; import com.opengamma.analytics.financial.model.finitedifference.PDEGrid1D; import com.opengamma.analytics.financial.model.finitedifference.PDEResults1D; import com.opengamma.analytics.financial.model.finitedifference.ThetaMethodFiniteDifference; import com.opengamma.analytics.financial.model.option.definition.Barrier; import com.opengamma.analytics.financial.model.option.definition.Barrier.BarrierType; import com.opengamma.analytics.financial.model.option.definition.Barrier.KnockType; import com.opengamma.analytics.financial.model.option.pricing.analytic.formula.BlackBarrierPriceFunction; import com.opengamma.analytics.financial.model.option.pricing.analytic.formula.EuropeanVanillaOption; import com.opengamma.analytics.financial.model.volatility.smile.fitting.interpolation.SurfaceArrayUtils; import com.opengamma.analytics.math.function.Function1D; /** * */ public class BarrierOptionPricer { private static final InitialConditionsProvider ICP = new InitialConditionsProvider(); private static final PDE1DCoefficientsProvider PDE = new PDE1DCoefficientsProvider(); private static final ThetaMethodFiniteDifference SOLVER = new ThetaMethodFiniteDifference(); private static final int DEFAULT_XNODES = 100; private static final int DEFAULT_TNODES = 50; private static final double DEFAULT_LAMBDA = 0.0; private static final double DEFAULT_BUNCHING = 1.0; private static final double DEFAULT_Z = 2.0; private final int _nTNodes; private final int _nXNodes; private final double _lambda; private final double _bunching; private final double _z = DEFAULT_Z; public BarrierOptionPricer() { _nTNodes = DEFAULT_TNODES; _nXNodes = DEFAULT_XNODES; _lambda = DEFAULT_LAMBDA; _bunching = DEFAULT_BUNCHING; } public BarrierOptionPricer(final int numXNodes, final int numTNodes, final double lambda, final double bunching) { _nXNodes = numXNodes; _nTNodes = numTNodes; _lambda = lambda; _bunching = bunching; } /** * Computes the price of a barrier option in the Black world. * @param option The underlying European vanilla option. * @param barrier The barrier. * @param rebate The rebate. This is paid <b>immediately</b> if the knock-out barrier is hit and at expiry if the knock-in barrier is not hit * @param spot The spot price. * @param costOfCarry The cost of carry (i.e. the forward = spot*exp(costOfCarry*T) ) * @param rate The interest rate. * @param sigma The Black volatility. * @return The price. */ public double getPrice(final EuropeanVanillaOption option, final Barrier barrier, final double rebate, final double spot, final double costOfCarry, final double rate, final double sigma) { Validate.notNull(option, "option"); Validate.notNull(barrier, "barrier"); final boolean isKnockIn = (barrier.getKnockType() == KnockType.IN); final boolean isDown = (barrier.getBarrierType() == BarrierType.DOWN); //in these pathological cases the barrier is hit immediately so the value is just from the rebate (knock-out) or a European option (knock-in) if (isDown && spot <= barrier.getBarrierLevel() || !isDown && spot >= barrier.getBarrierLevel()) { if (isKnockIn) { return blackPrice(spot, option.getStrike(), option.getTimeToExpiry(), rate, costOfCarry, sigma, option.isCall()); } return rebate; } if (isKnockIn) { return inBarrier(spot, barrier.getBarrierLevel(), option.getStrike(), option.getTimeToExpiry(), rate, costOfCarry, sigma, option.isCall(), rebate); } return outBarrier(spot, barrier.getBarrierLevel(), option.getStrike(), option.getTimeToExpiry(), rate, costOfCarry, sigma, option.isCall(), rebate); } /** * Computes the price of a one-touch out barrier option in the Black-Scholes world by solving the BS PDE on a finite difference grid. If a barrier is hit at any time before expiry, * the option is cancelled (knocked-out) and a rebate (which is often zero) is paid <b>immediately</b>. If the barrier is not hit, then a normal European option payment is made. <p> * If the barrier is above the spot it is assumed to be an up-and-out barrier (otherwise it would expire immediately) otherwise it is a down-and-out barrier * As there are exact formulae for this case (see {@link BlackBarrierPriceFunction}), this is purely for testing purposes. * @param spot The current (i.e. option price time) of the underlying * @param barrierLevel The barrier level * @param strike The strike of the European option * @param expiry The expiry of the option * @param rate The interest rate. * @param carry The cost-of-carry (i.e. the forward = spot*exp(costOfCarry*T) ) * @param vol The Black volatility. * @param isCall true for call * @param rebate The rebate amount. * @return The price. */ public double outBarrier(final double spot, final double barrierLevel, final double strike, final double expiry, final double rate, final double carry, final double vol, final boolean isCall, final double rebate) { final Function1D<Double, Double> intCon = ICP.getEuropeanPayoff(strike, isCall); final ConvectionDiffusionPDE1DStandardCoefficients pde = PDE.getBlackScholes(rate, rate - carry, vol); final boolean isUp = barrierLevel > spot; final double adj = 0.0; // _lambda == 0 ? ZETA * vol * Math.sqrt(expiry / (_nTNodes - 1)) : 0.0; double sMin; double sMax; BoundaryCondition lower; BoundaryCondition upper; if (isUp) { sMin = 0.0; sMax = barrierLevel * Math.exp(-adj); //bring the barrier DOWN slightly to adjust for discrete monitoring if (isCall) { lower = new DirichletBoundaryCondition(0.0, sMin); } else { final Function1D<Double, Double> lowerValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return Math.exp(-rate * tau) * strike; } }; lower = new DirichletBoundaryCondition(lowerValue, sMin); } upper = new DirichletBoundaryCondition(rebate, sMax); } else { sMin = barrierLevel * Math.exp(adj); //bring the barrier UP slightly to adjust for discrete monitoring sMax = spot * Math.exp(_z * Math.sqrt(expiry)); lower = new DirichletBoundaryCondition(rebate, sMin); if (isCall) { final Function1D<Double, Double> upperValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return Math.exp(-rate * tau) * (spot * Math.exp(carry * tau) - strike); } }; upper = new DirichletBoundaryCondition(upperValue, sMax); } else { upper = new DirichletBoundaryCondition(0.0, sMax); } } final MeshingFunction tMesh = new ExponentialMeshing(0, expiry, _nTNodes, _lambda); final MeshingFunction xMesh = new HyperbolicMeshing(sMin, sMax, spot, _nXNodes, _bunching); final PDEGrid1D grid = new PDEGrid1D(tMesh, xMesh); final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> pdeData = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>( pde, intCon, lower, upper, grid); final PDEResults1D res = SOLVER.solve(pdeData); //for now just do linear interpolation on price. TODO replace this with something more robust final double[] xNodes = grid.getSpaceNodes(); final int index = SurfaceArrayUtils.getLowerBoundIndex(xNodes, spot); final double w = (xNodes[index + 1] - spot) / (xNodes[index + 1] - xNodes[index]); return w * res.getFunctionValue(index) + (1 - w) * res.getFunctionValue(index + 1); } /** * Computes the price of a one-touch in barrier option in the Black-Scholes world by solving the BS PDE on a finite difference grid. If a barrier is hit at any time before expiry, * the option becomes a simple European (call or put). If the barrier is not hit a rebate is paid at the option expiry<p> * If the barrier is above the spot it is assumed to be an up-and-in barrier otherwise it is a down-and-in barrier * As there are exact formulae for this case (see {@link BlackBarrierPriceFunction}), this is purely for testing purposes. * @param spot The current (i.e. option price time) of the underlying * @param barrierLevel The barrier level * @param strike The strike of the European option * @param expiry The expiry of the option * @param rate The interest rate. * @param carry The cost-of-carry (i.e. the forward = spot*exp(costOfCarry*T) ) * @param vol The Black volatility. * @param isCall true for call * @param rebate The rebate amount. * @return The price. */ public double inBarrier(final double spot, final double barrierLevel, final double strike, final double expiry, final double rate, final double carry, final double vol, final boolean isCall, final double rebate) { final double outPrice = outBarrierSpecial(spot, barrierLevel, strike, expiry, rate, carry, vol, isCall, rebate); final double bsPrice = blackPrice(spot, strike, expiry, rate, carry, vol, isCall); return bsPrice + Math.exp(-rate * expiry) * rebate - outPrice; } /** * Computes the price of a one-touch out barrier option in the Black-Scholes world assuming the rebate is paid at the option expiry. This is NOT the case for a out barrier, * but is for an in, and as we must price an in barrier and a European option plus a bond (the rebate) minus an out, we need this special case. * @param spot The current (i.e. option price time) of the underlying * @param barrierLevel The barrier level * @param strike The strike of the European option * @param expiry The expiry of the option * @param rate The interest rate. * @param carry The cost-of-carry (i.e. the forward = spot*exp(costOfCarry*T) ) * @param vol The Black volatility. * @param isCall true for call * @param rebate The rebate amount. * @return The price. */ protected double outBarrierSpecial(final double spot, final double barrierLevel, final double strike, final double expiry, final double rate, final double carry, final double vol, final boolean isCall, final double rebate) { final Function1D<Double, Double> intCon = ICP.getEuropeanPayoff(strike, isCall); final ConvectionDiffusionPDE1DStandardCoefficients pde = PDE.getBlackScholes(rate, rate - carry, vol); final boolean isUp = barrierLevel > spot; final Function1D<Double, Double> rebateValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return rebate * Math.exp(-rate * tau); } }; double sMin; double sMax; BoundaryCondition lower; BoundaryCondition upper; if (isUp) { sMin = 0.0; sMax = barrierLevel; if (isCall) { lower = new DirichletBoundaryCondition(0.0, sMin); } else { final Function1D<Double, Double> lowerValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return Math.exp(-rate * tau) * strike; } }; lower = new DirichletBoundaryCondition(lowerValue, sMin); } upper = new DirichletBoundaryCondition(rebateValue, sMax); } else { sMin = barrierLevel; sMax = spot * Math.exp(_z * Math.sqrt(expiry)); lower = new DirichletBoundaryCondition(rebateValue, sMin); if (isCall) { final Function1D<Double, Double> upperValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return Math.exp(-rate * tau) * (spot * Math.exp(carry * tau) - strike); } }; upper = new DirichletBoundaryCondition(upperValue, sMax); } else { upper = new DirichletBoundaryCondition(0.0, sMax); } } final MeshingFunction tMesh = new ExponentialMeshing(0, expiry, _nTNodes, _lambda); final MeshingFunction xMesh = new HyperbolicMeshing(sMin, sMax, spot, _nXNodes, _bunching); final PDEGrid1D grid = new PDEGrid1D(tMesh, xMesh); final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> pdeData = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>( pde, intCon, lower, upper, grid); final PDEResults1D res = SOLVER.solve(pdeData); //for now just do linear interpolation on price. TODO replace this with something more robust final double[] xNodes = grid.getSpaceNodes(); final int index = SurfaceArrayUtils.getLowerBoundIndex(xNodes, spot); final double w = (xNodes[index + 1] - spot) / (xNodes[index + 1] - xNodes[index]); return w * res.getFunctionValue(index) + (1 - w) * res.getFunctionValue(index + 1); } protected double blackPrice(final double spot, final double strike, final double expiry, final double rate, final double carry, final double vol, final boolean isCall) { final Function1D<Double, Double> intCon = ICP.getEuropeanPayoff(strike, isCall); final ConvectionDiffusionPDE1DStandardCoefficients pde = PDE.getBlackScholes(rate, rate - carry, vol); final double sMin = 0.0; final double sMax = spot * Math.exp(_z * Math.sqrt(expiry)); BoundaryCondition lower; BoundaryCondition upper; if (isCall) { lower = new DirichletBoundaryCondition(0.0, sMin); final Function1D<Double, Double> upperValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return Math.exp(-rate * tau) * (spot * Math.exp(carry * tau) - strike); } }; upper = new DirichletBoundaryCondition(upperValue, sMax); } else { final Function1D<Double, Double> lowerValue = new Function1D<Double, Double>() { @Override public Double evaluate(final Double tau) { return Math.exp(-rate * tau) * strike; } }; lower = new DirichletBoundaryCondition(lowerValue, sMin); upper = new DirichletBoundaryCondition(0.0, sMax); } final MeshingFunction tMesh = new ExponentialMeshing(0, expiry, _nTNodes, _lambda); final MeshingFunction xMesh = new HyperbolicMeshing(sMin, sMax, spot, _nXNodes, _bunching); final PDEGrid1D grid = new PDEGrid1D(tMesh, xMesh); final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> pdeData = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>( pde, intCon, lower, upper, grid); final PDEResults1D res = SOLVER.solve(pdeData); //for now just do linear interpolation on price. TODO replace this with something more robust final double[] xNodes = grid.getSpaceNodes(); final int index = SurfaceArrayUtils.getLowerBoundIndex(xNodes, spot); final double w = (xNodes[index + 1] - spot) / (xNodes[index + 1] - xNodes[index]); return w * res.getFunctionValue(index) + (1 - w) * res.getFunctionValue(index + 1); } }