Estimating VGSA Parameters from Its Characteristic Function

# To estimate the parameters when only the characteristic function is known to us


# To estimate the parameters when only the characteristic function is known to us












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.

## Answer by Kevin (score 5, accepted)

https://quant.stackexchange.com/a/78557

### 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.