Recently I was working with a process named Variance Gamma with Stochastic Arrival (VGSA) and trying to fit this process on a given data.

To obtain VGSA, as explained in Carr et al. [2001], we take the VG process which is a homogeneous Lévy process and build in stochastic volatility by evaluating it at a continuous time change given by the integral of a Cox, Ingersoll, and Ross [1985] (CIR) process. Formally we define the CIR process $y(t)$ as the solution to the stochastic differential equation, $$ d y_t=\kappa\left(\eta-y_t\right) d t+\lambda \sqrt{y_t} d W_t $$ where $W_t$ is a Brownian motion, $\eta$ is the long-term rate of time change, $\kappa$ is the rate of mean reversion, and $\lambda$ is the volatility of the time change. The process $y(t)$ is the instantaneous rate of time change and so the time change is given by $Y(t)$ where $$ Y(t)=\int_0^t y(u) d u $$

The characteristic function for $Y(t)$ is given by $$ \begin{aligned} \mathbb{E}\left(e^{i u Y(t)}\right) & =\phi(u, t, y(0), \kappa, \eta, \lambda) \\ & =A(t, u) e^{B(t, u) y(0)} \end{aligned} $$ where $$ \begin{aligned} A(t, u) & =\frac{\exp \left(\frac{\kappa^2 \eta t}{\lambda^2}\right)}{\left(\cosh (\gamma t / 2)+\frac{\kappa}{\gamma} \sinh (\gamma t / 2)\right)^{2 \kappa \eta / \lambda^2}} \\ B(t, u) & =\frac{2 i u}{\kappa+\gamma \operatorname{coth}(\gamma t / 2)} \end{aligned} $$ with $$ \gamma=\sqrt{\kappa^2-2 \lambda^2 i u} $$

The stochastic volatility Lévy process, termed the VGSA process, is defined by $$ \begin{aligned} Z(t) & =X_{V G}(Y(t) ; \sigma, \nu, \theta) \\ & =\theta \gamma(Y(t) ; 1, \nu)+\sigma W(\gamma(Y(t) ; 1, \nu)) \end{aligned} $$

Thus $\sigma, \nu, \theta, \kappa, \eta$, and $\lambda$ are the six parameters defining the process. Its characteristic function is given by $$ \mathbb{E}\left(e^{i u Z_{V G S A}(t)}\right)=\phi\left(-i \Psi_{V G}(u), t, \frac{1}{\nu}, \kappa, \eta, \lambda\right) $$ where $\phi$ is the characteristic function of $Y(t)$ given in (6.68)and $\Psi_{V G}$ is the log characteristic function of the variance gamma process at unit time, namely, $$ \Psi_{V G}(u)=-\frac{1}{\nu} \log \left(1-i u \theta \nu+\sigma^2 \nu u^2 / 2\right) $$

We define the stock process at time $t$ by the random variable $$ S(t)=S(0) \frac{e^{(r-q) t+Z(t)}}{\mathbb{E}\left[e^{Z(t)}\right]} $$

We note that $$ \mathbb{E}\left[e^{Z(t)}\right]=\phi\left(-i \Psi_{V G}(-i), t, \frac{1}{\nu}, \kappa, \eta, \lambda\right) $$ which is equivalent to $e^{-\omega t}$ in the VG case. Therefore the characteristic function of the $\log$ of the stock price at time $t$ is given by $$ \mathbb{E}\left[e^{i u \log S_t}\right]=\exp \left(i u\left(\log S_0+(r-q) t\right)\right) \times \frac{\phi\left(-i \Psi_{V G}(u), t, \frac{1}{\nu}, \kappa, \eta, \lambda\right)}{\phi\left(-i \Psi_{V G}(-i), t, \frac{1}{\nu}, \kappa, \eta, \lambda\right)^{i u}} $$

For additional information, please refer to page 236 of the book titled "Computational Methods in Finance" by Ali Hirsa, specifically regarding the Variance-Gamma Stochastic Process (VGSA).

My approach involves initially inverse Fourier transforming the characteristic function and then utilizing the Log Likelihood method. Despite employing an initial guess and the Nelder–Mead method for minimization, the results obtained are often suboptimal, likely due to convergence to local minima.

I seek guidance on a more robust method to fit the VGSA model to my data for parameter estimation, as the current approach seems to be getting stuck at local minima. Any suggestions and code implementations in Python, R, or MATLAB would be highly appreciated.

  • $\begingroup$ Just to clarify: You have a time series of asset prices and try to fit the process to that history (ie, estimation under the real world measure $\mathbb{P}$). You don't want to calibrate the model to a set of observed option prices (under a risk-neutral measure $\mathbb{Q}$)? $\endgroup$
    – Kevin
    Commented Mar 1 at 17:45
  • $\begingroup$ @Kevin Sorry for the late reply. You are right from the Option pricing theoretic view, though my requirement is slightly different. I have a process here whose characteristic function is known to us but not the distribution. I want to estimate the parameters governing the process in this scenario from some real-world time series data. Please let me know if any further clarification is needed. $\endgroup$
    – Starlord22
    Commented Mar 2 at 18:18

1 Answer 1


General remarks

A difficult but very interesting problem. Some thoughts:

  • Use GMM instead of MLE. MLE is a special case of the generalised method of moments. The characteristic function gives you the moments via differentiating ($\mathbb{E}[X^n]=i^{-n}\varphi_X^{(n)}(0)$). Thus, you do not need to invert the characteristic function. As there are six parameters, you'd however need 6-7 moments which can be difficult to estimate from a time series.

  • Have you tried using fast Fourier transform? The increased speed might give you more capacity to focus on accuracy and stability? Carr, Geman, Madan, and Yor (2002, p. 320) explain the estimation of their CGMY model from historical data as follows

For each underlying asset, we formed the time series of daily log price relatives and then estimated the parameters of the Lévy density $C, G, M, Y,$ and $\eta$ from the mean-adjusted return data. Direct maximum likelihood estimation is computationally expensive, as it requires a Fourier inversion for each data point to evaluate the density, and these inversions must be nested into a gradient search optimization algorithm for the parameter estimation. The fast Fourier transform was used to invert the characteristic function once for each parameter setting. This method efficiently renders the level of the probability density at a prespecified set of values for returns. (...) With the density evaluated at these prespecified points, we binned the return series by counting the number of observations at each prespecified return point, assigning data observations to the closest prespecified return point. We then searched for parameter estimates that maximized the likelihood of this binned data. The reported estimates are thus for this binned maximum likelihood estimation using the fast Fourier transform.

  • Probably, using COS method instead of FFT would be even better.

  • There is actually a literature on how to better discipline the MLE estimation or using alternative estimators (eg MCMC). See here for a brief overview.

  • Not helpful at all but you might want to consider using a model with closed-form density (like Variance Gamma).

Details on COS method

I start with citing a well-known Lemma from analysis using the Fourier transform as an approximation of the original function.

Let $f:[a,b]\to\mathbb{R}$ be an integrable function with Fourier transform $\hat{f}$. Then, for $N\in\mathbb{N}$, we have \begin{align} f(x) &\approx \frac{2}{b-a} \sum_{n=0}^N \tilde{A_n} \cos\left( n\pi\frac{x-a}{b-a}\right), \\ \tilde{A_n} &= \begin{cases} \text{Re}\left(\hat{f}\left(\frac{n\pi}{b-a}\right)\cdot e^{-i\frac{n\pi a}{b-a}}\right) & \text{if } n\geq1, \\ \frac{1}{2} \text{Re}\left(\hat{f}\left(\frac{n\pi}{b-a}\right)\cdot e^{-i\frac{n\pi a}{b-a}}\right) & \text{if } n=0. \end{cases} \nonumber \end{align}

It's easy to prove. All credit to Fang and Oosterlee (2008). But you can already see where I am going. $\hat{f}$ is the known characteristic function, $f$ is the sought probability density function of $\ln(S_t)$. Critically, note that there is no Fourier inversion, no integral, no nothing. Just a sum. That means you can be very fast and focus on accuracy and stability.

A problem is the choice of $N$ and the bounds $a$ and $b$. Clearly, most density functions do not have compact support. However, even if they have fat tails, they can normally be safely truncated. If a log stock price is at $\ln(10)$, it's probably very unlikely you need to go all the way up to $\ln(1,000)$. Fang and Oosterlee have suggestions how you can choose $a$ and $b$ based on the cumulants of $\ln(S_t)$ but we don't know them without the parameters, so I would just take generously chosen values and do a sensitivity analysis. Pick some, estimate the parameters using MLE and then change $a$ and $b$ and restimate the parameters. If nothing changes, you're fine.

Below is a plot for the standard normal distribution with density $f:[-10,10]\to\mathbb{R},x\mapsto \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}$ and Fourier transform $\hat{f}(u)=e^{-\frac{1}{2}u^2}$. You see, it converges super quickly. For $N=20$, the plots are visually indistinguishable.

enter image description here

  • $\begingroup$ really thanks for this answer, I actually adopted the GMM method but due to the complicated nature of the function it become challenging for me to evaluate, please let me know if there is any guide/implementation to compute moments in such cases using the time series data computationally. I'ven't consider FFT (and COS), thank you for the guide, I'm currently going through that. I will really appreciate if you show some basic implementation FFT/COS in these scenario maybe with some simulated data in any language. $\endgroup$
    – Starlord22
    Commented Mar 5 at 12:10
  • $\begingroup$ @Starlord22 I added some details on COS. I think that's one of the fastest methods out there. It's not very demanding computationally and saves solving integrals numerically. Does that help? Please let me know if you have further questions! $\endgroup$
    – Kevin
    Commented Mar 6 at 9:07
  • $\begingroup$ thanks a lot for the answer i really helped me to understand the cos method, will work on it and if I need any further clarification will let you know. $\endgroup$
    – Starlord22
    Commented Mar 7 at 8:24
  • $\begingroup$ @Starlord22 happy to hear that! :) Let us know (in the comments) if you need more help and if it worked out - I'm curious to know ;) $\endgroup$
    – Kevin
    Commented Mar 7 at 8:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.