import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()Time-Varying Volatility
Modelling Stock Returns
A simple model of stock returns assumes that the daily return of a stock can be modeled as: \[ r_{t+1} = \mu + \sigma \varepsilon_{t+1} \] where \(\{\varepsilon_{t}\}_{t \geq 0}\) are independent and identically distributed (iid) random variables so that \(\varepsilon_{t} \sim N(0,1)\) for all \(t \geq 0.\)
The evidence, however, suggests that for many stocks and indices the variance moves over time. This implies that conditional on the information available up to time \(t\), the next period return might be better represented as: \[ r_{t+1} = \mu_{t} + \sigma_{t} \varepsilon_{t+1}. \]
If we denote by \(\ev_{t}\) the conditional expectation operator, that is, our best guess for the future conditional on the information that we have until time \(t,\) we have that: \[ \ev_{t}(r_{t+1}) = \mu_{t} + \sigma_{t} \ev_{t}(\varepsilon_{t+1}) = \mu_{t}, \] which shows that \(\mu_{t}\) represents the expected return conditional on the information available up to time \(t.\)
Estimating \(\mu_{t}\) is an important but very difficult question in finance. We still lack a good model of time-varying expected returns. Nevertheless, in this application we are interested in estimating \(\sigma_{t}\) which is not really affected by \(\mu_{t}\) if we use daily returns. Therefore, we will not pay too much attention to modelling \(\mu_{t}.\)
Similarly, if \(\var_{t}\) denotes the conditional variance we have that: \[ \var_{t}(r_{t+1}) = \sigma^{2}_{t} \var_{t}(\varepsilon_{t+1})= \sigma^{2}_{t}, \] implying that \(\sigma^{2}_{t}\) is the conditional variance of the daily stock returns.
A Simple Estimate of the Variance
Suppose we have collected \(n\) daily stock returns. We can then estimate \(\mu\) by using: \[ \hat{\mu} = \frac{1}{n} \sum_{i=1}^{n} r_{i}. \]
To estimate \(\sigma,\) we can first estimate the variance as: \[ \hat{\sigma}^{2} = \frac{1}{n-1} \sum_{i=1}^{n} (r_{i} - \hat{\mu})^{2}. \]
The standard deviation is then the square root of the variance. Note that in order to have an unbiased estimator of the variance we divide by \(n-1\) since we lost one degree of freedom when we estimated \(\mu\), which is used in our estimation of the variance. In practice, though, since we will be using datasets with a significant number of daily observations, dividing by \(n-1\) or \(n\) is virtually irrelevant.
Note that stocks do not trade on weekends and official holidays. Therefore, there are usually 252 trading observations per year. As a consequence, in order to annualize \(\mu\) or \(\sigma^{2}\) we will multiply by 252.
A First Look at the Data
Getting the Data
We will use the following libraries to estimate our model.
We first get data on the S&P 500 using yfinance. The variable start_date specifies the beginning of our dataset. The S&P 500 dataset is stored in sp500. We then compute daily returns using the Close price, drop the first observation since it is equal to NaN, and store the results in returns. Note that I multiply the returns by 100 to express them as a percentage.
start_date = '1993-01-01'
sp500 = yf.download('^GSPC', start=start_date, auto_adjust=False, progress=False, multi_level_index=False)
returns = 100 * sp500['Close'].pct_change().dropna()Plotting Daily Returns
Let us generate a plot of the data first.
plt.figure(figsize = (12,6))
plt.plot(returns)
plt.title('Daily S&P 500 Returns')
plt.xlabel('Date')
plt.ylabel('Return (%)')
plt.axhline(y=0, linewidth=1)
plt.xlim(returns.index.min(), returns.index.max())
plt.show()
Descriptive Statistics of Daily Returns
We see from the graph that daily returns are volatile and that volatility clusters over time — periods of large swings tend to be followed by more large swings. Let us now compute the mean and standard deviation of daily returns. Note that since the variance is annualized by multiplying by 252, the standard deviation is annualized by multiplying by the square root of 252. The numpy function sqrt() allows us to compute square roots.
Since returns is a Pandas dataframe, we can use the methods mean() and std() to compute the mean and the standard deviation.
avg = returns.mean()
sd = returns.std()
avg_ann = 252 * avg
sd_ann = np.sqrt(252) * sdOne way to present the results is to create a new Pandas dataframe.
pd.DataFrame([[avg, avg_ann], [sd, sd_ann]],
index=['Average Return (%)', 'Standard Deviation (%)'],
columns=['Daily', 'Annualized']).round(2)| Daily | Annualized | |
|---|---|---|
| Average Return (%) | 0.04 | 10.03 |
| Standard Deviation (%) | 1.16 | 18.38 |
Rolling Window Approach
One way to capture time-varying volatility is to use a rolling window when estimating the standard deviation of daily stock returns. In a rolling window approach, we estimate \(\sigma_{t}\) as: \[ \hat{\sigma}_{t}^{2} = \frac{1}{k-1} \sum_{i=1}^{k} (r_{t-i} - \hat{\mu}_{t})^{2}, \] where \[ \hat{\mu}_{t} = \frac{1}{k} \sum_{i=1}^{k} r_{t-i}. \]
In the example below, I use a window of 60 observations to estimate the standard deviation. Note that shift(1) is required to make sure that \(\sigma_{t}\) uses information up to time \(t-1\) in estimating the volatility.
k = 60
sd_rolling = returns.rolling(k).std().shift(1).dropna()
sd_rolling_ann = np.sqrt(252) * sd_rollingThe figure below plots the annualized rolling volatility estimate over time.
plt.figure(figsize = (12,6))
plt.plot(sd_rolling_ann)
plt.title('Annualized S&P 500 Volatility')
plt.xlabel('Date')
plt.ylabel('Standard Deviation (%)')
plt.axhline(y=sd_ann, linewidth=1)
plt.xlim(sd_rolling_ann.index.min(), sd_rolling_ann.index.max())
plt.show()
The figure shows that the rolling volatility reverts over time to the average level of volatility.
EWMA Volatility Estimation
A drawback of the rolling window approach is that it assigns equal weight to all \(k\) observations in the window and zero weight to all observations outside it. This can cause abrupt jumps in the volatility estimate when a large return enters or exits the window.
The Exponentially Weighted Moving Average (EWMA) model addresses this by assigning weights that decay geometrically as observations become older. The conditional variance is updated recursively as: \[ \sigma^{2}_{t} = \lambda \sigma^{2}_{t-1} + (1 - \lambda) r^{2}_{t-1}, \] where \(\lambda \in (0, 1)\) is the decay factor. A value of \(\lambda\) close to 1 makes the estimate slow to react to new information, while a smaller \(\lambda\) lets recent observations dominate. RiskMetrics, the industry standard developed by J.P. Morgan, recommends \(\lambda = 0.94\) for daily data.
Note that we can substitute the recursion backwards to express \(\sigma^{2}_{t}\) as a weighted sum of all past squared returns: \[ \sigma^{2}_{t} = (1 - \lambda) \sum_{i=1}^{\infty} \lambda^{i-1} r^{2}_{t-i}. \] The weights \((1-\lambda)\lambda^{i-1}\) sum to one and decline geometrically with \(i\), so recent returns receive the most weight.
Pandas implements this directly through ewm(). The alpha parameter in Pandas corresponds to \(1 - \lambda\), so for \(\lambda = 0.94\) we set alpha = 0.06. We use shift(1) to ensure \(\sigma_{t}\) only uses returns observed up to time \(t-1\).
lam = 0.94
var_ewma = returns.pow(2).ewm(alpha=1-lam, adjust=False).mean().shift(1).dropna()
sd_ewma = np.sqrt(var_ewma)
sd_ewma_ann = np.sqrt(252) * sd_ewmaThe figure below overlays the EWMA volatility estimate with the rolling window estimate for comparison.
plt.figure(figsize=(12, 6))
plt.plot(sd_rolling_ann, label='Rolling (k=60)', alpha=0.7)
plt.plot(sd_ewma_ann, label=r'EWMA ($\lambda=0.94$)', alpha=0.7)
plt.title('Annualized S&P 500 Volatility')
plt.xlabel('Date')
plt.ylabel('Standard Deviation (%)')
plt.axhline(y=sd_ann, linewidth=1)
plt.xlim(sd_ewma_ann.index.min(), sd_ewma_ann.index.max())
plt.legend()
plt.show()
Compared to the rolling window, the EWMA estimate avoids abrupt jumps when a large return enters or exits the window. Both approaches, however, treat \(\lambda\) and \(k\) as fixed parameters chosen by the analyst rather than estimated from the data. In the next notebook, we introduce the GARCH model, which estimates these parameters directly from the data and provides a more principled framework for modelling time-varying volatility.
Practice Problems
Problem 1 Generate a plot of the rolling volatility of NASDAQ Composite (^IXIC) over the same period.
Solution
start_date = '1993-01-01'
sp500 = yf.download('^IXIC', start=start_date, auto_adjust=False, progress=False, multi_level_index=False)
returns = 100 * sp500['Close'].pct_change().dropna()
sd_ann = np.sqrt(252) * returns.std()
k = 60
sd_rolling = returns.rolling(k).std().shift(1).dropna()
sd_rolling_ann = np.sqrt(252) * sd_rolling
plt.figure(figsize = (12,6))
plt.plot(sd_rolling_ann)
plt.title('Annualized Nasdaq Volatility')
plt.xlabel('Date')
plt.ylabel('Standard Deviation (%)')
plt.axhline(y=sd_ann, linewidth=1)
plt.xlim(sd_rolling_ann.index.min(), sd_rolling_ann.index.max())
plt.show()
Problem 2 Using the NASDAQ Composite data from the previous problem, compute the EWMA volatility with \(\lambda = 0.94\) and generate a plot comparing it with the rolling window estimate.
Solution
start_date = '1993-01-01'
sp500 = yf.download('^IXIC', start=start_date, auto_adjust=False, progress=False, multi_level_index=False)
returns = 100 * sp500['Close'].pct_change().dropna()
sd_ann = np.sqrt(252) * returns.std()
k = 60
sd_rolling = returns.rolling(k).std().shift(1).dropna()
sd_rolling_ann = np.sqrt(252) * sd_rolling
lam = 0.94
var_ewma = returns.pow(2).ewm(alpha=1-lam, adjust=False).mean().shift(1).dropna()
sd_ewma_ann = np.sqrt(252) * np.sqrt(var_ewma)
plt.figure(figsize=(12, 6))
plt.plot(sd_rolling_ann, label='Rolling (k=60)', alpha=0.7)
plt.plot(sd_ewma_ann, label=r'EWMA ($\lambda=0.94$)', alpha=0.7)
plt.title('Annualized Nasdaq Volatility')
plt.xlabel('Date')
plt.ylabel('Standard Deviation (%)')
plt.axhline(y=sd_ann, linewidth=1)
plt.xlim(sd_ewma_ann.index.min(), sd_ewma_ann.index.max())
plt.legend()
plt.show()