Estimating the parameters of a CEV Process

Reading time ~5 minutes

The CEV Process

In mathematical finance, the CEV or constant elasticity of variance model is a stochastic volatility model which was developed by John Cox in 1975. It captures the fact that the log returns may not have a constant volatility, which the usual geometric brownian motion assumes.

The equation of the CEV model is the following.

For we have the usual equation for the geometric Brownian motion:

Simulating processes

Using the discretization method, we can simulate these processes.

Geometric Brownian motion

The geometric Brownian motion can be simulated using the following class.

import numpy as np


class ProcessGBM:

    def __init__(self, mu, sigma):
        self._mu = mu
        self._sigma = sigma

    def Simulate(self, T=1, dt=0.001, S0=1.):
        n = round(T / dt)
        
        mu = self._mu
        sigma = self._sigma

        gaussian_increments = np.random.normal(size=n - 1)
        res = np.zeros(n)
        res[0] = S0
        S = S0
        sqrt_dt = dt ** 0.5
        for i in range(n - 1):
            S = S + S * mu * dt + sigma * \
                S * gaussian_increments[i] * sqrt_dt
            res[i + 1] = S

        return res

Following a similar logic, we can generate trajectories for a CEV process.

import numpy as np


class ProcessCEV:

    def __init__(self, mu, sigma, gamma):
        self._mu = mu
        self._sigma = sigma
        self._gamma = gamma

    def Simulate(self, T=1, dt=0.001, S0=1.):
        n = round(T / dt)
        
        mu = self._mu
        sigma = self._sigma
        gamma = self._gamma

        gaussian_increments = np.random.normal(size=n - 1)
        res = np.zeros(n)
        res[0] = S0
        S = S0
        sqrt_dt = dt ** 0.5
        for i in range(n - 1):
            S = S + S * mu * dt + sigma * \
                (S ** gamma) * gaussian_increments[i] * sqrt_dt
            res[i + 1] = S

        return res

Which can then be plotted:

import matplotlib.pyplot as plt

from ProcessGBM import ProcessGBM
from ProcessCEV import ProcessCEV

T = 20
dt = 0.01
plt.plot(ProcessGBM(0.05, 0.15).Simulate(T, dt), label="GBM")
plt.plot(ProcessCEV(0.05, 0.15, 0.5).Simulate(
    T, dt), label="CEV (gamma = 0.5)")
plt.plot(ProcessCEV(0.05, 0.15, 1.5).Simulate(
    T, dt), label="CEV (gamma = 1.5)")

plt.xlabel('Time index')
plt.ylabel('Value')

plt.title("Realization of stochastic processes")

plt.legend()

plt.show()

Simulation of various processes

Realisation of various stochastic processes with different parameters

And we can “see” the properties of the trajectories, a low value for yields trajectories which are slightly linear, whereas a higher than allows the volatility to increase for large values of the process (the trajectory in red).

Estimating the CEV parameters

Theory

A detailed method to estimate the various parameters can be found here: L. CHU, K & Yang, Hailiang & Yuen, Kam. (2001). Estimation in the Constant Elasticity of Variance Model. British Actuarial Journal. 7. 10.1017/S1357321700002233.

Note that they use slightly different notations, which will be kept in this section (the code for the estimator, however, translates everything back in terms of the original notations).

Calling

The paper recalls that is an estimate of , where is:

For a constant.

It looks like there is there a typo in the original paper, after a quick implementation (see below), I realized that the estimator does not work this this formula, but does so when I remove one of the . Therefore, in what follows, the following formula will be implemented:

Then, the variance of is minimized for a value of ( are negative constants)

And the authors suggest to start from a random value of and iterate over successive estimations of and . However, not being constant, this formula did not make a lot of sense to me. I chose a constant value of in the code, which provides good results.

Proof sketch : Itô’s lemma

In order to know where does comes from, we can apply Itô’s formula to :

Now:

Therefore, dividing by :

Finally:

And it works for any . To be honest I am quite uncomfortable with the meaning to give to every equation, but, at some level, it makes sens that “has to be” of the form above.

Python implementation

import numpy as np


class EstimatorCEV:

    def __init__(self, dt):
        self._dt = dt
        self._alpha0 = -5
    
    def Estimate(self, trajectory):
        sigma, gamma = self._evaluate_sigma_gamma(trajectory, self._alpha0)
        if sigma == None:
            return None, None, None
        else:
            mu = self._estimate_mu(trajectory)
            return (mu, sigma, gamma)

    def _log_increments(self, trajectory):
        return np.diff(trajectory) / trajectory[:-1]
    
    def _estimate_mu(self, trajectory):
        return np.mean(self._log_increments(trajectory)) / self._dt 

    def _log_increments_alpha(self, trajectory, alpha):
        mod_increments = self._log_increments(trajectory ** (1 + alpha))
        return mod_increments / (1 + alpha)

    def _evaluate_Vt(self, trajectory, alpha):
        lhs = self._log_increments_alpha(trajectory, alpha)
        rhs = self._log_increments(trajectory)
        center = 2 * (lhs - rhs) / (alpha * self._dt)
        return center

    def _evaluate_sigma_gamma(self, trajectory, alpha):
        if np.any(trajectory <= 0):
            return None, None
        
        Vts = self._evaluate_Vt(trajectory, alpha)
        if np.any(Vts <= 0):
            return None, None
        logVts = np.log(Vts)

        Sts = trajectory[:-1]  # removes the last term as in eq. (5)
        if np.any(Sts <= 0):
            return None, None
        logSts = np.log(Sts)

        ones = np.ones(Sts.shape[0])
        A = np.column_stack((ones, logSts))

        res = np.linalg.lstsq(A, logVts, rcond=None)[0]
        return (2 * np.exp(res[0] / 2), 0.5 * (res[1] + 2))


if __name__ == "__main__":

    from ProcessCEV import ProcessCEV

    def test(true_mu, true_sigma, true_gamma):
        dt = 0.001
        T = 10
        
        sample_mu = []
        sample_sigma = []
        sample_gamma = []

        for i in range(100):
            mu_est, sigma_est, gamma_est = EstimatorCEV(dt=dt).Estimate(ProcessCEV(
                true_mu, true_sigma, true_gamma).Simulate(T, dt=dt))
            
            if mu_est != None:
                sample_mu = [mu_est] + sample_mu
                sample_sigma = [sigma_est] + sample_sigma
                sample_gamma = [gamma_est ] + sample_gamma
    
        print("mu : " + str(true_mu) + " \t| est : " + str(np.mean(sample_mu)) + " \t| std : " + str(np.std(sample_mu)))
        print("sigma : " + str(true_sigma) + " \t| est : " + str(np.mean(sample_sigma)) + " \t| std : " + str(np.std(sample_sigma)))
        print("gamma : " + str(true_gamma) + " \t| est : " + str(np.mean(sample_gamma)) + " \t| std : " + str(np.std(sample_gamma)))
        print(10*"-")


    test(0.,0.5,0.8)
    test(0.2,0.5,0.8)
    test(0.2,0.5,1.2)
    test(0.,0.3,0.2)
    test(0.,0.5,2)

Running the tests…

mu : 0.0  | est : 0.08393998598504683   | std : 0.10950360432103644
sigma : 0.5   | est : 0.5290411001633833  | std : 0.010372415832707949
gamma : 0.8   | est : 0.7995912270171245  | std : 0.023552858445145052
----------
mu : 0.2  | est : 0.21126788898622506   | std : 0.1379314635624178
sigma : 0.5   | est : 0.5295722080080916  | std : 0.01067649238716792
gamma : 0.8   | est : 0.7981059837608774  | std : 0.019084913069205692
----------
mu : 0.2  | est : 0.20468181972280103   | std : 0.16048227152782532
sigma : 0.5   | est : 0.5299253491070517  | std : 0.008389773769199063
gamma : 1.2   | est : 1.2042365895022817  | std : 0.021453412919601324
----------
mu : 0.0  | est : 0.0743428420875073  | std : 0.07189306488901058
sigma : 0.3   | est : 0.31771919915814173   | std : 0.006063917451139664
gamma : 0.2   | est : 0.2080203482037684  | std : 0.039438575768452104
----------
mu : 0.0  | est : -0.009815876356783935   | std : 0.11178850054168586
sigma : 0.5   | est : 0.5311825321732557  | std : 0.013167973885022026
gamma : 2   | est : 2.0050079526000903  | std : 0.03448408469682434
----------

It worked !

Learning more

I used to study stochastic processes a while ago and reading: Pierre Henry-Labordère – Analysis, Geometry & Modeling in Finance (which I recommend if you are interested in stochastic processes and differential geometry and mathematical finance in the same time) and A Non-Gaussian Option Pricing Model with Skew L. Borland, J.P. Bouchaud led me to spend more time understanding them. Hence this small tuto, hope you liked it.

OCaml functor example : a progress bar

# IntroductionPython has an amazing progress bar library [tqdm](https://github.com/tqdm/tqdm). In a for loop or a map statement, ````tqdm...… Continue reading