Option Pricing with MATLAB

read

When considering some types of option, there sometimes exists a closed form solution which, under the Black-Scholes assumptions, delivers the ‘fair’ price of the option with respect to the various input parameters. For some of the more advanced option contracts, the closed form solution may not be so easy to obtain and thus a numerical method may be required. These numerical methods will be discussed in later parts, but first we present the closed form solution for one of the simplest types of option – the European vanilla call.

European call option

The Black-Scholes formula for the option price is given by

$C(S, \tau) = S e^{-q\tau} N(d_1) - Xe^{-r\tau}N(d_2),$ $d_1 = \frac{\ln(S/X)+(r-q+\sigma^2/2)\tau}{\sigma\sqrt{\tau}},$ $d_2 = d_1 - \sigma \sqrt{\tau},$

where the option parameters are

$N(.)$ : the cumulative distribution function of the standard normal distribution
$\tau=T-t$ : the time to maturity
$S$ : the current spot price of the underlying asset
$X$ : the strike price
$r$ : the risk free rate (annualized, continuously compounded)
$q$ : continuous dividend yield (annualized, continuously compounded)
$\sigma$ : the volatility of the underlying asset’s returns.

Below is the MATLAB implementation

% Black-Scholes formula for European vanilla call
% call syntax: c = BS_EurCall( S0, X, r, T, sigma, q)

function value = BS_EurCall( S0, X, r, T, sigma, q)

d1 = (log(S0/X) + (r - q + sigma^2/2)*T)/(sigma*sqrt(T));
d2 = d1 - sigma*sqrt(T);

value = S0.*exp(-q*T).*normcdf(d1) - X*exp(-r*T).*normcdf(d2);

end

However, assuming we didn’t know how to solve the Black-Scholes equation for the closed-form formula, we could use numerical methods such as binomial tree model or Monte Carlo simulation to obtain the option value. Below is the MATLAB implementation for each of these methods.

% Binomial tree model (BTM) program for European vanilla call options
% call syntax: c = btm_EurCall(S0,X,r,T,sigma,q,N)

function value = btm_EurCall(S0,X,r,T,sigma,q,N)
% set up lattice parameters
dt = T/N; dx = sigma*sqrt(dt);
u = exp(dx); d = 1/u;
df = exp(-r*dt);     % discount factor
p = (exp((r-q)*dt)-d)/(u-d);  % risk-neutral probability
% initialization
j = 1:1:N+1;  % range of index for price states
V = max(S0*u.^(2*j-N-2)-X,0);
% backward recursive through time
for n = N-1:-1:0   
j = 1:1:n+1;
V = df*(p*V(j+1)+(1-p)*V(j));
end
value = V(1);

end

% Monte-Carlo simulation program for European vanilla call options
% call syntax: c = MC_EurCall(S0,X,r,T,sigma,q,no_samples)

function value = MC_EurCall(S0,X,r,T,sigma,q,no_samples)
mu = r-q-sigma^2/2;
epsv = randn(no_samples,1);  % random standard normal numbers
ST = S0*exp(mu*T+epsv*sigma*sqrt(T));  % terminal prices in a vector
value = exp(-r*T)*mean(max(ST-X,0)); % discounted expectation

European put option

Given the call option formula, we can use the put-call parity to derive the price of the put option having the same underlying asset and strike price, which is given by

$\begin{align} P(S, \tau) &= C(S, \tau) + Xe^{-r\tau} - S \\ &= S[N(d_1)-1]+Xe^{-r\tau}[1-N(d_2)]\\ &= Xe^{-r\tau}N(-d_2) - Se^{-q\tau}N(-d_1) \end{align}$

Below is the MATLAB implementation

% Black-Scholes formula for European vanilla put
% call syntax: c = BS_EurPut( S0, X, r, T, sigma, q)

function value = BS_EurPut( S0, X, r, T, sigma, q)

d1 = (log(S0/X) + (r - q + sigma^2/2)*T)/(sigma*sqrt(T));
d2 = d1 - sigma*sqrt(T);

value = X*exp(-r*T).*normcdf(-d2) - S0.*exp(-q*T).*normcdf(-d1);

end

Similarly, the value of a put option can also be obtained using binomial tree method or Monte Carlo simulation. The codes are pretty much the same as in the case of call option, with a few minor changes.

% Binomial tree model (BTM) program for European vanilla call options
% call syntax: c = btm_EurPut(S0,X,r,T,sigma,q,N)

function value = btm_EurPut(S0,X,r,T,sigma,q,N)
% set up lattice parameters
dt = T/N; dx = sigma*sqrt(dt);
u = exp(dx); d = 1/u;
df = exp(-r*dt);     % discount factor
p = (exp((r-q)*dt)-d)/(u-d);  % risk-neutral probability
% initialization
j = 1:1:N+1;  % range of index for price states
V = max(X-S0*u.^(2*j-N-2),0);
% backward recursive through time
for n = N-1:-1:0   
j = 1:1:n+1;
V = df*(p*V(j+1)+(1-p)*V(j));
end
value = V(1);

end

% Monte-Carlo simulation program for European vanilla put options
% call syntax: c = MC_EurPut(S0,X,r,T,sigma,q,no_samples)

function value = MC_EurPut(S0,X,r,T,sigma,q,no_samples)
mu = r-q-sigma^2/2;
epsv = randn(no_samples,1);  % random standard normal numbers
ST = S0*exp(mu*T+epsv*sigma*sqrt(T));  % terminal prices in a vector
value = exp(-r*T)*mean(max(X-ST,0)); % discounted expectation

American call and put options

An American option can be exercised at any time, whereas a European option can only be exercised at the expiration date. This added flexibility of American options increases their value over European options in certain situations. Thus, we can say American Options = European Options + Premium where the Premium is greater than or equal to zero.

However it is suboptimal to exercise an American call option on a non-dividend-paying stock before the expiration date. This is because for a given movement in the stock price, the profit from holding an in the money call is equivalent to the profit from holding the stock. The call option, however, has the added benefit of protecting against the risk of a downward price movement below the strike price. Additionally, because of the time value of money, it costs more to exercise the option today at a fixed strike price K than in the future at K. Finally, there is an intrinsic time value of the option that would be lost by exercising the option prior to the expiration date. Hence, the price of an American and European call option without dividends should be the same.

The price of an American call option on an underlying asset that pays dividends, however, may diverge from its European counterpart. For an American call with dividends it may be beneficial to exercise the option prior to expiration. However a closed form formula for American option does not exist as the optimal exercise time is not specified but a variable to be determined itself. Thus we have to resort to numerical methods to obtain the option price.

Below is the MATLAB implementation for American call and put options

% BTM program for American vanilla call options
% call syntax: c = btm_AmeCall(S0,X,r,T,sigma,q,N)

function value = btm_AmeCall(S0,X,r,T,sigma,q,N)
% set up lattice parameters
dt = T/N; 
dx = sigma*sqrt(dt);
u = exp(dx); d=1/u;
df = exp(-r*dt);     % discount factor
p = (exp((r-q)*dt)-d)/(u-d);  % risk-neutral probability
% initialization
j = 1:1:N+1;  % range of index for price states
V = max(S0*u.^(2*j-N-2)-X,0);
% backward recursive through time
    for n = N-1:-1:0   
     j = 1:1:n+1;
     V = max(df*(p*V(j+1)+(1-p)*V(j)),S0*u.^(2*j-n-2)-X);
    end
value = V(1);

end

% BTM program for American vanilla put options
% call syntax: c = btm_AmePut(S0,X,r,T,sigma,q,N)

function value = btm_AmePut(S0,X,r,T,sigma,q,N)
% set up lattice parameters
dt = T/N; 
dx = sigma*sqrt(dt);
u = exp(dx); d=1/u;
df = exp(-r*dt);     % discount factor
p = (exp((r-q)*dt)-d)/(u-d);  % risk-neutral probability
% initialization
j = 1:1:N+1;  % range of index for price states
V = max(X-S0*u.^(2*j-N-2),0);
% backward recursive through time
for n = N-1:-1:0   
j = 1:1:n+1;
V = max(df*(p*V(j+1)+(1-p)*V(j)),X-S0*u.^(2*j-n-2));
end
value = V(1);

end

Chooser option

A standard chooser option gives its holder the right to choose, at a predermined time $T_c > t$ whether the T-maturity option is a standard European call or put with a common strike price $X$ for the remaining time to expiration $T - T_c$ . The payoff of the chooser option on the date of choice $T_c$ is

$V(S_{T_c},T_c) = \max(C(S_{T_c}, T-T_c; X),P(S_{T_c}, T-T_c; X)),$

where $T- T_c$ is the time to expiry in both call and put price formulas above, and $S_{T_c}$ is the asset price at time $T_c$ . For notational convenience, we take the current time $t = 0$ . Suppose the underlying asset pays a continuous dividend yield at rate $q$ . By the put-call parity relation, the above payoff function can be expressed as

$\begin{align} V(S_{T_c},T_c)&= \max(C,C+Xe^{r(T-T_c)}-S_{T_c}e^{-q(T-T_c)})\\ &= C + e^{-q(T-T_c)}\max(0,Xe^{-(r-q)(T-T_c)}-S_{T_c}). \end{align}$

Hence the chooser option can be viewed as the combination of one call with strike price $X$ and time to expiration $T$ and $e^{q(T-T_c)}$ units of put with strike price $Xe^{-(r-q)(T-T_c)}$ and time to expiration $T_c$ . Applying the Black-Scholes pricing approach, the value of the standard chooser option is given by

$\begin{align} V(S,0) &= Se^{-qT}N(x)-Xe^{-rT}N(x-\sigma\sqrt{T})+e^{-q(T-T_c)}\\ &\quad (Xe^{-(r-q)(T-T_c)}e^{-rT_c}N(-y+\sigma\sqrt{T_c})-Se^{-qT_c}N(-y))\\ &= Se^{-qT}N(x)-Xe^{-rT}N(x-\sigma\sqrt{T})+e^{-q(T-T_c)}\\ &\quad + Xe^{-rT}N(-y+\sigma\sqrt{T_c}) - Se^{-qT}N(-y), \end{align}$

where $S$ is the current asset price and

$x=\frac{\ln(S/X)+(r-q+\sigma^2/2)T}{\sigma\sqrt{T}}, \quad y=\frac{\ln(S/X)+(r-q)T+\sigma^2/2)T_c}{\sigma\sqrt{T_c}}$

Reference

Yue-Kuen Kwok. Mathematical Models of Financial Derivatives. Springer, Second edition.

Option Pricing with MATLAB - Part 1

European call option

European put option

American call and put options

Chooser option

Reference

Written by

Le Hoang Van

Supported by

Le Hoang Van