380 likes | 685 Views
The Black-Scholes Model & Basic Numerical Techniques. Chapter 2: ADVANCED OPTION PRICING MODEL. The Concepts Underlying Black-Scholes. The option price and the stock price depend on the same underlying source of uncertainty
E N D
The Black-Scholes Model & Basic Numerical Techniques Chapter 2: ADVANCED OPTION PRICING MODEL
The Concepts Underlying Black-Scholes • The option price and the stock price depend on the same underlying source of uncertainty • We can form a portfolio consisting of the stock and the option which eliminates this source of uncertainty • The portfolio is instantaneously riskless and must instantaneously earn the risk-free rate • This leads to the Black-Scholes differential equation
Risk-Neutral Valuation • The variable does not appear in the Black-Scholes equation • The equation is independent of all variables affected by risk preference • The solution to the differential equation is therefore the same in a risk-free world as it is in the real world • This leads to the principle of risk-neutral valuation
Applying Risk-Neutral Valuation • Assume that the expected return from the stock price is the risk-free rate • Calculate the expected payoff from the option • Discount at the risk-free rate
The Black-Scholes Formulas where N(x) is cumulative probability distribution function for a standardized normal distribution. It is the probability that a variable with a standard normal distribution will be less than x. = annualized standard deviation (volatility) of the continuously compounded return on the stock r = continuously compounded risk-free rate q = continuously compounded dividend yield
Matlab Program – Black-Scholes Formulas • function V = bsf ( S, E, r, q, sigma, t, T, CallPut ) %% BSF evaluates the Black-Scholes formula %% %% CallPut=1 for call -1 for put %% tau = T - t; if ( 0 < tau ) d1 = ( log ( S / E ) + ( r – q + 0.5 * sigma^2 ) * tau ) … / ( sigma * sqrt ( tau ) ); d2 = d1 - sigma * sqrt ( tau ); N1 = 0.5 * ( 1.0 + erf ( CallPut * d1 / sqrt ( 2.0 ) ) ); N2 = 0.5 * ( 1.0 + erf ( CallPut * d2 / sqrt ( 2.0 ) ) ); V = CallPut*(S * exp(-q*tau) * N1 - E * exp ( - r * tau ) * N2); else V = max ( CallPut*(S – E), 0.0 ); end
The Binomial Model S0u ƒu p S0 ƒ S0d ƒd (1– p ) • f=e-rT[pfu+(1-p)fd ]
The Binomial Model • The asset is assumed to only be allowed to either take an "up" step or a "down" step, where these steps are given as: for an up movement, and: for a down movement. • As the asset can only take an up or a down move, and not both simultaneously, the down movement can be simplified to:
The Binomial Model -- Cont’d • In a risk-neutral world the stock price grows at r-qrather than at rwhen there is a dividend yield at rate q • The probability, p, of an up movement must therefore satisfy pS0u+(1-p)S0d=S0e(r-q)T so that
Pricing options in CRR Binomial -- Program • function Binomial = CRR(CallPut, AssetP, Strike, RiskFree, Div, Time, Vol, nSteps) % Cox, Ross & Rubinstein (1979) Binomial Tree for European Call/Put Option Values based on % the following inputs: % CallPut = Call = 1, Put = -1; AssetP = Underlying Asset Price; Strike = Strike Price of Option % RiskFree = Risk Free rate of interest annualized eg. 0.05; Div = Dividend Yield of Underlying % Time = Time to Maturity in years; Vol = Volatility of the Underlying % nSteps = Number of Time Steps for Binomial Tree to take dt = Time / nSteps; RR = exp(RiskFree * dt); Up = exp(Vol * sqrt(dt)); Down = 1 / Up; Q_up = (exp((RiskFree - Div) * dt) - Down) / (Up - Down); Q_down = 1 - Q_up; Df = exp(-RiskFree * dt); %Df: Discount Factor
Pricing options in CRR Binomial – Program (Cont’d) %Populate all possible stock prices and option values on the end notes of the tree for i = 0:nSteps state = i + 1; St = AssetP * Up ^ i * Down ^ (nSteps - i); Value(state) = max(0, CallPut * (St - Strike)); end %Since value on the end nodes are known by above, % we start from nSteps-1 working backwards % double for loop: outer-every steps; inner-every nodes on each step for k = nSteps - 1 : -1 : 0 for i = 0:k state = i + 1; Value(state) = (Q_up * Value(state + 1) + Q_down * Value(state)) * Df; end end Binomial = Value(1);
Example: CRR - result • >> CRR(1,100,105,0.05,0,2,0.4,100) Binomial = 24.3440 • >> CRR(-1,100,105,0.05,0,2,0.4,100) Binomial = 19.3520
Finite Differences (FD) • As with any class of option, the price of the derivative is governed by solving the underlying partial differential equation. The use of finite difference methods allows us to solve these PDEs by means of an iterative procedure. • We can start by looking at the Black-Scholes partial differential equation: where dV is the change in the value of an option, dt is a small change in time. is the volatility of the underlying, S is the underlying price and is the carry (r-q). • By specifying initial and boundary conditions, one can attain numerical solutions to all the derivatives of the Black-Scholes PDE using a finite difference grid. The grid is typically set up so that partitions in two dimensions - space and time (in our case, we would be looking at the asset price and the change in time): • Once the grid is set up, there are three methods to evaluate the PDE at each time step. The difference between each of the three methods is contingent on the choice of difference used for time (i.e. forward, backward or central differences). Central differences is used for the space grid (S).
Explicit Finite Differences • Explicit FD uses forward differences at each time node t. By splitting the differential equation into the time element and space elements, we can apply forward differences to time as follows: if we substitute x = ln(S), the equation becomes: • Applying the finite differences method, the above equation can be broken down and approximated: becomes • For the space grid, we can apply central differences for all order of derivatives: becomes and becomes and becomes
Explicit Finite Differences • Combining the terms gives: • Which is the same as: where the probabilities of each of the nodes is: • This case is actually equivalent to the trinomial tree where probabilities can be assigned to the likelihood of an up move, a down move as well as no move. It can also be shown that the following approximation holds:
Explicit Finite-Difference • For stability and convergence, • As we can see, the result from this method is explicitly given because we know the value (the claim) at the boundary where the option expires. Then, we perform the calculation backwards in time until the valuation date. Since the time dependence (i) only depends on future dates (i+1) we can explicitly calculate the change, node by node backward in time.
Implicit Finite Differences • The implicit method takes backward differences for the time derivative but still using central differences for the space derivatives. • Although similar in nature to the explicit finite differences method, the implicit FD method is typically more stable and convergent than the explicit FD method - however, it is often more computationally intensive. • The approximation to the PDE under an implicit FD method is given by the following: Note that the main difference between the above equation and the one for the explicit FD method is in the selection of time step i. The subsequent simplification of the approximation and the associated probabilities is similar to that of the explicit FD method.
Crank-Nicholson Scheme • An improvement over the implicit FD method is the Crank-Nicolson Scheme which uses central differences for both time and space dimensions. The result is that over smaller time steps dt, the method is more accurate, stable and convergent than both implicit and explicit methods - however, like the implicit FD method (which requires evaluating equations at each time step) it is more computationally intensive than the explicit FD method. • The approximation for the PDE is given as:
Monte Carlo Simulations • In the field of mathematical finance, many problems, for instance the problem of finding the arbitrage-free value of a particular derivative, boil down to the computation of a particular integral. • When the number of dimensions (or degrees of freedom) in the problem is large, PDE's and numerical integrals become intractable, and in these cases Monte Carlo methods often give better results. For large dimensional integrals, Monte Carlo methods converge to the solution more quickly than numerical integration methods, require less memory and are easier to program.
Monte Carlo Simulations • if we suppose that our risk-neutral probability space is and that we have a derivative H that depends on a set of underlying instruments S1,...,Sn. Then given a sample ω from the probability space the value of the derivative is H(S1(ω),S2(ω),...,Sn(ω)) = :H(ω). Today's value of the derivative is found by taking the expectation over all possible samples and discounting at the risk-free rate. I.e. the derivative has value: where dfT is the discount factor corresponding to the risk-free rate to the final maturity date T years into the future. • Now suppose the integral is hard to compute. We can approximate the integral by generating sample paths and then taking an average. Suppose we generate N samples then
Monte Carlo Simulations • For example in the standard Black-Scholes model, the stock price evolves as dS = μ(t)Sdt + σ(t)SdWt. • To sample a path following this distribution from time 0 to T, we chop the time interval into M units of length δt, and approximate the Brownian motion over the interval dt by a single normal variable of mean 0 and variance δt. This leads to a sample path of Here each εis a draw from a standard normal distribution. • We obtain the Monte-Carlo value of this derivative by generating N lots of normal variables, creating N sample paths and so N values of H, and then taking the derivative. • The overall average Pmean is our Monte Carlo estimate of the option value. The computed variance, Pvar, may be used to give an approximate 95% confidence interval
Matlab Program -- Monte Carlo Simulations • function [Pmean, width] = MCbsf(S, K, r, q, sigma, T, N, CallPut) % MC Monte Carlo valuation for a European call % %%%%%%%%%% Problem and method parameters %%%%%%%%% %% S: Current price; K: Strike price, r: interest rate, sigma: volatility %% %% CallPut =1 for Call and -1 for Put %% %% T: Time to maturity, N: Number of sample paths %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% conf = [Pmean - width, Pmean + width] %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Svals = S*exp((r-q-0.5*sigma^2)*T + sigma*sqrt(T)*randn(N,1)); Pvals = exp(-r*T)*max(CallPut*(Svals-K),0); Pmean = mean(Pvals) width = 1.96*std(Pvals)/sqrt(N);
Newton’s Root Search Method • Solutions to nonlinear equations f(x)=0 • Can be found using Newton’s Method • Example: find the solution of f(x) = e-x – sin(x) • Given x0, iterate maxit times or until |xn-xn-1|tol.
Newton’s Method -- Program • function [f,f_prime] = myfun(x) % MYFUN– Evaluate f(x) = exp(x)-sin(x) % and its first derivative % [f,f_prime] = myfun(x) f=exp(-x)-sin(x); f_prime=-exp(-x)-cos(x); return • >> help feval • >> [f,f_prime]=feval(’myfun’,0);
Newton’s Method – Program (Cont’d) • function [x,n] = newtonf(fname,x0,tol,maxit) % NEWTON – Newton’s method for solving equations % [x,n] = NEWTON(fname,x0,tol,maxit) x = x0; n = 0; done=0; while ~done, n = n + 1; [f,f_prime]=feval(fname,x); x_new = x - f/f_prime; done=(n>maxit) | ( abs(x_new-x)<tol ); x=x_new; end • >> [x,n]=newtonf(’myfun’,0,1e-3,10) x = 0.588 n = 4
Monte Carlo Simulation: Option Pricing • We are going to price two types of option: • spread option: file MCSpreadOption.m. • This option can be priced in closed form solution but it is interesting to • see how we can simulate the paths of two stocks that are correlated. • Asian option: file Asian.m
Spread Option • Spread Option: A call option on the spread between two assets S1 and S2 has payoff with K a strike price. We can simulate the paths of the two stocks • for i = 0, 1, . . . , n − 1. Let assume that the two processes are correlated and ρ is the instantaneous correlation between Z1 and Z2.
Spread Option -- Codes • function [Pmean, width] = MCSpreadOption(S1, S2, K, r, v1, v2, rho, T, CallPut, nSimulations) Drift1 = (r - v1 ^ 2 / 2) * T; % Drift asset one Drift2 = (r - v2 ^ 2 / 2) * T; % Drift asset two v1Sqrdt = v1 * sqrt(T);v2Sqrdt = v2 * sqrt(T);sum=0; Epsilon1 = randn(1,nSimulations); % Random vector from a Stand.Norm.Distr. Epsilon2 = rho*Epsilon1+randn(1,nSimulations) * sqrt(1-rho^ 2); St1 = S1 * exp(Drift1 + v1Sqrdt * Epsilon1); St2 = S2 * exp(Drift2 + v2Sqrdt * Epsilon2); Pvals = exp(-r*T)*max(CallPut*(St1 - St2 – K),0); Pmean = mean(Pvals); width = 1.96*std(Pvals)/sqrt(nSimulations); • where the Levy representation of Geometric Brownian motions is used: CallPut =1 for calls and -1 for puts.
Spread Options -- Result • Let evaluate the option assuming that = 110, = 100, K = 1, r = .06, σ1 = .2, σ2 = .3, ρ = .5 and the number of simulations is equal to 10, 000. We assume that the evaluation date, t, is July 14th 2005 and the terminal date, T, is July 14th 2006. • >> MCSpreadOption Valuation Date Day = 14 Month = 07 Year = 2005 Terminal Date Day = 14 Month = 07 Year = 2006 • MCSpreadOption = 16.2346 Elapsed time is 1.336092 seconds.
Asian Options • Asian Option: It may be necessary to simulate paths of the underlying asset if the payoff of a derivative security depends on the value of the underlying asset at intermediate dates and not just the terminal value. Asian option is the simplest path-dependent option and the payoff is where for some dates t0< t1< ... < tm= T with T the date at which the payoff is received.
Asian Option -- Codes • function [Pmean, width] = Asian(S, K, r, q, v, T, nn, nSimulations, CallPut) dt = T/nn; Drift = (r - q - v ^ 2 / 2) * dt; vSqrdt = v * sqrt(dt); pathSt = zeros(nSimulations,nn); Epsilon = randn(nSimulations,nn); St = S*ones(nSimulations,1); % for each time step for j = 1:nn; St = St .* exp(Drift + vSqrdt * Epsilon(:,j)); pathSt(:,j)=St; end SS = cumsum(pathSt,2); Pvals = exp(-r*T) * max(CallPut * (SS(:,nn)/nn - K), 0); % Pvals dimension: nSimulations x 1 Pmean = mean(Pvals); width = 1.96*std(Pvals)/sqrt(nSimulations);
Asian Option -- Result • Let evaluate the option assuming that St= 100, K = 100, r = .06, σ = .2 and the number of simulations is equal to 5,000. We calculate the mean on daily basis and we assume that the evaluation date, t, is July 14th 2005, the terminal date, T, is July 14th 2006 and we count every day observation to calculate • >> Asian Valuation Date Day = 14 Month = 07 Year = 2005 Terminal Date Day = 14 Month = 07 Year = 2006 • Elapsed time is 115.923847 seconds. price = 6.1268