Modeling and Pricing Options for $VIX
Using a Regime-Switching Stochastic Process and Particle Expectation Maximization
One day I was bored, and while scrolling on X.com I saw a picture of the $VIX curve:
Its shape differs from that of a stock price curve, typically modeled using a geometric Brownian motion stochastic differential equation (GBM SDE):
\[dS_t=S_t(\mu dt+\sigma dW_t)\]
A very loose explanation:
An change in a stock price \(S_t\) on an infinitesimal time step \(dt\) is defined by drift \(\mu\), plus variance \(\sigma\) times the corresponding change of a Wiener process \(W_t\), all proportional to the stock price itself. A Wiener process \(W_t\) is a random process defined by the following:
- \(W_0=0\) almost surely.
- \(W\) has independent increments each following a standard normal distribution.
- \(W_t\) is almost surely continuous in \(t\).
$VIX is not valued by price — it is derived from the volatility of S&P 500 options. Specifically it is “a weighted average of its out-of-the-money call and put options” with expiry over the next 30 days:
\[V_{\text{VIX}}=\sqrt{\frac{2e^{r\tau}}\tau\left(\int_0^F\frac{P(K)}{K^2}dK+\int_F^\infty\frac{C(K)}{K^2}dK\right)}\]
- \(\tau\) — Days in a month (30)
- \(r\) — risk-free rate
- \(F\) — 30-day forward price of S&P 500
- \(P(K)\) and \(C(K)\) — put and call prices with strike \(K\) and maturity in 30 days
It can be thought of as the “stock market’s expectation of volatility based on S&P 500 index options.”
Relevant stochastic models
The stochastic process that provides the most accurate representation of $VIX is unclear — we can choose from several options.
I want to first discuss a generic mean-reverting Ornstein-Uhlenbeck process:
\[dx_t=\kappa(\theta-x_t)dt+\sigma dW_t\]
…where \(\kappa\) is reversion speed, \(\theta\) is the long-term mean, and \(\sigma\) is variance. An Ornstein-Uhlenbeck process’ drift \(\kappa(\theta-x_t)\) is strongly defined by its distance from the mean — for example \(x_t =\theta\) means drift equals 0.
Next I’ll bring up the Cox-Ingersoll-Ross (CIR) process (normally used as an interest short-rate model):
\[dx_t=\kappa(\theta -x_t)dt+\sigma\sqrt{x_t}dW_t\]
…the only difference in structure being the extra factor \(\sqrt{x_t}\). This produces larger variance at higher values, and low to vanishing variance at 0 with drift back to mean. A CIR process is guaranteed to be positive given the condition \(2\kappa\theta\geq\sigma^2\), known as the Feller condition.
In general $VIX will repeatedly “spike” and then subside back to some average, so the above models will be useful.
Derived $VIX modeling from the S&P 500
Considering $VIX is derived from GBM-following stock options, we can use a Heston model for pricing of S&P 500 options, within which a stochastic volatility process is derived. A Heston model consists of two SDEs:
- Price SDE: \( dS_t=\mu S_tdt+\sqrt{\nu_t}S_tdW_t^S \)
- Variance SDE: \( d\nu_t=\kappa(\theta-\nu_t)d+\xi\sqrt{\nu_t}dW_t^\nu \)
With initial variance \(\nu_0\), long-term average variance \(\theta\), reversion speed \(\kappa\), volatility of volatility \(\xi\), and the correlation between the two Wiener processes \(W_t^S\) and \(W_t^\nu\) equal to some \(\rho\). Notice the volatility \(\sqrt{\nu_t}\) is defined by a Cox-Ingersoll-Ross process itself.
As a higher-dimensional SDE, simulation is more computationally expensive.
The standard Bates model adds Poisson-distributed (randomly timed) jumps to the price evolution of a Heston model:
\[dS_t=\mu S_tdt+\sqrt{\nu_t}S_tdW_t^S+(J_t-1)S_{t^-}dN_t\]
…with Poisson process \(N_t(\lambda)\) and exponentially distributed intensity \(J_t\). It can capture the presence of surprise events that shock an asset price (say a Fed announcement boosting the S&P 500).
Direct modeling of $VIX using mean-reverting processes
Direct modeling of $VIX is simpler but carries a loss of information and degrades its arbitrage-free valuation. For fun and simplicity we will still look into it.
A generic Ornstein-Uhlenbeck process, as described above, is the simplest mean-reverting model. It allows for easy simulation and parameter estimation. That said, the process has normally distributed variation, does not fully capture $VIX’s “spike” behavior, and can go negative. We can potentially use it by applying a transformation to $VIX.
A Cox-Ingersoll-Ross process with Feller condition, as described above, is more accurate in distribution given behavior above and below mean, but incurs some cost in parameter estimation.
A CIR process with jump-diffusion is as follows:
\[dx_t=\kappa(\theta -x_t)dt+\sigma\sqrt{x_t}dW_t+(J_t-1)x_{t^-}dN_t\]
Of course we add the jump-diffusion term \((J_t-1)x_{t^-}dN_t\) to the aforementioned CIR process. This, although possessing even more complexity, seems like a pretty good direct model.
Finally we discuss regime-switching processes. In such a stochastic process we have \(K\) subprocesses each corresponding to a state \(S_t\in\{1,…,K\}\) that alternate in activity dependent on a Markov transition matrix \(P\). At each time step the conditional probability of transition from current state \(i\) to potential state \(j\) is defined by \(p_{ij}=P(S_{t+1}=j|S_t=i)\). In our case we can use 2 states we call “calm” and “crisis”.
I decided to try implementing a regime-switching process…here we go.
Implementing the RegimeProcess wrapper and Monte Carlo Engines
This project was my first or deepest interaction with several topics including stochastic process models, Monte Carlo pricing, the calibration and approximation methods involved, some important topics in C++ object-oriented programming, Python bindings, and more. It had multiple goals that changed as time progressed, overcame many issues and shortcomings, and was a good learning experience.
In general I aimed to establish:
- A C++ class RegimeProcess representing a generic set of stochastic processes that switch depending on the state of some Markov chain
- Compatibility with QuantLib’s ecosystem of data structures and algorithms — this potentially includes interest rate term structures, option payoff structures, pricing engines, etc.
- A European option pricing script runnable in a Jupyter notebook — leveraging Python’s flexibility and C++’s efficiency
Eventually I implemented the Monte Carlo pricing engine myself, and calibration infrastructure in Python — I’ll explain why in their respective sections.
RegimeProcess.cpp
The RegimeProcess is a stochastic process object that requires a set of subprocesses processes, a Markov transition matrix transitionMatrix, and a time step dt.
We write the definition:
RegimeProcess::RegimeProcess(
const std::vector<boost::shared_ptr<StochasticProcess>>& processes,
const Matrix& transitionMatrix,
Time dt)
: processes_(processes), P_(transitionMatrix), dt_(dt), K_(processes.size()) {
QL_REQUIRE(K_ > 0, "must provide at least one regime process");
QL_REQUIRE(
P_.rows() == static_cast<Size>(K_) && P_.columns() == static_cast<Size>(K_),
"transition matrix must be KxK"
);
// Check that all processes have the same dimensions
Size d = processes[0]->size();
for (const auto& p : processes) {
QL_REQUIRE(p->size() == d, "all processes must have the same dimension");
}
// Initialize RNG
rng_.seed(123456);
currentRegime_ = 0;
}We add a function to choose the regime of the next time step given the transition matrix:
Size RegimeProcess::sampleNextRegime(Size current) const {
std::vector<double> probs(K_);
for (Size j = 0; j < K_; ++j)
probs[j] = P_[current][j];
std::discrete_distribution<int> dist(probs.begin(), probs.end());
return static_cast<Size>(dist(rng_));
}A function to evolve the SDE based on the active regime:
Array RegimeProcess::evolve(
Time t0, const Array& x0, Time dt, const Array& dw
) const {
// Sample next regime index
Size next = sampleNextRegime(currentRegime_);
// Delegate evolution to the selected regime's process
Array x1 = processes_[next]->evolve(t0, x0, dt, dw);
// Update current regime
currentRegime_ = next;
return x1;
}Other relevant inherited functions/getters e.g.:
Matrix RegimeProcess::stdDeviation(Time t0, const Array& x0, Time dt) const {
return processes_[currentRegime_]->stdDeviation(t0, x0, dt);
}And a header file RegimeProcess.hpp structuring the class.
MonteCarloEngine.cpp
Monte Carlo simulation of random processes consists of generating numerous paths and analyzing the information that it reveals. In the case of option pricing, we repeatedly simulate the evolution of the underlying asset price from the present day, calculate the value of each path’s contract at maturity, average them, and discount the price based on the risk-free rate of return.
Monte Carlo simulation is used to value many types options with complex or nonexistent analytic solutions. It is computationally expensive. We will have to use it to value our regime-switching process options.
It's my understanding that the stateful property of the RegimeProcess is what warranted the explicit scripting of the Monte Carlo engines, rather than simple insertion into QuantLib’s built-in engines.
This project contains two engines; a generic RegimeSwitchingMCEngine, and a dedicated VIXRSMCEngine. The dedicated engine accounts for known factors such as 1 business day time steps, European/vanilla contracts and payoff structures, underlying valuation not measured in money, non-negativity, etc.
The generic engine is defined as follows:
RegimeSwitchingMCEngine::RegimeSwitchingMCEngine(
const boost::shared_ptr<RegimeProcess>& process,
Size numPaths,
Size timeSteps,
BigNatural seed)
: process_(process), numPaths_(numPaths), timeSteps_(timeSteps), seed_(seed) {
QL_REQUIRE(process_, "process cannot be null");
QL_REQUIRE(numPaths_ > 0, "number of paths must be positive");
QL_REQUIRE(timeSteps_ > 0, "number of time steps must be positive");
}And has a European pricing function with arbitrary payoff:
Real RegimeSwitchingMCEngine::priceEuropean(
Time maturity,
const Array& initialState,
Size initialRegime,
Rate riskFreeRate,
const boost::function<Real(const Array&)>& payoff) {
Time dt = maturity / timeSteps_;
Size dimensions = process_->size();
Real sumPayoffs = 0.0;
// Create uniform random sequence generator
MersenneTwisterUniformRng uniformGen(seed_);
RandomSequenceGenerator unifSeqGen(dimensions, uniformGen);
// Create inverse cumulative RSG for normal variates
InverseCumulativeRsg<RandomSequenceGenerator<
MersenneTwisterUniformRng>, InverseCumulativeNormal
> rng(unifSeqGen);
for (Size path = 0; path < numPaths_; ++path) {
process_->setRegime(initialRegime);
Array state = initialState;
Time t = 0.0;
for (Size step = 0; step < timeSteps_; ++step) {
// Get random normal vector
const std::vector<Real>& gaussian = rng.nextSequence().value;
Array dw(dimensions);
for (Size d = 0; d < dimensions; ++d) {
dw[d] = gaussian[d] * std::sqrt(dt);
}
// Evolve the process
state = process_->evolve(t, state, dt, dw);
t += dt;
}
// Calculate payoff at maturity
Real thisPayoff = payoff(state);
sumPayoffs += thisPayoff;
}
// Discount and average
Real discount = std::exp(-riskFreeRate * maturity);
return (sumPayoffs / numPaths_) * discount;
}The dedicated engine is defined similarly, but specifically has, for example, a European call pricing function:
Real VIXRSMCEngine::priceCall(
Real initialVIX,
Size initialRegime,
Real strike,
Size expiryDays,
Rate riskFreeRate) {
Time dt = 1.0 / 252.0; // Daily steps (business days convention)
Time maturity = expiryDays * dt;
Size dimensions = process_->size();
Real sumPayoffs = 0.0;
// Create uniform random sequence generator
MersenneTwisterUniformRng uniformGen(seed_);
RandomSequenceGenerator unifSeqGen(dimensions, uniformGen);
// Create inverse cumulative RSG for normal variates
InverseCumulativeRsg<RandomSequenceGenerator<
MersenneTwisterUniformRng>, InverseCumulativeNormal
> rng(unifSeqGen);
for (Size path = 0; path < numPaths_; ++path) {
process_->setRegime(initialRegime);
Array state = process_->initialValues();
state[0] = initialVIX / 100.0; // Convert to decimal
Time t = 0.0;
for (Size step = 0; step < expiryDays; ++step) {
const std::vector<Real>& gaussian = rng.nextSequence().value;
Array dw(dimensions);
for (Size d = 0; d < dimensions; ++d) {
dw[d] = gaussian[d] * std::sqrt(dt);
}
state = process_->evolve(t, state, dt, dw);
t += dt;
// Ensure non-negative VIX
if (state[0] < 0.0) state[0] = 0.0;
}
// Convert back to VIX percentage terms
Real finalVIX = state[0] * 100.0;
// Call payoff
Real payoff = std::max(finalVIX - strike, 0.0);
sumPayoffs += payoff;
}
// Discount and average
Real discount = std::exp(-riskFreeRate * maturity);
return (sumPayoffs / numPaths_) * discount;
}Finally there is an accompanying header file.
bindings.cpp and compilation
Usage of this infrastructure in Python and Jupyter notebooks requires Python bindings produced with PyBind11. We write some necessary functions create_cir_process, create_ou_process, numpy_to_array, array_to_numpy, and matrix_to_numpy. For example:
boost::shared_ptr<StochasticProcess> create_cir_process(
Real x0, Real kappa, Real theta, Real sigma) {
return boost::shared_ptr<StochasticProcess>(
new CoxIngersollRossProcess(kappa, theta, sigma, x0)
);
}Then we construct the module vixmodels in order to export to Python:
PYBIND11_MODULE(vixmodels, m) {
// ...
m.def("numpy_to_array", &numpy_to_array,
py::arg("arr"), "Convert numpy array to QuantLib Array");
// RegimeProcess wrapper
py::class_<RegimeProcess, boost::shared_ptr<RegimeProcess>>(m, "RegimeProcess")
.def(py::init<const std::vector<boost::shared_ptr<StochasticProcess>>&,
const Matrix&, Time>(),
py::arg("processes"), py::arg("transition_matrix"), py::arg("dt"))
.def("set_regime", &RegimeProcess::setRegime, py::arg("regime"))
.def("regime", &RegimeProcess::regime)
.def("size", &RegimeProcess::size)
.def("initial_values", &RegimeProcess::initialValues)
.def("drift", &RegimeProcess::drift, py::arg("t"), py::arg("x"))
.def("diffusion", &RegimeProcess::diffusion, py::arg("t"), py::arg("x"));
// ...
}Finally we have a simple build script build.sh and CMakeLists.txt.
Methods of calibration and pricing
We have the goal of calibrating a RegimeProcess to historical $VIX data, and by representing $VIX as a regime-switching pair of stochastic processes (of our choice), we will price European call and put options. I decided to set the pair of processes as a Cox-Ingersoll-Ross process and an Ornstein-Uhlenbeck process.
The properties and complexity of the process to be calibrated determines the calibration options available. While a single Ornstein-Uhlenbeck is suited for simple calibration, the usage of a state-dependent variance term \(\sigma\sqrt{x_t}dW_t\) in a CIR process rules out methods dependent on Gaussian distributions. More importantly the use of a pair of stochastic processes, each with unknown parameters and dependent on a hidden Markov model with an unknown transition matrix largely rules out analytical methods.
I am very unfamiliar with the following methods so I can only provide a brief explanation:
We use an iterative expectation-maximization (EM) algorithm — in particular a particle EM (pEM) algorithm. Given historical data, we estimate the state of the regime-switching process using a particle filter, itself a Monte Carlo algorithm. We leverage this to find the most likely parameters of each stochastic process given their log-likelihood functions, and repeat.
We receive the calibrated sets of parameters \(\Theta=\{\kappa_j,\theta_j,\sigma_j,P_{ij}\}\), and the probability of each regime’s activity given the day. With these parameters we build the RegimeProcess to perform Monte Carlo pricing.
IMPORTANT: \(P\)-measure and \(Q\)-measure
The fitted parameters are based on the \(P\)-measure — market prices subject to risk premiums and arbitrage opportunities. Accurate option pricing must be calibrated again, usually with an option surface, in order to produce \(Q\)-measure risk-neutral prices. Because of complexity and time constraints, this could not be added to the projects current form. Thus, the calculated option values are in “ballpark range”, but not accurate.
Implementing calibration and pricing infrastructure
We first create generic ParticleFilter and ParticleEM classes in particle_filter.py:
class ParticleFilter:
def __init__(
self,
num_particles: int, num_regimes: int, state_dim: int, seed: int = 42
):
self.N = num_particles
self.K = num_regimes
self.d = state_dim
self.rng = np.random.default_rng(seed)
def systematic_resampling(self, weights: np.ndarray) -> np.ndarray:
...
def filter(
self,
observations: np.ndarray, initial_state_sampler: Callable,
transition_fn: Callable, observation_fn: Callable,
observation_noise_std: float
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float]:
...
class ParticleEM:
def __init__(
self,
num_particles: int, num_regimes: int, state_dim: int, seed: int = 42
):
self.pf = ParticleFilter(num_particles, num_regimes, state_dim, seed)
self.K = num_regimes
self.d = state_dim
self.rng = np.random.default_rng(seed)
def backward_sampling(
self,
states: np.ndarray, regimes: np.ndarray, weights: np.ndarray,
transition_matrix: np.ndarray
) -> Tuple[np.ndarray, np.ndarray]:
...
def estimate_transition_matrix(
self,
regime_trajectories: List[np.ndarray]
) -> np.ndarray:
...
def fit(
self,
observations: np.ndarray, initial_params: Dict,
max_iterations: int = 50, tolerance: float = 1e-4,
num_trajectories: int = 10
) -> Tuple[Dict, List[float]]:
...
def _e_step(
self,
observations: np.ndarray, params: Dict
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float]:
# implemented by subclass only
def _m_step(
self,
observations: np.ndarray, trajectories_states: List[np.ndarray],
trajectories_regimes: List[np.ndarray], old_params: Dict
) -> Dict:
# implemented by subclass only<Now we build the subclass dedicated to the VIX CIR/OU process in calibration.py:
class VIXParticleEM(ParticleEM):
def __init__(
self,
num_particles: int = 1000, num_regimes: int = 2, seed: int = 42
):
super().__init__(num_particles, num_regimes, state_dim=1, seed=seed)
self.dt = 1.0 / 252.0 # Daily time step
self.rng = np.random.default_rng(seed)
def _cir_transition_vectorized(
self,
x: np.ndarray, kappa: float, theta: float, sigma: float, dt: float
) -> np.ndarray:
...
def _cir_transition_vectorized(
self,
x: np.ndarray, kappa: float, theta: float, sigma: float, dt: float
) -> np.ndarray:
...
def _sample_regimes_vectorized(
self,
current_regimes: np.ndarray, P: np.ndarray
) -> np.ndarray:
...
def _e_step(
self,
observations: np.ndarray, params: Dict
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float]:
P = params['transition_matrix']
process_params = params['process_params']
obs_noise_std = params['observation_noise_std']
def initial_state_sampler():
...
def transition_fn(states, regimes, t):
...
def observation_fn(states, regimes):
...
return self.pf.filter(
observations, initial_state_sampler, transition_fn,
observation_fn, obs_noise_std
)
def _m_step(
self,
observations: np.ndarray, trajectories_states: List[np.ndarray],
trajectories_regimes: List[np.ndarray], old_params: Dict
) -> Dict:
...
def initialize_vix_params(
observations: np.ndarray, num_regimes: int = 2
) -> Dict:
...
def calibrate_vix_model(
vix_data: np.ndarray, num_regimes: int = 2, num_particles: int = 1000,
max_iterations: int = 30, num_trajectories: int = 10, seed: int = 42
) -> Tuple[Dict, List[float]]:
...Finally we define the Python-wrapped pricing engines and script the remaining infrastructure in pricing.py:
...
try:
# vixmodels.cpython-313-darwin.so located in project/python
import vixmodels
except ImportError:
print("Error: vixmodels module not found. Did you build the C++ extension?")
print("Run: mkdir build && cd build && cmake .. && make && cd ..")
sys.exit(1)
class GenericMCPricer:
def __init__(
self,
regime_process, num_paths: int = 50000, time_steps: int = 252,
seed: int = 42
):
self.engine = vixmodels.RegimeSwitchingMCEngine(
regime_process, num_paths, time_steps, seed
)
def price_european(
self,
initial_val: float, initial_regime: int, strike: float,
maturity: float, option_type: str, risk_free_rate: float = 0.05
) -> float:
...
def vix_payoff_func(x):
value = float(x[0] * 100.0)
return float(payoff(value))
return self.engine.price_european(
maturity, initial_state, initial_regime, risk_free_rate,
vix_payoff_func
)
class VIXDedicatedMCPricer:
def __init__(
self,
regime_process, num_paths: int = 50000, seed: int = 42
):
self.engine = vixmodels.VIXRSMCEngine(regime_process, num_paths, seed)
def price_vix_option(
self,
initial_vix: float, initial_regime: int, strike: float, expiry_days: int,
option_type: str = 'call', risk_free_rate: float = 0.05
) -> float:
if option_type.lower() == 'call':
price = self.engine.price_call(
initial_vix,
initial_regime,
strike,
expiry_days,
risk_free_rate
)
elif option_type.lower() == 'put':
price = self.engine.price_put(
initial_vix,
initial_regime,
strike,
expiry_days,
risk_free_rate
)
else:
raise ValueError("option_type must be 'call' or 'put'")
return price
def load_vix_data(
filepath: str, date_col: str = 'Date', vix_col: str = 'VIX'
) -> pd.DataFrame:
...
def create_regime_process(
params: Dict, dt: float = 1.0/252.0
):
...
def infer_current_regime(
vix_data: np.ndarray, params: Dict, num_particles: int = 1000,
window_size: int = 20
) -> Tuple[int, np.ndarray]:
...
def initial_state_sampler():
...
def transition_fn(states, regimes, t):
...
def observation_fn(states, regimes):
...
...pricing_test.ipynb
Now we are able to build our RegimeProcess and Monte Carlo engine C++ objects through Jupyter Notebooks to perform pricing and calibration on historic data.
Calibration is performed as follows:
# load values
data_path = os.path.join('..', 'data', 'vix_data.csv')
vix_df = load_vix_data(data_path)
vix_values = vix_df['VIX'].values
calibrated_params, log_likelihoods = calibrate_vix_model(vix_values)
current_regime, regime_probs = infer_current_regime(vix_values, calibrated_params)
regime_process = create_regime_process(calibrated_params)
current_vix = vix_values[-1]
expiry = 28And we use a pricing engine as such:
pricer = GenericMCPricer(regime_process)
maturity = expiry / 252.0
for strike in range(19, 25, 2):
call_price = pricer.price_european(
initial_val=current_vix / 100.0,
initial_regime=current_regime,
strike=strike,
maturity=maturity,
option_type='call'
)
put_price = pricer.price_european(
initial_val=current_vix / 100.0,
initial_regime=current_regime,
strike=strike,
maturity=maturity,
option_type='put'
)We can also visualize, over time, the likelihood of each regime’s activity.
Results and conclusions
The dedicated Monte Carlo engine runs ~9x faster than the generic engine; it shaves off some Python overhead. Up to October 17th 2025, options with expiry 28 business days into the future displayed by Robinhood were as follows:
| Call Ask | Call Bid | Strike | Put Bid | Put Ask |
|---|---|---|---|---|
| $3.42 | $2.67 | 19 | $1.50 | $1.96 |
| $2.77 | $2.50 | 21 | $3.00 | $3.40 |
| $2.34 | $1.75 | 23 | $4.39 | $5.03 |
And behold our calculated prices:
| Generic Call | Dedicated Call | Strike | Generic Put. | Dedicated Put |
|---|---|---|---|---|
| $1.70 | $1.71 | 19 | $0.00 | $0.00 |
| $0.01 | $0..17 | 21 | $0.30 | $0.45 |
| $0.00 | $0.00 | 23 | $2.28 | $2.27 |
We are aiming for some value within the bid-ask spread. There are no alarming computational errors (earlier I was given something on the order of 10^98) but the remaining discrepancies are likely from the risk-adjusted valuation (\(Q\)-measure not accounted for). This would be a significant but necessary extra task for accurate valuation.
Visualization of regime probabilities also revealed something interesting. I first calibrated a RegimeProcess over ~20 months, followed by 72 months.
Short horizon:
Long horizon:
I never did clarify which of the CIR and OU processes were the “calm” and “crisis” regime. The difference in time horizon (and potentially some quirks of the initial values) has caused these definitions to flip (Regime 0 is CIR, Regime 1 is OU). It would help to have a more concrete definition of what we are looking for, a method of evaluation, or perhaps we can try many more combinations of processes.
I was also surprised by how strictly the RegimeProcess adhered to the calm regime, suggesting that times of crisis occur not during $VIX’s period spikes (which I suppose can be a property of a CIR process), but instead during genuine “newsworthy” crisis periods — Liberation Day and the beginning of the COVID-19 pandemic are clearly visible here, and the rest of the spikes can likely be researched easily. It also suggests that the rise of the $VIX level is a characterization of the crisis period, but not its subsidence.
This project, born by chance covered several topics, and can expand or be greatly refined in several directions. It did however come as part of an interest I had taken up in financial derivative pricing, and was a chance for me to develop some OOP and computational modeling skills. I hope you enjoyed reading.
GitHub link: https://github.com/zani-t/regime-process
References
- Doucet, A., & Johansen, A. M. (2009). A tutorial on particle filtering and smoothing.
- Cappé, O., Moulines, E., & Rydén, T. (2005). Inference in Hidden Markov Models.
- Cox, J. C., Ingersoll, J. E., & Ross, S. A. (1985). A theory of the term structure of interest rates.
- Hamilton, J. D. (1989). A new approach to the economic analysis of nonstationary time series.
- Ornstein, L. S., & Uhlenbeck, G. E. (1930). On the theory of the Brownian motion.
- Heston, S. L. (1993). A closed-form solution for options with stochastic volatility with applications to bond and currency options.
- Bates, D. S. (1996). Jumps and stochastic volatility: Exchange rate processes implicit in Deutsche Mark options.
- Chicago Board Options Exchange (CBOE). (2019). The Cboe Volatility Index – VIX: White Paper. Chicago Board Options Exchange.
- CenterPoint Securities. (n.d.). What Is the VIX Index and How Does It Work? Retrieved October 2025, from https://centerpointsecurities.com/vix-index/
Comments
Post a Comment