MATLAB Examples

Contents

% The Stochastic Volatility Inspired parameterization of the Implied
% Volatility Smile

raw SVI parameterization

a = -0.041;
b = 0.1331;
m = 0.3586;
rho = 0.3060;
sigma = 0.4153;
param_raw = [a; b; m; rho; sigma];

k = -1.5:0.01:1.5;
w_r = svi_raw(k,param_raw);

plot(k,w_r, 'ob');

natural SVI parameterization

param_natural = svi_convertparameters( param_raw, 'raw', 'natural');
w_n = svi_natural(k,param_natural);

plot(k,w_r, 'ob');
hold on
plot(k,w_n, 'xr');
legend('Raw parameterizatio','Natural parameterization');
hold off

% check that results are equivalent
display(any(abs(w_n-w_r)>1e-10));
ans =

     0

Jump-Wings parameterization

tau = 1;
param_jw = svi_convertparameters( param_raw, 'raw', 'jumpwing', tau);
w_j = svi_jumpwing(k,param_jw, tau);

plot(k,w_r, 'ob');
hold on
plot(k,w_j, 'xr');
legend('Raw parameterizatio','Jump-Wing parameterization');
hold off

% check that results are equivalent
display(any(abs(w_j-w_r)>1e-10));
ans =

     0

Test conversion in other direction

v = 0.01742625;
psi = 0.1752111;
p = 0.6997381;
c = 1.316798;
vt = 0.0116249;
param_jw = [v; psi; p; c; vt];
param_raw = svi_convertparameters( param_jw, 'jumpwing', 'raw', tau);
param_natural = svi_convertparameters( param_jw, 'jumpwing', 'natural', tau);

w_r = svi_raw(k,param_raw);
w_n = svi_natural(k,param_natural);
w_j = svi_jumpwing(k,param_jw, tau);

plot(k, w_r, 'ob')
hold on
plot(k, w_n, 'xr')
plot(k, w_j, 'xg')
hold off

% check that results are equivalent
display(any(abs(w_n-w_r)>1e-10));
display(any(abs(w_j-w_r)>1e-10));
display(any(abs(w_j-w_n)>1e-10));
ans =

     0


ans =

     0


ans =

     0

Arbitrage-free Surface SVI

theta = 0.1:0.1:1;
rho = 0.3;
phifun = 'heston_like';
phi_param = (1+abs(rho))/4 + 1;

[totalvariance] = svi_surface(k, theta, rho, phifun, phi_param);
plot(k, totalvariance)

phifun = @(theta, lambda) 1./(lambda*theta) .* (1 - (1 - exp(-lambda*theta))./(lambda*theta));
param_natural = [0; 0; rho; theta(4); phifun(theta(4),phi_param)];
totalvariance_nat = svi_natural(k,param_natural);

plot(k, totalvariance_nat, 'ob')
hold on
plot(k, totalvariance(:,4), 'xr');
hold off

display(any(abs(totalvariance_nat-totalvariance(:,4)) > 1e-10));
ans =

     0

calibrate SSVI to Sample Data

clear
load sampledata
% the sample data set consists of bid and ask prices of options written on
% the S&P 500 index on 15-Sep-2011, the day before triple witching.

% remove non-usable data: non-positive bid-ask spread, bid below 3/8, itm options
forward = close.*exp((interest_rate-dividend_yield).*maturity);
pos_remove = ask <= bid | bid < 3/8 | (strike < forward & put==0) | (forward <= strike & put==1);
ask = ask(~pos_remove);
bid = bid(~pos_remove);
dividend_yield = dividend_yield(~pos_remove);
forward = forward(~pos_remove);
interest_rate = interest_rate(~pos_remove);
maturity = maturity(~pos_remove);
put = put(~pos_remove);
strike = strike(~pos_remove);
clear close pos_remove

% prepare data
bid(put==1) = bid(put==1) + exp(-interest_rate(put==1).*maturity(put==1)).*(forward(put==1)-strike(put==1));
ask(put==1) = ask(put==1) + exp(-interest_rate(put==1).*maturity(put==1)).*(forward(put==1)-strike(put==1));
mid = (bid+ask)/2;
log_moneyness = log(strike./forward);
clear put

% estimate implied volatility
implied_volatility_bid = blsimpv(forward.*exp(-interest_rate.*maturity), strike, interest_rate, maturity, bid);
implied_volatility_ask = blsimpv(forward.*exp(-interest_rate.*maturity), strike, interest_rate, maturity, ask);
implied_volatility     = blsimpv(forward.*exp(-interest_rate.*maturity), strike, interest_rate, maturity, mid);

% plot total implied volatility
total_implied_volatility = implied_volatility.^2 .* maturity;
maturities = sort(unique(maturity));
T = length(maturities);
figure
hold on
for t = 1:T
    pos = maturity==maturities(t);
    [x, idx] = sort(log_moneyness(pos));
    y = total_implied_volatility(pos);
    y = y(idx);
    plot(x,y);
end
hold off

% fit SSVI to data
phifun = 'power_law';
parameters = fit_svi_surface(implied_volatility, maturity, log_moneyness, phifun);


% get fitted total implied variances
total_implied_variance = implied_volatility.^2 .* maturity;
model_total_implied_variance = zeros(size(total_implied_variance));
model_implied_volatility = zeros(size(total_implied_variance));
for t = 1:T
    pos = maturity==maturities(t);
    [model_total_implied_variance(pos), model_implied_volatility(pos)] = ...
        svi_jumpwing(log_moneyness(pos), parameters(:,t), maturities(t));
end

% plot fitted data versus market data
for t = 1:T
    pos = maturity==maturities(t);

    subplot(3,5,t);
    plot(log_moneyness(pos), implied_volatility_ask(pos), 'xb');
    hold on
    plot(log_moneyness(pos), implied_volatility_bid(pos), 'xb');
    plot(log_moneyness(pos), model_implied_volatility(pos), 'xr');
    hold off
end

inter-/ extrapolation of SVI fit

clear
load sampledata
% the sample data set consists of bid and ask prices of options written on
% the S&P 500 index on 15-Sep-2011, the day before triple witching.

% remove non-usable data: non-positive bid-ask spread, bid below 3/8, itm options
forward = close.*exp((interest_rate-dividend_yield).*maturity);
pos_remove = ask <= bid | bid < 3/8 | (strike < forward & put==0) | (forward <= strike & put==1);
ask = ask(~pos_remove);
bid = bid(~pos_remove);
dividend_yield = dividend_yield(~pos_remove);
forward = forward(~pos_remove);
interest_rate = interest_rate(~pos_remove);
maturity = maturity(~pos_remove);
put = put(~pos_remove);
strike = strike(~pos_remove);
clear close pos_remove

% prepare data
bid(put==1) = bid(put==1) + exp(-interest_rate(put==1).*maturity(put==1)).*(forward(put==1)-strike(put==1));
ask(put==1) = ask(put==1) + exp(-interest_rate(put==1).*maturity(put==1)).*(forward(put==1)-strike(put==1));
mid = (bid+ask)/2;
log_moneyness = log(strike./forward);
clear ask bid put

% estimate implied volatility
implied_volatility     = blsimpv(forward.*exp(-interest_rate.*maturity), strike, interest_rate, maturity, mid);

% fit SSVI to data
phifun = 'power_law';
[parameters, theta, maturities] = fit_svi_surface(implied_volatility, maturity, log_moneyness, phifun);


log_moneyness = -1.5:0.1:1.5;

[~, idx] = ismember(maturities, maturity);
forward_theta = forward(idx);
interestrate_theta = interest_rate(idx);

tau_interp = 30/365.25;
forward_interp = interp1(maturities, forward_theta, tau_interp, 'linear');
interest_interp = interp1(zcb_days, zcb_rate, tau_interp*365.25, 'pchip');


[call_price, implied_volatility, total_implied_variance] = svi_interpolation(log_moneyness, tau_interp, ...
    forward_interp, interest_interp, parameters, theta, maturities, forward_theta, interestrate_theta);

strike = forward_theta*exp(log_moneyness);
plot(strike, call_price);