Overview

The objective of this report is to price several European call options on the S&P 500 Index with a maturity of 1 year (250 trading days) and strike prices from 100 to 2000 (intervals of 100) USD. This procedure is repeated under different correlation scenarios and parametric assumptions. By the end, it will be shown how stochastic volatility processes can generate model prices that exhibit implied volatility smiles similar to the Black-Scholes model.

Additional assumptions in this space include:

A monte carlo simulation method is implemented to model option prices that follow a continuous-time stochastic process. Additionally, the “GUIDE”, “RND”, “RQuantLib”, and “fOptions” CRAN R-packages are leveraged to assist in the required numerical computations.

Background Assumptions and Stochastic Processes

Continuous-Time Stochastic Processes

A continuous-time stochastic process meausures the uncertain variability in a random variable at any point in time. Moreso, the random variable we are modeling - the S&P 500 index - is a continuous variable. We also assume that the innovations, or price changes, in the random variable can only be practically forecasted given the current realized value, i.e. observations realized in the past are irrelevant concerning prediction tasks. However, a concession must be made. Past history of a random variable, specifically dealing with financial data, is useful in determining the statistical properties of the random variable’s stochastic process.

In the space of Markov stochastic processes, changes in random variables are independent normally distributed, which leads to the property that incremental changes in variances across any unit of time are additive. Wiener processes are a special form of Markov stochastic processess. They model the innovations of a random variable as independent and following a standard normal distribution, \(N(\mu=0,\sigma^2=1)\). The variance of this process equals \(\Delta t\), i.e. the variabiliy in a Wiener process is a function of time. \[\Delta z = \epsilon\sqrt{\Delta t} \]

Generalized Wiener Processes

Generalized Wiener processes relaxe the assumption that changes in a random variable follow a Wiener process with a drift and variance rate of 0 and 1, respectively. This model can be defined as

\[\Delta x = a\Delta t + b\Delta z\]

This property is important because it models changes in a random variable as normally distributed with a drift rate \(a\) and variace rate \(b^2\) with the following characteristics.

\[\mu(\Delta x) = a\Delta t\] \[\sigma(\Delta x) = b\sqrt{\Delta t}\] \[\sigma^2(\Delta x) = b^2 \Delta t\]

Ito Processes

The generalized Wiener process describd above assumes the drift and variance coefficients are constant. An Ito process (Ito, 1951) expresses the change in a random variable, which is in this context an option contract, as a function of the underlying in addition to time. The drift and variance rate change with respect to time.

\[\Delta x = a(x,t)\Delta t + b(x,t)\Delta z\]

Correlated Processes in Shocks

In this report, we model a stochastic process of a single random variable, the S&P500, in which the variance exhibits serial-correlated properties. This assumption allows us to model leverage effects and volatility clustering, of which both are empirical realities in financial time series data. This can be accomplished by drawing random uncorrelated samples from a standard normal distribution and, following Cholesky’s decomposition (factorization), deriving a correlation structure over \(\Delta t\):

\[\epsilon_1=\epsilon_1\] \[\epsilon_2=\rho\epsilon_1+\epsilon_2\sqrt(1-\rho^2)\]

Such a process provides additional practical usecases when modeling financial time series and pricing derivative contract on underlyings that exhibit this type of correlated stochastic process.

The Black-Scholes-Merton Option Pricing Model (BS)

Fundamental to BS is the property of risk-neutral valuation, which relies on the assumption that investors are indifferent, or netural, toward risk such that risk premiums are not an increasing function of risk levels. Also,the value of a derivative is independent of the underlying’s expected return. this leads us to assume that, in a risk-neutral world, the expected return on the underlying, and all available investment opportunities, is the risk-free rate. Subsequently, the expected payoff is discounted at the risk-free rate. Though deriving the BS differential equation is beyond the scope of this assignment, it is worth noting that an investor’s risk aversion is not a factor.

One of the principles that underline BS is the notion that the payoff of an EU option can be replicated by holding a defined amount of shares of the underlying, \(\Delta\), and borrowing at the risk-free rate over the life of the option. It follows that the cost to set up this replicating portfolio is equivalent to the payoff of the option, a property consistent with the Law of One Price (LOP). \(\Delta\) shares in the underlying within the replicating portfolio, along with price of the underlying, constantly ajust, ensuring the payoffs match.

Some additional assumptions include:

This report focuses on pricing EU call options. The Black-Scholes-Merton formula that solves for the price of a EU call (and put) option is expressed below. \(d_1\) and \(d_2\) represent the cumulative probability distribution function (CDF) of a standard normal distribution. Note that \(Ke^{-rT}N(d_2)\) reflects the cost of borrowing to construct the replicating portfolio.

\[d_1 = \frac{\ln(\frac{S}{K}) + (r + \frac{\sigma^2}{2})T]} {\sigma\sqrt T}\]

\[d_2 = d_1 - \sigma\sqrt T\]

The BS factors in the price of the underlying security to price derivative contracts. We leverage this closed-form solution throughout the remainder of the report to price EU calls and compare results with stochastic volatility models (SVM).

Benchmark Pricing of EU Call Option under the BS

We price 20 EU call options on a sequence of strikes from 100 to 2000, which increase by 100 USD. The volatility parameter is set equal to the long-run annualized volatility, 20%, which is assumed to be constant over the life of the option. The “RND” package is leverage to compute the call prices.

Strike Price (K) vs. Black-Scholes Call Prices
K BS_prices
100 966.2422200
200 867.4844399
300 768.7266599
400 669.9688884
500 571.2131446
600 472.5301653
700 374.6289189
800 280.3798672
900 195.5139026
1000 126.1912139
1100 75.4348810
1200 42.0061523
1300 21.9726266
1400 10.8948348
1500 5.1655563
1600 2.3604392
1700 1.0467287
1800 0.4531036
1900 0.1924187
2000 0.0805017

For the remainder of the report, these prices will remain the same. The output is consistent what we know about EU call options: the price of a call option is higher at lower strike prices, particularly as the strike price decreases relative to the spot price of 1,065 USD.

Monte Carlo Simulation of Stochastic Volatility

The following procedure applies a recursive Monte Carlo simulation method to simulate 1000 sample paths of the S&P 500 index over the course of a year on a daily frequency (250 days). The life of the option is 1 year, so the final index level for each path is used in the payoff calculation. Also, the payoff of the call option is dependent on the path followed by the underlying.

As discussed above, we are modeling the value of call option as a stochastic process with respect to volatility. Unlike the BS, we relax the assumption of constant volatility and assume volatility is stochastic. As mentioned previously, in the risk-neutral world, where \(\mu\) can be replaced with the risk-free rate, \(\r\), we can express price changes in the underlying index as following a Wiener process. Following Ito’s lemma,

\[d\ln(S) = (\mu - \frac {\sigma^2}{2})dt + \sigma dz \] \[\ln (S_{t+\Delta t})-\ln S_t = (\mu - \frac{\sigma^2}{2})\Delta t +\sigma\epsilon\sqrt{\Delta t}\] or alternatively, the path of the underlying under the risk-neutral probability can be expressed as,

\[S_{t+\Delta t} = S_te^{(r- \frac {\sigma^2}{2})\Delta t + \sqrt{\Delta t V_t}\epsilon^1_{t+\Delta t}}\], where \(V_t\) is the model’s variance rate. \(V_t\) has a drift coefficient, \(a\), and operates similarly to mean-reversion; it pulls the underlying’s variance rate back to the the long-run level of variance, \(V_L\), by a rate equivalently expressed by \(a\).

In addition to simulating a path for the underlying, we need to create a path for a second stochastic variable: the index’s volatility. This stochastic process can be represented as,

\[V_{(t+\Delta t)}= V_t + a(V_L-V_t)\Delta t + c\sqrt{\Delta t |V_t|}*[\rho\epsilon^1_{t+\Delta t}+\sqrt{(1-\rho^2)} \epsilon^2_{t+\Delta t}]\], where \(\epsilon^1\) and \(\epsilon^2\) are sampled from a standard normal distribution.

Let the payoff of a EU call option be expressed as \(max(S_T-K_i,0)\). The price of each call option is derived as follows:

\[S^{path,n}_{250} - K_i\]

,where \(n = 1,2,...,1000\) and \(i = 100,200,...,2000\).

-Take the average of the option payoffs on the 250th day.

\[\mu_{Payoff} = \frac {\sum_{i=1}^{1000}S^{250}_{i} - K_n}{1000}\]

\[ C_{EU} = \mu_{Payoff}e^{-rT}\]

This procedure is repeated 20 times for each respective strike price. In instances where \(K\) is very high, say 2000 USD, the probability of observing negative payoffs at the end of the simulation period is quite high. In that case, we specify a logical boolean statement that, in instances where the payoff on the 250th day is less than zero, the payoff is set equal to zero.

Scenario I: Correlation Structure of -30%

The Monte Carlo simulation is illustrated below. An important distinction to keep in mind is that \(\rho = -30 \%\).

Notice the extreme (blue) path that deviates from the rest by a significant margin. As history shows, realizing an exteme event like this is quite possible, but standard long-normal models of asset prices follow a normal distribution and underestimate the liklihood of observing extreme events in the tails of the distribution.

As derived above, we compute the prices of the call options.

Strike Price (K) vs. Black-Scholes Call Prices
K SVM_prices
100 959.4995794
200 860.7417993
300 762.0175948
400 663.7107706
500 566.1387606
600 469.3213866
700 373.8533130
800 280.8848923
900 193.2482861
1000 116.9839275
1100 60.7556898
1200 29.5385240
1300 14.3754594
1400 7.6408402
1500 4.7372014
1600 3.2268591
1700 2.3259874
1800 1.7113608
1900 1.2948267
2000 0.9985533

Next, we can compare the price differences between the log-normal BS model and the SVM by computing and plotting the implied volatilities of both models.

To derive more information, we compute the implied distribution of the S&P 500 from the SVM and BS prices.

Scenario II: Stochastic Volatility Model with Zero Correlation.

The same procedure in scenario I is repeated, with the exception that the correlation structure between the shock’s is changed to zero. It has been shown that (Hull & White, 1988) an EU option with a stochastic volatility process that is uncorrelated with the underlying’s price, has the same price as the BS price.

Strike Price (K) vs. SVM Call Prices
K SVM_prices
100 953.558182
200 854.800402
300 756.042622
400 657.458472
500 559.567490
600 462.412412
700 366.604888
800 273.398353
900 186.662513
1000 112.545038
1100 60.396447
1200 32.383240
1300 17.863174
1400 10.707046
1500 7.556885
1600 5.883333
1700 4.799302
1800 3.910482
1900 3.150889
2000 2.543513

Scenario III: Stochastic Volatility Model with Altered Drift Coefficients

In this last scenario, the original correlation structure of -30% is restored, but the \(c\) coefficient is adjusted from 0.85 to 0.15.

Analysis and Economic Interpretations

A Word About Implied Volatility

Implied Volatility (IV) can be described as the level (or magnitude) of volatility the market expects the underlying to realize in the future. As such, implied volatility is derived from option prices and is leveraged by traders and other investment professionals to monitor the market’s expectations of volatility on the underlying. This is particularly useful if the underlying is a broad market index.

The term “implied” can be attributed to the BS model. By inputing the implied volatility into the pricing model, we reach the market price of the option. As such, the market price of options are based on the BS model, and the implied volatility is the level of volatility at which the market price and BS price are equivalent, i.e.

\[P_{BS}=P_{MKT}\]

Furthermore, we can think of volatility smiles as the relationship between implied volatility levels and the strike price. In the figures above, we see that implied volatility is higher for options ITM or OTM as compared to options ATM. More importantly, the rate of change in IV rises sharply as the strike price decreases and the call option falls further ITM. This can be viewed a number of ways, though it is possible that the IV smile / smirk suggests the market expects downside volatility as more probable than upside volatility.

In this report, IV depends only on the strike price, although alternative methods allow IV to depend on the strike price, as well as the time to maturity (etc.).

Volatility Puzzle Explained

Under the Black-Scholes model, volatility is assumed to be constant; we should observe a straight horizontal line. However, as we have seen up until now, IV under BS is not constant; rather, volaitlity increases as the option moves away from the ATM strike price.

  • As the strike price, K, decreases, the price of an EU call option increases under the assumption that a lower strike price provides more margin and likelihood for the call to expire ITM. Increasing payoffs, under the rules of arbitrage pricing theory, must be met with higher permiums.

  • At a fixed maturity, as stirke prices decline, the volatility parameter must increase as well in order to maintain equality, \(P_{BS}=P_{MKT}\). What this tells us is that as the strike price declines, for options with equity underlyings, IV increases. The market’s expectations of higher future volatility continue to rise as K declines.

There are several possible expanations for this. One of the main arguments is the empirical observation of negative correlations between equity prices and volatility. Let it be known that the S&P 500 is an index of 500 individual stock (equity) positions. Nonetheless, as prices rise, volatility levels fall. The revese is also true, volatility levels rise as prices decline. This is known as the leverage effect. Another important argument that sheds light on this puzzle is the volatility feedback effect and crashphobia. According to Hull (2018),

“As volatility increases(decreases) because of external factors, investors require a higher (lower) return and as a result the stock price declines (increases).”

—Hull, 2018 Options, Futures, and Other Derivatives, 10th Edition

Crashphobia refers to the fear of huge market selloffs; “crashes”. As such, investors price in that fear into option prices, as seen in the implied volatility smile (smirk). Again, from Hull,

“Traders are concerned about the possibility of another crash similar to October, 1987, and they price options accordingly. There is some empirical support for this explanation. Declines in the S&P 500 tend to be accompanied by a steepening of the volatility skew. When the S&P increases, the skew tends to become less steep.”

—Hull, 2018 Options, Futures, and Other Derivatives, 10th Edition

Alternative methods of computing IV are used, including dividing the strike price by the spot price of the underlying- \(\frac {K}{S_0}\).

Stochastic Volatility Revisited

The central takeaway from implied volatility is that it is used to realistically gauge market volatility in a nonconstant, nonlognormal space. Under the assumptions of a stochastic volatility model, we realistically assume some degree of correlation between asset prices and volatility. SVM’s allow investment professional to measure the market’s future expectations of underlying volatility with a relativley greater degree of precision and realism. In other words, the biases of BS are mitigated with a SVM.

Refer to Scenario I:

  • Implied Volatility

    • Comparing the Implied Volatility charts, we find that the SVM exhibits higher levels of implied volatility as K increases from K 1000 on, whereas the BS shows IV relatively flat. IV is modeled similarly between both models as K decreases below 1000.

  • Implied Distribution

    • Let the implied distribution of the underlying asset be a function of the option price. The option price tells us something about the distribution of the undelrlying’s price over the maturity of the option. Using the derived price data, we can derive a probability density function (PDF) of the undelrying. We use the following equation to derive the PDF of the S&P 500.

\[\pi(k_n)=\frac{\frac{c(k_{n+1})-c(k_n)}{k_{n+1}-k_n}-\frac{c(k_n)-c(k_{n-1})}{k_n-k_{n-1}}}{k_n-k_{n-1}}\]

  • Based on the implied distributions for the BS and SVM, the SVM indicates the probability of extreme left-tail events are higher (heavier) compared to the implied PDF derived by BS. This is visible by the area under the left-side of the curve. Under the assumption that volatility and underlying prices are negatively correlated, we would expect to see more density in the left-tail, as the probability of price decreases is higher under increasingly-volatile environments.

Refer to Scenario II:

Under Scenario II, the volatility in the undelrying is a stochastic process. The change made is that the correlation between prices and volatility is set to zero.

  • Implied Volatility

    • Compared to Scenario I, we see further evidence of the benefits of implementing a SVM to price options. The IV increases as the option moves away from the ATM strike. However, in this scenario, the IV at higher stirkes is greater than in scenario I. This appears logical considering that in the first scenario, volatility was negatively correlated with asset prices. Increases in asset prices are met with less volatility as reflected in the low IV levels at high strikes. Here, IV is relatively higher and slightly increasing with the strike price.

  • Implied Distribution

    • The implied distributions are fairly similar, with the exception that their is more density in the right-tail of scenario II’s implied distribution compared to Scenario I. We are limited in the conclusions we can draw due to the small number of strike prices. Their are, however, more advanced numerical procedures for interpolating density values at other strike prices.

Refer to Scenario III:

Lastly, under the third scenario, we assume the original negative correlation structure of -30% in scenario I, but adjust the drift coefficient \(c\) from 0.85 (under scenario I) to 0.15. \(c\) represents the volatility of the variance process by which the SVM follows. Therefore, it is expected that implied volatility levels will be lower and the density in the tails of the implied distribution will be lower compared to scenario I.

  • Implied Volatility

    • It appears that the initial strike prices to the left/right of the ATM strike have higher IV’s in scenario I than scenario III. Additionally, IV is higher as strikes increase in scenario I. As the strike price reaches the lower bound, however, both scenarios appear to be similar.

  • Implied Distribution

    • Judging from the two distributions, it appears that under scenario I, where \(c\) = 0.85, their is more density outside the center of the PDF, particularly in the left-tail, compared to scenario III. The implied distribution in scenario III shows more density near the central tendency and in the right tail. This result can be attributed to the different values of the volatilty-of-variance coefficient. At higher levels, the variance process will exhibit more volatility. Under a negative correlation scenario, this is reflected in heavier left-tails, as in scenario I. At lower levels, the volatility process is more stable, as indicated by heavier (lighter) right (left) tails in scenario III.

Concluding Remarks

Stochastic volatility models assist investment professionals in pricing derivative securities with the advantage of relaxing assumptions of constant volatility and lognormal processes. Prices can be derived under realistic and empirical conditions of nonconstant volatility. We have derived and shown the implications of varying assumptions and parameters by implementing a monte carlo simulation method.

Appendix: R Code

# set.seed(6360)
# S <- 1065
# r <- 0.0125
# vol <- 0.2
# a <- 0.95
# c <- 0.15
# p <- -0.3
# K <- seq(100,2000,100)
# TTM <- 1
# t <- 1/250
# shock_1 <- matrix(0, ncol = 1000, nrow = 250)
# for (i in 1:ncol(shock_1)) {
#  shock_1[,i] <- rnorm(250, mean = 0, sd = 1)
# } 
# shock_2 <- matrix(0, ncol = 1000, nrow = 250)
# for (i in 1:ncol(shock_2)) {
#  shock_2[,i] <- rnorm(250, mean = 0, sd = 1)
# }

### MC Simulation of Underlying

  ## Model Stochastic Volatility Process

    vol_matrix <- matrix(0, ncol = 1000, nrow = 250)
    colnames(vol_matrix) <- paste("Sim", 1:ncol(vol_matrix))
    rownames(vol_matrix) <- paste("Day", 1:250)
    vol_null <- vol^2
    for (j in 1:ncol(vol_matrix)) {
  
      vol_matrix[1,j] <- vol_null+c*sqrt(t*abs(vol_null))*(p*shock_1[1,j]+sqrt(1-p^2)*shock_2[1,j])
  
      for (i in 2:nrow(vol_matrix)) {
        vol_matrix[i,j] <- vol_matrix[i-1,j] + a*(vol_null-vol_matrix[i-1,j])*t+c*
        sqrt(t*abs(vol_matrix[i-1,j]))*(p*shock_1[i,j] + sqrt(1-p^2)*shock_2[i,j])
      }
    }
  

  ## Model Price Changes in Underlying Index 

    price_matrix <- matrix(0, ncol = 1000, nrow = 250)
    colnames(price_matrix) <- paste("Sim", 1:ncol(price_matrix))
    rownames(price_matrix) <- paste(" Day", 1:250)
    price_null <- S    
    for (j in 1:ncol(price_matrix)) {
  
      price_matrix[1,j] <-price_null*exp((r-0.5*vol_null)*t+sqrt(t*vol_null)*shock_1[1,j])
  
      for (i in 2:nrow(price_matrix)) {
        price_matrix[i,j] <- price_matrix[i-1,j]*exp((r-0.5*abs(vol_matrix[i-1,j]))*t + 
                             sqrt(t*abs(vol_matrix[i-1,j]))*shock_1[i,j])
      }
    }

#    matplot(price_matrix, type = 'l',pch = NULL, 
#        main = '1000 Simulations of Index Level over 250 Days Following a 
#            Stochastic Volatility Process', xlab = "Days", 
#        ylab = 'Index Price Level', sub = "c parameter = 0.15")

### Price the Benchmark - EU Call Options using Black-Scholes OPM

  ## Black-Scholes Model
    BS_prices <- as.vector(rep(0,20))
    for (i in 1:length(BS_prices)) {
      BS_prices[i] <- price.bsm.option(s0=S, k=K[i], r=r, te=1, sigma=vol, y=0)$call
    }
    
    BS_df <- cbind(K, BS_prices)
    BS_df <- as.data.frame(BS_df)

### Price the MC Simulated Stochastic Volatility Model EU Call Options    

  ## Derive and Discount Average Option Payoff Per K - Yield 20 EU Call Option Prices
    SVM_prices <- rep(0, 20)
    payoff <- matrix(0, nrow = 1000, ncol = 1)
    for (i in 1:20) {
      payoff <- price_matrix[250,] - K[i]
      for (k in 1:length(payoff)) {
        if (payoff[k] < 0) {
          payoff[k] = 0
        }
      }
      SVM_prices[i] <- as.vector(mean(payoff) * exp(-r * 1))
    } 
    
    SVM_df <- cbind(K, SVM_prices)
    SVM_df <- as.data.frame(SVM_df)

### Compute Implied Volatility
  
  ## Black Scholes Model
    
    IV_BS <- as.vector(rep(0,20))
    
    for (i in 1:20) {
      
      IV_BS[i] <- GBSVolatility(price = BS_df[i,2], TypeFlag = "c", S = S, 
                                 X = BS_df[i,1], Time = 1, r = r, b = 0)
    }
    
  ## Stochastic Volatility Model
    
    IV_SVM <- as.vector(rep(0,20))
    
    for (i in 1:20) {
      
      IV_SVM[i] <- GBSVolatility(price = SVM_df[i,2], TypeFlag = "c", S = S, 
                                 X = SVM_df[i,1], Time = 1, r = r, b = 0)
    }

### Charts
 

#  plot.default(K, BS_prices, type = 'l', lwd = 4, 
#               main = "EU Call Price vs.Strike Price (BS)",
#               ylab = "Call Price", xlab = "Strike Price (K)", col = "cyan")
#
#  plot.default(K, SVM_prices, type = 'l', lwd = 4, 
#               main ="EU Call Price vs.Strike Price (SVM)", 
#               ylab = "Call Price", xlab = "Strike Price (K)", col = "black", 
#              sub = "Modeled via Monte Carlo Simulation of Stochastic Volatility Process;c parameter = 0.15")
#  
#  plot.default(K, IV_BS, type = 'l', lwd = 4, 
#               main = "Implied Volatility vs. Strike Price (BS)",
#               ylab = "Implied Volatility", xlab = "Strike Price (K)" , col = 'cyan')
#  
#  plot.default(K, IV_SVM, type = 'l', lwd = 4,
#               main = "Implied Volatility vs. Strike Price (SVM",
#               ylab = "Implied Volatility", xlab = "Strike Price (K)", col = 'black',
#               sub = "Modeled via Monte Carlo Simulation of Stochastic Volatility Proces;c parameter = 0.15")

### Compute Implied Distributions

  ## SVM Implied Distribution
    SVM_ID <- as.vector(rep(0,20))
      x <- SVM_prices
      k <- K
    for (i in 2:19) {
      SVM_ID[i] <- ((x[i+1]-x[i])/(k[i+1]-k[i])-((x[i]-x[i-1])/(k[i]-k[i-1])))/(k[i]-k[i-1])
    }
  
  ## BS Implied Distribution
    BS_ID <- as.vector(rep(0,20))
      x <- BS_prices
      k <- K
    for (i in 2:19) {
      BS_ID[i] <- ((x[i+1]-x[i])/(k[i+1]-k[i])-((x[i]-x[i-1])/(k[i]-k[i-1])))/(k[i]-k[i-1])
    }
    
#  plot.default(K, SVM_ID, type = 'l', lwd = 4, col = 'black', xlab = "Strike", 
#         ylab = "Density",
#         main = "Implied Distribution from EU Call Prices & Strikes (SVM) vs. Strikes",
#         sub = "Derived from a Stochastic Volatility Model; c parameter = 0.15")