simulate

Consider a European up-and-in call option on a single underlying stock. The evolution of this stock's price is governed by a Geometric Brownian Motion (GBM) model with constant parameters:

Assume the following characteristics:

The goal is to simulate various paths of daily stock prices, and calculate the price of the barrier option as the risk-neutral sample average of the discounted terminal option payoff. Since this is a barrier option, you must also determine if and when the barrier is crossed.

This example performs antithetic sampling by explicitly setting theAntitheticflag totrue, and then specifies an end-of-period processing function to record the maximum and terminal stock prices on a path-by-path basis.

Create a GBM model usinggbm.

barrier = 120;% barrierstrike = 100;% exercise pricerate = 0.05;%年无风险的老鼠esigma = 0.3;% annualized volatilitynPeriods = 63;% 63 trading daysdt = 1 / 252;% time increment = 252 daysTime = nPeriods * dt;% expiration time = 0.25 yearsobj = gbm(rate, sigma,'StartState', 105);

Perform a small-scale simulation that explicitly returns two simulated paths.

rng('default')% make output reproducible[X, T] = obj.simBySolution(nPeriods,'DeltaTime', dt,...'nTrials', 2,'Antithetic', true);

Perform antithetic sampling such that all primary and antithetic paths are simulated and stored in successive matching pairs. Odd paths (1,3,5,...) correspond to the primary Gaussian paths. Even paths (2,4,6,...) are the matching antithetic paths of each pair, derived by negating the Gaussian draws of the corresponding primary (odd) path. Verify this by examining the matching paths of the primary/antithetic pair.

plot(T, X(:,:,1),'blue', T, X(:,:,2),'red') xlabel('Time (Years)'), ylabel('Stock Price'),...title('Antithetic Sampling') legend({'Primary Path''Antithetic Path'},...'Location','Best')

To price the European barrier option, specify an end-of-period processing function to record the maximum and terminal stock prices. This processing function is accessible by time and state, and is implemented as a nested function with access to shared information that allows the option price and corresponding standard error to be calculated. For more information on using an end-of-period processing function, seePricing Equity Options.

Simulate 200 paths using the processing function method.

rng('default')% make output reproduciblebarrier = 120;% barrierstrike = 100;% exercise pricerate = 0.05;%年无风险的老鼠esigma = 0.3;% annualized volatilitynPeriods = 63;% 63 trading daysdt = 1 / 252;% time increment = 252 daysTime = nPeriods * dt;% expiration time = 0.25 yearsobj = gbm(rate, sigma,'StartState', 105); nPaths = 200;% # of paths = 100 sets of pairsf = Example_BarrierOption(nPeriods, nPaths); simulate(obj, nPeriods,'DeltaTime', dt,...'nTrials', nPaths,'Antithetic', true,...'Processes', f.SaveMaxLast);

Approximate the option price with a 95% confidence interval.

optionPrice = f.OptionPrice(strike, rate, barrier); standardError = f.StandardError(strike, rate, barrier,...true); lowerBound = optionPrice - 1.96 * standardError; upperBound = optionPrice + 1.96 * standardError; displaySummary(optionPrice, standardError, lowerBound, upperBound);

Up-and-In Barrier Option Price: 6.6572 Standard Error of Price: 0.7292 Confidence Interval Lower Bound: 5.2280 Confidence Interval Upper Bound: 8.0864

functiondisplaySummary(optionPrice, standardError, lowerBound, upperBound) fprintf(' Up-and-In Barrier Option Price: %8.4f\n',...optionPrice); fprintf(' Standard Error of Price: %8.4f\n',...standardError); fprintf(' Confidence Interval Lower Bound: %8.4f\n',...lowerBound); fprintf(' Confidence Interval Upper Bound: %8.4f\n',...upperBound);end

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (SDEMRD): $d{X}_t=S(t)[L(t)-X_t]dt+D(t,X_t^{\frac{1}{2}})V(t)dW$

where $$D$$is a diagonal matrix whose elements are the square root of the corresponding element of the state vector.

Create acirobject to represent the model: $d{X}_t=0.2(0.1-X_t)dt+0.05X_t^{\frac{1}{2}}dW$.

cir_obj = cir(0.2, 0.1, 0.05)% (Speed, Level, Sigma)

cir_obj = Class CIR: Cox-Ingersoll-Ross ---------------------------------------- Dimensions: State = 1, Brownian = 1 ---------------------------------------- StartTime: 0 StartState: 1 Correlation: 1 Drift: drift rate function F(t,X(t)) Diffusion: diffusion rate function G(t,X(t)) Simulation: simulation method/function simByEuler Sigma: 0.05 Level: 0.1 Speed: 0.2

Use the optional name-value argument for'Scheme'to specify a scheme to simulate the sample paths.

[paths,times] = simulate(cir_obj,10,'ntrials',4096,'scheme',“quadratic-exponential”);

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (SDEMRD): $d{X}_t=S(t)[L(t)-X_t]dt+D(t,X_t^{\frac{1}{2}})V(t)dW$

where $$D$$is a diagonal matrix whose elements are the square root of the corresponding element of the state vector.

Create acirobject to represent the model: $d{X}_t=0.2(0.1-X_t)dt+0.05X_t^{\frac{1}{2}}dW$.

cir_obj = cir(0.2, 0.1, 0.05)% (Speed, Level, Sigma)

cir_obj = Class CIR: Cox-Ingersoll-Ross ---------------------------------------- Dimensions: State = 1, Brownian = 1 ---------------------------------------- StartTime: 0 StartState: 1 Correlation: 1 Drift: drift rate function F(t,X(t)) Diffusion: diffusion rate function G(t,X(t)) Simulation: simulation method/function simByEuler Sigma: 0.05 Level: 0.1 Speed: 0.2

Use optional name-value inputs for thesimByEulermethod that you can call throughsimulateinterface using thesimulate'Scheme'for“quadratic-exponential”. The optional inputs forsimByEulerdefine a quasi-Monte Carlo simulation using the name-value arguments for'MonteCarloMethod','QuasiSequence', and'BrownianMotionMethod'.

[paths,time] = simulate(cir_obj,10,'ntrials',4096,'montecarlomethod','quasi','quasisequence','sobol','BrownianMotionMethod','brownian-bridge','scheme',“quadratic-exponential”);