A weak approximation for the Wiener–Hopf factorization

Abstract

The Wiener–Hopf factorization plays a crucial role in studying various mathematical problems. Unfortunately, in many situations, the Wiener–Hopf factorization cannot provide closed form solutions and one has to employ some approximation techniques to find its solutions. This article provides several weak, approximation for a given Wiener–Hopf factorization problem. Application of our finding in spectral factorization and Lévy processes have been given.

Keywords:

• principal value integral
• Hölder condition
• Padé approximant
• continued fraction
• Fourier transform
• Hilbert transform
• Shannon sampling theorem
• spectral factorization
• Lévy processes

AMS Subject Classifications:

• 30E25
• 11A55
• 42A38
• 60G51
• 60j50
• 60E10

Public Interest Statement

This article provides several weak approximation for a given Wiener–Hopf factorization problem. Application of our finding in spectral factorization and Lévy processes have been given.

1. Introduction

Roughly speaking, the Wiener–Hopf factorization problem is a technique to find a single complex-valued function $\mathrm{\Phi }$ in which its radial limits, say ${\mathrm{\Phi }}^{\pm }$, are respectively analytic and bounded separately in the upper and lower complex half planes (i.e. $C_{+}:=\{\mathit{\lambda }\in C:\mathfrak{I}(\mathit{\lambda })\ge 0\}$ and $C_{-}:=\{\mathit{\lambda }\in C:\mathfrak{I}(\mathit{\lambda })\le 0\}$) and satisfy ${\mathrm{\Phi }}^{+}(\mathit{\omega }){\mathrm{\Phi }}^{-}(\mathit{\omega })=g(\mathit{\omega })$, where $\mathit{\omega }\in R$ and $g$ is a zero index function which satisfies the Hölder condition.

The Wiener–Hopf factorization has proved remarkably useful in solving an enormous variety of model problems in a wide range of branches of physics, mathematics, and engineering. Subjects for which the problem is applicable range from neutron transport (Noble, 1988), geophysical fluid dynamics (Davis, 1987; Kaoullas & Johnson, 2010), diffraction theory (Noble, 1988), fracture mechanics (Freund, 1998; Shakib, Akhgarian, & Ghaderi, 2015), non-destructive evaluation of materials (Achenbach, 2012), a wide class of integral equations (Payandeh Najafabadi & Kucerovsky, 2009,2014b), acoustics (Abrahams & Wickham, 1990), elasticity (Norris & Achenbach, 1984; Ogilat, 2013), electromagnetics (Daniele, 2014; Sautbekov & Nilsson, 2009), water wave phenomena (Chakrabarti & George, 1994; Kim, Schiavone, & Ru, 2011), geophysics (Davis, 1987), financial mathematics (Beheshti, Payandeh Najafabadi, & Farnoosh, 2013; Fusai, Abrahams, & Sgarra, 2006), distribution of extrema in a wide class of Lévy processes (Payandeh Najafabadi & Kucerovsky, 2011,in press), statistical decision problems (Kucerovsky, Marchand, Payandeh Najafabadi, & Strawderman, 2009), etc.

The key steps to solve a Wiener–Hopf factorization is decomposing of the kernel $g$ into a product of two terms, $g^{+}$ and $g^{-}$, where $g^{+}$ and $g^{-}$ are analytic and bounded in the upper and the lower complex half planes, respectively. Such decomposition can be expressed in terms of a Sokhotski–Plemelj integral (see Equation, 1), but this form presents some difficulties in numerical work due to slow evaluation and numerical problems caused by singularities near the integral contour (see Kucerovsky & Payandeh Najafabadi, 2009, for more details). To overcome these problems, several approximation methods have been considered (see Abrahams, (2000); Kudryavtsev & Levendorskiĭ, 2009; Kuznetsov, 2010; Rawlins, 2012 among others). But, as far as we known, ( i) none of them provides any estimation bound for their approximation methods; (ii) most of them need uniform convergence, which is usually hard to achieve.

This article studies the problem of solving a Wiener–Hopf factorization problem, approximately. Then, it provides (i) an $L_p(R)$, $1<p\le 2$, approximation for a Wiener–Hopf factorization problem; (ii) estimation bounds for such approximation technique; (iii) application of our findings in spectral factorization and Lévy processes. This article has been developed as the following. Section 2 collects some useful elements which are used later. The main contribution of this article on approximating solutions of a given a Wiener–Hopf factorization problem has been given in Section 3. Application of our findings has been given in Section 4. Concluding remarks has been given in Section 5.

2. Preliminaries

Now, we collect some lemmas which are used later.

Definition 1

A function $f$ in $L_1(R)\cap L_2(R)$ is said to be an exponential-type T function on the domain $D$ if there are positive constants M and T such that $|f(\mathit{\omega })|\le Mexp\{T|\mathit{\omega }|\}$, for $\mathit{\omega }\in D$.

The well-known Paley–Wiener theorem states that the Fourier transform of an $L_2(R)$ function vanishes outside of an interval $[-T,T]$, if and only if the function is of exponential-type T (see Dym & McKean, 1972, p. 158, for more details). The exponential-type functions are continuous functions which are infinitely differentiable everywhere and have a Taylor series expansion over every interval (see Champeney, 1987, p. 77; Walnut, 2002, p. 81). These functions are also called band-limited functions (see Bracewell, 2000, p. 119, for more details on bandlimited functions) (which are equivalent to exponential-type functions by the above stated Paley–Wiener theorem). The index of a complex-valued function $f$ on a smooth oriented curve $\mathrm{\Gamma }$, such that $f(\mathrm{\Gamma })$ is closed and compact, is defined to be the winding number of $f(\mathrm{\Gamma })$ about the origin (see Payandeh Najafabadi, 2007, §1, for more technical details). Computing the index of a function is usually a key step to determine the existence and number of solutions of a Wiener–Hopf factorization problem. The Sokhotski–Plemelj integral of a function s which satisfies the Hölder condition and it is defined by a principal value integral, as follows.(1) $$\begin{array}{c}\hfill {\mathit{\varphi }}_s(\mathit{\lambda }):=\frac{1}{2\mathit{\pi }i}{ʄ}_R\frac{s(x)}{x-\mathit{\lambda }}dx,\text{for}\mathit{\lambda }\in C.\end{array}$$(1)

The following are some well-known properties of the Sokhotski-Plemelj integral, proofs can be found in Ablowitz and Fokas (1990, §7), Gakhov (1990, §2), and Pandey (1996, §4), among others. The radial limit of the Sokhotski-Plemelj integral of s, given by ${\mathit{\varphi }}_s^{\pm }(\mathit{\omega })={lim}_{\mathit{\lambda }\to \mathit{\omega }+i{0}^{\pm }}{\mathit{\varphi }}_s(\mathit{\lambda })$ can be represented as the jump formula. i.e. ${\mathit{\varphi }}_s^{\pm }(\mathit{\omega })=\pm s(\mathit{\omega })/2+{\mathit{\varphi }}_s(\mathit{\omega })$, (or ${\mathit{\varphi }}_s^{\pm }(\mathit{\omega })=\pm s(\mathit{\omega })/2+H_s(\mathit{\omega })/(2i)$) where $H_s(\mathit{\omega })$ is the Hilbert transform of s and $\mathit{\omega }\in R$.

The Hausdorff-Young theorem states that: If s is a function in $L_p(R)$. Then, its Fourier transform, say $\widehat{s}$, is an $L_{p^{\ast }}(R)$ function that satisfies $||\widehat{s}{||}_{p^{\ast }}\le {(2\mathit{\pi }/p)}^{-1/p}/{p^{\ast }}^{1/p^{\ast }}{||s||}_p$, where $1<p\le 2$ and $1/p+1/p^{\ast }=1$, see Pandey (1996) for more details. From the Hausdorff–Young Theorem, one can observe that if $\{s_n\}$ is a sequence of functions converging in $L_p(R)$, $1<p\le 2$, to s. Then, the Fourier transforms of $s_n$ converge in $L_{p^{\ast }}(R)$, to the Fourier transform of s, whenever $1/p+1/p^{\ast }=1$. Using the Hausdorff–Young theorem, Payandeh Najafabadi and Kucerovsky (2014a) established that Hilbert transform of an $L_p(R)$, $1<p\le 2$ function s, say $H_s$, satisfies(2) $$\begin{array}{c}\hfill ||H_s{||}_p\le {||s||}_p.\end{array}$$(2)

Form this observation, one may conclude that “if $\{f_n\}$, $n\ge 1$, is a sequence of functions which converge in $L_p(R)$, $1<p\le 2$, to $f$. Then, the Hilbert transforms of $f_n$’s also converge in $L_p(R)$ to the Hilbert transform of $f$”.

The following, from Kucerovsky and Payandeh Najafabadi (2009), recalls some further useful properties of functions in $L_p(R)$ space.

Lemma 1

Suppose s and r are functions in $L_p(R)$, and suppose that $|s|$ and $|r|$ are bounded above by a. Then,

(i)	$||\sqrt{s}-\sqrt{r}{||}_p\le \frac{1}{2\sqrt{a}}||s-{r||}_p;$
(ii)	$||lns-{lnr||}_p\le a^{-1}||s-{r||}_p$, whenever s and r are positive-valued functions;
(iii)	$||e^{-is/2}-e^{-ir/2}{||}_p\le \frac{1}{2}||s-{r||}_p$, whenever s and r are real-valued functions;
(iv)	$||1/s-1/{r||}_p\le a^{-2}||s-{r||}_p$.

The followings recall definition and some useful properties on a mixture-gamma distribution, which plays an important role for the next sections (see Bracewell, 2000 for more details).

Definition 2

(mixture-gamma family of distributions) A non-negative random variable $X$ is said to be distributed according to a mixture-gamma distribution if its density function is given by(3) $$\begin{array}{c}\hfill p(x)=\sum _{k=1}^{\mathit{\nu }}\sum _{j=1}^{n_{\mathit{\nu }}}{c}_{kj}\frac{{\mathit{\alpha }}_k^j{x}^{j-1}}{(j-1)!}{e}^{-{\mathit{\alpha }}_{k}x},x\ge 0,\end{array}$$(3)

where $c_{k_j}$ and ${\mathit{\alpha }}_k$ are positive value which satisfy ${\sum }_{k=1}^{\mathit{\nu }}{\sum }_{j=1}^{n_{\mathit{\nu }}}{c}_{k_j}=1$.

Lemma 2

The characteristic function of a distribution (or equivalently the Fourier transform of its density function), say $\widehat{p}$, has the following properties:

(i)	$\widehat{p}$ is a rational function if and only if the density function belongs to the mixture-gamma family given by ;
(ii)	$\widehat{p}(0)=1;$ and the norm of $\widehat{p}(\mathit{\omega })$ bounded by 1.

3. Main results

Definition 3

The Wiener–Hopf factorization is the problem of finding a sectionally analytic function $\mathrm{\Phi }$ whose upper and lower radial limits at the real line, say ${\mathrm{\Phi }}^{\pm }$, satisfy(4) $$\begin{array}{c}\hfill {\mathrm{\Phi }}^{+}(\mathit{\omega }){\mathrm{\Phi }}^{-}(\mathit{\omega })=g(\mathit{\omega }),\text{for}w\in R,\end{array}$$(4)

where $g$ is a given continuous function satisfying a Hölder condition on $R$. Moreover, $g$ is assumed to have zero index, to be non-vanishing on $R$, and bounded above by 1.

Payandeh Najafabadi and Kucerovsky (2011) established that sectionally analytic functions ${\mathrm{\Phi }}^{\pm }$ satisfying a zero index Wiener–Hopf factorization (3) can be found by$$\begin{array}{c}\hfill {\mathrm{\Phi }}^{\pm }(\mathit{\lambda })=exp\{\pm {\mathit{\varphi }}_{lng}(\mathit{\lambda })\mp {\mathit{\varphi }}_{lng}(0)\},\mathit{\lambda }\in C,\end{array}$$

where ${\mathit{\varphi }}_{lng}$ stands for the Sokhotski–Plemelj integration of $lng$. Using the jump formula, the above ${\mathrm{\Phi }}^{\pm }$ can be re-stated as(5) $$\begin{array}{c}\hfill {\mathrm{\Phi }}^{\pm }(\mathit{\omega })=\sqrt{g(\mathit{\omega })}exp\{\pm \frac{i}{2}(H_{lng}(0)-H_{lng}(\mathit{\omega }))\},\end{array}$$(5)

where $H_{lng}$ stands for the Hilbert transform of $lng$.

In many situations, the Wiener–Hopf factorization problems (3) cannot be solved explicitly and has to be solved approximately (see Kucerovsky & Payandeh Najafabadi, 2009, for more details). The following develops an approximate technique to solve Equation 3.

Theorem 1

Suppose $g$ in the the Wiener-Hopf factorization problem (3) is a given, bounded (above by a), zero index function, satisfies the Hölder condition and $g(0)=1$. Moreover, suppose that there is a sequence of sectionally analytic functions ${\mathrm{\Phi }}_n^{\pm }(\mathit{\omega })=\sqrt{g(\mathit{\omega })}exp\{\pm \frac{i}{2}(H_{lng}(0)-H_{lng}(\mathit{\omega }))$, where $g_n$-s converge (in $L_p(R)$, $1<p\le 2$, sense) to $g$. Then, sectionally analytical solution of the Wiener–Hopf factorization problem (3), say ${\mathrm{\Phi }}^{\pm }$, can be approximated by ${\mathrm{\Phi }}_n^{\pm }$ and the error estimate satisfies$$\begin{array}{c}\hfill ||{\mathrm{\Phi }}_n^{\pm }-{\mathrm{\Phi }}^{\pm }{||}_p\le \frac{1}{2a\sqrt{a}}||g_n-{g||}_p^2+\frac{3}{2\sqrt{a}}||g_n-g{||}_p.\end{array}$$

Proof

Set $k(\mathit{\omega }):=-H_{lng}(\mathit{\omega })+H_{lng}(0)$ and $k_n(\mathit{\omega }):=-H_{ln{g}_n}(\mathit{\omega })+H_{ln{g}_n}(0)$. Now, from Equation 5, Equation 2, and Lemma 1, observe that$$\begin{array}{cc}\hfill ||{\mathrm{\Phi }}_n^{\pm }-{\mathrm{\Phi }}^{\pm }{||}_p& =||\sqrt{g_n}{e}^{\pm i{k}_n/2}-\sqrt{g}{e}^{\pm ik/2}{||}_p\hfill \\ \hfill & \le \left[||\sqrt{g_n}-\sqrt{g}{||}_p+||\sqrt{g}{||}_p\right]||e^{\pm i{k}_n/2}-e^{ik/2}{||}_p+|e^{\pm ik/2}|||\sqrt{g_n}-\sqrt{g}{||}_p\hfill \\ \hfill & \le \frac{1}{2}\left[||\sqrt{g_n}-\sqrt{g}{||}_p+||\sqrt{g}{||}_p\right]||-H_{ln{g}_n}(\mathit{\omega })+H_{ln{g}_n}(0)+H_{lng}(\mathit{\omega })-H_{lng}(0){||}_p\hfill \\ \hfill & +|e^{\pm ik/2}|||\sqrt{g_n}-\sqrt{g}{||}_p\hfill \\ \hfill & \le \left[||\sqrt{g_n}-\sqrt{g}{||}_p+||\sqrt{g}{||}_p\right]||H_{ln{g}_n}-H_{lng}{||}_p+||\sqrt{g_n}-\sqrt{g}{||}_p\hfill \\ \hfill & \text{since}k\text{and}{k}_n\text{are real-valued functions}\hfill \\ \hfill & \le \left[||\sqrt{g_n}-\sqrt{g}{||}_p+||\sqrt{g}{||}_p\right]||ln(g_n)-ln(g){||}_p+||\sqrt{g_n}-\sqrt{g}{||}_p\hfill \\ \hfill & \le \left[\frac{1}{2\sqrt{a}}||g_n-g{||}_p+\sqrt{a}\right]\frac{1}{a}||g_n-{g||}_p+\frac{1}{2\sqrt{a}}||g_n-g{||}_p\hfill \\ \hfill & =\frac{1}{2a\sqrt{a}}||g_n-{g||}_p^2+\frac{3}{2\sqrt{a}}||g_n-g{||}_p.\hfill \end{array}$$$\square$

In our belief, the most favorable situation is to approximate $g$ by a sequence of rational functions, which are obtained from a Padé approximant or a continued fraction expansion.

Using the Shannon sampling theorem, the following develops an elegant scheme that allows to provide an explicit solution for a Wiener–Hopf factorization Equation 3.

Theorem 2

Suppose $g$ in the Wiener–Hopf factorization problem (Equation 3) is a given, $L_1(R)\cap L_2(R)$, bounded, zero index function, satisfies the Hölder condition. Moreover, suppose that there is a sequence of $ln(g)$ is an exponential-type T function, then unique solutions of the Wiener–Hopf factorization (Equation 3) can be explicitly determined by$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\pm }(\mathit{\omega })=exp\left\{\pm \sum _{n=-\infty }^{\infty }ln(g(\frac{2n}{T}))\frac{e^{\pm i\mathit{\pi }(T\mathit{\omega }-2n)}-1}{2i\mathit{\pi }(T\mathit{\omega }-2n)}\right\},\end{array}$$

Proof

Using the fact that $ln(g)$ is an exponential-type T function, one can decompose $ln(g(\mathit{\omega }))$ as $ln(g(\mathit{\omega }))=K_{+}(\mathit{\omega })+K_{-}(\mathit{\omega })$, where$$\begin{array}{cc}\hfill K_{\pm }(\mathit{\omega })& =\pm \sum _{n=-\infty }^{\infty }ln(g(\frac{2n}{T}))\frac{exp\{\pm i\mathit{\pi }(T\mathit{\omega }-2n)\}-1}{2i\mathit{\pi }(T\mathit{\omega }-2n)}.\hfill \end{array}$$

Sectionally analytical properties of $K_{\pm }(\mathit{\omega })$ in $C_{\pm }$ has been established by Kucerovsky and Payandeh Najafabadi (2009). $\square$

The following theorem provides the error bound for approximate solution arrives from the Shannon sampling theorem.

Theorem 3

Suppose $g$ in the Wiener–Hopf factorization problem () is a given, $L_1(R)\cap L_2(R)$, bounded, zero index function, satisfies the Hölder condition and $ln(g)$ is an exponential-type T function. Moreover, suppose that there is a sequence of $g_m$ in $L_1(R)\cap L_2(R)$ and $ln(g_m)$ are exponential-type T functions. Then, approximate solutions of the Wiener–Hopf factorization problem () can be determined by$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\pm }^{(m)}(\mathit{\omega })=\pm exp\left\{\pm \sum _{n=-\infty }^{\infty }ln(g_m)(\frac{2n}{T}))\frac{e^{\pm i\mathit{\pi }(T\mathit{\omega }-2n)}-1}{2i\mathit{\pi }(T\mathit{\omega }-2n)}\right\}\end{array}$$

and the error bound satisfies$$\begin{array}{c}\hfill |{\mathrm{\Phi }}_{\pm }^{(m)}-{\mathrm{\Phi }}_{\pm }|\le \Vert ln(g_m)-ln(g)\Vert ,\end{array}$$

where the norm is defined by $\Vert M\Vert :={sup}_{ij}{\left\{{\int }_{-\infty }^{\infty }{|M_{ij}(x)|}^{2}dx\right\}}^{1/2}$.

Proof

The proof is straightforward by a double application of Lemma (3.2, Payandeh Najafabadi & Kucerovsky, 2014a) and Theorem 2. $\square$

4. Applications

This section provides the application of the above results in two different contexts. The first subsection considers the problem of finding spectral factorization, whenever its corresponding spectral density function has been given. The second subsection derives density/probability functions of extrema random variables of a given Lévy processes.

4.1. Application to spectral factorization

Definition 4

Suppose $\mathit{\rho }(x)$ is an autocorrelation of the stochastically stationary process $X(t)$, i.e. $\mathit{\rho }(x):=Corr(X)(t),X(t-x))$. Then, the spectral function $S(\mathit{\lambda })$ is defined by $\mathit{\rho }(\mathit{\tau })={\int }_{-\infty }^{\infty }{e}^{i\mathit{\lambda }\mathit{\tau }}S(\mathit{\lambda })d\mathit{\mu }(\mathit{\lambda })$, where $\mathit{\mu }$ is a given measure.

The spectral density function $S$ has properties that: (1) for stochastically stationary processes, it can be understood as the Fourier transform of the autocorrelation $\mathit{\rho }(·)$ [The Wiener–Khintchine’s theorem: Reinsel (1997, p. 219)]; (2) it is a Hermitian function (Koopmans, 1995, p. 122); (3) it defines almost everywhere (a.e.) on the interval $[-\mathit{\pi },-\mathit{\pi }]$ (Wilson, 1972); and (4) $S$ is integrable and has a Fourier series expansion $S(\mathit{\theta })={\sum }_{k=-\infty }^{\infty }{\mathit{\gamma }}_k{e}^{ik\mathit{\theta }}$, where ${\mathit{\gamma }}_k={\int }_{-\mathit{\pi }}^{\mathit{\pi }}\mathbf{S}(\mathit{\theta })e^{-ik\mathit{\theta }}d\mathit{\theta }$ (Wilson, 1972).

The spectral factorization plays a crucial role in a wide range of scientific fields, such as communications (Magesacher & Cioffi, 2011), system theory (Janashia, Lagvilava, & Ephremidze, 2011), optimal control (Johannesson, Rantzer, & Bernhardsson, 2011), filtering theory (Anderson & Moore, 2005), network theory (Belevitch, 1968; Ivrlac & Nossek, 2014) prediction theory of stationary processes (Brockwell & Davis, 2002), deriving forward expression from a backward one in the ARIMA processes (Brockwell & Davis, 2009). In such applications, the spectral density factorization is the most difficult step. Since Wiener’s seminal efforts create a computational method for spectral factorization, several authors have developed different methods to do so, but none of them have provided a method that has an essential superiority over all others (Janashia et al., 2011). On the other hand, most of these methods impose some extra restriction on the spectral density (e.g. to be real or rational function, or to be non-singular on the boundary). The following represents an explicit solution for the problem of spectral factorization for spectral densities.

Proposition 1

Suppose $S$ is a spectral density function. Then, the right and left spectral factorizations $L_{\pm }$ are given by$$\begin{array}{c}\hfill L_{\pm }(\mathit{\omega })=exp\left\{\pm \sum _{n=-\infty }^{\infty }ln(S(\frac{2n}{T}))\frac{e^{\pm i\mathit{\pi }(T\mathit{\omega }-2n)}-1}{2i\mathit{\pi }(T\mathit{\omega }-2n)}\right\}.\end{array}$$

Proof

The exponential-type T condition of $ln(S)$ arrives from the fact that the spectral density can be understood as an extension version of the inverse Fourier transform on the autocorrelation function which is bounded function. The rest of proof arrives from an application of Theorem 2 along the fact that $S$ is Hermitian function which its index is zero (Voronin, 2010,2011). $\square$

4.2. Application to Lévy processes

Suppose $X_t$ is a one-dimensional real-valued Lévy process started from $X_0=0$ and defined by a triple $(\mathit{\mu },\mathit{\sigma },\mathit{\nu }):$ the drift $\mathit{\mu }\in R$, volatility $\mathit{\sigma }\ge 0$, and the jumps measure $\mathit{\nu }$ which is given by a non-negative function defined on $R\backslash \{0\}$ satisfying ${\int }_{R}min\{1,x^2\}\mathit{\nu }(dx)<\infty$. Moreover, suppose that random stopping time $\mathit{\tau }(q)$ has either a geometric (with parameter $q\in (0,1)$) or an exponential distribution (with parameter $q>0$) and independent of the Lévy process $X_t$ which $\mathit{\tau }(0)=\infty$. The Lévy–Khintchine formula states that the characteristic exponent $\mathit{\psi }$ (i.e. $\mathit{\psi }(\mathit{\omega })=ln(E(exp(i\mathit{\omega }{X}_1))),\mathit{\omega }\in R$) can be represented by(6) $$\begin{array}{c}\hfill \mathit{\psi }(\mathit{\omega })=i\mathit{\mu }\mathit{\omega }-\frac{1}{2}{\mathit{\sigma }}^2{\mathit{\omega }}^2+{\int }_R(e^{i\mathit{\omega }x}-1-i\mathit{\omega }x{I}_{[-1,1]}(x))\mathit{\nu }(dx),\mathit{\omega }\in R.\end{array}$$(6)

The extrema of the Lévy process $X_t$ are given by(7) $$\begin{array}{cc}& M_q=sup\{X_s:s\le \mathit{\tau }(q)\}\hfill \\ \hfill & I_q=inf\{X_s:s\le \mathit{\tau }(q)\}.\hfill \end{array}$$(7)

The Wiener–Hopf factorization is a well-known technique to study the characteristic functions of the extrema random variables (Bertoin, 1996). Namely, the Wiener–Hopf factorization states that: (i) random variables $M_q$ and $I_q$ are independent; (ii) product of their characteristic functions equal to the characteristic function of Lévy process $X_t;$ (iii) random variable $M_q$ ($I_q$) is infinitely divisible, positive (negative), and has zero drift.

In the cases that the characteristic function of Lévy process $X_t$ either a rational function or can be decomposed as a product of two sectionally analytic functions in the closed upper and lower half complex planes $C^{+}$ and $C^{-}$. The characteristic functions of random variables $M_q$ and $I_q$ can be determined explicitly. Lewis and Mordecki (2005) considered a Lévy process $X_t$ which its negative jumps distributed according to a mixture-gamma family of distributions (given by Definition 2) and an arbitrary positive jumps measure. They established such process has the characteristic function which can decompose as a product of a rational function and an arbitrary function, which are analytic in $C^{+}$ and $C^{-}$, respectively. Moreover, they provided an analog result for a Lévy process whose its corresponding positive jumps measure follows from a mixture-gamma family of distributions while its negative jumps measure is an arbitrary one, more details can be found in Lewis and Mordecki (2008).

Unfortunately, in the most situations, the characteristic function of the process neither is a rational function nor can be decomposed as a product of two analytic functions in $C^{+}$ and $C^{-}$. Therefore, the characteristic functions of $M_q$ and $I_q$ should be expressed in terms of a Sokhotski-Plemelj integral (see Equation 1). But, this form, also, presents some difficulties in numerical work due to slow evaluation and numerical problems caused by singularities near the integral contour. To overcome these difficulties, approximation methods have to be considered.

It is well known that a Lévy process $X_t$ which its jumps distribution follows from the phase-type distribution has a rational characteristic function (Doney, 1987). Kuznetsov (2010) utilized this fact and approximated a jumps measure $\mathit{\nu }$ a ten-parameter family of Lévy processes (named $\mathit{\beta }-$family of Lévy process) by a sequence of the phase-type measures. Then, he determined the characteristic functions of random variables $M_q$ and $I_q$, approximately.

The following theorem represents an estimation bound for approximated extrema’s density/probability functions of a Lévy process.

Theorem 4

Suppose $X_t$ is a Lévy process defined by a triple $(\mathit{\mu },\mathit{\sigma },\mathit{\nu })$ and its random stopping time $\mathit{\tau }(q)$ has been distributed according to either a geometric or an exponential distribution. Moreover suppose that there is a sequence of jumps measure ${\mathit{\nu }}_n$ which satisfies the following two conditions:

A1)

They converge in $L_p(R),1<p\le 2$, to $\mathit{\nu }$ and ${\int }_{-1}^{1}x{\mathit{\nu }}_n(dx)={\int }_{-1}^{1}x\mathit{\nu }(dx);$

A2)

their corresponding characteristic exponents ${\mathit{\psi }}_n$ (arrived by the Lévy-Khintchine Formula 6(6) $$\begin{array}{c}\hfill \mathit{\psi }(\mathit{\omega })=i\mathit{\mu }\mathit{\omega }-\frac{1}{2}{\mathit{\sigma }}^2{\mathit{\omega }}^2+{\int }_R(e^{i\mathit{\omega }x}-1-i\mathit{\omega }x{I}_{[-1,1]}(x))\mathit{\nu }(dx),\mathit{\omega }\in R.\end{array}$$(6) ) as well as the characteristic exponents $\mathit{\psi }$ (correspondence with jumps measure $\mathit{\nu }$) are bounded above by M.

Then, density/probability of supremum and infimum of Lévy process $X_t$, say respectively $f_q^{+}$ and $f_q^{-}$, can be approximated by sequence of density/probability functions $f_{q,n}^{+}$ and $f_{q,n}^{-}$ which has the following error bound.

(i)	For exponentially distributed stopping time $\mathit{\tau }(q)$,$$||f_q^{\pm }-f_{q,n}^{\pm }{||}_p\le \frac{q^2}{2M^4{(2\mathit{\pi })}^{1/p}}||{\mathit{\nu }}_n-{\mathit{\nu }||}_p^2+\frac{3q}{2M^2}||{\mathit{\nu }}_n-\mathit{\nu }{||}_p;$$
(ii)	For geometric stopping time $\mathit{\tau }(q)$,$$||f_q^{\pm }-f_{q,n}^{\pm }{||}_p\le \frac{{(1-q)}^2}{2M^4{(2\mathit{\pi })}^{1/p}}||{\mathit{\nu }}_n-{\mathit{\nu }||}_p^2+\frac{3(1-q)}{2M^2}||{\mathit{\nu }}_n-\mathit{\nu }{||}_p.$$

Proof

From Bertoin (1996), one can observe that the Fourier transform of $M_q$ and $I_q$ density functions, say respectively ${\mathrm{\Phi }}^{+}$ and ${\mathrm{\Phi }}^{-}$, satisfy either the Wiener–Hopf factorization ${\mathrm{\Phi }}^{+}(\mathit{\omega }){\mathrm{\Phi }}^{-}(\mathit{\omega })=q/(q-\mathit{\psi }(\mathit{\omega })$, where $\mathit{\omega }\in R$ (for exponentially distributed stopping time) or the Wiener–Hopf factorization ${\mathrm{\Phi }}^{+}(\mathit{\omega }){\mathrm{\Phi }}^{-}(\mathit{\omega })=(1-q)/(1-q\mathit{\psi }(\mathit{\omega })$, where $\mathit{\omega }\in R$ (for geometric stopping time). Now, from the fact that expressions $q{(q-\mathit{\psi }(·))}^{-1}$ and $(1-q){(1-q\mathit{\psi }(·))}^{-1}$ are the characteristic function of Lévy process $X_t$ respectively for exponential and geometric stopping time, observe that both expressions are bounded by 1 (property of the characteristic function given by Lemma 2, part ii). For part (i), from Theorem 1 observe that$$\begin{array}{cc}\hfill ||{\mathrm{\Phi }}_n^{\pm }-{\mathrm{\Phi }}^{\pm }{||}_{p^{\ast }}& \le \frac{1}{2}||\frac{q}{q-{\mathit{\psi }}_n}-\frac{q}{q-\mathit{\psi }}{||}_{p^{\ast }}^2+\frac{3}{2}||\frac{q}{q-{\mathit{\psi }}_n}-\frac{q}{q-\mathit{\psi }}{||}_{p^{\ast }}\hfill \\ \hfill & \le \frac{q^2}{2M^4}||{\mathit{\psi }}_n-{\mathit{\psi }||}_{p^{\ast }}^2+\frac{3q}{2M^2}||{\mathit{\psi }}_n-\mathit{\psi }{||}_{p^{\ast }}\hfill \\ \hfill & \le \frac{q^2}{2M^4{(2\mathit{\pi })}^{2/p}}||{\mathit{\nu }}_n-{\mathit{\nu }||}_p^2+\frac{3q}{2M^2{(2\mathit{\pi })}^{1/p}}||{\mathit{\nu }}_n-\mathit{\nu }{||}_p,\hfill \end{array}$$

where $1/p^{\ast }+1/p=1$. The second inequality arrives from part (iv) of Lemma 1, while the third inequality obtains from Equation 6 along with conditions $A_2$ and an application of Hausdorff-Young theorem. The rest of proof arrives from an application of the Hausdorff-?Young theorem. Proof of part (ii) is quite similar. $\square$

It would worthwhile to mention that, in the case of ${\int }_{-1}^{1}x{\mathit{\nu }}_n(dx)=:c_n\ne d:={\int }_{-1}^{1}x\mathit{\nu }(dx)$. One may obtain sequence ${\mathit{\xi }}_n=\frac{d}{c_n}{\mathit{\nu }}_n$ which satisfy the desire condition.

The following utilizes result of the above theorem and provides a procedure to find the extrema’s density/probability of a wide class of Lévy processes, approximately.

Procedure 1

Suppose $X_t$ is a Lévy process with bounded characteristic exponents $\mathit{\psi }$. Moreover, suppose that random stopping time $\mathit{\tau }(q)$ has been distributed according to either a geometric or an exponential distribution. Then, by the following steps, one can approximate (in $L_p(R),1<p\le 2$, sense) density/probability functions of the extrema random variables $M_q$ and $I_q$.

Step 1.	Approximate jumps measure $\mathit{\nu }$ with either a phase-type or a mixture-gamma density function, say ${\mathit{\nu }}^{\ast }$, where ${\int }_{-1}^{1}x{\mathit{\nu }}^{\ast }(dx)={\int }_{-1}^{1}x\mathit{\nu }(dx)$ and $||\mathit{\nu }-{\mathit{\nu }}^{\ast }{||}_p\le \mathit{\epsilon };$
Step 2.	Decompose rational function $q/(q-\mathit{\psi }(\mathit{\omega }))$ (or $(1-q)/(1-q\mathit{\psi }(\mathit{\omega }))$) into product of two rational and sectionally analytic functions, say $g^{\pm }$, in $C^{\pm }$, respectively;
Step 3.	Obtain, approximate, density/probability functions of $M_q$ and $I_q$ by the inverse Fourier transform of $g^{+}$ and $g^{-}$, respectively.

In many situations, it is more convenient to approximate the characteristic exponents $\mathit{\psi }$, rather than the jumps measure $\mathit{\nu }$. The following extends results of Theorem 4 to such situations.

Corollary 1

Suppose $X_t$ is a Lévy process defined by a triple $(\mathit{\mu },\mathit{\sigma },\mathit{\nu })$ and its random stopping time $\mathit{\tau }(q)$ has been distributed according to either a geometric or an exponential distribution. Moreover suppose that there is a sequence of bounded characteristic exponents ${\mathit{\psi }}_n$ (i.e. $|{\mathit{\psi }}_n|\le M$) which converge in $L_p(R),1<p\le 2$, to the characteristic exponent of the process $\mathit{\psi }$, and their corresponding jumps measure ${\mathit{\nu }}_n$ satisfies ${\int }_{-1}^{1}x{\mathit{\nu }}_n(dx)={\int }_{-1}^{1}x\mathit{\nu }(dx)$.

Then, density/probability function of the supremum and the infimum of Lévy process $X_t$, say respectively $f_q^{+}$ and $f_q^{-}$, can be approximated by sequence of density/probability functions $f_{q,n}^{+}$ and $f_{q,n}^{-}$ which has the following error bound.

(i)	For exponentially distributed stopping time $\mathit{\tau }(q)$,$$||f_q^{\pm }-f_{q,n}^{\pm }{||}_p\le \frac{q^2{(2\mathit{\pi })}^{1/p}}{2M^4}||{\mathit{\psi }}_n-{\mathit{\psi }||}_{p^{\ast }}^2+\frac{3q{(2\mathit{\pi })}^{1/p}}{2M^2}||{\mathit{\psi }}_n-\mathit{\psi }{||}_{p^{\ast }};$$
(ii)	For discrete stopping time $\mathit{\tau }(q)$,$$||f_q^{\pm }-f_{q,n}^{\pm }{||}_p\le \frac{{(1-q)}^2{(2\mathit{\pi })}^{1/p}}{2M^4}||{\mathit{\psi }}_n-{\mathit{\psi }||}_{p^{\ast }}^2+\frac{3(1-q){(2\mathit{\pi })}^{1/p}}{2M^2}||{\mathit{\psi }}_n-\mathit{\psi }{||}_{p^{\ast }}.$$

Using the fact that, the characteristic exponent $\mathit{\psi }(i\mathit{\omega }),\mathit{\omega }\in R$, is a real-valued function, Bertoin (1996). One can suggest that following procedure to generate approximation density/probability functions for $M_q$ and $I_q$.

Procedure 2

Suppose $X_t$ is a Lévy process with bounded characteristic exponents $\mathit{\psi }$. Moreover, suppose that random stopping time $\mathit{\tau }(q)$ has been distributed according to either a geometric or an exponential distribution. Then, by the following steps, one can approximate (in $L_p(R),1<p\le 2$, sense) density/probability functions of the extrema random variables $M_q$ and $I_q$.

Step 1.	Approximate (in $L_p(R),1<p\le 2$, sense) the characteristic exponent $\mathit{\psi }(i\mathit{\omega })$ by a rational function, generated by the Padé approximant or the continued fraction, say ${\mathit{\psi }}^{\ast }(i\mathit{\xi });$
Step 2.	Decompose rational function $q/(q-{\mathit{\psi }}^{\ast }(\mathit{\omega }))$ (or $(1-q)/(1-q{\mathit{\psi }}^{\ast }(\mathit{\omega }))$) into product of two rational and sectionally analytic functions, say $g^{\pm }$, in $C^{\pm }$, respectively;
Step 3.	Obtain, approximate, density/probability functions of $M_q$ and $I_q$ by the inverse Fourier transform of $g^{+}$ and $g^{-}$, respectively.

Using Lemma 2, one can readily, conclude that the above two procedures approximate density/probability functions of $M_q$ and $I_q$ by the mixture-gamma density functions.

Now, we provide several examples.

Example 1

Consider a 1-stable Lévy process with a jumps measure $\mathit{\nu }(dx)={c|x|}^{-2}{I}_{(-\infty ,\infty )}(x)dx$. One can, readily, show that the characteristic exponent of such process is $\mathit{\psi }(\mathit{\omega })=(c_1+c_2)|\mathit{\omega }|\{1+i\mathit{\beta }\text{sgn}\mathit{\omega }\frac{2}{\mathit{\pi }}ln(|\mathit{\omega }|)\}+i\mathit{\omega }\mathit{\eta }$, where $\mathit{\eta }$ is a normalized real-valued. The natural logarithm $ln(|\mathit{\lambda }|)$ has the continued fraction$$ln(|\mathit{\lambda }|)={\frac{\mathit{\lambda }-1}{1}}_{+}{\frac{1^2(\mathit{\lambda }-1)}{2}}_{+}{\frac{1^2(\mathit{\lambda }-1)}{3}}_{+}{\frac{2^2(\mathit{\lambda }-1)}{4}}_{+}{\frac{2^2(\mathit{\lambda }-1)}{5}}_{+}{\frac{3^2(\mathit{\lambda }-1)}{7}}_{+}{\frac{3^2(\mathit{\lambda }-1)}{9}}_{+}\dots ,\mathit{\lambda }\in C.$$(see Jones & Thron, 1931 or Cuyt, Petersen, Verdonk, Waadeland, & Jones, 2008, among others). In Practice, in the above-continued fraction, it has to cut off somewhere and obtain a rational function for $ln(·)$. Consequently, an expression $q{(q-\mathit{\psi }(·))}^{-1}$ can be approximated by rational function $P(\mathit{\lambda })/Q(\mathit{\lambda })$. Now, the characteristic functions for the extrema can be obtained after decomposing $P(\mathit{\lambda })/Q(\mathit{\lambda })$ into product of two rational, analytic, and bounded functions, in $C^{+}$ and $C^{-}$.

Example 2

Suppose $X_t$ is a Lévy process with independent and continuous $\mathit{\tau }(q)$ and a jumps measure $\mathit{\nu }(dx)=exp\{\mathit{\alpha }x\}cosec{h}^2(x/2)dx$. The characteristic exponent for such Lévy process is given by$$\begin{array}{c}\hfill \mathit{\psi }(\mathit{\omega })=\frac{{\mathit{\sigma }}^2{\mathit{\omega }}^2}{2}+i\mathit{\rho }\mathit{\omega }+4\mathit{\pi }(\mathit{\omega }-i\mathit{\alpha })coth(\mathit{\pi }(\mathit{\omega }-i\mathit{\alpha }))-4\mathit{\gamma },\end{array}$$

where $\mathit{\gamma }=\mathit{\pi }\mathit{\alpha }cot(\mathit{\pi }\mathit{\alpha })$, $\mathit{\rho }=4{\mathit{\pi }}^2\mathit{\alpha }+\frac{4\mathit{\gamma }(\mathit{\gamma }-1)}{\mathit{\alpha }}-\mathit{\mu }$, $\mathit{\omega }\in R$, and $\mathit{\alpha },\mathit{\mu },\text{and}\mathit{\sigma }$ are given. The continued fraction for $tanh(\mathit{\lambda })$ is given by$$\begin{array}{c}\hfill tanh(\mathit{\lambda })=\frac{\mathit{\lambda }}{1}+{\frac{{\mathit{\lambda }}^2}{3}}_{+}{\frac{{\mathit{\lambda }}^2}{5}}_{+}{\frac{{\mathit{\lambda }}^2}{7}}_{+}{\frac{{\mathit{\lambda }}^2}{9}}_{+}\cdots ,\mathit{\lambda }\in C,\end{array}$$(Cuyt et al., 2008). After cutting off the above-continued fraction at the $n^{\text{th}}$ term, one may be obtained a rational function, say ${\mathit{\rho }}_n$, for $tanh(·)$. Substituting the ${\mathit{\rho }}_n$ in an expression $q(q-\mathit{\psi }(·))$, we may obtain a rational function, say $P(\mathit{\lambda })/Q(\mathit{\lambda })$. Therefore, the characteristic functions for the extrema can be found by decomposing a rational function $P(\mathit{\lambda })/Q(\mathit{\lambda })$ as a product of two rational, analytic, and bounded functions in $C^{+}$ and $C^{-}$. This observation verifies Kuznetsov’s (2010) result.

Example 3

Metron model is a Lévy process with a jumps measure $\mathit{\nu }(x)=a{(\mathit{\delta }\sqrt{2\mathit{\pi }})}^{-1}exp\{-{(x-\mathit{\mu })}^2/(2{\mathit{\delta }}^2)\}$ and characteristic exponent $\mathit{\psi }(\mathit{\omega })=i\mathit{\mu }\mathit{\omega }-{\mathit{\sigma }}^2{\mathit{\omega }}^2/2+a\{e^{-{\mathit{\delta }}^2{\mathit{\omega }}^2/2+i\mathit{\mu }\mathit{\omega }}-1\}$, where a, $\mathit{\delta }$, $\mathit{\mu }$, and $\mathit{\sigma }$ are given and $\mathit{\omega }\in R$. Metron model has tail behaviors heavier than the Gaussian but all exponential moments are finite. The continued fraction for the exponential function $exp(·)$ is given by$$\begin{array}{c}\hfill e^{\mathit{\lambda }}={\frac{1}{1}}_{-}{\frac{\mathit{\lambda }}{1}}_{+}{\frac{\mathit{\lambda }}{2}}_{-}{\frac{\mathit{\lambda }}{3}}_{+}{\frac{\mathit{\lambda }}{2}}_{-}{\frac{\mathit{\lambda }}{5}}_{+}{\frac{\mathit{\lambda }}{2}}_{-}\cdots ,\mathit{\lambda }\in C,\end{array}$$(Cuyt et al., 2008). After cutting off the above-continued fraction somewhere, and substituting the arrived rational function in an expression $q(q-\mathit{\psi }(\mathit{\lambda }))$, we may obtain a rational function, say $P(\mathit{\lambda })/Q(\mathit{\lambda })$. The characteristic functions for the extrema can be found by decomposing $P(\mathit{\lambda })/Q(\mathit{\lambda })$ as a product of two rational, analytic, and bounded functions in $C^{+}$ and $C^{-}$.

5. Conclusion and suggestion

This paper provides two techniques to solve a Wiener–Hopf factorization problem. Application of our findings, in two different contexts, has been given. It would be worthwhile mentioning that: (1) In the situation where given function $ln(g)$ is not of exponential-type function. We suggest to approximate such function with an exponential-type function which pointwise converges to such function (see Kucerovsky and Payandeh Najafabadi (2009) for more details); (2) Other applications our findings can be developed to other situations where the Wiener–Hopf factorization is applicable, such as finding first/last passage time and the overshoot, the last time the extrema was archived, several kind of option pricing, etc.

Acknowledgements

Thanks to an anonymous reviewer for his/her constructive comments. The author Dan Z. Kucerovsky acknowledges Natural Sciences and Engineering Research Council (NSERC).

Additional information

Funding

This work was supported by Natural Sciences and Engineering Research Council (NSERC) of Canada.

Notes on contributors

Amir T. Payandeh Najafabadi

Amir T. Payandeh Najafabadi is an Associate Professor in Department of Mathematics sciences at Shahid Behashti University, Tehran, Evin (Email address: [email protected]). He was born on Sep 3, 1973. He received his PhD from University of New Brunswick, Canada in 2006. He has published 28 papers and was co-author of two books. His major research interests are: Statistical Decision Theory, Lévy processes, Risk theory, Riemann-Hilbert problem, & integral equations.