Computational method for probability distribution on recursive relationships in financial applications

TL;DR

This paper introduces a recursive computational method to accurately determine probability distributions of functions of observed financial data, applicable to various financial instruments and statistical tests.

Contribution

It presents a novel recursive approach for calculating distributions in financial applications, improving accuracy over existing methods.

Findings

High accuracy demonstrated in numerical experiments

Applicable to diverse financial functions and instruments

Efficient recursive computation method

Abstract

In quantitative finance, it is often necessary to analyze the distribution of the sum of specific functions of observed values at discrete points of an underlying process. Examples include the probability density function, the hedging error, the Asian option, and statistical hypothesis testing. We propose a method to calculate such a distribution, utilizing a recursive method, and examine it using various examples. The results of the numerical experiment show that our proposed method has high accuracy.

Equations

$$Y=\sum _{i=1}^{N}h(X_{i-1},X_i)$$

$$f_{Y|X_0}(y|x_0)=\frac{{\mathrm{d}}^2\mathbb{E}[(Y-y{)}^{+}|X_0=x_0]}{\mathrm{d}{y}^2}.$$

$$\mathbb{E}[(Y-y{)}^{+}]=\mathbb{E}\left[{\left(\sum _{i=1}^{N}h(X_{i-1},X_i)-y\right)}^{+}\right],$$

$$g_n(y|x_n)=\mathbb{E}\left[{\left(\sum _{i=n+1}^{N}h(X_{i-1},X_i)-y\right)}^{+}|X_n=x_n\right].$$

$$g_n(y|x_n)={\int }_{\mathbb{R}}{g}_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})f(x_{n+1}|x_n)\mathrm{d}{x}_{n+1}$$

$$g_{N-1}(y|x_{N-1})$$

$$={\int }_{\mathbb{R}}(h(X_{N-1},X_N)-y{)}^{+}f(x_N|x_{N-1})\mathrm{d}{x}_N.$$

$$g_0(y|x_0)=\mathbb{E}[(Y-y{)}^{+}|X_0=x_0],$$

$$f_{Y|X_0}(y|x_0)=\frac{{\partial }^2{g}_0(y|x_0)}{\partial y^2}.$$

$${\bar{F}}_n(y|x_n)={\int }_{\mathbb{R}}{\bar{F}}_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})f(x_{n+1}|x_n)\mathrm{d}{x}_{n+1}$$

$${\bar{F}}_{N-1}(y|x_{N-1})={\int }_{\mathbb{R}}-{\text{mathbbm}1}_{\{y<h(x_{N-1},x_N)\}}f(x_N|x_{N-1})\mathrm{d}{x}_N$$

$$F_n(y|x_n)={\int }_{\mathbb{R}}{F}_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})f(x_{n+1}|x_n)\mathrm{d}{x}_{n+1}$$

$$F_{N-1}(y|x_{N-1})={\int }_{\mathbb{R}}\text{mathbbm}{1}_{\{y\ge h(x_{N-1},x_N)\}}f(x_N|x_{N-1})\mathrm{d}{x}_N.$$

$$f_n(y|x_n)={\int }_{\mathbb{R}}{f}_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})f(x_{n+1}|x_n)\mathrm{d}{x}_{n+1}$$

$$f_{N-1}(y|x_{N-1})={\int }_{\mathbb{R}}\frac{\partial }{\partial y}\text{mathbbm}{1}_{\{y\ge h(x_{N-1},x_N)\}}f(x_N|x_{N-1})\mathrm{d}{x}_N$$

$$g_{n+1}(y|x_{n+1})=x_{n+1}g_{n+1}\left(\frac{y}{x_{n+1}}\Big{|}1\right)=x_{n+1}\bar{g}_{n+1}\left(\frac{y}{x_{n+1}}\right)$$

$$g_n(y|x_n)$$

$$={\int }_{\mathbb{R}}{x}_{n+1}{\bar{g}}_{n+1}\left(\frac{y-h(x_n,x_{n+1})}{x_{n+1}}\right)f(x_{n+1}|x_n)\mathrm{d}{x}_{n+1}.$$

$${\bar{g}}_n(y)={\int }_{\mathbb{R}}{x}_{n+1}{\bar{g}}_{n+1}\left(\frac{y-h(1,x_{n+1})}{x_{n+1}}\right)f(x_{n+1}|1)\mathrm{d}{x}_{n+1}$$

$${\bar{g}}_{N-1}(y)=\mathbb{E}[(h(1,X_N)-y{)}^{+}|X_{N-1}=1].$$

$$F_n(y)={\int }_{\mathbb{R}}{F}_{n+1}\left(\frac{y-h(1,x_{n+1})}{x_{n+1}}\right)f(x_{n+1}|1)\mathrm{d}{x}_{n+1}$$

$$F_{N-1}(y)={\int }_{\mathbb{R}}\text{mathbbm}{1}_{\{y\ge h(1,x_N)\}}f(x_N|1)\mathrm{d}{x}_N.$$

$$g_0(y|x_n)=\mathbb{E}\left[{\left(\sum _{i=1}^N{X}_i-y\right)}^{+}|X_0=x_0\right],$$

$${\int }_{\mathbb{R}}{\mathrm{e}}^{-\mathrm{i}\nu y}{f}_{N-1}(y|x_{N-1})\mathrm{d}y$$

$${\int }_{\mathbb{R}}{\mathrm{e}}^{-\mathrm{i}\nu y}{f}_n(y|x_n)\mathrm{d}y$$

$$[x_{\text{min}},x_{\text{max}}]\times [y_{\text{min}},y_{\text{max}}]$$

$$F_n(y|x_n)={\int }_{\mathbb{R}}{F}_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})f(x_{n+1}|x_n)\mathrm{d}{x}_{n+1}.$$

$$F_{n+1}(y-h(x_n,x_{n+1})|\cdot )=1$$

$$F_{n+1}(y-h(x_n,x_{n+1})|\cdot )=0$$

$$F_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})=F_{n+1}(y-h(x_n,x_{n+1})|x_{\text{max}})$$

$$F_{n+1}(y-h(x_n,x_{n+1})|x_{n+1})=F_{n+1}(y-h(x_n,x_{n+1})|x_{\text{min}})$$

Full text

1 Introduction

This paper introduces a recursive method to compute interesting quantities related to probability distributions in various financial applications. The method is versatile, and hence, with slight modifications, it is easy to apply the basic framework to various applications. More precisely, the method is based on a convolution-like formula, applied to compute the distribution of the sum of values of one-dimensional processes observed at discrete points. Financial applications include numerical densities of asset price or volatility models, hedging error distributions, arithmetic Asian option prices, and statistical hypothesis tests.

Various kinds of stochastic processes are used in quantitative finance, such as the Cox-Ingersoll-Ross model (CIR), the constant elasticity of variance model (CEV), stochastic volatility, and GARCH models. The probability distributions of the processes in financial models can be used for risk management, asset pricing, hedging analysis, parameter estimation, and statistical hypothesis testing. In many cases, the closed form formulas for the density function of stochastic models are not known, and it is advantageous to develop a numerical procedure to compute the probability distributions or density functions.

When trading a financial option, the investor usually performs a hedging procedure to reduce risk. In general, continuous models of asset price movements assume a continuous hedging process. However, in practice, because continuous trading is not applicable, a discrete time hedging strategy is applied, and hence, the discrete time hedging error occurs even in the complete market model. Many studies examine the discrete time hedging error in financial options. Sepp, (2012) derived a numerically approximated distribution of the delta-hedging error based on the characteristic function for a jump diffusion model. Park et al., (2016) computed the delta-hedging error based on a recursive method for a jump diffusion model, which is the same framework this study proposes. This study extends this to a Lévy model and shows it is possible to easily adapt the method to not only delta-hedging processes but also other trading strategies, such as minimum variance hedging.

The arithmetic Asian option is a financial derivative whose payoff is the arithmetic average of the underlying asset prices observed at future times. Asian options are safer with respect to the manipulation of underlying asset prices that may occur when they are close to maturity than European options and financial instruments suitable for less frequently traded assets (Musiela and Rutkowski,, 2006). Since the closed-form formula is not available for the Asian option price, we study numerical approximation and simulation methods (Kemna and Vorst,, 1990; Věcěr,, 2002). Our example is consistent with Lee, (2014), which computes European option based Asian option prices; however, we directly apply risk-neutral probability density in this study.

We also examine an example of the statistical hypothesis test. In general, a parametric statistical test depends on the probability distribution of a test statistic. For typical sample mean tests, the test statistics are generally approximated by a t distribution; however, if the corresponding random variable is far from the normal distribution, it might undermine the accuracy of the test. Therefore, a more exact distribution will be helpful in performing a more reliable test. We provide an example of a skewness test where the recursive method is applied. In general, financial asset return distributions are negatively skewed (Fama,, 1965; French et al.,, 1987; Cont,, 2001) and the third moment of financial asset distribution has been extensively studied (Kraus and Litzenberger,, 1976; Harvey and Siddique,, 2000; Christoffersen et al.,, 2006; Choe and Lee,, 2014; Lee,, 2016). Our example demonstrates a method to compute critical values and statistical power.

The rest of the paper is organized as follows: Section 2 explains the basic recursive method and the process of applying the numerical procedure to compute the probability density functions. Section 3 applies the proposed method to the examples. Section 4 concludes the paper.

2 Basic method

2.1 Derivation

Let $X$ be a continuous stochastic process defined on the time horizon $[0,T]$ or a discrete stochastic process defined on time indexes $0=t_{0},t_{1},\cdots,t_{N}=T$. If $X$ is continuous, we are particularly interested in the behaviors of $X_{t_{i}}$ for the discrete observation times $0=t_{0},t_{1},\cdots,t_{N}=T$. When necessary, we can introduce a complete filtered probability space $(\Omega,\mathcal{F},\mathbb{P})$ over $[0,T]$, with filtration $\{\mathcal{F}_{t}\}_{t\in[0,T]}$.

In numerous financial applications, it is advantageous to examine the distribution of some finite summation of

[TABLE]

for some function $h$, where $X_{i}=X_{t_{i}}$. A subsequent section will explain specific examples in financial practice, including a numerical probability density function, arithmetic Asian option price, the distribution of the realized variance, the hedging cost distribution, and statistical hypothesis testing.

To compute the distribution of $Y$, we propose a numerical scheme based on a recursive relationship. One method to represent the conditional probability density function of $Y$ given $X_{0}=x_{0}$ is based on the second derivative of the expectation of the European option payoff or the rectified unit linear function:

[TABLE]

This approach is in line with the method introduced in Breeden and Litzenberger, (1978), which derived the state price density function, which is similar to the risk-neutral density function, based on European option prices. When the expectation is under a risk-neutral measure, the density function $f_{Y|X_{0}}$ is also considered to be the risk-neutral density function. We are interested in both physical and risk-neutral probabilities. This study examines the numerical method for computing the conditional probability density function $f_{Y|X_{0}}$ based on Eq. (1).

To calculate

[TABLE]

consider the following relationship. Define

[TABLE]

The above formula is similar to the $\mathcal{F}_{n}$-conditional expectation of the European option payoff, by regarding $\sum_{i=n+1}^{N}h(X_{i-1},X_{i})$ as an asset price and $y$ as a strike price. Let $f_{X_{n+1}|X_{n}}(x_{n+1}|x_{n})$ or simply $f(x_{n+1}|x_{n})$, be the transition probability density function from $X_{n}=x_{n}$ to $X_{n+1}=x_{n+1}$. Then, we derive the following successive relationships for $0\leq n<N-1$, as follows:

[TABLE]

and

[TABLE]

Since

[TABLE]

we have

[TABLE]

Or, by differentiating both sides of Eq. (2) for $y$, we can express

[TABLE]

where $\bar{F}_{n}(y|x_{n})=\frac{\partial g_{n}(y|x_{n})}{\partial y}$, and similarly,

[TABLE]

where $\mathbbm{1}$ denotes the indicator function. In the above equations, $\bar{F}_{n}=-1+F_{n}$, where $F_{n}(y|x_{n})$ is a conditional cumulative distribution function of $Y_{n}:=\sum_{i=n+1}^{N}h(X_{i-1},X_{i})$, given $X_{n}=x_{n}$. Based on the discussion so far, we define the following.

Definition 1**.**

We define

[TABLE]

and

[TABLE]

In addition, by differentiating both sides of (4) for $y$, we obtain

[TABLE]

and

[TABLE]

where $f_{n}(y|x_{n})=\frac{\partial F_{n}(y|x_{n})}{\partial y}$, and the derivative is distributional.

Remark 1**.**

If $X$ is scale invariant, that is, $g_{n}(y|x_{n})=x_{n}g_{n}\left(\frac{y}{x_{n}}|1\right)$, and let \bar{g}_{n}(y)=g_{n}(y\big{|}1), we obtain

[TABLE]

and

[TABLE]

By setting $x_{n}=1$,

[TABLE]

and

[TABLE]

Alternatively, we can use the cumulative distribution function to express

[TABLE]

and

[TABLE]

The above method applies to the one-dimensional function, and hence, the computational cost is much less.

One straightforward example is an arithmetic Asian option price under geometric Brownian motion. In this case,

[TABLE]

which can be interpreted as the Asian option price at time [math], with underlying price $X_{0}=x_{0}$, and strike price $y$. Note that the price is equal to $x_{0}$ times the Asian option price with underlying price $X_{0}=1$ and strike price $y/x_{0}$. In other words, $g_{0}(y|x_{0})=x_{0}g_{0}\left(\frac{y}{x_{0}}|1\right)$, and this is also applied to every $n$.

Remark 2**.**

We obtain the intuitive forms by applying the Fourier transform to Eqs. (6) and (7), and changing the order of integrals. By definition, the Fourier transform of the left-hand side of Eq. (7) is the expectation of $\exp(-\mathrm{i}\nu Y_{N-1})$.

[TABLE]

Similarly, the Fourier transform of the left side of Eq. (6) is the expectation of $\exp(-\mathrm{i}\nu Y_{n})$.

[TABLE]

In the second equality, we substitute $z$ for $y-h(x_{n},x_{n+1})$.

2.2 Numerical procedure

This subsection describes the numerical algorithm when applying the recursive method in computing a probability distribution. There are specific considerations in applying the numerical method for every application; however, in this subsection, we examine the common concerns in applying the computational procedure.

Selection of object function

For the numerical procedure, we determine whether to use function $F_{n}$ or $f_{n}$, although this distinction has no significant effect on the results. Both $F_{n}$ and $f_{n}$ have nice properties that stabilize the numerical procedure. Since $F_{n}$ is theoretically bounded between 0 and 1, it can be easily corrected, even if $F_{n}$ is outside the bounded region, owing to numerical error. The density function $f_{n}$ converges to 0 as $y$ goes $\infty$ or $-\infty$; we can assume that integration with large $y$ in absolute value is almost zero. This property of $f_{n}$ tends to make implementation easier when conducting numerical procedures; therefore, many examples in this study were based on $f_{n}$. Meanwhile, in Subsection 3.2, a singular point in the density function, such as in Dirac measure, makes it easier to use $F$ in the numerical procedure.

Adaptive meshing

For $F_{n}(y|x_{n})$ or $f_{n}(y|x_{n})$, the numerical domain is bounded by the region

[TABLE]

where $x_{\textrm{min}}$ and $x_{\textrm{max}}$ are generally fixed, whereas $y_{\textrm{min}},y_{\textrm{max}}$ are usually dynamic throughout the iteration. The adaptive change in the domain of $y$ during the numerical procedure is due to the change in the reasonable numerical support of the conditional distribution of $Y$ throughout the iteration.

For example, if the values of $y$ tend to increase as the numerical procedure proceeds, the numerical domain of $y$ changes accordingly. Let us suppose our interest lies in the distribution of $Y=\sum_{i=1}^{N}X_{i}$. The sufficient numerical domain for $\sum_{i=n}^{N}X_{i}$ is generally larger than $\sum_{i=n+1}^{N}X_{i}$. Therefore, it is natural to expand the domain of $[y_{\textrm{min}},y_{\textrm{max}}]$ for $F_{n}$, as $n$ proceeds from $N$ to 1.

If the numerical domain of $y$ is expanding, the number of intervals that divide the domain could become too large. Instead of increasing the number of intervals over $y$, we fix the total number of the discretized points over $[y_{\textrm{min}},y_{\textrm{max}}]$ to prevent the grid size from becoming too large. The dynamic allocation algorithm is straightforward. Since $F_{n}$ or $f_{n}$ converges to 0 or 1, as $y$ approaches $\infty$ or $-\infty$ for all $x_{n}$, with given threshold $\epsilon$, we increase $y_{\max}$ or decrease $y_{\min}$ as $F_{n}$ or $f_{n}$ reaches the convergence criteria; for example, $f(y_{\max})<\epsilon$. Since we fix the number of discretized points over $[y_{\textrm{min}},y_{\textrm{max}}]$, say $M$, the step size in $y$, $\Delta y=(y_{\textrm{max}}-y_{\textrm{min}})/M$ this too changes, as $y_{\textrm{min}},y_{\textrm{max}}$ change.

Referencing previous function

Another concern is referencing previous function values of $F_{n+1}$ or $f_{n+1}$ at step $n$, as we compute

[TABLE]

When query point $(y-h(x_{n},x_{n+1}),x_{n+1})$ falls within the grid $[y_{\textrm{min}},y_{\textrm{max}}]\times[x_{\textrm{min}},x_{\textrm{max}}]$ defined at time $t_{n+1}$, we simply retrieve values such as $F_{n+1}(y-h(x_{n},x_{n+1})|x_{n+1})$ or $f_{n+1}(y-h(x_{n},x_{n+1})|x_{n+1})$ using interpolations, such as linear or piecewise cubic Hermite interpolations. This works well even using the nearest value.

When $y-h(x_{n},x_{n+1})$ is outside $[y_{\textrm{min}},y_{\textrm{max}}]$, we assign the extrapolated value depending on the function, $F_{n+1}$ or $f_{n+1}$. Since $F_{n+1}$ is a cumulative distribution function, $F_{n+1}(y|\cdot)$ approaches 1 when $y$ approaches $\infty$ and $F_{n+1}(y|\cdot)$ approaches [math] when $y$ approaches $-\infty$. Therefore, it is natural to assign

[TABLE]

when $y-h(x_{n},x_{n+1})>y_{\textrm{max}}$ and

[TABLE]

when $y-h(x_{n},x_{n+1})<y_{\textrm{min}}$. Similarly, for $f_{n+1}$, we assign $f_{n+1}(y-h(x_{n},x_{n+1})|\cdot)=0$ when $y-h(x_{n},x_{n+1})>y_{\textrm{max}}$ or $y-h(x_{n},x_{n+1})<y_{\textrm{min}}$.

When $x_{n+1}$ is outside $[x_{\textrm{min}},x_{\textrm{max}}]$ and $y-h(x_{n},x_{n+1})\in[y_{\textrm{min}},y_{\textrm{max}}]$, we assign

[TABLE]

when $x_{n+1}>x_{\textrm{max}}$ and

[TABLE]

when $x_{n+1}<x_{\textrm{min}}$. In this case, $x_{n+1}$ is far from $x_{n}$ and the transition probability $f(x_{n+1}|x_{n})$ is relatively small, and hence, $F_{n+1}(y-h(x_{n},x_{n+1})|x_{n+1})f(x_{n+1}|x_{n})$ is close to zero. Therefore, error due to discrepancy between the true value of $F_{n+1}(y-h(x_{n},x_{n+1})|x_{n+1})$ and its approximation is small, relative to the overall integration. We summarize this explanation in Figure 1 and present the algorithm in Algorithm 1333For an example of Matlab code, see https://github.com/ksublee/Recursive_method.

2.3 Computational cost

The computational cost depends on the number of time steps and the size of grid $[x_{\textrm{min}},x_{\textrm{max}}]\times[y_{\textrm{min}},y_{\textrm{max}}]$. To compute $f_{n}(y|x)$, we perform numerical integration for every point $(x,y)$ in the grid. Let the numbers of $x$ and $y$ in the grid be $N_{x}$ and $N_{y}$, respectively, and the number of steps for numerical integration in Eqs (4) or (6) be $N_{z}$. Then, the computational complexity of each time step is proportional to $N_{x}N_{y}N_{z}$.

This is similar to an alternative method such as the Fourier transform described in Remark LABEL:remark:Fourier. More precisely, let $\xi_{n}=\mathbb{E}[\mathrm{e}^{-\mathrm{i}\nu Y}|\mathcal{F}_{n}]$. Then, the equation corresponding to the recursive relationship is

[TABLE]

, and hence, we perform numerical integration with respect to $x_{n}$ on every point of $(x_{n-1},\nu)$ over a grid $[x_{\textrm{min}},x_{\textrm{max}}]\times[\nu_{\textrm{min}},\nu_{\textrm{max}}]$, for every time step. Therefore, essentially, the Fourier transform and the recursive methods have the same time complexity. One advantage of our method over the Fourier method is that we do not have to apply Fourier transform in the final step to retrieve the distribution or density function.

3 Application

In this section, we apply our proposed method to several examples. As we introduce various distinct examples, please note that each subsection uses different notations.

3.1 Numerical density for diffusion models

This subsection demonstrates the computation of probability densities, or likelihood functions of various diffusion models numerically based on the recursive method.

3.1.1 CIR model

Consider a square-root process $X$, also known as the Cox-Ingersoll-Ross model (Cox et al.,, 1985), defined by:

[TABLE]

The density function of the transition probability from $X_{0}=x_{0}$ to $X_{t}=x$ of the square root process is given by

[TABLE]

where

[TABLE]

and $I_{q}$ denotes the modified Bessel function of the first kind of order $q$.

Although the closed-form solution of the density function is available, for illustrative purposes, we examine the approximation method based on recursive relation and discretization to compute the probability density function of $X_{N}$ for some $t_{N}$. The approximate distribution of $\Delta X_{n}=X_{n}-X_{n-1}$ with $\Delta t=t_{n}-t_{n-1}$ by normal distribution in typical Monte Carlo simulations is as follows:

[TABLE]

Let

[TABLE]

then

[TABLE]

and the recursive method can be applied to compute the distribution of $Y$, that is, $X_{N}-X_{0}$. For the integrand of the recursive method, we use the conditional probability density function $f_{n}(y|x_{n})$. We obtain similar results using the conditional distribution function $F_{n}(y|x_{n})$.

Figure 2 compares the numerically computed density function, closed-form formula, and simulation results. The parameter settings for the square root process are $\kappa=11,\theta=0.2,\gamma=1.5$ and for the numerical procedure $\Delta t=1/1250,N=100,[x_{\min},x_{\max}]=[0,0.6]$ with $\Delta x=0.002$; The number of intervals for the $y$ domain are 1,000. Time is annualized, and hence, $t_{N}=20$ days. On the left side of Figure 2, the numerically computed density function based on the recursive method is very close to the closed-form formula and simulation histogram. On the right side of Figure 2, the conditional density $f_{Y|X_{0}}$ is plotted as a function of $X_{0}=x_{0}$ and $Y=X_{N}-X_{0}=y$. The global error is decreasing with the increasing number of intervals for the $y$-axis as plotted in Fig 3.

3.1.2 CEV model

The constant elasticity of variance model (CEV) identifies the leverage effect of volatility being negatively correlated with asset prices (Cox and Ross,, 1976) with the following formula in the stock price process:

[TABLE]

The closed-form formula of the stock price distribution is not known and it is worthwhile to compute the density function with the numerical method. The conditional probability density function is computed using the same method as the square root process in the previous subsection.

The conditional probability density functions of $Y$, $f_{Y|X_{0}}(y|x_{0})$ with various $x_{0}$ (right) and the density function of $X_{t}$ with $X_{0}=1$ is presented in Figure 4. Comparing the numerical probability density function and the simulation histogram shows that the recursive method generates a more precise density function. The parameter setting for the CEV model example is $\mu=0.05,\sigma=0.2,\gamma=0.7$. For the numerical procedure, $\Delta t=1/1250,N=200,[x_{\min},x_{\max}]=[0.5,1.5]$ with $\Delta x=0.005$, the number of intervals on the $y$-axis are 200 and the tolerance for the dynamic allocation is $10^{-8}$.

3.1.3 Stochastic volatility model and integrated variance

This subsection computes the probability density function based on the numerical likelihood method for various stochastic volatility models. Consider a stochastic volatility model such that

[TABLE]

with parameters of $a\in\{0,1\}$ and $b\in\{1/2,1,3/2\}$. This classification is from Christoffersen et al., (2010). When $a=0$ and $b=1/2$, the square root process is used for stochastic volatility, as in Heston, (1993).

Other than $a$ and $b$, we fix the parameter setting $\kappa=11,\theta=0.2,\gamma=0.8,V_{0}=0.2$, and for the numerical procedure, $\Delta t=1/1250,N=100$. Figure 5 presents the numerically computed various probability density functions. For stochastic volatility models, the numerical probability density functions are close to the simulation histograms.

The integrated variance is defined by

[TABLE]

It is known that the realized variances in the stochastic volatility models converge to the integrated variances (Barndorff-Nielsen,, 2002). However, the closed-form formulas of the unconditional distribution for integrated variances of stochastic volatility are generally not known. For discussion on special cases of affine models, see Broadie and Kaya, (2006). Since the variance process is unobservable, sometimes the integrated variance is more interesting, as it can be approximated by the quadratic process of underlying return.

Using the recursive method with $h(x_{n},x_{n+1})=x_{n}$, and approximating

[TABLE]

the probability density functions of the integrated variance under various stochastic models are presented in Figure 6. The basic setting is the same as in previous cases, except for the formula of $h$. Compared with the simulation results, the numerical density functions are precise for all stochastic volatility models.

3.2 The GARCH model

We can calculate the numerical density function for not only continuous models but also discrete time models. Consider the GARCH model (Bollerslev,, 1986) for the variance $\sigma^{2}$ and log-return $\varepsilon$:

[TABLE]

We are interested in the distributions of both time $N$ variance, $\sigma^{2}_{N}$, and total return, $\sum_{i=1}^{N}\varepsilon_{i}$.

The conditional distribution of $\alpha\varepsilon^{2}_{n}$ with given $\sigma_{n}$ is represented by the gamma distribution, such that

[TABLE]

where the first argument of $\Gamma(,)$ is the shape parameter and the second argument is the scale parameter in the gamma distribution. Therefore, the transition probability from $\sigma^{2}_{n}$ to $\sigma^{2}_{n+1}$ can be represented by the shifted gamma probability density function such that

[TABLE]

First, let us examine the distribution of $\sigma^{2}_{N}$ with given $\sigma_{0}$. Since the closed form of the transition probability density $f_{\sigma^{2}_{n+1}|\sigma^{2}_{n}}$ is available, by setting $h(x_{n},x_{n+1})=x_{n+1}-x_{n}$, and hence, $Y=\sigma^{2}_{N}-\sigma^{2}_{0}$, the theory is based on the recursive method.

Since the probability density function contains a non-smooth point, it is preferable to use the cumulative distribution function, $F_{n}$, for the recursive procedure. Figure 7 shows the histogram of the simulated GARCH process and numerically computed density function of $\sigma^{2}_{N}$ with $N=20$. For the parameter settings, we utilize $\omega=0.001$, $\alpha=0.05$ and $\beta=0.9$. For the numerical procedure, $x$-grid is set to be $[x_{\textrm{min}},x_{\textrm{max}}]=[0.01,0.05]$ with $\Delta x=0.0002$ and the number of intervals for the $y$-axis is 150. The tolerance level for $y$ is $10^{-8}$. The right side of the figure presents the typical shape of $F_{0}$. For any starting point $x_{0}$, along the $x$-axis, the conditional cumulative distribution function of $Y$ is retrieved along the $y$-axis.

Second, we set $h(x_{n},x_{n+1})=x_{n}/N$, $Y=\sum_{i=1}^{N}\sigma^{2}_{n}/N$, and compute $\sum_{i=1}^{N}\varepsilon_{i}$. Since

[TABLE]

follows the conditional normal distribution, we can compute the total log-return from time $t_{0}$ to time $t_{N}$. Figure 8 presents the results. On the left side is the probability density function of $\sum_{i=1}^{N}\sigma_{n}^{2}/N$ and the probability density function of $\sum_{i=1}^{N}\varepsilon_{i}$ is on the right side. The figures show that the numerically computed GARCH return distribution is close to the simulation result.

3.3 Hedging error

Many researchers such as Sepp, (2012), have studied errors occurring from hedging strategies of an option under time-discretization. We consider hedging errors in a framework in which the underlying process follows the exponential Lévy model, as in Madan et al., (1998). This subsection can be regarded as an extension of Park et al., (2016).

Let $\gamma_{t}$ be a gamma process, which is a Lévy process with independent and gamma distributed increments, with a mean rate parameter of 1 and variance parameter $\nu$. In other words, the Lévy measure of $\gamma$ is represented by $\nu z^{-1}\exp(-z/\nu)$ for jump size $z$. Assume that the price of the underlying asset under a risk-neutral measure follows the exponential variance gamma model:

[TABLE]

where $X_{t}$ is a variance gamma process represented by a time changed Brownian motion

[TABLE]

and $\omega=(1/\nu)\log(1-\theta\nu-\frac{1}{2}\sigma^{2}\nu)$. The Lévy measure of the variance gamma process is represented by

[TABLE]

Madan et al., (1998) derived the probability density function of $S_{t}$, with $S_{0}=1$, as follows:

[TABLE]

From the above formula, the transition probability density function from $S(t_{n})=x_{n}$ to $S(t_{n+1})=x_{n+1}$ can be computed by numerical integration.

There are several practical methods to compute the European call option price under the variance gamma process. We use the fast Fourier transform method based on the dampened option price as follows:

[TABLE]

for some constant, $\alpha$, as explained in Carr and Madan, (1999). With this setting, the European call option price with $S_{0}=1$ and maturity, $T$, is represented by

[TABLE]

where $\psi_{T}(u)=\mathbb{E}^{\mathbb{Q}}[\mathrm{e}^{\mathrm{i}u\log S_{T}}]$.

We test two ways of hedging the European call option $C$: delta hedging and minimal variance hedging, as proposed in Föllmer and Sondermann, (1986). The minimal variance hedging ratio of the European call option, $C$, under the variance gamma model is represented by

[TABLE]

For a more detailed explanation, consult Cont et al., (2007).

Consider an investor with a short position in the call option and hedging strategy with long position in the underlying asset. Under the trading strategy $\phi$, the trading error between time $t_{i}$ and $t_{i+1}$ is defined by the difference between $t_{i+1}$-realized value of the hedged portfolio and the price of the risk-free asset as follows:

[TABLE]

The total error is

[TABLE]

Similarly, for the delta hedging strategy, the total error is

[TABLE]

where $D$ denotes the delta hedging ratio. Since the terms inside the summation in Eqs. (9) and (10) are functions of the underlying process, we can apply the numerical recursive method to compute the distribution of the hedging errors.

Figure 9 compares the numerically computed probability density functions of delta hedging errors (left) and minimum variance hedging errors (right) of European call options with the simulation histograms. The European call options have strike price $K=0.9,1,1.05$ with $S_{0}=1$ listed from top to bottom in Figure 9. The parameter setting is $\sigma=0.2,\theta=1.2,\nu=0.001,r=0.02,\delta t=1/250$.

3.4 Arithmetic Asian option

As mentioned, an Asian option is a financial derivative that is more robust to the manipulation of underlying asset prices than the European option. Since the closed form formula for arithmetic Asian option prices is not known, simulation, approximation, or computational methods are generally used. The recursive method can also be applied to compute the arithmetic Asian option prices.

Let $S$ be an underlying asset price process. The arithmetic Asian option price over observation points $t_{1},\cdots,t_{N}=T$ with strike price $K$ is represented by

[TABLE]

where the expectation is in terms of risk-neutral probability. To compute the above expectation, we need to compute the distribution of $Y=\frac{1}{N}\sum_{i=1}^{N}S_{i}$. Using the recursive method, setting $h(x_{i-1},x_{i})=x_{i}/N$, or $h(x_{i-1},x_{i})=x_{i}$, and subsequently, rescaling distribution, we can compute the risk-neutral distribution of $Y$. The Asian option price is determined using the recursive method whenever the risk-neutral transition probability is available. The following can be considered an alternative method to that proposed in Lee, (2014).

We assume that the stock price process follows a variance Gamma process, as in Eq (8). The parameter settings are $\sigma=0.2,\theta=1.2,\nu=0.001,r=0.02$, and $S_{0}=1$. The overall method is the same as in subsection 3.3, except the function form of $h$. The comparison between the numerical method and simulation is presented in Figure 10 and the two results are quite similar.

3.5 Skewness test

In this section, we examine how the proposed recursive method can be used in hypothesis testing of the third moment of the asset return. The distributions of asset returns tend to be skewed to the left. Although it is generally not easy to measure the exact third moment of the return distribution, the importance of the third moment has been acknowledged and extensively studied (Kraus and Litzenberger,, 1976; Harvey and Siddique,, 2000).

Let $R$ be a stationary return process. For the hypothesis testing, the null hypothesis is $\mathrm{H}_{0}:\mathbb{E}[(\Delta R)^{3}]=0$ and the alternative hypothesis is $\mathrm{H}_{1}:\mathbb{E}[(\Delta R)^{3}]<0$. For example, consider a jump diffusion model as follows:

[TABLE]

where $\mu$ is drift, $\sigma$ is volatility, $N$ is a Poisson process with intensity $\lambda$, and $J$ follows a normal distribution with mean $\mu_{J}$ and standard deviation $\sigma_{J}$. For simplicity, we set $\mu=-\lambda\mu_{J}$, such that $R$ becomes a martingale (with respect to suitable filtration). Under this assumption, the statistical hypothesis can be modified as follows: $\mathrm{H}_{0}:\mu_{J}=0$ and $\mathrm{H}_{1}:\mu_{J}<0$. The test is similar to a simple t-test; however, the distribution of $(\Delta R)^{3}$ does not follow the normal distribution, and hence, it is advantageous to compute the exact distribution of $(\Delta R)^{3}$ using the recursive method.

We compute the critical values that determine the rejection of the null hypothesis for sample sizes with given significance level, $\alpha=0.05$ under the null hypothesis (see Figure 11). The null hypothesis is rejected when the test statistic; the sample mean of $(\Delta R)^{3}$ is less than the corresponding critical value with given sample size. To compute the critical values, the numerical probability density function of $(\Delta R)^{3}$ is computed to a $h(x_{n},x_{n+1})=(x_{n+1}-x_{n})^{3}$ and the presumed parameter settings of $\lambda=10,\sigma_{J}=0.01,\mu_{J}=0$ and $\sigma=0.1975$.

The statistical power, typically denoted by $1-\beta$, is the probability that the test correctly rejects the null hypothesis when the null hypothesis is invalid. We examine a power curve versus sample size where $\mu_{J}$ is presumed to be $-0.05$, to imply negative skewness, and the other parameters are the same as in the previous case. The right side of Figure 11 shows the increase in the power curve with increasing sample size. The curve implies that if we seek 90% statistical power, this model will need approximately 90 samples. This section presents an example of the third moments test, and the recursive method is deemed to be applicable for various statistical tests.

4 Conclusion

This study proposed a recursive formula for the distribution of specific functions and detailed the application of the numerical procedure. Various examples, including the numerical density function, the hedging error distribution, the arithmetic Asian option pricing, and statistical hypothesis testing, showed that the proposed method is quite precise. The method is versatile, and we expect that more applications will become available not only in finance but also in various probabilistic analysis. This study applied the method to a one-dimensional model, and future studies can extend the method to the two-dimensional process model.