Genius Mathematics Consultants

Model validation is a strict regulatory requirement for many financial services businesses, essential for risk management and front office profitability. It’s also critical for many other industries as well. In the medical device, aviation and mining industries, model failure could not only lead to financial loss, but cost lives. In so many industries, models must be thoroughly validated both before deployment, and at regular intervals after deployment. They must also be stress tested on unusual scenarios to ensure they will remain robust.

In order to speed up and broaden the scope of model validation, many organisations are considering building automated model validation tools. In this article, we’ll look at how automated model validation differs from manual validation, and how both are an essential part of a validation ecosystem.

What is model validation?

Model validation is a process to regularly test and monitor mathematical models. This spans the conceptual soundness of the model, the correctness of the coding implementation, and the interactions between systems such as data interfaces. The context could be option pricing or risk models in finance, machine learning models to detect mineral deposits in mining, automated trading systems, or medical monitoring devices. Increasingly, AI and machine learning models need to be validated. And of course, the wide-spread use of impressive but error-prone AI tools like ChatGPT necessitates a whole new world of model validation. The process typically involves:

• Reviewing the correctness of the mathematical methodology and documenting it
• Assessing the any model limitations or boundaries of validity
• Building an independent model, either of the same methodology or a different one, and ensuring the output of both models are within some acceptable tolerance
• Monitoring the ongoing appropriateness of the model considering changes in downstream and upstream systems, and changes in the environment the model operates in
• Stress testing the model to ensure robustness under unusual scenarios (eg financial downturns)
• Documenting all findings

What is automated model validation?

Automated model validation is the use of software tools and frameworks to autonomously and systematically test and monitor models. Typically these tools:

• Automatically pull in and format data from source systems (where a human validator might need to examine a GUI one number at a time)
• Verify huge numbers of model outputs by comparing against independent model implementations

Benefits of automated model validation

A downside of automated model validation systems is they may require a significant initial investment to build. So what are the advantages?

• Once set up, automated model validation is far less time intensive than model validation. An automated system could examine thousands of trades in a trading book, where a manual validator has time to check a dozen.
• Automated systems make it easy to quickly retest the models after system updates, or set up regular periodic validations or continuous monitoring.
• Scalability: automated systems can have a scale beyond the scope of human validators
• Automated systems may free up human validators to focus more on conceptual model review
• To the extent that less manual work is required, they may reduce costs.

Benefits of manual model validation

• Manual validation involves reviewing the assumptions and methodology of the model, often in a changing environment. This is something which doesn’t happen at all with automated model validation.
• How can you be sure the independent implementation used by the automated validation system is correct? To an extent, automated validation actually results in two models that need to be validated. And where both models have similar assumptions or methodologies, it’s entirely possible for the models to agree and both be wrong.
• Manual model validators also provide another very important function – the preparation of detailed validation reports clarifying the methodology, assumptions and limitations of the model where documentation is usually sorely lacking.
• Human validators can improve model design, not just test model outputs.
• Credibility with the regulator: in the case of regulated industries, human expert analysis may be considered essential.

Why a hybrid approach is best practice

The best approach to model validation is a hybrid approach, combining manual model validation with automated model validation. This is where manual model validators carefully check model assumptions, work to improve and maintain the accuracy of the automated model, and produce high quality documentation about the models. At the same time, the automated system allows for much more frequent and large scale testing and monitoring of trades and models. Having found serious errors in models that had already been validated several times, this author can tell you that the need for expert validators will never go away.

We offer model validation consulting services, including both manual model validation and the design of automated model validation systems. Contact us to learn more.

Choosing between a vendor trading system and an in-house solution is one of the most important technology decisions banks, trading firms, and hedge funds make. This article explores the trade-offs and helps decision-makers navigate this complex choice – build or buy?

Should you build your own quantitative trading system in-house, or to pay ongoing licensing fees for a vendor trading platform such as Murex or Calypso? This decision has profound implications for the cost, flexibility, and competitiveness of your trading operations, and it’s not a decision that can be easily changed in the future.

Vendor trading systems will likely be faster to market, with faster access to extensive functionality, and models that have undergone many cycles of validation already. But any functionality you don’t need simply adds unnecessary complexity. Systems like Murex have a significant learning curve to set up and configure, with Murex technical support staff required to be hired at your firm. This is in additional to technical support staff at the vendor that you will be billed for. Adding your own models and functionality may be harder, more expensive or not possible at all. Model development and bug fixing is likely to be much faster with all source code in-house, leading to a more agile business. Model documentation can be maintained in-house, instead of being at the mercy of how conscientious the vendor has been with their documentation.

For regulated entities, vendor systems can make model validation and regulatory reporting more complex. Without source code access, determining what models are doing or how they are failing becomes a complex business problem, and demonstrating compliance to auditors and regulators more time consuming, which increases cost and operational risk.

In terms of cost, developing your own system might cost more initially, but vendor fees and hidden charges can easily be underestimated. These will have no expiry and will eventually accumulate to more than the cost of building your own system. Once you have signed on, moving away from the vendor setup would be so expensive that you’re at their mercy in terms of fees.

The interests of the vendor and your own firm are never entirely aligned. Vendors are juggling hundreds of clients. When engaging in model validation of the vendor systems, I’ve noticed that vendors are often reluctant to admit their systems have a fault. In a particularly egregious case, I saw a vendor deliberately obstruct model validation staff in order to avoid admitting there was a serious problem with their model. After months of fighting with the vendor, they finally fixed the fault and whole book got revalued.

Vendor Trading Systems (e.g., Murex, Calypso)

Pros:

• Speed to market: Because vendor systems are already built, you will likely be operational sooner.
• Broad functionality: Immediate access to a wide range of models and configuration options, including derivative pricing models and market risk models.
• Model confidence: The models of very large vendors have been repeated validated by their clients, reducing the likelihood of faults.

Cons:

• Very high licencing and maintenance costs: The cost of signing on is only the beginning, as vendors charge further fees for support staff, unlocking more advanced models or modules, and even fixing faults.
• Vendor lock-in: once you have heavily invested in the vendor’s systems, it’s too expensive to change to another option. This gives you little negotiating power with the vendor.
• Less flexibility: You’re largely limited to the existing functionality the system and any inherent design limitations.
• Integration: The vendor systems must be integrated with your own
• Unused Complexity: the broad functionality needed for the vendor’s many clients makes configuration and use more complicated, and increases the risk of model misconfiguration.
• Model documentation: Although vendors provide documentation, the quality of this can vary. With no direct access to the source code, it can be challenging to determine what the models are doing.
• Politics and misaligned interests: The vendor may be more focused on extracting money, juggling many other clients, and avoiding admitting to any faults with their systems, than in helping your business.

Building an In-House Trading System

Pros:

• Fully customizable: An in-house system is tailored for exactly your requirements
• Streamlined: No unused functionality adding complexity
• Unique or cutting-edge: Instead of using the system everyone else is using, you can innovate for a competitive advantage.
• Agility: Rapid implementation of fixes and new features with fewer layers of bureaucracy and go-betweens.
• Cheaper: Likely cheaper in the long run, but with significant upfront costs.
• Model documentation: you have direct access to the source code when questions arise about the behaviour of the models

Cons:

• Longer time to market: The vendor solution can likely be operational sooner.
• Required expertise: You may require more quantitative expertise in-house. However, you could alternatively partner with a suitable consulting firm.

At Genius Mathematics Consultants, we specialize in providing quant support to help firms design and implement in-house trading systems. Our consulting services span quantitative and algorithmic trading, derivative pricing models, and risk model design such as Market VaR calculations. We ensure that your technology not only runs efficiently but also meets regulatory standards. We bring the ability to translate cutting-edge research into reliable production code. By partnering with us, you gain access to flexible quant expertise that supports your team in building systems tailored to your strategies—without the cost and complexity of maintaining in-house quant staff. To learn more about how your business can partner with our consultants, Contact Us.

I previously wrote an article for Dorian Trader which explored using ridge regression to develop a trading strategy. The article showed how ridge regression helped prevent overfitting as compared to regular regression, and was often uncannily good at predicting the peaks and troughs of trading data. See the original article here.

I also created a follow up article for Crypto News testing the same strategy on cryptocurrencies.

In those articles, the strategy was limited to only using a single asset. That is, the signals used as inputs were calculated only from the price history of the same asset that we were trading. In this article, we’ll explore whether the price history of closely related assets can contain valuable predictive information. This could be, for example, stocks from a similar industry, or different crypto currencies, both of which tend to move together. This is usually called statistical arbitrage or pairs trading (but note that despite the name, one could do it with any number of related assets, not necessarily two).

The basic mechanism

Pairs trading or statistical arbitrage makes use of the tendency for two assets to move together. It may mean that when one asset moves, the other is more likely to follow suit. It may also mean that the further apart the two assets become, the more likely they are to converge together again.

When we applied ridge regression to the single asset case, we used moving averages and regression lines of various lengths as our main signals. Here, we will do the same, except that we will have an additional set of signals generated from the second asset. We will use three moving averages and three regression lines of different lengths, so the regression has a choice of using whichever is more predictive. You could also use more than three, especially initially, if you are unsure of how much time might elapse between divergence and convergence of the related assets.

The mechanism here is that by including moving averages of different lengths, a linear model like ridge regression can implicitly include their differences as well, so that it can capture how much each asset has diverged (up or down) compared to its earlier values. It also implicitly includes the difference between the divergence of the first asset and the divergence of the second asset, so that it can capture how much one asset has diverged relative to the other.

If we believe that there is a lag between divergence and convergence, we could consider introducing a lag between the signal calculation time and the time at which the ridge regression tries to predict the price. However, by including signals like moving averages calculated over a variety of time windows, the model will be able to automatically calculate a few lagged moving averages anyway (as these are a linear combination of existing moving averages).

An additional issue springs to mind in the two asset case, is that the two assets may have very different scales, and the ratio between the two could drift over time. For this reason, one might consider taking various ratios of the signals and including them as additional input signals:

• For each asset, we can generate additional signals by taking the ratios of moving averages of different lengths, and the ratios of regression lines of different lengths.
• We’ll also add the ratios of the same signal between the two assets. This will allow the ridge regression to track how much the the assets have diverged in terms of their ratio.

The reality is, if the data window is not too large, the scales of the assets won’t change too much over the data set, and the relevant ratios can be assumed constant and will be automatically captured by the coefficients chosen by the ridge regression. So one might make a case that it’s not necessary to include the ratios when fitting over a relatively small time window (we will use one week). However, I’ve chosen to include all the mentioned ratios anyway. If nothing else, being able to see the coefficients for these signals provides more insight into how the ridge regression is working.

At the moment the model trains only on the most recent week of data. Although more data might seem better, more recent data is more relevant and allows the model to adapt to recent market conditions.

Cryptocurrencies

We use one week of data ending on the 14th of August 2025, with the test set being approximately the final three hours of data.

It’s most likely that movements in the big coins are going to presage movements in the smaller coins, rather than the other way round. For this reason it’s most sensible to use BTC as an input when trying to predict ETH. Below we show the coefficients that the ridge regression has assigned to each of our signals.

RegLine5 = -0.707
RegLine20 = 0.509
RegLine50 = -0.403
MovAvg5 = -0.323
MovAvg20 = 0.185
MovAvg50 = 1.127
MovAvg5/20 = -0.442
MovAvg5/50 = 1.025
MovAvg20/50 = -0.674
RegLine5/20 = 0.055
RegLine5/50 = -0.205
RegLine20/50 = 0.129
Close_BTC-USD = 0.014
RegLine5_BTC-USD = -0.046
RegLine20_BTC-USD = 0.542
RegLine50_BTC-USD = -0.142
MovAvg5_BTC-USD = -0.051
MovAvg20_BTC-USD = -0.192
MovAvg50_BTC-USD = -0.235
MovAvg5/20_BTC-USD = 0.249
MovAvg5/50_BTC-USD = -0.496
MovAvg20/50_BTC-USD = 0.291
RegLine5/20_BTC-USD = 0.096
RegLine5/50_BTC-USD = -0.18
RegLine20/50_BTC-USD = 0.117
MovAvg5/5 = -0.555
MovAvg20/20 = 0.139
MovAvg50/50 = 0.116
Strategy profit = 107.596
Buy Hold profit = 19.161

The first 12 signals are all calculated from ETH price data. They consist of the linear regression lines (across three windows of different lengths), moving averages, ratios of the moving averages and ratios of the linear regression predictions. The following 12 signals are identical except calculated using BTC price data. The final three signals are the ratios of the the ETH moving averages to the BTC moving averages.

Since ridge regression tries to make the coefficients small when it doesn’t harm the predictive power too much, and the signals are scaled prior to fitting, one can to an extent judge the significance of each signal by the size of the coefficient.

• The most significant ETH signals are the 50 minute moving average, and the ratio between the 5 and 50 minute moving averages. The other moving average ratios, and the regression lines, are also significant.
• The BTC signals are less significant than the ETH signals with smaller coefficients, as you would expect. The most significant are the 20 minute regression line, and the ratio between the 5 and 50 minute moving average.
• The ratio between the 5 minute BTC and 5 minute ETH moving average also makes a contribution.

Looking at the plot, we see what we often see with ridge regression – it can be very good at placing sell markers at the peaks, and buy markers at the troughs. This gives a visual representation of why the strategy profit is 107.6 while the buy and hold profit is only 19.2.

But how much impact is the second asset, BTC, having here?

Isolating the impact of the second asset

We naturally want to know how much predictive power is coming from the second asset (BTC), and how much is coming from ETH’s own price history. We can see from the size of the coefficients that the ETH price history is more significant than the BTC price history.

To further evaluate this question, we can omit the BTC signals and see how much the prediction / profitability degrades.

Now when we do this exercise, it’s important that we make the test and fit sets the same. This is because when the test set is outside of the data seen during fitting, it can be totally arbitrary and could favour any one of the models simply by chance. This makes comparison of different models very difficult. This might be mitigated by using a very large test set, which would have to be assumed to be representative of all future possibilities. However, this requires that a much larger amount of data be obtained, and the fact that market regimes change over time raises it’s own challenges. By using overlapping test and fit sets, we can test which of the three models can match the data most successfully, albeit without considering potential overfitting. In this case, I’ve used the entire week’s dataset for both testing and fitting (which is why the numbers are different from the previous section). For reference, the buy and hold profit was 945.3 in this case.

• Fit and backtest the model using only ETH signals – Strategy profit = 1343.0, Rsquared = 0.011
• Fit and backtest the model using both – Strategy profit =1556.5, Rsquared = 0.014

These numbers make a possible case that including the BTC signals improves the strategy. However, the case isn’t water tight due to the fact that including additional signals increases the opportunity for overfitting.

More than two assets

Statistical arbitrage, or the principles behind pairs trading, can work with any number of assets. For example, one could use a whole basket of stocks from the same industry. In the case of crypto trading, we can consider using other coins in our prediction. Movements in smaller coins are unlikely to have as big an impact as Bitcoin on the price of ETH. Let’s try the exercise again including the next four largest coins by market capitalization: XRP, BNB, SOL and DOGE.

Using five coins as additional signals: Strategy profit = 2051.7, Rsquared = 0.015

The strategy profit is distinctly higher than when using only ETH, and when using ETH and BTC as inputs. All of the coins had signals with decent sized coefficients, suggesting that the ridge regression found them significant. However, as mentioned already, overfitting is a concern here.

Risks of multi-asset strategies

It’s important to consider that a strategy involving the signals of other assets may be more dangerous than one involving only one asset. This is because although a group of assets may often move similarly, there’s the possibility for their behaviour to significantly diverge for periods of time. For example, two crypto coins may often move similarly as they are both affected by general sentiment about crypto. But if one coin were to experience significant growth or decline based on events pertaining only to that asset, the relationship may cease to hold. Safeguards need to be added to the strategy to enable it to disregard the signals of other assets under certain circumstances.

I recently wrote an article for Dorian Trader on using a neural network to replicate option pricing models. The idea behind using neural networks in option pricing, is that while the original numerical model might be slow, the neural network (once trained) is lightning fast at pricing options.

In that article, I fitted a neural network to the vanilla Black-Scholes formulae as a proof-of-concept. I used ChatGPT to generate an initial piece of code, before tweaking it myself to improve the accuracy and configure the output how I wanted it. So this exercise is also an interesting demonstration of AI-assisted coding.

You can read the full article here

Below I’ve attached the python code that I used. The code trains the neural network using your Nvidia GPU if you have one, else it uses the CPU.

There’s many parameters you can alter, including the number of options to use in fitting, number and size of the layers in the neural network, and number of epochs to train for.

An interesting exercise would be to try to modify the code to work with more complex options, such as barrier options or American options.

Python Code

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import TensorDataset, DataLoader

# Black-Scholes formula for European call options
def bs_price(S, K, T, r, sigma):
    from torch.distributions import Normal
    sqrtT = torch.sqrt(T)
    d1 = (torch.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * sqrtT)
    d2 = d1 - sigma * sqrtT
    norm = Normal(0., 1.)
    return S * norm.cdf(d1) - K * torch.exp(-r * T) * norm.cdf(d2)

# Generate dataset
N = 600_000
S = torch.rand(N, 1) * 100 + 1
K = torch.rand(N, 1) * 100 + 1
T = torch.rand(N, 1) * 2 + 0.01
r = torch.rand(N, 1) * 0.1
sigma = torch.rand(N, 1) * 0.5 + 0.01
prices = bs_price(S, K, T, r, sigma)

# Transform the data so S and T are one
K2 = K / S
T2 = torch.ones_like(K)
sigma2 = sigma * torch.sqrt(T)
DF = torch.exp(-r * T)
prices2 = prices / S

X = torch.cat([K2, T, r, sigma2, DF], dim=1)
y = prices2

train_size = int(0.8 * N)
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

batch_size = 512
dataset = TensorDataset(X_train, y_train)
loader = DataLoader(dataset, batch_size=batch_size, shuffle=True)

# Network configuration
hidden_size = 256
num_layers = 6

def build_model(input_dim, hidden_size, num_layers):
    layers = [nn.Linear(input_dim, hidden_size), nn.ReLU()]
    for _ in range(num_layers - 1):
        layers += [nn.Linear(hidden_size, hidden_size), nn.ReLU()]
    layers.append(nn.Linear(hidden_size, 1))
    return nn.Sequential(*layers)

class BSNet(nn.Module):
    def __init__(self):
        super(BSNet, self).__init__()
        self.model = build_model(5, hidden_size, num_layers)
    def forward(self, x):
        return self.model(x)

# Initialize
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
net = BSNet().to(device)
optimizer = optim.Adam(net.parameters(), lr=1e-3, weight_decay=1e-6)
criterion = nn.MSELoss()

# Training loop
epochs = 200
for epoch in range(1, epochs + 1):
    net.train()
    running_loss = 0.0
    for xb, yb in loader:
        xb, yb = xb.to(device), yb.to(device)
        preds = net(xb)
        loss = criterion(preds, yb)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        running_loss += loss.item() * xb.size(0)
    epoch_loss = running_loss / train_size
    if epoch % 50 == 0:
        print(f"Epoch {epoch}/{epochs}, Training MSE: {epoch_loss:.10f}")

def eval_test():
    net.eval()
    with torch.no_grad():
        X_device = X_test.to(device)
        preds = net(X_device)
        return criterion(preds, y_test.to(device)).item()

test_loss = eval_test()
print(f"\nTest MSE after {epochs} epochs: {test_loss:.10f}")

def price_option(S_val, K_val, T_val, r_val, sigma_val):
    net.eval()
    with torch.no_grad():
        valsbs = torch.tensor([[S_val, K_val, T_val, r_val, sigma_val]], device=device)
        vals = torch.tensor([[K_val/S_val, T_val, r_val, sigma_val*torch.sqrt(torch.tensor(T_val)), torch.exp(torch.tensor(-r_val*T_val))]], device=device)
        nn_p = net(vals).cpu().item()*S_val
        bs_p = bs_price(valsbs[:,0:1], valsbs[:,1:2], valsbs[:,2:3], valsbs[:,3:4], valsbs[:,4:5]).cpu().item()
    print(nn_p, bs_p)

print("\nTesting on some options:")
price_option(55, 60, 1.2, .03, .12)
price_option(340, 330, 1.4, .02, .15)
price_option(131, 131, 0.6, .02, .1)

Is measure theory really needed in finance?

Measure theory is an area of mathematics that was created to generalise the theory of integration to more exotic kinds of functions. It has applications in pure mathematics, particularly to the theory of differential equations – see, for example, Sobolev spaces. At some point, finance academics started using the language of measure theory when writing papers and textbooks about financial derivative pricing. However, the field of measure theory isn’t really necessary in finance because in almost all cases the functions that arise in derivative pricing in practice are nice, smooth functions that don’t require considerations from measure theory like sigma algebras. Couching derivative pricing in the language of measure theory simply makes it difficult for non-mathematicians to understand things which they would otherwise be able to understand.

The main concept from measure theory which is useful in derivative pricing is that of a “measure”. However, we only need a greatly simplified concept of it. In finance, a measure is basically just a smooth function which is used to weight or assign probabilities to different outcomes.

Futures convexity corrections

In this section, we discuss the convexity correction which accounts for the difference between forwards and futures prices. Later, we will identify another convexity correction that is required for late or early payment of interest (such as for swaps in arrears).

In this discussion, everything is with respect to the risk neutral measure, which takes into account investor’s risk preferences so that by definition expectations agree with market prices.

Consider modelling the short rate (instantaneous interest rate) \(r\) by a stochastic differential equation with time dependent drift \(\theta(t)\):

\[dr = \theta(t)dt + \sigma dz.\]

The forward rate \(F(t,T,T + \tau)\) between \(T\) and \(T+\tau\), as seen from time \(t\), is some function of the values of \(r(t)\) at each time \(t\). You’ll remember from calculus that when you want to find the rate of change of a function of a function, you use the chain rule, which generates an additional term. In Stochastic calculus, the chain rule is known as Ito’s Lemma, and it actually has three terms – it has an additional second order (or convexity) term that you don’t get in ordinary calculus. So, a function of a stochastic variable obeys its own stochastic differential equation which we can generate using Ito’s lemma. Because of the new third term, this new stochastic differential equation can have a drift term even if the original equation does not. In other words, even if \(\theta(t) = 0\), the forward rate equation can have a drift term, so that it’s expected value at time \(T\) is not the same as it’s value now at time \(t\). The precise equation for the expected future value is derived in this addendum to Hull.

We’ve seen that the forward rate is not a martingale under the risk neutral measure. In other words, it is not equal to its future expected value because the equation describing it has a drift term. By contrast, futures are martingales and are equal to their future expected values. The reason for this is simple. Because Futures are margined daily, their market value is by definition zero. So, since they cost nothing to buy, their expected future value, with respect to the risk neutral measure, must be zero as well. This means that the expected value tomorrow, or any future day, must be equal to today’s value.

This means that the equation derived in Hull for the expected future value of the forward, is also equal to the convexity correction between forwards and futures prices.

Change of numéraire

Changing the numéraire means changing the unit of measure used to value financial instruments. For example, changing to a different currency, or changing to the value of the same currency at some future time after it has appreciated due to the time value of money. We will refer to the wikipedia article on numeraires.

We consider the value of an asset \(S\) in terms of a numéraire \(M\). We assume the existence of a so-called risk neutral measure \(Q\) under which asset prices are martingales, that is the expected values at time \(T\) are the same as their current values:

\[\frac{S(t)}{M(t)} = E_Q \left[ \frac{S(T)}{M(T)} \right]\]

We define a new measure \(Q^N\) at time \(T\) by weighting all probabilities by the factor \(\frac{N(T)}{M(T)}\). Then we have the formula

\[E_{Q^N} \left[ \frac{S(T)}{N(T)} \right] = E_Q \left[\frac{N(T)}{M(T)} \frac{S(T)}{N(T)} \right] / E_Q \left[ \frac{N(T)}{M(T)} \right]\]

To understand this, note that when we switch back to calculating the expectation with respect to \(Q\), we have to add the weighting (or scaling) factor inside the expectation. In addition, because on the LHS everything is in terms of the numeraire \(N\), we have to have the division on the RHS so both sides are in the same numeraire (imagine a currency change). Note also that the expectation in the denominator is simply equal to the same expression evaluated at time \(t\), using the martingale property above.

The forward measure

Now we consider a particular example that we will need shortly. The time value of money means that a dollar today is worth more than a dollar in the future, leading to the use of discounting in pricing financial instruments. We let \(P(t,T)\) be the value at time \(t\) of a zero coupon bond which pays a cashflow of \(1\) at time \(T\) (effectively, this is the discount factor between \(t\) and \(T\)). Let’s choose a change of numeraire given by \(N(t) = P(t,T)\) and \(M(t) = P(t,T+\tau)\). Then we have

\[ \frac{N(t)}{M(t)} = \frac{1}{P(t,T,T+\tau)} = E_Q \left[ \frac{N(T)}{M(T)} \right], \]

where \(P(t,T,T+\tau)\) represents the forward rate as seen from time \(t\), and we have used the martingale property described above.

The new measure associated with this numeraire change is called the \(T + \tau\)-forward measure, and expectations with respect to this measure are denoted by \(E^{T + \tau}\).

Vanilla swap pricing and the forward measure

For a vanilla swap, the rate fixes at the start of each period, but the payment is settled at the end of each period.

Pricing a vanilla swap requires no complex modelling, as the valuation just involves computing the interest payments as if the future interest rates were equal to their current forward values. This is possible because the expected future value of the forward rate (with respect to the forward measure) is equal to its present value. In technical language, the forward rate is a “martingale”.

Consider a swap period starting at time \(T\) and ending at time \(T + \tau\), as seen from the current time \(t\). Let \(F(t,T,T + \tau)\) denote the forward interest rate over this period. The Martingale property means that

\[E^{T + \tau}\{F(T,T,T+\tau|\mathcal{F}_t\}=F(t,T,T+\tau).\]

The “filtration” \(\mathcal{F}_t\) simply indicates that the expectation is as seen from time \(t\), and we neglect it in what follows. The notation \(E^{T + \tau}\) means that the expectation is taken with respect to the (\(T + \tau\))-forward measure.

A clear explanation of this fact can be found in Brigo & Mecurio [1]. Rearranging the definition of a (simply compounded) forward rate slightly, we have

\[F(t,T,T + \tau)P(t,T + \tau) = \frac{1}{Y(T, T+\tau)}(P(t,T) – P(t,T+\tau)),\]

where \(Y\) represents the year fraction.

We now consider taking the expectation of both sides. Because of the discount factor on the LHS which scales all possible forward values by a constant factor, we are taking the expectation of the forward rate with respect to the (\(T + \tau\))-forward measure. On the other hand, the RHS is just the difference between the today-prices of two zero coupon bonds of different maturities. Since the latter are real assets trading in the market, their expected value must be equal to their current value. It follows that the same is true for the left hand side, i.e. the forward rate scaled by the discount factor.

Swap in arrears pricing

Swaps in arrears differ from vanilla swaps in that the both fixing and settlement occur at the start of each period. In this case, the expected future interest rate is no longer equal to the current forward value, necessitating a so-called convexity correction. The value of a swap in arrears is given by

\[V(t) = P(t,T) E^T \left[ F(t,T,T + \tau) \right].\]

For swaps in arrears, the expectation is taken at \(T\) instead of \(T + \tau\), so we can’t use the martingale property from the previous section. Using the change of numéraire formulae in a previous section, we can change to the (\(T + \tau\))-forward measure which introduces a discount factor inside the expectation as follows:

\[V(t) = P(t,T) E^{T+\tau} \left[ \frac{F(t,T,T + \tau)}{P(T,T+\tau)} \right] P(t,T,T+\tau)\]

\[= P(t,T+\tau) E^{T+\tau} \left[ \frac{F(t,T,T + \tau)}{P(T,T+\tau)} \right].\]

Because of the factor of \(1/P(T,T+\tau)\) inside the expectation, we can’t use the martingale property of forward rates here.

Attempting the quantify the magnitude of this correction is known as the “convexity correction”. It’s an amount that gets added to the forward rates before using them to price swaps, and calculating it requires that you first make a choice of stochastic interest rate model.

Why Forward Rate Agreements (FRAs) don’t require a convexity adjustment.

A FRA is similar to a single period swap, but not quite. A key difference is that the rate is fixed in advance, at the start of the accural period, just like a swap in arrears. However, unlike a swap in arrears, FRAs do not require a convexity adjustment, and one naturally wonders why. The answer lies in the second key difference between a FRA and a swap – for a FRA, the floating leg is discounted by the forward floating rate between the start and end dates of the accrual period. This discount factor, as seen from time \(t\), is

\[ P(t,T,T+\tau) = P(t,T+\tau)/P(t,T).\]

The value of a FRA is given by

\[V(t) = P(t,T) E^T \left[ F(t,T,T + \tau) P(T,T+\tau) \right],\]

which differs from a swap in arrears by the discount factor inside the expectation.

Note that if we now change to the (\(T + \tau\))-forward measure, as we did above for swaps in arrears, the new discount factor \(P(T,T+\tau)\) would exactly cancel with the change of measure \(1/P(T,T+\tau)\). Thus, FRAs don’t require a convexity correction.

Girsanov’s theorem

We’ve seen above how forwards are martingales under the forward measure, but have a drift term under the risk neutral measure. So it seems like one could get rid of a drift term by an appropriate measure, or conversely, change measure at the expense of adding a drift term. Girsanov’s theorem is the formalization of this idea.

[1] Brigo & Mecurio, Interest Rate Models – Theory and Practice, Springer Finance 2007

What are perpetual options?

Everyone familiar with crypto trading will have heard of perpetual futures. These products are defined by regular settlements, say once per day at times \(t_1, t_2, t_3,…\), where the holder pays the amount \(V(t,S(t)) – P(S(t))\). Here, \(P\) is the payoff of the future, given by \(P(S) = S\). The function \(V\) is the current price at which the perpetual is trading on the market.

Perpetual options come about by simply replacing the payoff P in the above construct with the payoff off an option. That is, \(Max(S-K, 0)\) for a call option and \(Max(K-S, 0)\) for a put option.

For each daily settlement, you can view \(P\) as the cashflow you receive, and \(V\) as the cashflow you have to pay. The function \(V\), known as the funding fee, is often viewed as a penalty which forces the value of the perpetual to converge towards the payoff as the time approaches each daily expiry. Otherwise, the holders of the perpetuals are penalized by having to pay more than what they receive. As we’ll see, another way to view it in the case of a perpetual option, is as the premium of a portfolio of options.

One of the arguments used to motivate perpetuals is that they prevent the available liquidity from being divided among options of many different expiries. Another is that they avoid the inconvenience and spread fees of having to manually roll contracts as each option expires.

Perpetual options were originally discussed in this short paper.

Perpetual options as a portfolio of options

Let’s just consider for a moment the second part of the settlement only. Since at each time \(t_i\) we receive the payoff \(P_i = P(S(t_i))\), we can view this collection of cashflows as exactly that from a portfolio of options with expiries \(t_1, t_2, t_3,…\).

It makes sense therefore that the payments \(V_i = V(t_i,S_i)\) can be conceived of as the premium payments corresponding to the option portfolio. Now pay attention to the following fact: since at each time \(t_i\) the holder of the perpetual pays the current value of the option portfolio, this means that the notionals of all remaining options must double.

about the payment  as the option premium payment, since the investor pays the current value of the portfolio, it’s clear that the notional of the options must double at each payment time (excluding the one that’s expiring at the same time, of course). Since the kth option notional will be doubled one more time than the (k-1)th option, each successive option in the portfolio will have half the notional than the last one. In other words, the notionals must follow a geometric series with ratio ½.

Continuous funding

The notional halfing each time you move to the option with an expiry one day longer means that the notionals obey an exponential relationship. Instead of using a base of one half, we can use a base of \(e\) as long as we introduce a constant out the front. Thus, assuming our perpetual is a call, we could write the price of our perpetual (option portfolio) as

\[V = A\sum_{i=1}^{\infty} e^{-t_i/T} C_{t_i},\]

for some constant \(A\), where \(C_{t_i}\) is the price of a call option with expiry \(t_i\), and \(T\) is the interval between option expiries (assumed to be one day in our discussion above). If we now imagine decreasing the spacing of the options so that the sum becomes an integral (while leaving \(T\) unchanged), we get what’s known as the continuous funding case:

\[V = A\int_{0}^{\infty} e^{-t/T} C_{t}dt\]

Pricing models for perpetual options

Since the product can be conceived of as a (discrete or continuous) collection of options, the problem of pricing reduces to the problem of pricing options of different maturities. Assuming a Black-Scholes framework, this reduces to the problem of  fitting an arbitrage free volatility surface to whatever implied volatilities exist in the market for vanilla options. Algorithms exist that can find the closest (least squares) fit subject to no arbitrage conditions. If there is limited data, a parametric approach could be considered, as these have far fewer degrees of freedom than a surface obtained by interpolating between options which exist in the market. At the more complex end of the scale, local and stochastic vol models are both parametric approaches but could be overkill.

As usual, volatilities for options whose strikes or expiries fall in between those in the market are generated by interpolating the volatility surface. For example, interpolating between time 0 and the first option can be done by assuming the volatility surface is linear in variance (\(V^2T\) interpolation).

Continuously funded everlasting options would need to be priced as a numerical integral, with options of different expiries being assigned different volatilities from the surface.

Variants where the strike is a moving average of recent prices can also be handled using the corresponding option pricing formulae for these products.

In other words, as long as you have the volatility data and model required to price regular options, perpetual prices follow without too much difficulty.

In terms of evaluating the integral above, the exponential weighting term means that we should only need to consider about 10 funding intervals (each of length \(T\)) before the contribution of additional terms becomes negligible. The integral can then be carried out by approximating it by a finite number of steps as usual.

How do the price and greeks of perpetual options compare to vanilla options?

The price/greeks of an everlasting are going to be a weighted sum of the prices/greeks of options of increasing expiries. It should be noted that the price of an everlasting will depend on how the volatility of an option depends on expiry according to the chosen vol surface. For simplicity, I’ve assumed the volatility surface is flat. If the funding period is short such as a single day, assuming volatility has no term structure may be reasonable.

To see how the price and greeks of a perpetual option relates to that of a single vanilla option with a single expiry, I’ve used python code to produce some graphs.

Focusing on calls, I’ve plotted the prices of an 8 hour, strike 100 call. I’ve set the volatility to 80% to accentuate the difference between the two prices which is otherwise hard to see.

The two prices are surprisingly close. The vanilla is slightly more valuable near at the money, and slightly less valuable in the wings. This is seen more clearly when plotting the difference between the two prices:

Since the price graphs are very similar, the Greeks will be very similar as well. We plot the delta difference:

The gamma difference:

And the vega difference:

The patterns are interesting, but their small magnitude wouldn’t have much impact on a standard market risk or margin amount (VaR) calculation, for which approximate values of greeks would be sufficient. Thus the risk profile is very similar to that of a dated option with maturity equal to the length of the funding period.

Hedging perpetual options

Viewing it as a portfolio of regular options, you can hedge it using a portfolio of whatever instruments you would normally use for the regular options. However, unless you only wish to hedge the nearest few expiries, this might require purchasing many instruments. Conceptually, purchasing spot in order to delta hedge could be accomplished by either combining the spot amounts required for each of the deltas of the individual options, or purchasing an amount of spot based on the delta of the everlasting option.

Alternatively, depending on the hedge, it might be possible to hedge all expiries simultaneously using another perpetual product, if it exists.

Arbitrage with vanilla options

In the standard options market one can look for various kinds of arbitrage, including violations of put call parity and inconsistencies with the volatility surface. This includes calendar arbitrage, where an option with closer maturity is assigned an implied volatility that leads it to be more valuable than an option with further maturity, and butterfly arbitrage, where the option price as a function of strike is not convex.

If the everlasting option is viewed as a finite portfolio of options, and one of those options is mispriced compared to an equivalent dated option with the same expiry, then one simply buys one and sells the other. In someone believes that the perpetual is mispriced at many expiries, they can buy/sell a portfolio of dated options alongside buying or selling the everlasting option.

Python code

Below is the python code used for this article.

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

T = 1/365/3
K = 100
r = 0
V = 0.8

def Call(V, S, K, T, r):
    if T > 0:
        d1 = (np.log(S/K) + (r + 0.5*V**2)*T)/(V*np.sqrt(T))
        d2 = (np.log(S/K) + (r - 0.5*V**2)*T)/(V*np.sqrt(T))
        call = S*norm.cdf(d1,0.0,1.0) - K*np.exp(-r*T)*norm.cdf(d2,0.0,1.0)
        return call
    else:
        return S - K

def ECall(V, S, K, T, r):
    
    u = np.sqrt(1 + 8/(V*V*T))
    
    intrinsic = max(S-K,0)
    
    if S >= K:
        timevalue = (K/u)*(S/K)**(-(u-1)/2)
    else:
        timevalue = (K/u)*(S/K)**((u+1)/2)
    
    return intrinsic + timevalue

def Delta(V, S, K, T, r, Pricefunction):
    step = 0.0001*S
    return(Pricefunction(V,S+step,K,T,r) - Pricefunction(V,S-step,K,T,r)) / (2*step)

def Gamma(V, S, K, T, r, Pricefunction):
    step = 0.0001*S
    return (Pricefunction(V,S+step,K,T,r) -2*Pricefunction(V,S,K,T,r) + Pricefunction(V,S-step,K,T,r))/(step**2)

def Vega(V, S, K, T, r, Pricefunction):
    step = 0.0001*V
    return (Pricefunction(V+step,S,K,T,r) - Pricefunction(V,S,K,T,r))/step

# Price plot
plt.figure(1)
plt.xlabel('Spot', fontsize=12)
plt.ylabel('Call Price', fontsize=12)
xlist = np.linspace(90,110,300)

ylistvan = [Call(V, S, K, T, r) for S in xlist]
plt.plot(xlist,ylistvan,label='Vanilla', linewidth=1)

ylistever = [ECall(V, S, K, T, r) for S in xlist]
plt.plot(xlist,ylistever,label='Everlasting', linewidth=1)

#plt.plot(xlist[:-1],ylist,label='Everlasting')
plt.legend(loc='upper right')

# Price difference plot
plt.figure(2)
plt.xlabel('Spot', fontsize=12)
plt.ylabel('Everlasting - vanilla price', fontsize=12)
ylistdiff = [ylistever[i] - ylistvan[i] for i in range(len(ylistever))]
plt.plot(xlist,ylistdiff,label='Everlasting - vanilla', linewidth=1)

# Delta difference plot
plt.figure(3)
plt.xlabel('Spot', fontsize=12)
plt.ylabel('Everlasting - vanilla delta', fontsize=12)
ylistdeltavan = [Delta(V, S, K, T, r, Call) for S in xlist]
ylistdeltaever = [Delta(V, S, K, T, r, ECall) for S in xlist]

ylistdeltadiff = [ylistdeltaever[i] - ylistdeltavan[i] for i in range(len(ylistever))]
plt.plot(xlist,ylistdeltadiff,label='Everlasting - vanilla', linewidth=1)

# Gamma difference plot
plt.figure(4)
plt.xlabel('Spot', fontsize=12)
plt.ylabel('Everlasting - vanilla gamma', fontsize=12)
ylistgammavan = [Gamma(V, S, K, T, r, Call) for S in xlist]
ylistgammaever = [Gamma(V, S, K, T, r, ECall) for S in xlist]

ylistdeltadiff = [ylistgammaever[i] - ylistgammavan[i] for i in range(len(ylistever))]
plt.plot(xlist,ylistdeltadiff,label='Everlasting - vanilla', linewidth=1)

# Vega difference plot
plt.figure(5)
plt.xlabel('Spot', fontsize=12)
plt.ylabel('Everlasting - vanilla vega', fontsize=12)
ylistvegavan = [Vega(V, S, K, T, r, Call) for S in xlist]
ylistvegaever = [Vega(V, S, K, T, r, ECall) for S in xlist]

ylistdeltadiff = [ylistvegaever[i] - ylistvegavan[i] for i in range(len(ylistever))]
plt.plot(xlist,ylistdeltadiff,label='Everlasting - vanilla', linewidth=1)

To mitigate credit risk in trading, one or both parties may be required to post both initial and variation margin. While initial margin is posted by both counterparties at the inception of the trade, variation margin is periodically exchanged (daily or even intra-daily) during the life of the trade due to changes in the mark to market value of the position.

The margin protects against a counterparty defaulting if the trade doesn’t move in their favour. In the case of centrally cleared derivatives, an intermediary manages the margining requirements of the counterparties. For OTC derivatives which are reasonably standardized, central clearing may still be a regulatory requirement. However, for OTC derivatives that are bespoke and illiquid, determining appropriate margining requirements may be more challenging due to the difficulty in regular valuation of the contracts. Margining requirements have become more stringent since the GFC.

While one does not want to be inadequately protected against non-payment by a counterparty, excessive margining also imposes costs on the derivatives markets. The immediate questions is, how does one estimate the amount of margin required by a counterparty to ensure that, at some chosen statistical confidence level, they can meet their obligations? This question is the subject of this article.

While margining is often a regulatory requirement, financial services firms including proprietary trading firms and brokers/exchanges, may wish to conduct their own margin or collateral modelling to ensure they are protected against default by counterparties due to sudden changes in valuation of their trades. An example of this in the context of cryptocurrency exchanges is provided later in the article.

Margin modelling methodologies

Approaches to margin modelling may have much in common with methodologies used for market risk calculations. In particular, value at risk (VAR) and expected shortfall methods. The VAR method assume a distribution type with which to model risk factor (market data) shifts (typically a normal distribution), calibrates the distributions to the market somehow (typically historical simulation, which analyses the shifts in market variables over some specified time interval), and then calculates the 99% quantile worst market move (or some other specified confidence level) of the derivative or portfolio valuation over some given time horizon. Expected shortfall methods are a slight modification on this approach which look at the mean of all outcomes above the 99% quantile.

More frequent margin calls will reduce the credit risk (due to the shorter time horizon available for the market to move), at the expense of more operational inconvenience.

These methods require careful choice of a historical window for calibration. If only the most recent data is used, rather than a historically stressed period, it may cause over leveraging during more stable market periods, leading to higher exposure when the markets become more tumultuous.

Another question that arises is how to calculate margin requirements for an entire portfolio instead of a single trade. A portfolio may have offsetting trades, and a decision must be made about what correlations to assume and what level of diversification benefit to allow.

The ISDA standard initial margin model (SIMM)

To avoid disputes over margining amounts, it’s important to have a single, transparent and standardized approach that all parties must adhere to. The International Swaps and Derivatives Association prescribes an approach known as the standard initial margin model. The ISDA documentation for this approach can be found here.

The ISDA includes four product classes: Interest rates and FX (ratesFX), credit, equity and commodity.

Six “risk classes” are specified. These are market data categories on which trade valuations depend. They are interest rate, credit (qualifying), credit (non qualifying), equity, commodity, and FX. Each trade depends on one or more “risk factors” from each category, such as interest rate curves, equity prices and exchange rates. The model postulates that the moves in any given market data or “risk factor” are normally distributed, which is also the approach commonly taken for market risk. The changes in trade valuations are also assumed to be normally distributed, something which is not actually true for nonlinear products like options even when the risk factor shifts are normal, and which becomes less true the larger the shifts are. The valuation shifts depend on the changes in the risk factors through trade sensitivities.

We first discuss how to calculate the initial margin for an individual trade. The approach is to calculate the delta, vega and interest rate sensitivities of the trade. The ISDA has calculated numbers representing how much they believe each risk factor will move at 99% confidence, which then multiply the sensitivities. The process by which the user finds the correct number is fairly intricate, and we won’t attempt to reproduce every detail here. The approach also takes into account the ISDA’s view on correlations between risk factors.

For options, a “curvature” term is also added. It is actually the vega contribution scaled by a “scaling factor”, which represent’s the ISDA’s view on how gamma and vega are related for vanilla options.

Having calculated the initial margin for an individual trade, we now discuss how to aggregate trades to calculate the initial margin of the entire portfolio. A fantastic simplifying feature of normal distributions, is that the sum of two normal distributions is another normal distribution with some new standard deviation. The new standard deviation is a function of the old ones, calculated using what I will call the “square root rule”: \(\)

\[ \sigma_{agg} = \sqrt{\sum_i \sigma_i + \sum_r \sum_{s \neq r} \psi_{rs} \sigma_r \sigma_s},\]

where \(\psi\) is a matrix of correlations between the distributions. One can aggregate the normally distributed valuation shifts for different risk factors for any given trade, giving a total normal distribution for each trade. One can then aggregate all these normal distributions across all trades in the same way, yielding a new normal distribution corresponding to the total valuation shift in the whole portfolio..

The ISDA methodology has three differences from the procedure described above. Firstly, it applies the square root rule to the initial margins, rather than the standard deviations of the distributions as described above. Since an initial margin is just some quantile of a normal distribution distribution, it is proportional to the standard deviation, so that the square root rule applies to the initial margins just as like it does to the standard deviations.

Secondly, It’s easy to verify that for the square root rule, the order in which we aggregate all the normal distributions doesn’t matter. The ISDA approach aggregates the normal distributions in a specific order. For each product class, and within each risk class \(X\) within the product class, it first aggregates the delta margin across all trades, and similarly the vega margin, curvature margin and base correlation margin. It then aggregates these four distributions as a simple sum:

\[ IM_X = DeltaMargin_X + VegaMargin_X + CurvatureMargin_X + BaseCorrMargin_X \]

This is the third difference. This sum is more conservative than that obtained by using the square root rule, and assumes that there is no diversification benefit between these four types of moves. The ISDA may believe that in stressed conditions, these different classes of risk factors may move together.

The ISDA then uses the square root rule to aggregate these initial margins across all risk classes within the product class. Just like with the four margin types, the ISDA stipulates that no diversification benefit may be assumed between the product classes, so that the total amount of initial margin required is the sum of those for each of the four product classes:

\[ SIMM_{Product} = SIMM_{RatesFX}+SIMM_{Credit}+SIMM_{Equity}+SIMM_{Commodity}\]

That is, above the product class level, we don’t apply the square root rule but a straight sum. I’m not sure whether this is based in part on some analysis the ISDA has done concerning correlations between product classes, or is simply conservatism.

In any event, we’ve now calculated the total initial margin across our entire portfolio, as stipulated by the ISDA’s standard initial margin model.

FX sensitivities

Consider an FX trade with underlying \(S = CCY1CCY2\) and present value \(PV(S)\) expressed in CCY2. When calculating sensitivities, SIMM considers a 1% shift \(dS=.01S\). The sensitivity with respect to \(CCY1\) is defined to be the corresponding shift

\[sen_{CCY1} = dPV.\]

The sensitivity with respect to \(CCY2\) is defined to be

\[sen_{CCY1} = -Sd(PV/S).\]

In other words, we convert \(PV\) to \(CCY1\), find the shift, and then convert back. Using the chain rule for differentiation, it’s easy to show that

\[sen_{CCY1} + sen_{CCY2} = .01PV.\]

Once one sensitivity is calculated, the other can be easily found using this formula.

SIMM contribution of option premium settlement

For options, the option premium is typically settled two business days after the trade date. When a SIMM calculation is performed for an option trade before the option premium has settled, the premium cashflow is included in the calculation of the sensitivities. Since the option premium is equal and opposite to the usual PV of the trade, this means that the “total PV” including the premium payment is zero at the inception of the trade. Using the above formula, this means that the two sensitivities are equal and opposite.

The contribution that the premium payment makes to \(sen_{CCY1}\) and \(sen_{CCY2}\) depends on which currency settlement occurs in. Either way, since it is a fixed and known payment, its sensitivity to one of the currencies will be zero, and its sensitivity to the other currency will be \(.01Premium\).

Bitmex’s insurance fund

Given the high volatilities of cryptocurrencies, market risk and margin modelling is particularly challenging. When some major cryptocurrency exchanges combine this with very high leverages, managing their risk becomes particularly challenging.

On a traditional exchange, if the loss on a trade exceeds the posted margin, the exchange will insist that the trader put forward additional margin amount to cover the loss, and chase the debt if necessary.

Bitmex have taken a margining approach different to that used in traditional finance. They call this their insurance fund. In essence it means that Bitmex will liquidate your position at half margin, and appropriate the amount remaining after liquidating your position at the best available price. This money is then used to compensate for other trades by other participants which lost more than the margin amount. This means that traders which lost more than their margin are compensated, at the expense of traders who didn’t lose all of their margin.

To new participants, it seems surprising that the exchange could “steal” half your margin and keep it. However, Bitmex point out that, on the upside, traders are guaranteed to never lose more than their margin amount.

A curious feature of quantitative finance is the way pricing models are calibrated to the market. For bonds, one adds a Z-spread to the rates which matches market quotes. The Z-spread allows the model to incorporate the market’s views about credit default risk and future interest rates. For options priced using the Black-Scholes model, one prices using the implied volatilities that reproduce market prices. In this way, pricing models are really just interpolation tools to interpolate between the prices of products trading in the market in order to price different contracts with different parameters which are not currently represented in the market.

For options priced using the constant volatility Black-Scholes model, this creates the strange situation where options with different strikes and expiries are imagined to be based on underlying assets with different volatilities – even though they all have the same underlying which can only realise a single volatility.

Another issue concerns the pricing of more exotic derivatives like Asian options and barrier options. When adding a barrier to a European option, it’s no longer clear whether it’s appropriate to use implied volatilities backed out from European options. Unfortunately, barrier options may not themselves be sufficiently liquid in the market to generate their own implied volatilities (such an implied volatility surface would now have three dimensions – strike, expiry and barrier level). Furthermore, unlike vanilla options, barrier option valuation depends on volatility term structure, so that replacing a variable volatility with it’s constant average value isn’t appropriate. Similarly, for Asian options, valuation depends on the path taken by the underlying and by its volatility, not just on their average values.

We know that volatilities are not, in reality, constant. And the assumption that they are leads to difficulties in reproducing the market prices of European options, and meaningfully pricing less liquid exotic derivatives. Furthermore, empirical analysis of stock data may indicate that volatility has, statistically speaking, a dependence on the current underlying level. This naturally leads us to consider various volatility option pricing models where the volatility depends on both time and the underlying level.

Local Volatility Models

In the local volatility model (LV), the constant volatility of Black-Scholes is replaced by an instantaneous volatility \(\sigma(S_t,t)\) which depends on both the current underlying value and time. The underlying evolves under the stochastic differential equation

\[dS_t = (r_t – d_t)S_t dt + \sigma(S_t,t) S_t dW_t,\]

where \(r\) is the risk free rate, \(d\) is the dividend yield, and \(W_t\) is a Wiener process (i.e. a brownian motion or random walk).

As long as the volatility surface is arbitrage free, local volatility models can be calibrated to exactly reproduce the market option prices and volatility smile at each pillar maturity trading in the market. Calibration can be done computationally efficiently using Dupire’s formula. Dupire’s formula gives the volatility \(\sigma(K,T)\) at some strike \(K\) and expiry \(T\), in terms of the price \(C(K,T)\) of a call option. The formula is

\[\sigma(K,T)^2 = \frac{\frac{\partial C}{\partial K} – (r-q)(C – K\frac{\partial C}{\partial K})}{\frac{1}{2}K^2\frac{\partial^2 C}{\partial K^2}},\]

where \(r\) is the risk free rate and \(q\) is the dividend yield. Slightly confusingly, this function seems to define the local volatility function \(\sigma(S,t)\) at a strike \(K\) instead of at some underlying level \(S\). However, the assumption is simply that, in order to calculate the local volatility at \(S\), you set \(K=S\) in this formula.

Thus, the problem of calibrating the model to market prices (i.e. obtaining the function \(\sigma(S,t)\)), is reduced to the problem of generating the function of call prices \(C(K,T)\), for every strike and expiry. This of course must be done starting with market prices for only a finite number of strikes and expiries. There are a variety of methods proposed for constructing an arbitrage free surface of call prices, and we refer to Rasmussen for a survey.

Once the local volatility function has been obtained, one proceeds to price derivatives numerically. For example, one can use the Monte Carlo method to generate a large number of possible paths using the stochastic differential equation above. Alternatively, one might use PDE finite differencing methods to price the derivative.

A significant downside of local vol models is that \(dV/dS\) has the opposite from that observed in the market. When spot increases, a local vol model predicts that volatility should decrease, and vice versa. Since the observed market behaviour is the opposite, this makes the model unsuitable for hedging. Also, smile flattens with maturity.

Another issue which arises when pricing certain types of derivatives, is the flattening of the forward smile. Local volatility models are fitted to the prices of European options. It’s known that the vol smile of European options becomes flatter for longer maturities (in other words, the implied vol is more constant, having less dependence on strike). This means that when you consider the forward volatility between \(T_1\) and \(T_2\), you find that it is also quite flat. This is not what’s observed in practice, where volatility smiles in the future are not in general less flat than those observed today.

To address these issues, we have stochastic volatility models.

Stochastic Volatility Models

Local volatility models are deterministic in that, given future values for \(t\) and \(S_t\), we know the volatility will be precisely that prescribed by the local volatility function \(\sigma(S,T)\). Stochastic volatility (SV) models allow the volatility to evolve according to its own stochastic differential equation with a random term, just like the underlying \(S\). Stochastic volatility models are of the form

\[dS_t = (r_t – d_t)S_t dt + \sigma_t S_t dW_t,\]

\[d \sigma_t^2 = \alpha(\sigma_t^2,t) dt + \beta(\sigma_t^2,t) dB_t,\]

where \(B\) is a second Wiener process, correlated to \(W\) with some correlation \(\rho\).

Stochastic volatility models are typically more challenging to calibrate to market data than a local volatility model. With few exceptions (such as the Hagen formula for the SABR model), stochastic volatility models are solved numerically, and calibrated to the market numerically. Standard optimization routines can be used to solve the model, check the level of agreement with market prices, adjust the model parameters, and repeat until a sufficiently close fit is obtained. For example, a least squares objective function may be used to minimise the difference between the prices generated by the model and those seen in the market.

A prominent example of a stochastic volatility model which is the market standard for vanilla interest rate derivatives, see our main article on the SABR model. Another stochastic volatility model, which is commonly used for equities, is the Heston model, which we discuss in the next section.

Heston Model

The Heston model is a particular example of a stochastic volatility model which is popular for equity derivatives. The second equation above takes the form

\[d \sigma_t^2 = \theta (\omega – \sigma_t^2) dt + \varepsilon \sigma_t dB_t,\]

where \(\theta\), \(\omega\) and \(\varepsilon\) are constants. Note that the first term causes the variance to tend to mean revert to the value \(\omega\).

Local Stochastic Volatility Models

Finally, we mention models which are a kind of hybrid between local volatility models and stochastic volatility models. These are sometimes known as “local stochastic volatililty” or LSV models. They aim to incorporate the best characteristics of both models. See for example the following papers:

• Calibrating and pricing with embedded local volatility models
• Option Pricing under Hybrid Stochastic and Local Volatility

The SABR model, introduced by Hagen et al, is a stochastic volatility model which is commonly used for valuing interest rate derivatives. It models the time evolution of some forward rate \(F_t\), of some maturity \(T\) and tenor \(L\), as follows:

\[ dF_t = \alpha_t F_t^\beta dW_t,\]

\[ d\alpha_t = v \alpha_t dZ_t,\]

where \(W_t\) and \(Z_t\) are two Wiener processes with constant correlation \(\rho\).

Unlike local volatility models, the model reproduces the correct smile dynamics, so it can be used for hedging. The SABR model is popular because it is simpler than some other stochastic volatility models, and allows an approximate pricing formula, in terms of the SABR model parameters, which gives the price in terms of the implied Black volatility. That is, the implied volatility that when inserted into Black’s model reproduces the price generated by the SABR model. This makes the model much easier to calibrate, and in general it fits the market well. The Hagen formula is straight forward but long, so we don’t reproduce it here.

The SABR model quickly became the market standard for vanilla interest rate derivatives like caps, floors and swaptions.

The first equation above shows that the forward \(F_t\) will remain positive under the evolution of the equation. Thus, the SABR model needs to be modified if zero and negative rates are considered possible.

Fitting the model

The SABR model is fitted for a single forward rate, with some expiry T and some tenor L. As mentioned already, a direct formula, in terms of the SABR model parameters, gives you a good approximation for the Black implied vol for any given strike. Therefore, using a set of strikes for which there exist market prices (and therefore implied vols), one can do a least squares fitting which minimises the deviation of the SABR approximate Black vols from the market implied vols. That is, one wishes to find SABR parameters which minimise

\[ \sum_{i} (\sigma_{K_i}^{SABR} – \sigma_{K_i}^{BS}),\]

where the sum is over the set of strikes for a given expiry. The SABR volatilities are computed from the SABR parameters using the Hagen implied volatility formula, avoiding the need to numerically simulate the SABR model. Since \(\beta\) and \(\rho\) both affect the slope of the volatility smile, it’s common to fix \(\beta\) according to some view of the market, and fit only the remaining three parameters using this method.

This procedure can be repeated for all required expiries and tenors for which there are market prices.

Finally, if one wishes to price an option on a forward rate of an expiry and tenor for which there is no analogue in the market, one can construct a set of SABR parameters by interpolating between the SABR parameters of nearby expiries and tenors. Alternatively, one can interpolate the implied volatilities directly.

Swaptions and caps/floors

For swaptions, one uses the swap rate as the underlying. That is, the rate that makes the swap value equal to zero. For caps and floors, one must first do a bootstrapping process to obtain the individual caplet/floorlet vols from the trading caps and floors.

The parameters alpha, beta and rho.

The variable \(\alpha_t\) represents the volatility of the forward process. The first SABR parameter \(\alpha\) is a non-negative parameter which is defined as the initial value of \(\alpha_t\) at time \(t=0\). It’s effect is approximately a parallel shift of the volatility smile vertically.

The parameter \(\beta\) satisfies \(0 \leq \beta \leq 1\). As \(\beta\) increases, the slope of the volatility smile decreases. Note that for \(\beta\) close to zero, the SABR model is close to a normal, rather than lognormal, model.

As mentioned already, \(\rho\) represents the correlation between changes in the underlying and changes in it’s volatility. Like \(\beta\), the effect of \(\rho\) is to change the slope of the volatility smile.

The parameter \(v\) is the volatility of the volatility. A higher value for \(v\) gives a greater curvature to the volatility smile.

Limitations of the SABR model

The SABR model has some limitations worth noting. Firstly, as it only models a single forward rate, it should not be used to price derivatives whose valuation depends on more than one forward rate. This is because forward rates with nearby maturities and tenors are not independent. However, note that swaptions are nonetheless commonly priced using the SABR model by modelling the swap rate, instead of the constituent swaplets. Secondly, the equation lacks a mean reversion term, limiting its ability to accurately model interest rates which tend to be eventually pulled back to normal values. Finally, the convenient Black implied vol formula can become inaccurate under certain circumstances.

For pricing non-vanilla derivatives which depend on multiple forward rates (of various expiries and tenors), the SABR model must be modified so that forward rates of nearby tenors and expiries are not independent. See the SABR libor market model.

This is some code I wrote to generate a volatility surface in C++ using options and futures data, and use it to calculate price and Greeks for an ATM option. To simulate live trading, the code updates the volatility surface and option price as new data becomes available.

Note that the code is currently configured to execute for the first 150,000 lines of market data. This can be changed on by altering max_remaining_lines on line 62

Simplifications

Currently, the code creates only a single volatility smile for a single maturity. Of course, the code can quickly be adapted to generate a whole volatility surface, as this is just a number of smiles at different expiries with a suitable method of interpolating between them.

I used the future as a proxy for the forward, just based on the data I had. This is not strictly correct, as they can differ by a convexity correction. Also, rather than incorporate interest rate data, interest rates are currently assumed to be zero. Since the expiry used was short, this was acceptable.

Calculating Implied Volatilities

Implied volatilities are calculated using Newton’s method.

I’ve found that Newton’s method can fail to converge when converting volatility pillars from moneyness to delta, particularly for problematic vol surface data such as data containing arbitrage. However, these issues may not arise for a single implied volatility calculation. But convergence guarantees, and appropriate behaviour for all input data, should be further investigated for live trading code.

It’s also worth checking out our article on removing arbitrage from volatility surfaces.

Volatility Surface Construction

I took ATM to be the strike equal to the Futures value at expiry (that is, ATMF). Another possible convention is to use the strike for which a call and put have a delta of the same absolute magnitude.

For simplicity, I adopted the sticky strike convention to avoid having to convert volatility surfaces to delta. That is, the pillars of the vol surface are fixed strikes. I also adopted the common market convention of parameterising the vol surface using calls for strikes above ATM, and puts for strikes below ATM.

While cubic spline is market standard, this would have required either using a third party library or implementing it manually. Given my time requirements for this project, I chose linear interpolation between strike pillars. Flat extrapolation was chosen before the first and after the last pillar.

Volatility interpolation will not proceed if one of the two required pillar vols is missing. If the additional complexity was worthwhile, one could proceed to use the next pillar over if available. One could also attempt to use put in place of call (or vice versa) if the intended market price was not available. Note that the volatility surface output displays “-1” when volatility data is not yet available. As long as a valid volatility has been calculated at least once in the past, this will continue to be used. A possible modification would involve removing volatilities based on market prices that are either too old or no longer existing in the market.

Data Quality

The code uses the mid price if both bid and ask quantities are non-zero. If only one is non-zero, it will use that value. If both quantities are zero, the data is ignored and the volatility surface is not updated. That is, the row is interpreted as missing data.

Some data rows in the options input data seemed erroneous. In this case, Newton’s method could not find a volatility solution. I imposed a max_iterations of 100, and discarded such rows. Similarly, a volatility > 10 (1000%) was assumed to be due to erroneous data.

Greeks

The greeks delta, gamma, vega and theta are calculated using finite differences and a reasonable choice of step size. Although infinitesimal greeks are available via the Black-Scholes formula, traders may be more interested in the consequences of small, finite movements rather than infinitesimal movements. The greeks can be quite sensitive to the chosen step size.

Because I used Black’s formula for pricing which takes the forward as input, and spot data was not available, I outputed forward delta and forward gamma (that is, I shifted the forward value instead of the spot value). Delta will differ from spot by forward/spot and gamma by that amount squared. Theta is chosen to be a 1-day theta.

The decision was made to output updated option price and greeks whenever option call price changed by a tolerance of at least 0.00001. This condition could be modified depending on what exactly the reason is for outputting changed prices.

Code Structure

The primary auxillary functions are:

• “Black” to price calls and puts
• “Implied Vol” to use Newton’s method to calculate a volatility from a market price
• “PriceAndGreeks” to return a 5-vector price, delta, gamma, vega, theta.
• “InterpolateVol” to use linear interpolation and flat extrapolation to extract a volatility from the vol surface data.
• “UpdateVol” to update a volatility at a particular strike if either the corresponding market price has changed, or the forward has changed.

Presently, the remainder of the code, including file input/output and the primary loop through market data lines, occurs in the “main” function. Given more time, the main function could be broken down into multiple functions for cleaner, more reader and more reusable code (although the value of this depends on the ultimate purpose and future direction of the project).

Code Testing

The code begins with a few simple unit tests to test pricing, calculation of implied volatilities and volatility interpolation. Basic reasonableness examinations have been done on the output data, however more extensive testing should be done before putting such code into production.

C++ code

#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <limits>

using namespace std;

vector<vector<string>> getData(string fileName, string delimiter);
vector<string> getDataSingleLine(string fileName, string delimiter);
void writeData(vector<vector<string>> Z, string fileName, string delimiter);
vector<string> split(const string& str, const string& delim);
double Black(double V, double F, double K, double T, double r, char callput);
double Ncdf(double x);
double ImpliedVol(double initialVolGuess, double targetPrice, double tol, double F, double K, double T, double r, char callput);
vector<double> PriceAndGreeks(double V, double F, double K, double T, double r, char callput);
double interpolateVol(double strike, double F, vector<double> CallStrikes, vector<double> CallVols, vector<double> PutVols);
void UpdateVol(int updated_symbol, double FuturePrice, double T, vector<double>& CallStrikes, vector<double>& CallVols, vector<double>& PutStrikes, vector<double>& PutVols, vector<double>& CallPrices, vector<double>& PutPrices);

int main()
{
    // Unit tests
    cout << "Unit tests" << '\n';
    cout << "Call price should be 4.26894926: " << Black(0.1,100,99,1,0.05,'C') << '\n';
    cout << "Call implied vol should be 0.1: " << ImpliedVol(0.05, 4.26894926, 0.00001, 100, 99, 1, 0.05, 'C') << '\n';
    cout << "Put price should be 3.31771983: " << Black(0.1,100,99,1,0.05,'P') << '\n';
    cout << "Put implied vol should be 0.1: " << ImpliedVol(0.15, 3.31771983, 0.00001, 100, 99, 1, 0.05, 'P') << '\n';
    vector<double> teststrikes{1,2,3,4,5};
    vector<double> test_vols_call{.10,.12,.14,.16,.18};
    vector<double> test_vols_put{.102,.122,.142,.162,.182};
    cout << "Interpolated vol should be .102: " << interpolateVol(0.9, 3, teststrikes, test_vols_call, test_vols_put) << '\n';
    cout << "Interpolated vol should be .18: " << interpolateVol(5.1, 3, teststrikes, test_vols_call, test_vols_put) << '\n';
    cout << "Interpolated vol should be .131: " << interpolateVol(2.5, 3, teststrikes, test_vols_call, test_vols_put) << '\n';

    // Load data
    string market_data_filename = "market_data.csv";
    string static_info_filename = "static_info.csv";

    // Output files
    ofstream fittingfile("fitting_output.csv");
    ofstream optionfile("option_output.csv");

    // Format static data.
    // This code assumes order of headings and that calls are labelled 1 to 145 and puts 146 to 290 with none missing.
    // Code would need to be modified to handle more general input files or potential missing data.
    vector<vector<string>> static_info = getData(static_info_filename, ",");
    double expiry_time = stod(static_info[1][5]); // Assumes all options have the same expiry

    vector<double> CallStrikes(145);
    vector<double> PutStrikes(145);

    for(int i = 1; i < static_info.size(); i++)
    {
        string symbol = static_info[i][0].substr(4, symbol.size() - 4);
        if (static_info[i][3].compare("C") == 0) {CallStrikes[stoi(symbol)-1] = stod(static_info[i][4]);}
        else if (static_info[i][3].compare("P") == 0) {PutStrikes[stoi(symbol)-146] = stod(static_info[i][4]);}
    }

    // Load market data
    // To simulate the arrival of live trading data, data is loaded one line at a time, with volatility surface calculations updated every line.

    ifstream file(market_data_filename);
    int max_remaining_lines = 150000; // Used to run code for a limited number of rows only
	string line = "";

	double current_time;
	int updated_symbol; // symbol that is updated this data row
	double T; // time to expiry
	getline(file, line); // First row are headings

	// Most recent price update for each strike. -1 represents no data yet.
    vector<double> CallPrices(145,-1);
    vector<double> PutPrices(145,-1);
    vector<double> CallVols(145,-1);
    vector<double> PutVols(145,-1);
    vector<double> OptionCallPrice(5,-1);

    double FuturePrice = -1;

    // Output files
    vector<string> heading_row;
    heading_row.push_back("Timestamp");
    for(int i = 1; i < CallStrikes.size(); i++)
    {
        heading_row.push_back(to_string(CallStrikes[i])); // Assuming call strikes and put strikes are the same, which seems to be the case
    }
    for(int l = 0; l < heading_row.size(); l++)
    {
        fittingfile << heading_row[l] << ",";
    }
    fittingfile << "\n";

    vector<string> heading_row2 {"Timestamp","Vol","Call price","Call delta","Call gamma","Call vega","Call theta","Put price","Put delta","Put gamma","Put vega","Put theta"};
    for(int l = 0; l < heading_row2.size(); l++)
    {
        optionfile << heading_row2[l] << ",";
    }
    optionfile << "\n";

    vector<vector<string>> pricing_output;

    // Main loop. Reads each line in market data input and updates futures price, vol surface and option price/greeks where appropriate.
	while (getline(file, line) && max_remaining_lines > 0)
	{
		vector<string> column = split(line, ",");

		max_remaining_lines = max_remaining_lines - 1;

		current_time = stod(column[1]);
        T = (expiry_time - current_time)*pow(10,-9)/31536000; // convert nanoseconds to years

        updated_symbol = stoi(column[0].substr(4, column[0].size() - 4));

        // Code to find most appropriate price to use for current symbol
        double current_bid = stod(column[2]);
        int bid_quantity = stod(column[3]);
        double current_ask = stod(column[4]);
        int ask_quantity = stod(column[5]);
        if (bid_quantity == 0 && ask_quantity == 0) {continue;} // Interpreted as no data update
        if (current_bid == 0 && current_ask == 0) {continue;} // Interpreted as no data update
        double mid = 0;
        // if one has quantity zero, use the other
        if (bid_quantity == 0 and ask_quantity > 0){mid = current_ask;}
        else if (bid_quantity > 0 and ask_quantity == 0){mid = mid + current_bid;}
        else {mid = (current_bid + current_ask)/2;}

        // Update prices
        if (updated_symbol == 0) {FuturePrice = mid;}
        else if (updated_symbol <= 145)
        {
            CallPrices[updated_symbol-1] = mid;
        }
        else
        {
            PutPrices[updated_symbol-146] = mid;
        }

        // Update volatilities
        if (FuturePrice < 0) {continue;} // unable to calculate volatilities without a forward

        if (updated_symbol == 0) // If Future price has changed, all vols must be updated
        {
             for(int l = 1; l < CallVols.size()+ PutVols.size() + 1; l++)
             {
                 UpdateVol(l, FuturePrice, T, CallStrikes, CallVols, PutStrikes, PutVols, CallPrices, PutPrices);
             }
        }
        else{UpdateVol(updated_symbol, FuturePrice, T, CallStrikes, CallVols, PutStrikes, PutVols, CallPrices, PutPrices);}

        // If V or F have changed, output ATM option price and greeks.
        double impliedATMvol = interpolateVol(FuturePrice, FuturePrice, CallStrikes, CallVols, PutVols);

        // If required vol pillars are still missing from vol surface data, no option price can be output
        if (impliedATMvol > 0)
        {
            double oldprice = OptionCallPrice[0];
            OptionCallPrice = PriceAndGreeks(impliedATMvol, FuturePrice, FuturePrice, T, 0, 'C');
            if(abs(OptionCallPrice[0] - oldprice) > 0.00001) // only write a row to file if price has changed.
            {
                vector<double> OptionPutPrice = PriceAndGreeks(impliedATMvol, FuturePrice, FuturePrice, T, 0, 'P');
                optionfile << to_string(current_time) << ",";
                optionfile << to_string(impliedATMvol) << ",";
                for(int l = 0; l < OptionCallPrice.size(); l++)
                {
                    optionfile << to_string(OptionCallPrice[l]) << ",";
                }
                for(int l = 0; l < OptionPutPrice.size(); l++)
                {
                    optionfile << to_string(OptionPutPrice[l]) << ",";
                }
                optionfile << "\n";
            }
        }

        // Output fitting csv
        // epoch nanoseconds, fitted volatility at each strike
        vector<string> fitting_row;
        fitting_row.push_back(to_string(current_time));
        for(int i = 1; i < CallVols.size(); i++)
        {
            if (CallStrikes[i] < FuturePrice) {fitting_row.push_back(to_string(PutVols[i]));}
            else {fitting_row.push_back(to_string(CallVols[i]));}
        }

        for(int l = 0; l < fitting_row.size(); l++)
        {
            fittingfile << fitting_row[l] << ",";
        }
        fittingfile << "\n";

    }

	file.close();
	fittingfile.close();
	optionfile.close();

    return 0;
}

// Update the vol at a given strike
void UpdateVol(int updated_symbol, double FuturePrice, double T, vector<double>& CallStrikes, vector<double>& CallVols, vector<double>& PutStrikes, vector<double>& PutVols, vector<double>& CallPrices, vector<double>& PutPrices)
{
        double tolerance = 0.0000001;
        double r = 0; // Assumption of zero interest rates
        double initialVolGuess;

        if (updated_symbol <= 145)
        {
            if(CallPrices[updated_symbol-1] < 0) {return;}

            double mid = CallPrices[updated_symbol-1];
            if (CallVols[updated_symbol-1] < 0) {initialVolGuess = 0.15;}
            else {initialVolGuess = CallVols[updated_symbol-1];}
            double impvol = ImpliedVol(initialVolGuess, mid, tolerance, FuturePrice, CallStrikes[updated_symbol-1], T, 0, 'C');
            if (impvol > 0 && impvol < 10)
            {
                CallVols[updated_symbol-1] = impvol; // ImpliedVol returns -1 if Newton fails to converge
            }
        }
        else
        {
            if(PutPrices[updated_symbol-146] < 0) {return;}

            double mid = PutPrices[updated_symbol-1];
            if (PutVols[updated_symbol-146] < 0) {initialVolGuess = 0.15;}
            else {initialVolGuess = PutVols[updated_symbol-146];}
            double impvol = ImpliedVol(initialVolGuess, mid, tolerance, FuturePrice, PutStrikes[updated_symbol-146], T, 0, 'P');
            if (impvol > 0 && impvol < 10)
            {
                PutVols[updated_symbol-146] = impvol;
            }
        }

}

// Linearly interpolate the vol smile for some strike
double interpolateVol(double strike, double F, vector<double> CallStrikes, vector<double> CallVols, vector<double> PutVols)
{
    // Assume callstrike list is identical to putstrike list, and callstrike list is ordered

    // Flat extrapolation
    if (strike <= CallStrikes[0]) {return PutVols[0];}
    if (strike >= CallStrikes[CallStrikes.size() - 1]) {return CallVols[CallVols.size() - 1];}

    int rightindex = 0;
    for(int i = 0; i < CallStrikes.size(); i++)
    {
        if (strike < CallStrikes[i]){rightindex = i; break;}
    }

    double right_strike = CallStrikes[rightindex];
    double left_strike = CallStrikes[rightindex - 1];
    double right_vol;
    double left_vol;

    // Use put for strike below ATM and call otherwise
    if (left_strike < F) {left_vol = PutVols[rightindex - 1];}
    else {left_vol = CallVols[rightindex - 1];}

    // If vols don't exist yet return -1
    if (left_vol < 0 or right_vol < 0) {return -1;}

    if (right_strike < F) {right_vol = PutVols[rightindex];}
    else {right_vol = CallVols[rightindex];}

    return left_vol + (right_vol - left_vol)*(strike - left_strike)/(right_strike - left_strike);
}

vector<double> PriceAndGreeks(double V, double F, double K, double T, double r, char callput)
{
    double deltastep = 0.0001*F;
    double gammastep = 0.01*F;
    double vegastep = 0.0001*V;
    double thetastep = min(1.0/365,T);

    vector<double> result(5,0);
    double price = Black(V, F, K, T, r, callput);
    result[0] = price;
    result[1] = (Black(V, F + deltastep/2, K, T, r, callput) - Black(V, F - deltastep/2, K, T, r, callput))/(deltastep/2);
    result[2] = (Black(V, F + gammastep, K, T, r, callput) + Black(V, F - gammastep, K, T, r, callput) - 2*price) / pow(gammastep, 2);
    result[3] = (Black(V + vegastep/2, F, K, T, r, callput) - Black(V - vegastep/2, F, K, T, r, callput))/(vegastep/2);
    result[4] = (Black(V, F, K, T - thetastep, r, callput) - price)/thetastep;

    return result; // price, delta, gamma, vega, theta
}

// Uses Newton's method to calculate implied volatility from market price
double ImpliedVol(double initialVolGuess, double targetPrice, double tol, double F, double K, double T, double r, char callput)
{
    double derivativeStep = min(0.00001, tol/100);
    double currentVol = initialVolGuess;
    double functionValue = Black(currentVol, F, K, T, r, callput) - targetPrice;
    double derivativeValue;
    int max_iterations = 100;

    while (abs(functionValue) > tol){
        derivativeValue = (Black(currentVol + derivativeStep, F, K, T, r, callput) - functionValue - targetPrice)/derivativeStep;
        currentVol = currentVol - functionValue / derivativeValue;
        functionValue = Black(currentVol, F, K, T, r, callput) - targetPrice;
        max_iterations = max_iterations - 1;
        if (max_iterations == 0) {return -1;}
    }

    return currentVol;
}

// Prices a put or call option using Black's formula and the forward value
double Black(double V, double F, double K, double T, double r, char callput)
{
    constexpr double lowest_double = std::numeric_limits<double>::lowest();
    if (T < lowest_double) {return F-K;} // avoid divide by 0
    double d1 = (log(F/K) + pow(V,2)*T/2)/(V*sqrt(T));
    double d2 = d1 - V*sqrt(T);
    if (callput == 'C')
    {
        return exp(-r*T)*(F*Ncdf(d1)-K*(Ncdf(d2)));
    }
    else if (callput == 'P')
    {
        return exp(-r*T)*(-F*Ncdf(-d1)+K*(Ncdf(-d2)));
    }
}

// Normal cumulative cdf function
double Ncdf(double x)
{
    return erfc(-x / sqrt(2))/2.0;
}

void writeData(vector<vector<string>> Z, string fileName, string delimiter)
{
	ofstream file(fileName);

    for(int k = 0; k < Z.size(); k++)
    {
        for(int l = 0; l < Z[0].size()-1; l++)
        {
            file << Z[k][l] << delimiter;
        }
        file << Z[k][Z[0].size()-1] << "\n";
    }

	file.close();
}

vector<vector<string>> getData(string fileName, string delimiter)
{
	ifstream file(fileName);

	vector<vector<string> > Z;

	string line = "";

	while (getline(file, line))
	{
		vector<string> vec = split(line, delimiter);
		Z.push_back(vec);
	}

	file.close();

	return Z;
}

vector<string> getDataSingleLine(string fileName, string delimiter)
{
	ifstream file(fileName);

	vector<string> data;

	string line = "";

    getline(file, line);
    data = split(line, delimiter);

	file.close();

	return data;
}

vector<string> split(const string& str, const string& delim)
{
    vector<string> tokens;
    size_t prev = 0, pos = 0;
    do
    {
        pos = str.find(delim, prev);
        if (pos == string::npos) pos = str.length();
        string token = str.substr(prev, pos-prev);
        tokens.push_back(token);
        prev = pos + delim.length();
    }
    while (pos < str.length() && prev < str.length());
    return tokens;
}