Physics and Derivatives: Effective-Potential Path-Integral Approximations of Arrow-Debreu Densities

TL;DR

This paper introduces an effective-potential path-integral method, inspired by physics, to accurately approximate transition probabilities and Arrow-Debreu densities for diffusions, offering a computationally efficient alternative for derivatives pricing.

Contribution

The paper develops a semi-analytical approximation method based on path-integral techniques for diffusions, demonstrating high accuracy and efficiency in derivatives pricing.

Findings

Accurately approximates transition probabilities for complex models.

Performs well even in high volatility regimes.

Offers a computationally efficient alternative to numerical schemes.

Abstract

We show how effective-potential path-integrals methods, stemming on a simple and nice idea originally due to Feynman and successfully employed in Physics for a variety of quantum thermodynamics applications, can be used to develop an accurate and easy-to-compute semi-analytical approximation of transition probabilities and Arrow-Debreu densities for arbitrary diffusions. We illustrate the accuracy of the method by presenting results for the Black-Karasinski and the GARCH linear models, for which the proposed approximation provides remarkably accurate results, even in regimes of high volatility, and for multi-year time horizons. The accuracy and the computational efficiency of the proposed approximation makes it a viable alternative to fully numerical schemes for a variety of derivatives pricing applications.

Figures

Tables

Table 1. Table 1: Black-Karasinski T 𝑇 T maturity zero-coupon bonds obtained with the GTFK approximation, the Exponent Expansion (EE) of Ref.  Stehlíková and Capriotti ( 2014 ) , the Karhunen-Loéve (KL) expansion of Ref.  Daniluk and Muchorski ( 2016 ) to first and second order, and by solving numerically the associated PDE. The parameters of the BK process are: mean-reversion speed a = 0.1 𝑎 0.1 a=0.1 , level b = ln ⁡ 0.04 𝑏 0.04 b=\ln 0.04 , volatility σ = 0.85 𝜎 0.85 \sigma=0.85 , and initial rate r 0 = 0.06 subscript 𝑟 0 0.06 r_{0}=0.06 .

$T$	$\mathrm{EE}$	KL(1)	KL(2)	GTFK	PDE
0.1	0.9939 (0.00%)	0.9939 (0.00%)	0.9939 (0.00%)	0.9939 (0.00%)	0.9939
0.5	0.9681 (0.00%)	0.9681 (0.00%)	0.9681 (0.00%)	0.9681 (0.00%)	0.9681
1.0	0.9331 (0.00%)	0.9331 (0.00%)	0.9331 (0.00%)	0.9331 (0.00%)	0.9331
2.0	0.8581 (0.01%)	0.8580 (0.02%)	0.8581 (0.01%)	0.8582 (0.00%)	0.8582
3.0	0.7845 (0.01%)	0.7842 (0.05%)	0.7844 (0.02%)	0.7847 (0.01%)	0.7846
5.0	0.6595 (0.04%)	0.6582 (0.24%)	0.6593 (0.08%)	0.6602 (0.06%)	0.6598
10.0	-	0.4545 (1.69%)	0.4601 (0.48%)	0.4628 (0.10%)	0.4623
20.0	-	0.2440 (9.06%)	0.2592 (3.38%)	0.2672 (0.41%)	0.2683

Table 2. Table 2: GARCH linear SDE T 𝑇 T maturity zero-coupon bonds obtained with the GTFK approximation, the Exponent Expansion (EE) of Ref.  Capriotti et al. ( 2019 ) , and by solving numerically the associated PDE. The parameters of the process are: mean-reversion level a = 0.1 𝑎 0.1 a=0.1 , level b = 0.04 𝑏 0.04 b=0.04 , volatility σ = 0.6 𝜎 0.6 \sigma=0.6 , and initial rate y 0 = 0.06 subscript 𝑦 0 0.06 y_{0}=0.06 .

$T$	$\mathrm{EE}$	GTFK	PDE
0.1	0.9940 (0.00%)	0.9940 (0.00%)	0.9940
0.5	0.9707 (0.00%)	0.9707 (0.00%)	0.9707
1.0	0.9429 (0.00%)	0.9429 (0.00%)	0.9429
2.0	0.8914 (0.03%)	0.8920 (0.03%)	0.8917
3.0	0.8459 (0.08%)	0.8472 (0.07%)	0.8466
5.0	0.7834 (1.40%)	0.7717 (0.12%)	0.7726
7.5	-	0.6923 (1.45%)	0.7025
10.0	-	0.6223 (3.92%)	0.6477

Table 3. Table 3: GARCH linear SDE T 𝑇 T maturity zero-coupon bonds obtained with the GTFK approximation, the Exponent Expansion (EE) of Ref.  Shaimerdenova ( 2015 ) , and by solving numerically the associated PDE. The parameters of the process are: mean-reversion level a = 0.1 𝑎 0.1 a=0.1 , level b = 0.02 𝑏 0.02 b=0.02 , volatility σ = 0.5 𝜎 0.5 \sigma=0.5 , and initial rate y 0 = 0.01 subscript 𝑦 0 0.01 y_{0}=0.01 .

$T$	$\mathrm{EE}$	GTFK	PDE
0.1	0.9990 (0.00%)	0.9990 (0.00%)	0.9990
0.25	0.9975 (0.00%)	0.9975 (0.00%)	0.9975
0.5	0.9949 (0.00%)	0.9949 (0.00%)	0.9949
1.0	0.9896 (0.00%)	0.9896 (0.00%)	0.9896
2.5	0.9728 (0.02%)	0.9723 (0.03%)	0.9726
5.0	0.9359 (0.62%)	0.9403 (0.15%)	0.9417
10.0	0.8315 (5.10%)	0.8709 (0.60%)	0.8762

Equations

$$d{Y}_t={\mu }_y(Y_t)dt+{\sigma }_y(Y_t)d{W}_t,$$

$$\psi^{Y}_{\lambda}(y_{T},y_{0},T)=\mathbb{E}\Big{[}\delta(Y_{T}-y_{T})e^{-\lambda\int_{0}^{T}du\,r_{u}}\Big{|}\,Y_{0}=y_{0}\Big{]},$$

$${\int }_{A}d{y}_T{\psi }_0^Y(y_T,y_0,T)\equiv \mathbb{P}\left[Y_T\in A|Y_0=y_0\right].$$

$$V(0)=\mathbb{E}\Big{[}e^{-\int_{0}^{T}du\,r_{u}}P(r_{T})\Big{|}\,Y_{0}=y_{0}\Big{]},$$

$$V(0)=\int d{y}_T{\psi }_1^Y(y_T,y_0,T)P(y_T),$$

$$Z_{\lambda }(r_0,T)=\int d{y}_T{\psi }_{\lambda }^Y(y_T,y_0,T),$$

$$\hat{\mathcal{H}}=\frac{{\hat{p}}^2}{2m}+V(\hat{x})$$

$$\hat{\rho }=e^{-\beta \hat{\mathcal{H}}}$$

$$\rho (x_T,x_0,T)$$

$${\int }_{x(0)=x_0}^{x(T)=x_T}\mathcal{D}[x(t)]e^{S[x(t)]},$$

$$S[x(t)]=-\frac{1}{\hslash }{\int }_0^{T}dt\left[\frac{m}{2}{\dot{x}}^2(t)+V(x(t))\right],$$

$$\displaystyle\left(\frac{m}{{2\pi\hbar\Delta t}}\right)^{N/2}\int\ldots\int\prod_{i=1}^{N-1}dx_{i}\exp{\big{[}S(x_{i},x_{i-1})\big{]}}~{},$$

$$S(x_i,x_{i-1})=$$

$$-\frac{\Delta t}{\hslash }\left[\frac{m}{2}\frac{(x_i-x_{i-1}{)}^2}{\Delta t}+V((x_{i-1}+x_i)/2)\right].$$

$$\bar{x}[x(t)]=\frac{1}{T}{\int }_0^T\text{-}\text{-}dtx(t),$$

$$\rho (x_T,x_0,T)=\int d\bar{x}{\rho }_{\bar{x}}(x_T,x_0,T),$$

$${\rho }_{\bar{x}}(x_T,x_0,T)=$$

$${\int }_{x(0)=x_0}^{x(T)=x_T}\mathcal{D}[x(t)]\delta \left(\bar{x}-\frac{1}{T}{\int }_0^{T}dtx(t)\right)e^{S[x(t)]},$$

$$V_{\bar{x}}(x)=w(\bar{x})+\frac{m}{2}{\omega }^2(\bar{x})(x-\bar{x}{)}^2,$$

$${\bar{\rho }}_{\bar{x}}(x_T,x_0,T)$$

$$={\int }_{x(0)=x_0}^{x(T)=x}\mathcal{D}[x(t)]\delta \left(\bar{x}-\frac{1}{T}{\int }_0^{T}dtx(t)\right)e^{S_{\bar{x}}[x(t)]},$$

$$S_{\bar{x}}[x(t)]=-\frac{1}{\hslash }{\int }_0^{T}dt\left[\frac{m}{2}{\dot{x}}^2(t)+V_{\bar{x}}(x(t))\right],$$

$${\bar{\rho }}_{\bar{x}}(x_T,x_0,T)=\sqrt{\frac{m}{2\pi \beta {\hslash }^2}}{e}^{-\beta w(\bar{x})}\frac{f}{\sinh f}\times$$

$$\frac{1}{\sqrt{2\pi \alpha }}\exp \left[-\frac{{\xi }^2}{2\alpha }-\frac{m\omega \coth f}{4\hslash }(x_T-x_0{)}^2\right],$$

$$\alpha (\bar{x})=\frac{\hslash }{2m\omega (\bar{x})}\left(\coth f(\bar{x})-\frac{1}{f(\bar{x})}\right).$$

$${\bar{\rho }}_{\bar{x}}(x_0,x_0,T)=\sqrt{\frac{m}{2\pi \beta {\hslash }^2}}{e}^{-\beta w(\bar{x})}\frac{f}{\sinh f}\times$$

$$\frac{1}{\sqrt{2\pi \alpha }}\exp \left[-\frac{{\xi }^2}{2\alpha }\right],$$

$$\mathcal{Z}\equiv$$

$$\sqrt{\frac{m}{2\pi \beta {\hslash }^2}}\int d\bar{x}{e}^{-\beta V_{\mathrm{eff}}(\bar{x})},$$

$$V_{\mathrm{eff}}(\bar{x})=w(\bar{x})+\frac{1}{\beta }\ln \frac{\sinh f(\bar{x})}{f(\bar{x})}.$$

$$⟨⟨V(\bar{x}+\xi )⟩⟩$$

$$=w(\bar{x})+\frac{m}{2}{\omega }^2(\bar{x})\alpha (\bar{x}),$$

$$⟨⟨V^{\prime \prime }(\bar{x}+\xi )⟩⟩$$

$$⟨⟨F(\bar{x}+\xi )⟩⟩$$

$$=e^{\frac{\alpha (\bar{x})}{2}{\partial }_x^2}F(\bar{x}),$$

Full text

Path integrals; Semi-classical methods; Stochastic processes; Arrow-Debreu pricing; Derivative pricing; Zero-coupon bonds; Black-Karasinski model; Inhomogeneous Brownian Motion; GARCH.

Introduction

Path integrals Feynman et al. (2010), also known as Wiener integrals in stochastic calculus Wiener (1921a, b); Kac (1966), are a well-established mathematical formalism which has been used for a long time in Physics to develop accurate approximations and efficient computational techniques Kleinert (2009).

Among these, so-called semi-classical methods Kleinert (2009) play a central role. These approximations can be developed in several ways which, while sharing the same limiting behavior, lead to genuinely different results. The renowned Wentzel-Kramers-Brillouin approximation Wentzel (1926); Kramers (1926); Brillouin (1926), which is equivalento to a saddle-point approximation of the path integral Kleinert (2009); Rajaraman (1975); Kakushadze (2015), and the Wigner-Kirkwood expansion Wigner (1932); Kirkwood (1933); Fujiwara et al. (1982); Hillery et al. (1984), are well-known theoretical devices in this context.

A prominent role among semi-classical approximations is played by so-called effective potential methods Feynman et al. (2010); Feynman (1998) based, borrowing renormalization group ideas, on ‘integrating out’ the fluctuations around a ‘classical’ trajectory. Although exact in principle, the calculation can be performed only at some level of approximation, using a perturbation scheme in which the choice of the unperturbed system plays a crucial role in the quality of the approximation.

A particularly successful effective potential approximation is the one stemming on a simple and nice idea originally due to Feynman Feynman et al. (2010) and independently developed by Giachetti and Tognetti Giachetti and Tognetti (1985) and Feynman and Kleinert Feynman and Kleinert (1986) (GTFK), which is based on a self-consistent (non-local) harmonic approximation of the effective potential in a sense that will become clear in the following sections.

Basically, the GTFK effective potential is employed within the usual classical formalism, but accounts for the quantum nature of a system through suitable renormalization parameters it contains; hence, the approximation does not immediately lead to final results, but reduces a quantum-mechanical problem to a classical one, to be treated by any known method. Physicists know that this amounts to an enormous simplification.

The most appealing aspect is that the classical behavior is fully accounted for by the GTFK potential, so it opened the way to face challenging quantum systems whose classical analogues were known to be characterized by peculiar nonlinear excitations, e.g., those dubbed solitons in 1D or vortices in 2D. The latter are the ‘engine’ of a topological phase transition, for the study Kosterlitz and Thouless (1973) of which Michael Kosterlitz and David Thouless (KT) earned the 2016 Nobel prize. By the GTFK method it has been possible to establish that some real magnetic compounds do show a KT transition.

Other quantum systems that were succesfully treated by (suitable generalizations of) the same method, are frustrated antiferromagnets, e.g., the so-called two-dimensional (2D) $J_{1}$-$J_{2}$ model Capriotti et al. (2004), and 2D Josephson-junction arrays, which can be artificially fabricated, also with the inclusion of resistors; in the latter case, the effective potential could be naturally extended to account for the related dissipative coupling with the environment Cuccoli et al. (1997).

The connection between the so-called euclidean path integrals Feynman et al. (2010); Kleinert (2009), namely those employed to describe the thermodynamics of quantum systems, and the formalism of derivatives pricing has also been known since the seminal papers of Linetsky (1997) and Bennati et al. (1999) (see also the recent review Kakushadze (2015)). In particular, it is a known fact that a variable following a non-linear diffusion process can be described by the same formalism used to model the finite-temperature properties of a quantum particle in a potential which is linked to the drift of the diffusion, where the role of the mass is played by the inverse of the volatility squared, that of the temperature by the inverse of time and that of quantum fluctuations by the Brownian noise Bennati et al. (1999). The interest in financial engineering for the path-integral formalism mainly stems from the possibility of developing accurate approximation schemes, that are not otherwise available, or known, in traditional formulations of stochastic calculus Bennati et al. (1999); Kakushadze (2015); Capriotti (2006).

In this paper, we will consider the application of the GTFK method to generalized short-rate models of the form $r_{t}=r(Y_{t})$ with $Y_{t}$ following the non-linear diffusion process specified by the following stochastic differential equation (SDE)

[TABLE]

for $t>0$, where $\mu(Y_{t})$ and $\sigma(Y_{t})$ are the drift and volatility functions, respectively, $Y_{0}=y_{0}$, and $W_{t}$ is a standard Brownian motion.

Short-rate models are of paramount importance in financial modeling, providing the foundation of many approaches used for pricing of both interest rate and credit derivatives (Andersen and Piterbarg, 2010; O’Kane, 2010). In particular, celebrated affine models (Duffie et al., 2000) like those of (Vasicek, 1977), (Hull and White, 1990) and (Cox et al., 1985), play a prominent role. This is mainly due to their analytical tractability allowing one to derive closed-form expressions for fundamental building blocks like zero-coupon bonds or, in the context of default intensity models (O’Kane, 2010), survival probabilities.

Unfortunately, the availability of closed-form solutions comes often at the price of less than realistic properties of the underlying rates. For instance, Gaussian models such as those of Vasicek (1977) and Hull and White (1990), when calibrated to financial data, typically imply that rates can assume negative values with sizable probabilities. While this can be possibly not a problem for interest-rate models, especially in a low interest-rate environment, it is not consistent with absence of arbitrage in the context of default intensity models (O’Kane, 2010). On the other hand, square-root diffusions such as that of Cox et al. (1985) - while guaranteed to be non-negative - may give rise to distributions of the par swap rate, see Andersen and Piterbarg (2010); Li et al. (2018), that do not admit values below a finite threshold and may be considered therefore unrealistic.

Unfortunately, more realistic models lacks the same degree of analytical tractability as that shown by affine models. As a result, although widely used in practice, their implementations rely on computationally intensive partial differential equations (PDE) or Monte Carlo (MC) methods for the calculation of bond prices or survival probabilities. This is particularly onerous in the context of multi factor problems, notably the ones involving the calculation of valuation adjustments (XVA), cf. Gregory (2010), that are currently very prominent in financial engineering. Indeed, these applications require Monte Carlo simulations and, e.g., the valuation of conditional bond prices or survival probabilities at different points of the simulated paths, which are expensive to compute for models that lack closed-form solutions for these quantities. In this context, reliable analytical approximations are particularly important to reduce the numerical burden associated with these computations.

More specifically, in this paper we will focus on developing approximations of the so-called (generalized) Arrow-Debreu (AD) densities, see (Andersen and Piterbarg, 2010; Karatzas and Shreve, 1991), also known as Green’s functions, which are the fundamental building blocks for pricing contingent claims. These are defined, in this setting, as

[TABLE]

where $\lambda$ is a real number, and $\delta(\cdot)$ is the standard Dirac’s delta function. This, for $\lambda=0$, gives the transition density, of paramount importance for maximum-likelihood estimations in econometrics Aït Sahalia (1999), such that

[TABLE]

The price at time $t=0$ of a European option with expiry $T$ and payout of the form $P(r_{T})$,

[TABLE]

can be obtained by integrating the product of the payout function and the ($\lambda=1$) AD density over all the possible values of the short rate at time $T$, namely

[TABLE]

where the integration is performed over the range of the function $y_{T}=r^{-1}(r_{T})$. In particular, the moment generating function for the random process $\int_{0}^{T}\-du\,r_{u}$ can be obtained for $P\equiv 1$,

[TABLE]

which, for $\lambda=1$, gives the value at time $t=0$ of a zero-coupon bond with maturity $T$ Andersen and Piterbarg (2010). In the context of default intensity models, where the default of a firm is modeled by the first arrival of a Poisson process with time-dependent intensity $r_{t}$, O’Kane (2010), Eq. (6) for $\lambda=1$ represents the survival probability up to time $T$, conditional on survival up to time $t=0$. This is the fundamental building block for the evaluation of cash flows that are contingent on survival or default, see O’Kane (2010).

The structure of the paper is as follows. We start by reviewing the formalism of the GTFK effective potential method in the context of the path-integral formulation of quantum statistical mechanics. We then make the connection between the formalism used in quantum Physics and the one used in finance by reviewing the path-integral formulation of AD densities for non-linear diffusion and we show how the GTFK approximation can be used in the mathematical setting of stochastic calculus in order to develop a semi-analytical approximation for the generalized AD densities (2), and zero-coupon bonds (6) for non-linear diffusion of the form (1). Remarkably, the GTFK method, yielding exact results in the limit of zero volatility and time to maturity as any semi-classical approximation, is also exact whenever the drift potential is quadratic, which means it is exact, as we will recall, for the Vasicek Vasicek (1977) and quadratic model Kakushadze (2015). We finally illustrate the remarkable accuracy of the GTFK method for models for which an analytical solution is not available via the application to the so-called Black-Karasinski (BK) model Black and Karasinski (1991) and the so-called GARCH linear stochastic differential equation (SDE) Li et al. (2018); Capriotti et al. (2019), both of particular relevance for the valuation of credit derivatives.

Effective Potential Approximation in Quantum Statistical Mechanics

We start by recalling the path-integral formalism of quantum thermodynamics for a non-relativistic particle of mass $m$ described by the standard Hamiltonian

[TABLE]

where $\hat{x}$ and $\hat{p}$ are the canonical coordinate and momentum operators such that $[\hat{x},\hat{p}]=i\hbar$, with $\hbar$ the reduced Planck’s constant, and where $V(\hat{x})$ is the potential the particle is subject to.

The quantum thermodynamical properties of the particle at temperature ${\cal T}$ can be described by the density matrix Feynman et al. (2010),

[TABLE]

where $\beta=1/k_{B}{\cal T}$, with $k_{B}$ the Boltzmann’s constant. The elements of the density matrix, in the coordinate representation, can be expressed in terms of Feynman’s path integral Feynman et al. (2010) as

[TABLE]

where the path integration is defined over all paths $x(t)$ such that $x(0)=x_{0}$ and $x(T)=x_{T}$, with $T=\beta\hbar$ the so-called euclidean time and the functional

[TABLE]

is the euclidean action. The functional integration in Eq. (9), is formally defined as the limit for $N\to\infty$ of the expression

[TABLE]

with $\Delta t=T/N$, $x_{N}\equiv x_{T}$ and

[TABLE]

Although the evaluation of the path integral in Eq. (9) is possible just in a few cases for simple potentials, the formalism allows for new kinds of approximations. In particular, here we pursue an approximation stemming on an idea originally due to Feynman, that consists in classifying the paths according to an equivalence relation, and consequently decompose the integral into a first sum over all paths belonging to the same class, and a second one over all the equivalence classes. In particular, equivalent paths are those who share the average point, defined as the functional

[TABLE]

so that each equivalence class is labelled by a real number $\bar{x}$ representing the common average point and we can factor out in Eq. (9) an ordinary integral over $\bar{x}$, namely

[TABLE]

where the reduced density matrix

[TABLE]

represents the contribution to the path integral in Eq. (9) that comes from those paths that have $\bar{x}$ as average point.

As the path integration has been reduced to paths belonging to the same class, we can develop a specialized approximation for each class. In particular, the GTFK method approximates the potential in the action Eq. (10) with a quadratic potential in the displacement from the average point $\bar{x}$

[TABLE]

where the parameters $w(\bar{x})$ and $\omega^{2}(\bar{x})$ are to be optimized so that the trial reduced density matrix

[TABLE]

with the action given by

[TABLE]

best approximates the reduced density matrix in Eq. (Effective Potential Approximation in Quantum Statistical Mechanics). Note that one does not need to include a linear term in the trial potential (16), since it would give a vanishing contribution to the trial action (18), due to the very definition of $\bar{x}$.

The path integral in Eq. (Effective Potential Approximation in Quantum Statistical Mechanics), corresponding to the harmonic action (18) can be worked out analytically Cuccoli et al. (1995a), giving

[TABLE]

where $\xi=(x_{T}+x_{0})/2-\bar{x}$, $f=\beta\hbar\omega(\bar{x})/2$ and

[TABLE]

The diagonal elements of the reduced density matrix read in particular

[TABLE]

taking a suggestive form in terms of a Gaussian distribution with mean $\bar{x}$ and variance $\alpha(\bar{x})$, describing the fluctuations around the average point. In particular, the so-called partition function, ${\cal Z}$, Feynman (1998) assumes the classical form

[TABLE]

where the GTFK effective potential reads:

[TABLE]

In order to close the approximation we still need to devise an optimization scheme for the parameters $w(\bar{x})$ and $\omega(\bar{x})$ in Eq. (16). For example, we could simply identify the trial potential (16) with the expansion of $V(\bar{x})$ up to second order by setting $w(\bar{x})=V(\bar{x})$ and $\omega(\bar{x})=V^{\prime\prime}(\bar{x})$ for any $\bar{x}$. However, this approximation has limitations. For instance, it can happen that $V^{\prime\prime}(\bar{x})$ is negative: in this case, writing $f=\beta\hbar\omega/2$ as $f=i\phi$, $\alpha$ can be analytically continued as $\alpha=(\beta\hbar^{2}/4m)(1/\phi^{2}-\cot\phi/\phi)$, which diverges to $+\infty$ for $\phi\to\pi^{-}$ (or $f^{2}\to-\pi^{2}$) and is negative for $\phi>\pi$ ($f^{2}<-\pi^{2}$). As a consequence, if $\omega^{2}(\bar{x})$ is negative, for sufficiently large time horizons $T$ we have $f^{2}<-\pi^{2}$ and $\alpha(\bar{x})<0$. In this situation, the reduced density matrix (Effective Potential Approximation in Quantum Statistical Mechanics) is not well defined and the approximation breaks down.

A more robust approximation can be devised by observing that the Gaussian density $\bar{\rho}_{\bar{x}}(x_{0},x_{0},T)$ has to be close to $\rho_{\bar{x}}(x_{0},x_{0},T)$, so that $V_{\bar{x}}(x)$ must approximate $V(x)$ not only at $\bar{x}$: this is accomplished by requiring the equality of the Gaussian averages of the true and the trial potentials, and of their derivatives up to the second one

[TABLE]

with the short-hand notation

[TABLE]

and $\alpha(\bar{x})$ given by Eq. (20). The equations above impose that the expectation value according to the Gaussian probability distribution in Eq. (Effective Potential Approximation in Quantum Statistical Mechanics) of the potential and of its second order expansion are in agreement with each other, for every value of $\bar{x}$. Under the GTFK approximation the quantum effects are embedded in the notion of the effective potential (23) which is a renormalized version of the potential $V(x)$ where $\alpha(\bar{x})\equiv\langle\!\langle\xi^{2}\rangle\!\rangle$ – representing the average quadratic fluctuations around $\bar{x}$ due to the quantum effects – is the renormalization parameter. Note that Eq. (25) is self consistent, meaning that its solution $\omega^{2}(\bar{x})$ in turn determines the variance $\eqref{eq.alpha}$.

It can be shown that the above determination of the parameters $w({\bar{x}})$ and $\omega({\bar{x}})$ satisfies a variational principle based on the so-called Jensen-Feynman inequality, ${\cal Z}\geq{\cal Z}_{0}\,e^{\langle S-S_{0}\rangle_{0}}$, where the functional average is taken with whatever trial action $S_{0}$, ${\cal{Z}}_{0}$ being the corresponding partition function. Indeed, taking $S_{0}=S_{\bar{x}}$ and maximizing the right-hand side of the inequality one just finds Eqs. (24) and (25).

The GTFK method, becomes exact in both limits of high-temperature $\beta\to 0$ and vanishing quantum effects $\hbar/m\to 0$, for which the parameter $\alpha$ vanishes as $\beta\hbar^{2}/12m$ and the effective potential (23), coincides with the exact classical potential:

[TABLE]

so that the partition function in Eq. (Effective Potential Approximation in Quantum Statistical Mechanics) coincides with the well-known exact classical result Feynman (1998).

The effective potential can be compared to the semiclassical effective potential introduced by Wigner and Kirkwood Wigner (1932); Kirkwood (1933); Fujiwara et al. (1982); Hillery et al. (1984) (WK), that was substantially found as an expansion in $\beta$ and $\hbar$ of the exact classical effective potential $V_{\rm ex}$, defined such that the quantum density bears the classical form

[TABLE]

The WK expansion is in principle exact, but only the first few terms are practically affordable, and while lowering the temperature all terms soon diverge. One has indeed Jizba and Zatloukal (2014)

[TABLE]

This apparently disagrees from the expansion (27) because the comparison is a little subtle: indeed, $V_{\rm eff}$ has not to be directly compared with $V_{\rm{ex}}$, because, in order to obtain $\rho(x_{0},x_{0},T)$ one cannot integrate over $x_{0}$ as made in Eq. (Effective Potential Approximation in Quantum Statistical Mechanics), but rather over $\bar{x}$. Accounting for this Kleinert (1986), the WK and the GTFK effective potentials do agree Vaia and Tognetti (1990); Cuccoli et al. (1992). Similarly, GTFK is distinct from the exponential power series expansion of Makri and Miller (1989), previously applied successfully in the financial context Capriotti (2006); Stehlíková and Capriotti (2014); Capriotti et al. (2019), and which we will use as one of the benchmarks when discussing our numerical results.

With respect to these approaches, the GTFK method has a strong advantage: it still gives a meaningful representation of the thermodynamics down to zero temperature, where it is equivalent to the so-called self-consistent harmonic approximation Koehler (1966a, b), that was initially applied to quantum crystal lattices. Therefore, increasing temperature from zero the accuracy increases more and more, because the renormalization parameter $\alpha(\bar{x})$ decreases. The price to pay is that one still has to solve the classical problem with the effective potential, but this is nevertheless a huge simplification, especially in view of the plenty of methods that have been developed to treat classical systems. In particular, thanks to the fact that the nonlinear character of the potential is kept, the GTFK approach allows for studying quantum systems whose classical counterpart is characterized by nonlinear excitations (solitons, vortices) and constitutes a much simpler and clearly interpretable alternative to heavy numerical approaches, such as Quantum Monte Carlo.

The GTFK approach is also distinct from other semi-classical path-integral approximations, like the Wentzel-Kramers-Brillouin (WKB) Wentzel (1926); Kramers (1926); Brillouin (1926) or the equivalent saddle-point approximations Kleinert (2009); Rajaraman (1975); Kakushadze (2015), which are based on a power-series expansion of the action around the classical trajectory $x_{\rm{c}}(t)$ rather than around the average point, i.e., the density matrix, Eqs. (9) and (10), is expressed as

[TABLE]

where $x_{\rm{c}}(t)$ obeys the classical equation of motion $\delta{S}/\delta{x(t)}=0$ and satisfies the boundary conditions $x_{\rm{c}}(0)=x_{0}$ and $x_{\rm{c}}(T)=x_{T}$, while the path summation is over closed paths $\tilde{x}(t)=x(t){-}x_{\rm{c}}(t)$ with the expanded action

[TABLE]

The WKB approximation is exact for a quadratic potential, and, the first term being of order $\hbar^{-1}$, it can include the effect of tunneling (for instance, in a double-well potential) at variance with the GTFK; however, one has to consider that it is not crucial to account for tunneling effects, as they are soon overwhelmed by quantum thermal fluctuations and are practically absent in many-body systems; moreover, beyond few relatively simple cases, the evaluation of the path integral (31) is generally hard, mainly due to the dependence of $\tilde{S}$ upon the classical path. On the other hand, the non-local nature of the GTFK approximation yields the possibility of tuning two families of parameters, $w(\bar{x})$ and $\omega(\bar{x})$, allowing one to look for the best approximation of the true action in a richer space, while preserving the property of being exact in the classical limit and for harmonic actions. By ‘richer space’ we mean that the trial action, thank to its dependence on the average-point functional, is much more general than the local actions corresponding to physical potentials. The GTFK can also be systematically improved, at least in principle, without suffering from the divergencies appearing instead in most perturbative approaches Kleinert (2009).

The generalizations of the GTFK approach to many degrees of freedom, as well as to Hamiltonian systems Cuccoli et al. (1992); Cuccoli et al. (1995a), have found numerous applications in Physics and Physical Chemistry. Besides the tests on simple models with one degree of freedom Feynman and Kleinert (1986); Janke and Kleinert (1986, 1987); Vaia and Tognetti (1990), it is noteworthy that the very first paper regarded the 1D sine-Gordon model Giachetti and Tognetti (1985); Giachetti et al. (1988a), whose classical version is characterized by the existence of topological nonlinear excitations, the solitons, that determine an anomaly of thermodynamic quantities like the specific heat: the GTFK method allowed for the first time to quantify the same anomaly for the quantum system, and was shown to agree with the outcomes of hard Quantum Monte Carlo calculations Giachetti et al. (1988a) and to admit a renormalized continuum limit in agreement with exact ‘Bethe Ansatz’ calculations Giachetti et al. (1988b).

Among many accomplishments, one should mention the quantitative explanation Cuccoli et al. (1991) of experimental data regarding a quasi-1D magnet CsNiF3, that behaves similarly to the sine-Gordon model, while a major one has been the study of 2D quantum anisotropic magnets Cuccoli et al. (1995b); Cuccoli et al. (1998), whose classical counterpart shows the topological phase transition studied Kosterlitz and Thouless (1973) by Kosterlitz and Thouless (KT); the GTFK approach allowed also to quantitatively characterize Cuccoli et al. (2006) earlier experiments, showing that magnetic and calorimetric measurements performed in 1983 were the first known experimental observation of KT behavior in a real magnet; a further success in the magnetic realm was providing a consistent picture of the elusive Ising phase transition in a frustrated model such as the 2D quantum $J_{1}$-$J_{2}$ Heisenberg antiferromagnet Capriotti et al. (2004).

2D Josephson-junction arrays are also typical KT systems: the effective potential was extended to include the dissipative effect of resistive shunts among the junctions used in experiments, getting quantitative accuracy for the phase diagram Cuccoli et al. (2000). The versatility of the GTFK potential is witnessed also by recent applications in the theoretical interpretation of thermal expansion measurements obtained by $x$-ray absorption spectroscopy in alloys Yokoyama and Eguchi (2013); Yokoyama et al. (2018).

Path-Integral formulation of Stochastic Calculus

In this section, we briefly review how the formalism of stochastic calculus can be recast in the language of path-integrals in Euclidean time, focussing for simplicity on the case of a single SDE as in Eq. (1). As a first step, in order to simplify the derivation, it is convenient to transform the original process into an auxiliary one, $X_{t}$, with constant volatility $\sigma$. Following Aït Sahalia (1999), this can be achieved in general through the so-called Lamperti’s transform

[TABLE]

A straightforward application of Ito’s Lemma gives the stochastic differential equation satisfied by $X_{t}$ for $t\geq 0$:

[TABLE]

where

[TABLE]

Here, $y=\gamma^{-1}(x)$ is the inverse of the transformation (32). The generalized AD density (2) for the processes $X_{t}$ and $Y_{t}$ are related by the Jacobian associated with (32) giving

[TABLE]

It is well known, see e.g., Andersen and Piterbarg (2010); Karatzas and Shreve (1991), that the generalized AD density (2) for the process (33) satisfies the following conjugate forward (Fokker-Planck) partial differential equation (PDE)

[TABLE]

with the initial condition $\psi_{\lambda}(x_{t},x_{0},0)=\delta(x_{0}-x_{t})$.

A path-integral representation of the AD density can be constructed Bennati et al. (1999) starting from the Euler approximation, correct up to $O(\Delta t)$, for the solution of the Fokker-Planck PDE (Path-Integral formulation of Stochastic Calculus)

[TABLE]

Using the Markov property, the equation above gives a prescription to write the solution of the Fokker-Planck equation in the form of a convolution product of short-time AD densities as:

[TABLE]

with $\Delta t=T/N$, $x_{N}\equiv x_{T}$ and

[TABLE]

where the term

[TABLE]

arises, at order $O(\Delta t)$, from using the analytically convenient Stratonovich mid-point discretization Bennati et al. (1999). As a result, the limit $N\to\infty$ of Eq. (Path-Integral formulation of Stochastic Calculus) can be formally written as

[TABLE]

where

[TABLE]

has the same form of the density matrix in Eq. (9), the functional

[TABLE]

has the same form of the euclidean action in Eq.(10),

[TABLE]

can be called drift potential and we have defined

[TABLE]

in order to give Eq. (43) a suggestive Lagrangian structure as in Eq. (10).

The key observation is that the path integral in Eq. (42) is formally equivalent to density matrix in Eq. (9) describing the quantum termodynamics of a particle of mass $m=\hbar/\sigma^{2}$ in a potential $\hbar V(x)$, at temperature ${\cal T}=\hbar/k_{B}T$ (such that $\beta\hbar=T$).

The GTFK can be therefore applied straightforwardly and here for convenience we restate the results with the notation of stochastic calculus:

[TABLE]

where $\xi=(x_{T}+x_{0})/2-\bar{x}$, $f=\omega(\bar{x})T/2$ and

[TABLE]

with $w(\bar{x})$ and $\omega(\bar{x})$ solutions of the self-consistent equations:

[TABLE]

The GTFK method, becomes exact in the limit of short time to maturity $T\to 0$ and vanishing volatility $\sigma\to 0$ for which the parameter $\alpha$ vanishes as $\sigma^{2}T/12$. Furthermore, given the form of the chosen trial potential, for harmonic actions, the GTFK approximation is, in fact, exact. This is for instance the case for the Vasicek model Vasicek (1977) as it will be illustrated in the next section.

Numerical Results

In this section we illustrate the effectiveness of the GFTK approach by discussing its application to a few diffusions processes of the form (1), starting from two cases in which the method gives exact results, namely the Vasicek and the so-called quadratic short-rate model. We then discuss the Black-Karasinski (BK) Black and Karasinski (1991) and GARCH linear SDE model Li et al. (2018); Capriotti et al. (2019) - for which the AD density (2) or zero-coupon bonds (6) are not know analytically - by presenting the comparison of the GTFK results with those obtained by solving numerically the relevant PDEs and by employing other approximations.

.1 Vasicek model

The Vasicek model Vasicek (1977) is a simple example of affine process Duffie et al. (2000)

[TABLE]

where $a$ is the mean-reversion speed, $b$ the mean-reversion level, $\sigma$ the volatility, and $r(X_{t})=X_{t}$. The drift potential (44) is given by the quadratic form

[TABLE]

The path integral for quadratic potentials is known to be analytically tractable and corresponds in quantum Physics to the so-called harmonic oscillator Feynman et al. (2010). In this case, the GTFK self-consistent conditions (48) and (49) read:

[TABLE]

and the reduced density matrix (Path-Integral formulation of Stochastic Calculus) reads:

[TABLE]

with $\alpha=\sigma^{2}/2a(\coth f-1/f)$, $f=aT/2$, both independent of $\bar{x}$. The integral over $\bar{x}$ in Eq. (14) can then be performed analytically giving, after a somewhat tedious but straightforward calculation,

[TABLE]

where $\bar{\sigma}^{2}=\sigma^{2}(1-\exp(-2aT))/2a$, in agreement with the known result Jamshidian (1989).

.2 Quadratic Short Rate Model

In the quadratic short rate model, the short rate is defined as

[TABLE]

with $X_{t}$ following the OU diffusion (50), which is positive definite for $\beta>0$ and $\gamma^{2}<4\beta$. In this case, the drift potential (44) reads

[TABLE]

while the GTFK conditions, (48) and (49), can be determined as

[TABLE]

which, as in the Vasicek model discussed above give a frequency $\omega$ that is not dependent on the average point and a function $w(\bar{x})$ which is quadratic in $\bar{x}$. Also in this case the Gaussian integration can be performed analytically leading to the exact result.

.3 Black-Karasinki Model

The BK Black and Karasinski (1991) model is a conspicuous example of a diffusion that is particularly suitable for financial applications because the short rate at any time horizon follows an intuitive lognormal distribution. Unfortunately, it lacks the same degree of analytical tractability as that shown by affine models. As a result, although widely used in practice, BK implementations rely on computationally intensive numerical simulations based on PDE or Monte Carlo Andersen and Piterbarg (2010).

The short rate in the BK model is defined as

[TABLE]

with $X_{t}$ following the OU diffusion (50). In this case, the drift potential (44) reads

[TABLE]

while the GTFK conditions, (48) and (49), can be determined with some straightforward algebra as

[TABLE]

with the second to be solved self-consistently with the renormalization parameter in Eq. (47).

In Fig. 1 we plot the GTFK self-consistent parameters $\omega^{2}(\bar{x})$, and $\alpha(\bar{x})$ and the diagonal trial reduced density matrix $\bar{\rho}_{\bar{x}}(x_{0},x_{0},T)$ in Eq. (Path-Integral formulation of Stochastic Calculus) as a function of the average point $\bar{x}$ for different strength of the diffusive effects, namely of the time to maturity and volatility. For weak diffusive effects, the parameter $\alpha(\bar{x})$ is relative small and the trial reduced density matrix has a sharp peak around $x_{0}$. In this region, both $\alpha(\bar{x})$ and $\omega^{2}(\bar{x})$ display a weak dependence on $\bar{x}$ which signals the adequacy of a local harmonic approximation to capture the purely diffusive effects in the problem. However, as the diffusive effects increase, with larger volatility and/or time to maturity, the renormalization parameter $\alpha(\bar{x})$ increases, the trial density broadens and both $\alpha(\bar{x})$ and $\omega^{2}(\bar{x})$ display a more marked dependency on the average point $\bar{x}$, signaling that a non-local approximation is needed to best capture the diffusive effects given an harmonic ansatz of the effective potential.

An illustration of the accuracy of the BK AD densities (2) obtained with the GTFK approximation is displayed for a high volatility case in Fig. 2, for different values of time to maturity, by comparing with a numerical solution of the Fokker-Planck equation (Path-Integral formulation of Stochastic Calculus). Here we observe that the GTFK approximation is hardly distinguishable from the PDE result up to $T=5$, and remains very accurate even for large time horizons.

This is also confirmed by the results for zero-coupon bonds (6) reported in Table 1 illustrating how the GTFK method compares favorably with the results obtained with recently proposed semi-analytical approximations, namely the Exponent Expansion (EE) Stehlíková and Capriotti (2014), and the Karhunen-Loéve (KL) expansions Daniluk and Muchorski (2016) when benchmarked agains a numerical solution of the associated PDE. In particular, for short time horizons, the GTFK approximation has comparable accuracy with the EE. For larger time horizons, the GTFK compares better and better and remains very accurate even when the EE, which has a finite convergence ratio in $T$, eventually breaks down. Similarly, the GTFK method has better accuracy than the first order KL expansion, and comparable accuracy with the second order KL expansion for short time horizons, while it has significantly better accuracy for large time horizons. Even for time horizons as large as 20 years the GTFK approximation produces zero-coupon bond prices within 50 basis points from the exact result, as also illustrated in Fig. 3. Similar conclusions can also be drawn when comparing with other recently proposed approaches as those in Refs. Tourrucôo et al. (2007); Antonov and Spector (2011).

.4 GARCH Linear SDE

As an example of a more challenging application, we then consider the GARCH linear SDE or Inhomogenous Geometric Brownian Motion Li et al. (2018); Capriotti et al. (2019) model, which is a special case of the so-called Continuous Elasticity of Variance (CEV) diffusion Cox and Ross (1976), namely

[TABLE]

with $r(Y_{t})=Y_{t}$.

The process defined by the SDE in Eq. (62) can be shown to be strictly positive Kloeden and Platen (1992). As a result, like the BK model, it is well suited to represent default intensities. It can be also shown to have probability density profiles which are more intuitive than those generated by the widely used square-root processes Cox et al. (1985); Li et al. (2018). Unfortunately, even if it can be solved exactly Kloeden and Platen (1992) it does not admit a closed form for the (generalized) AD prices (2).

Under the Lamperti’s transformation (32) for this process, namely $X_{t}=\log Y_{t}$, Eq. (62) reads

[TABLE]

with

[TABLE]

The drift potential (44) associated with the SDE (63) reads therefore

[TABLE]

which is related to the so-called Morse potential Bentaïba et al. (1994). The GTFK conditions, (48) and (49), can be determined with some straightforward algebra as

[TABLE]

Examples of AD densities (2) obtained with the GTFK approximation for the GARCH linear SDE are displayed in Fig. 4, for different values of the diffusion parameters, with a comparison with a numerical solution of the Fokker-Planck equation (Path-Integral formulation of Stochastic Calculus). Here we observe that the GTFK approximation, as in the BK case, is difficult to distinguish from the PDE result up to several years maturity, and for large enough volatilities. As in the BK case, the accuracy of the approximations depends on the chosen model parameters, and the maturity being considered. The approximation becomes less accurate for larger maturities $T$ and volatility. The behaviour with respect to the mean-reversion speed $a$ is instead less clear-cut as this parameter affects both the variance of the process and the non-linearity of the drift potential (.4).

The accuracy of the GTFK method for the GARCH linear SDE is also illustrated for zero-coupon bonds (6) in Table 2 and 3 for two sets of model parameters, showing how the GTFK method compares favorably with the results obtained with recently proposed semi-analytical approximations, namely the EE Capriotti et al. (2019), when benchmarked agains a numerical solution of the associated PDE. In general, although less accurate than in the BK case, due to the more complex form of the drift potential (.4), the approximation produces satisfactory results for maturities up to several years even in regimes of high volatility.

Conclusions

An effective-potential path-integral formalism of quantum statistical mechanics – dubbed GTFK after the authors Giachetti and Tognetti (1985); Feynman and Kleinert (1986) who originally introduced it – has been widely utilized in Physics for the study of the quantum thermodynamics of condensed matter systems. The method is based on a self-consistent harmonic approximation of the pure-quantum contributions to the thermodynamics, while fully accounting for the classical behaviour of the system Cuccoli et al. (1995a). As a semiclassical approach, it is exact in the high-temperature and zero-quantum fluctuations limits but, remarkably, it also gives a meaningful representation in the zero-temperature limit, where it is equivalent to a self-consistent harmonic approximation of the potential.

By exploiting the path-integral formulation of stochastic calculus, we have shown how the GTFK approach can be used to develop an accurate semi-analytical approximation of (generalized) Arrow-Debreu densities, and zero-coupon bonds for non-linear diffusions. The method is exact in the limit of zero volatility, zero time to maturity, and for Ornstein-Ulhenbeck diffusions.

The GTFK provides remarkably accurate results for the Black-Karasinski and GARCH linear SDE for interest rates or default intensities, even for high volatilities and long time horizons, with results that compare favorably with previously presented approximation schemes Tourrucôo et al. (2007); Antonov and Spector (2011); Daniluk and Muchorski (2016); Stehlíková and Capriotti (2014); Capriotti et al. (2019), with expressions that are more compact and easier to compute, and less severe limitations arising from a finite convergence radius in the time to maturity or volatility. Similarly to the approach in Capriotti (2006), the range of application of the expansion can be further extended to even larger time horizons by means of a fast numerical convolution Bennati et al. (1999).

The GTFK approximation can be potentially improved in one of two ways: by pursuing higher-order corrections as in the so-called variational perturbation theory Kleinert (2009) or by its generalization to Hamiltonian systems Cuccoli et al. (1992); Cuccoli et al. (1995a) that would allow avoiding the non-linearities in the potential introduced (e.g., as for the GARCH linear SDE) via the Lamperti’s transformation (32).

The accuracy and ease of computation of the GTFK method makes it a computationally efficient alternative to fully numerical schemes such as binomial trees, PDE or Monte Carlo for the calculation of transition densities – whether for the maximization of classical likelihoods or the computation of posterior distributions – and for the evaluation of European-style derivatives. This is of practical utility e.g., for econometric applications Aït Sahalia (1999), for speeding up pricing or calibration routines for valuation of derivatives Andersen and Piterbarg (2010) or in the context of time consuming multi-factor simulations that are common place in financial engineering in a variety of applications Hull (2017).

Acknowledgements.

It is a pleasure to acknowledge Jim Gatheral, Tao-Ho Wang and Mehdi Sonthonnax for useful discussions. The authors are grateful to Prof. Valerio Tognetti for igniting in them the passion for Path Integrals, and for his warm support throughout the years.