<jupyter_start><jupyter_text>Listed Volatility and Variance Derivatives**Wiley Finance (2017)**Dr. <NAME> | The Python Quants GmbHhttp://tpq.io | [@dyjh](http://twitter.com/dyjh) | http://books.tpq.io  Realized Variance and Variance Swaps  Introdution This chapter discusses basic notions and concepts needed in the context of variance swaps and futures. It covers among others the following topics:* **realized variance**: the basic measure on which variance swaps and variance futures are defined* **variance swap**: the definition of a variance swap and some numerical examples* **mark-to-market**: the mark-to-market valuation approach for a variance swap* **variance swap on EURO STOXX 50**: simple re-calculation of a variance swap given historical data* **variance vs. volatility**: major differences between the two measures  Realized Variance Historical or realized variance $\sigma^2$ generally is defined as $$\sigma^2 \equiv \frac{252}{N} \cdot \sum^{N}_{n=1}R_n^2$$where, for a time series $S_n, n=0, 1, ..., N$, the log returns are given by$$R_n \equiv \log \frac{S_{n}}{S_{n-1}}$$Here, it is assumed that there are 252 trading days per year and that the average daily return is zero. The simple application of these definitions yields values as decimals. Scaling by a factor of $100^2=10000$ gives values in percent.$$\sigma^2 \equiv 10000 \cdot \frac{252}{N} \cdot \sum^{N}_{n=1}R_n^2$$To simplify notation, we use the notation $\sigma^2$ instead of $\hat{\sigma}^2$ for the realized variance from here on. The concept of realized variance is easily illustrated by the use of historical data for the EURO STOXX 50 stock index. To this end, we read data from the index provider's Web site http://www.stoxx.com with Python and the pandas library. For details on using the pandas library for interactive financial analytics see chapter _Python Introduction_ or refer to Hilpisch (2018): _Python for Finance_.  As usual, we start with some Python library imports.<jupyter_code>import math
import numpy as np
import pandas as pd<jupyter_output><empty_output><jupyter_text>First, we need the complete URL of the data set.<jupyter_code>url = 'https://hilpisch.com/lvvd_eikon_eod_data.csv'<jupyter_output><empty_output><jupyter_text>Second, we read the data with the pandas library from that source.<jupyter_code>es = pd.read_csv(url, index_col=0, parse_dates=True)<jupyter_output><empty_output><jupyter_text>Let us inspect the most final five data rows.<jupyter_code>es.tail()<jupyter_output><empty_output><jupyter_text>Third, we select the EURO STOXX 50 index data from the just downloaded and imported data set, i.e. the data sub-set for the symbol SX5E. Using this sub-set, we generate a new pandas ``DataFrame`` object to store the data. The historical time series of daily closing levels of the EURO STOXX 50 can then easily be inspected by a call of the ``plot`` method. The following figure shows the graphical output.<jupyter_code>from pylab import mpl, plt
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
data = pd.DataFrame({'SX5E': es['.STOXX50E']})
data.plot(figsize=(10, 6));<jupyter_output><empty_output><jupyter_text>Historical index levels of EURO STOXX 50 index. Fourth, the log returns are calculated (in vectorized fashion, i.e. simultaneously over the whole time series) and stored as a new column in the pandas ``DataFrame`` object.<jupyter_code>data['R_n'] = np.log(data['SX5E'] / data['SX5E'].shift(1))<jupyter_output><empty_output><jupyter_text>Let us inspect the last five data rows of this new DataFrame object.<jupyter_code>data.tail()<jupyter_output><empty_output><jupyter_text>In the fifth step, we calculate the realized variance, again in vectorized fashion. With the following code we calculate the realized variance for every single date of the time series.<jupyter_code>## np.cumsum calculates the element-wise cumulative sum of an array/time series
## np.arange(N) gives an array of the form [0, 1, ..., N-1]
data['sigma**2'] = 10000 * 252 * (np.cumsum(data['R_n'] ** 2) 
                                  / np.arange(len(data)))<jupyter_output><empty_output><jupyter_text>The third column of the ``DataFrame`` object now contains the realized variance.<jupyter_code>data.tail()<jupyter_output><empty_output><jupyter_text>In the sixth and final step, one can now compare the index level time series with the realized variance over time graphically &mdash; see the following figure.<jupyter_code>data[['SX5E', 'sigma**2']].plot(subplots=True,
                                 figsize=(10, 8),
                                 grid=True);<jupyter_output><empty_output><jupyter_text>Historical index levels of EURO STOXX 50 index and realized variance (long-term). Now let us implement the same approach for a shorter, recent period of time, i.e. the second half of the year 2015. The realized variance has to be re-calculated since there is now a new starting date.<jupyter_code>## select time series data with date later/earlier than given dates
short = data[['SX5E', 'R_n']][(data.index > '2015-7-1')
                             & (data.index <= '2015-12-31')]

## calculate the realized variance in percent values
short['sigma**2'] = 10000 * 252 * (np.cumsum(short['R_n'] ** 2)
                                   / np.arange(len(short)))<jupyter_output><empty_output><jupyter_text>The first five rows of the new ``DataFrame`` object:<jupyter_code>short.head()<jupyter_output><empty_output><jupyter_text>A graphical comparison of the EUROS STOXX 50 time series data with its realized variance for the shorter time frame is displayed in the following figure.<jupyter_code>short[['SX5E', 'sigma**2']].plot(subplots=True,
                                 figsize=(10, 8),
                                 grid=True);<jupyter_output><empty_output><jupyter_text>Historical index levels of EURO STOXX 50 index and realized variance (short-term).  Variance Swaps Nowadays, variance swaps are popular financial instruments for volatility/variance trading and hedging purposes. See, for instance, the paper Bossu et al. (2005) for an overview of the features and characteristics of variance swaps.  Definition of a Variance Swap A variance swap is a financial instrument that allows investors to trade future realized variance against current implied volatility (the "strike"). The characteristics and payoff of a variance swap are more like those of a forward contract than those of a typical swap on interest rates, currencies, equities, etc.The payoff $h_T$ of a variance swap maturing at some future date $T$ is$$h_T = \sigma_{0,T}^2 - \sigma_K^2$$with $\sigma_K^2$ being the variance strike and $\sigma_K$ the volatility strike. At inception, i.e. at $t=0$, the volatility strike is set such that the value of the variance swap is zero. This implies that the volatility strike is set equal to the implied volatility $\sigma_i(0,T)$ for the maturity $T$.  Numerical Example Consider a Black-Scholes-Merton (1973) world with a geometric Brownian motion driving uncertainty for the index level of relevance (cf. Black and Scholes (1973) and Merton (1973)). The risk-neutral stochastic differential equation (SDE) in this model (without dividends) is given by: $$dS_t = r S_t dt + \sigma S_t dZ_t$$$S_t$ is the index level at time $t$, $r$ the constant risk-less short rate, $\sigma$ the instantaneous volatility and $Z_t$ a standard Brownian motion. For a comprehensive treatment of this and other continuous time financial models refer, for example, to Björk (2009). Instantaneous volatility (and variance) in this model world is constant which makes implied volatility also constant, say $\sigma_t =\sigma_i = \sigma = 0.2$. Given, for example, a Monte Carlo simulation of this model, realized variance might deviate from $\sigma^2 = 0.2^2 = 0.04$. Let us implement such a Monte Carlo simulation for the model. An Euler discretization scheme for the above SDE is given for $t \geq \Delta t$ by:$$S_t = S_{t-\Delta t} \exp \left( \left( r -\frac{\sigma^2}{2}\right)  \Delta t + \sigma \sqrt{\Delta t} z_t  \right)$$with  $\Delta t$ being the fixed time interval used for the discretization and $z_t$ a standard normally distributed random variable. A Python implementation might look like follows (refer to Hilpisch (2018) for details on Monte Carlo simulation with Python) .<jupyter_code>import sys
sys.path.append('./scripts')
import variance_swaps
variance_swaps.generate_path??<jupyter_output><empty_output><jupyter_text>Using this function and providing numerical parameters returns a pandas ``DataFrame`` with a single simulated path for the model. A sample path is shown in the following figure.<jupyter_code>S0 = 100  # initial index level
r = 0.005  # risk-less short rate
sigma = 0.2  # instantaneous volatility
T = 1.0 # maturity date
M = 50  # number of time intervals

data = variance_swaps.generate_path(S0, r, sigma, T, M, seed=100000)
data.plot(figsize=(10, 5));<jupyter_output><empty_output><jupyter_text>Sample path based on geometric Brownian motion. Given such a simulated path, one can calculate realized variance over time in the same fashion as above for the EURO STOXX 50 index.<jupyter_code>data['R_t'] = np.log(data['index'] / data['index'].shift(1))
## scaling now by M / T since returns are not necessarily daily returns
data['sigma**2'] = 10000 * M / T * (np.cumsum(data['R_t'] ** 2)
                                    / np.arange(len(data)))
data.tail()<jupyter_output><empty_output><jupyter_text>The following figure shows the results graphically.<jupyter_code>data[['index', 'sigma**2']].plot(subplots=True,
                                 figsize=(10, 8),
                                 grid=True);<jupyter_output><empty_output><jupyter_text>Geometric Brownian motion sample path with realized variance. In this case, the payoff $h_T$ of the variance swap at maturity is:<jupyter_code>data['sigma**2'].iloc[-1] - 20 ** 2<jupyter_output><empty_output><jupyter_text>In general, variance swaps have a notional that differs from a value of 1, i.e. the above value would have to be multiplied by a notional not equal to 1. It is market practice to define the variance swap notional in volatility terms:$$Notional = \frac{VegaNotional}{2 \cdot Strike}$$This can be done consistently due to the following relationship for a derivative instrument $f$ depending on some underlying $S$ with volatility $\sigma$ (and satisfying further technical assumptions):\begin{eqnarray*}\frac{\partial f}{\partial \sigma} &=& \frac{\partial f}{\partial (\sigma^2)} \cdot 2\sigma \\\Leftrightarrow \frac{\partial f}{\partial (\sigma^2)} &=& \frac{\frac{\partial f}{\partial \sigma}}{2\sigma}\end{eqnarray*} Say, we want a Vega notional of 100,000 currency units, i.e. we want a payoff of 100,000 currency units per volatility point difference (e.g. when realized volatility is 1 percentage point above the volatility strike). The variance notional then is:$$Notional = \frac{100000}{2 \cdot 20} = 2500$$<jupyter_code>Notional = 100000. / (2 * 20)
Notional<jupyter_output><empty_output><jupyter_text>Given this value for the variance notional, the payoff of the variance swap in the above numerical example would be:<jupyter_code>Notional * (data['sigma**2'].iloc[-1] - 20 ** 2)<jupyter_output><empty_output><jupyter_text>Mark-to-Market What about the value of a variance swap over time? A major advantage of working with variance (instead of volatility) is that variance is additive over time (when a mean of about zero is assumed). This gives rise to the following present value of a variance swap at time $t$, for a constant short rate $r$.$$V_t = Notional \cdot e^{-r(T-t)} \cdot \left( \frac{t \cdot \sigma_{0,t}^2 + (T-t) \cdot \sigma_i^2(t,T)}{T}-\sigma_K^2 \right) $$ The major component of the mark-to-market value of the variance swap is the time weighted average of realized variance $\sigma_{0,t}^2$ up until time $t$ and implied variance $\sigma_i^2(t,T)$ for the remaining life time from $t$ onwards.In the model economy, $\sigma_i^2(t,T)=\sigma_K^2=\sigma^2=400$. Therefore:$$V_t = Notional \cdot e^{-r(T-t)} \cdot \left( \frac{t \cdot \sigma_{0,t}^2 + (T-t) \cdot 400}{T} - 400 \right) $$For $t=0$, this obviously gives $V_t = 0$ as desired. This is readily implemented in Python given that we already have realized variance in a pandas ``DataFrame`` object. We calculate it again in vectorized fashion for $t=0, \Delta t, 2\Delta t..., T$.<jupyter_code>dt = T / M
t = np.arange(M + 1) * dt
t
sigma_K = 20
data['V_t'] = Notional * np.exp(-r * (T - t)) * ((t * data['sigma**2']
                        + (T - t) * sigma_K ** 2) / T - sigma_K ** 2) 
data.tail()<jupyter_output><empty_output><jupyter_text>Graphically, we get the result as presented in the following figure.<jupyter_code>data[['index', 'sigma**2', 'V_t']].plot(subplots=True,
                                 figsize=(10, 8),
                                 grid=True);<jupyter_output><empty_output><jupyter_text>Geometric Brownian motion sample path with realized variance and variance swap mark-to-market values.  Vega Sensitivity What is the sensitivity of the mark-to-market value of a variance swap with regard to implied volatility? Recall that the value itself is given by$$V_t = Notional \cdot e^{-r(T-t)} \cdot \left( \frac{t \cdot \sigma_{0,t}^2 + (T-t) \cdot \sigma_i^2(t,T)}{T}-\sigma_K^2 \right)$$Differentiation with respect to $\sigma_i$ gives a Vega of$$Vega_t = \frac{\partial V_t}{\partial\sigma_i} = Notional \cdot e^{-r(T-t)} \cdot \frac{T-t}{T} \cdot2 \sigma_i(t,T) $$At inception of the variance swap, we have a Vega of \begin{eqnarray*}Vega_0 &=& \frac{\partial V_0}{\partial\sigma_i} \\       &=& \frac{VegaNotional}{2 \cdot \sigma_i(0,T)} \cdot e^{-r(T-t)} \cdot \frac{T-t}{T} \cdot2 \sigma_i(0,T) \\       &=& e^{-rT} \cdot VegaNotional\end{eqnarray*}In this case, Vega equals the discounted Vega notional. For general $t$, we get\begin{eqnarray*}Vega_t &=& \frac{\partial V_t}{\partial\sigma_i} \\       &=& \frac{VegaNotional}{2 \cdot \sigma_i(t,T)} \cdot e^{-r(T-t)} \cdot \frac{T-t}{T} \cdot2 \sigma_i(t,T) \\       &=& e^{-r(T-t)} \cdot VegaNotional \cdot \frac{T-t}{T}\end{eqnarray*}This illustrates that Vega sensitivity diminishes over time and that it is proportional to the time-to-maturity.  Variance Swap on the EURO STOXX 50 We are now ready to do a historical re-calculation of a variance swap on the EURO STOXX 50. We will re-calculate a variance swap during June 2015. To this end we also use VSTOXX sub-index data for the shortest maturity available which provides us with a time series for the correct implied volatilities. EURO STOXX 50 data is already available.<jupyter_code>es.info()<jupyter_output><empty_output><jupyter_text>The VSTOXX data can be read from the same source (see chapter _Data Analysis and Strategies_).<jupyter_code>url = 'https://hilpisch.com/vstoxx_eikon_eod_data.csv'
vs = pd.read_csv(url, index_col=0, parse_dates=True)
vs.info()<jupyter_output><empty_output><jupyter_text>The data column `.V6I1` contains the index values (= implied volatility) for the neareast option series maturity available (i.e within a maximum of one month). For example, on the 1st of June 2015, the index values represent implied volatilities for the maturity on the third Friday in June 2015, i.e. the 19th of June. Maturity $T$ then is:<jupyter_code>## 15 trading days
T = 15.<jupyter_output><empty_output><jupyter_text>The variance swap we want to re-calculate should start on 1st of June 2015 and shall have a maturity until the 19th of June. It shall have a Vega notional of 100,000 EUR. First, let us select and collect the data needed from the available data sets.<jupyter_code>data = pd.DataFrame(es['.STOXX50E'][(es.index > '2015-5-31')
                        & (es.index < '2015-6-20')])
data.columns = ['SX5E']
data['V6I1'] = vs['.V6I1'][(vs.index > '2015-5-31')
                        & (vs.index < '2015-6-20')]<jupyter_output><empty_output><jupyter_text>The new data set looks as follows. Note that the VSTOXX sub-index is only available up until two days before the maturity date.<jupyter_code>data<jupyter_output><empty_output><jupyter_text>We forward fill the `NaN` values since for the vectorized calculations to follow we want to use these data points but they will have a negligible or zero influence anyway.<jupyter_code>data = data.fillna(method='ffill')<jupyter_output><empty_output><jupyter_text>We save the data set for later re-use.<jupyter_code>h5 = pd.HDFStore('./data/SX5E_V6I1.h5')
h5['SX5E_V6I1'] = data
h5.close()<jupyter_output><empty_output><jupyter_text>The implied volatility on 1st of June was 25.871%. This gives rise to a variance swap strike of $\sigma_K^2=25.871^2=669.31$. For a Vega notional of 100,000, the variance notional therefore is$$Notional=\frac{100000}{2 \cdot 25.871} = 1932.67$$<jupyter_code>data['V6I1'][0]
sigma_K = data['V6I1'][0]
Notional = 100000 / (2. * sigma_K)
Notional<jupyter_output><empty_output><jupyter_text>Three time series have to be calculated now:* log returns of the EURO STOXX 50 index* realized variance of the EURO STOXX 50 index* mark-to-market values of the variance swapFirst, the log returns.<jupyter_code>data['R_t'] = np.log(data['SX5E'] / data['SX5E'].shift(1))<jupyter_output><empty_output><jupyter_text>Second, the realized variance which we scale by a factor of 10,000 to end up with percent values and not decimal values.<jupyter_code>data['sigma**2'] = 10000 * 252 * (np.cumsum(data['R_t'] ** 2)
                                  / np.arange(len(data)))<jupyter_output><empty_output><jupyter_text>Third, the mark-to-market values. We start with the array of elapsed days.<jupyter_code>t = np.arange(1, 16)
t<jupyter_output><empty_output><jupyter_text>We assume a fixed short rate of 0.1%.<jupyter_code>r = 0.001
data['V_t'] = np.exp(-r * (T - t) / 365.) * ((t * data['sigma**2']
                + (T - t) * data['V6I1'] ** 2) / T - sigma_K ** 2)<jupyter_output><empty_output><jupyter_text>The initial value of the variance swap is zero.<jupyter_code>data['V_t'].loc['2015-06-01'] = 0.0<jupyter_output><empty_output><jupyter_text>The complete results data set is given below.<jupyter_code>data<jupyter_output><empty_output><jupyter_text>The payoff of the variance swap at maturity given the variance notional then is:<jupyter_code>Notional * data['V_t'][-1]<jupyter_output><empty_output><jupyter_text>Finally, a plot of the major results is presented in the following figure.<jupyter_code>data[['SX5E', 'sigma**2', 'V_t']].plot(subplots=True,
                                       figsize=(10, 8));<jupyter_output><empty_output><jupyter_text>EURO STOXX 50 historical index levels with realized variance and futures prices. We save the data for use in the next chapter.<jupyter_code>h5 = pd.HDFStore('./data/var_data.h5', 'a')
h5['var_swap'] = data
h5.close()<jupyter_output><empty_output><jupyter_text>Variance vs. Volatility Both variance and volatility are tradable asset classes. The following sub-sections discuss some differences between the two measures of variability and the two asset classes, respectively. See also Bennett and Gil (2012) on this topic.  Squared Variations Squared variations are in many application scenarios the better measure for variability compared to simple variations. By squaring variations, one makes sure that variations do not cancel each other out. Since volatility is generally defined as the square root of variance, both measures avoid the cancelling of positive and negative variations.  Additivity in Time Although both volatility and variance avoid the cancelling out of variations, there is a major difference between both when it comes to additivity. While variance is additive (linear) in time, volatility is convex (non-linear) in time.Assume we have $N$ return observations and assume $0<M<N$. We then have: \begin{eqnarray*}\sigma^2 &\equiv& \frac{252}{N} \cdot \sum^{N}_{n=1}R_n^2 \\         &=& \frac{252}{N} \cdot \left( \sum^{M}_{n=1}R_n^2 + \sum^{N}_{n=M+1}R_n^2 \right) \\         &=& \frac{252}{N} \cdot \sum^{M}_{n=1}R_n^2 + \frac{252}{N} \cdot \sum^{N}_{n=M+1}R_n^2  \\         &\equiv& \sigma_1^2 + \sigma_2^2\end{eqnarray*} Here, $\sigma_1^2$ is the  variance for the first part and $\sigma_2^2$ for the second part of the return observations. Note that one needs to keep the weighting factor constant at $\frac{252}{N}$ in order to retain additivity.This aspect can be illustrated by a simple numerical example. Consider first a function to calculate realized variance that we can re-use.<jupyter_code>## function to calculate the realized variance
rv = lambda ret_dat: 10000 * 252. / N * np.sum(ret_dat ** 2)<jupyter_output><empty_output><jupyter_text>Second, a simple example data set ...<jupyter_code>data = np.array([0.01, 0.02, 0.03, 0.04, 0.05])<jupyter_output><empty_output><jupyter_text>... of length $N=5$.<jupyter_code>N = len(data)
N<jupyter_output><empty_output><jupyter_text>Then, we easily see additivity.<jupyter_code>rv(data[:2]) + rv(data[2:])
rv(data)<jupyter_output><empty_output><jupyter_text>Next, we use the EURO STOXX 50 index data from before. Let us have a look at the year 2013 and the two halfs of the year.<jupyter_code>data = pd.DataFrame(es['.STOXX50E'][(es.index > '31-12-2012')
                             & (es.index < '01-01-2014')])
data.columns = ['SX5E']
## we need log returns
data['R_t'] = np.log(data['SX5E'] / data['SX5E'].shift(1))<jupyter_output><empty_output><jupyter_text>We have 256 index level observations and 255 return observations.<jupyter_code>N = len(data) - 1
N
var_1st = rv(data['R_t'][data.index < '2013-07-01'])
var_1st
var_2nd = rv(data['R_t'][data.index > '2013-06-30']) 
var_2nd<jupyter_output><empty_output><jupyter_text>Again, additivity is given for the realized variance.<jupyter_code>var_1st + var_2nd
var_full = rv(data['R_t'])
var_full<jupyter_output><empty_output><jupyter_text>Obviously, this is different when considering realized volatility instead of variance.<jupyter_code>vol_1st = math.sqrt(rv(data['R_t'][data.index < '2013-07-01']))
vol_1st
vol_2nd = math.sqrt(rv(data['R_t'][data.index > '2013-06-30']))
vol_2nd
vol_1st + vol_2nd
vol_full = math.sqrt(rv(data['R_t']))
vol_full<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1149642423424639, 'fraction_numerical': 0.02914586425583592, 'mean_word_length': 4.626012145748988, 'pattern_counts': {'":': 0, '<': 198, '<?xml version=': 0, '>': 196, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '47', 'dir': 'jupyter-structured-clean-dedup', 'id': '576544', 'license': 'CNRI-Python', 'n_tokens_mistral': 7435, 'n_tokens_neox': 6871, 'n_words': 2487, 'path': 'code/09_variance_swaps.ipynb', 'repo_name': 'yhilpisch/lvvd', 'size': '22243'}