Abstract
The solution of the Kolmogorov backward equation is expressed as a functional integral by means of the Feynman–Kac formula. The expectation value is approximated as a mean over trajectories. In order to reduce the variance of the estimate, importance sampling is utilized. From the optimal importance density, a modified drift function is derived which is used to simulate optimal trajectories from an Itô equation. The method is applied to option pricing and the simulation of transition densities and likelihoods for diffusion processes. The results are compared to known exact solutions and results obtained by numerical integration of the path integral using Euler transition kernels. The importance sampling leads to strong variance reduction, even if the unknown solution appearing in the drift is replaced by known reference solutions. In models with low-dimensional state space, the numerical integration method is more efficient, but in higher dimensions it soon becomes infeasible, whereas the Monte Carlo method still works.
Similar content being viewed by others
Notes
Notation In the following, the components of vectors and matrices are denoted by Greek letters, e.g. \(f_{\alpha }(x), \alpha =1,\ldots ,p\), and partial derivatives by commas or subscripts, i.e. \(f_{\alpha ,\beta }:=\partial f_{\alpha }/\partial x_{\beta }= \partial _{\beta } f_{\alpha }=(f_{x})_{\alpha \beta }\). The Jacobian matrix \(\partial f/\partial x\) is denoted as \(f_{x}\). The derivative \(\phi _{x}=\phi _{,\alpha }\) of a scalar function \(\phi \) is defined as a row vector. For a bilinear form \(f'Ag\) one obtains \((f'Ag)_{x}= f_{\alpha }A_{\alpha \beta }g_{\beta ,\gamma }+ f_{\alpha ,\gamma }A_{\alpha \beta }g_{\beta }+ f_{\alpha }A_{\alpha \beta ,\gamma }g_{\beta } =f'Ag_{x}+g'A'f_{x}+ (f\otimes g)'A_{x}\). A sum convention is used for the Greek indices, i.e. \(f_{\alpha }g_{\alpha }=\sum _{\alpha }f_{\alpha }g_{\alpha }\). Latin indices denote time, e.g. \(f_{j \alpha }=f_{\alpha }(x(t_{j}))\).
The time ordering is necessary since \(H(t)=L+v\) and \(H(s)\) do not commute in general, i.e. \(H(t) H(s)-H(s) H(t)\ne 0\). The time ordering operator \(\, \mathop {T}\limits ^{\rightarrow }[H(t) H(s)]=H(s) H(t); s<t\) interchanges the order; cf. Abrikosov et al. (1963, p. 48).
\(C(x_{0},\tau _{0}):=C(\tau _{T},x_{0},\tau _{0},\delta \tau )\) is defined in Eq. (11). A variant of \(C\) fulfils \(C_{\tau _{0}}+L_{1} C+v C+O(\delta \tau )=0\), but then \(\beta (\varvec{x},\varvec{\tau })\) in (18) must be replaced by \(\exp \{\int _{\tau _{0}}^{t}v(x(\tau ),\tau )\,\mathrm{d} \tau \}\). See Appendix B, Eq. (45).
Alternatively, one can use the backward equation \(\mu _{s}(t,x,s)=\int y\; p_{s}(y,t|x,s)\,\mathrm{d}y = -L(x,s)\) \(\mu (t,x,s)\) and expand the solution \(\mu (t,x,s)=U(t,s)x\) (Dieterich et al. 1977; Aït-Sahalia 2002).In the time-invariant case one has the short time expansion \(\mu _{L}(t,x,s)=\sum _{l=0}^{L}\frac{(s-t)^{l}}{l!}L(x)^{l}x.\) The variance is given by \(\Sigma (t,x,s)=U(t,s)x x'-\mu \mu '.\)
In a narrow sense, this terminology is only used for \(v=0\), but see Karatzas and Shreve 1991, p. 369.
References
Abrikosov, A., Gorkov, L., Dzyaloshinsky, I.: Methods of Quantum Field Theory in Statistical Physics, p. 48. Dover, New York (1963)
Aït-Sahalia, Y.: Maximum likelihood estimation of discretely sampled diffusions: a closed-form approximation approach. Econometrica 70(1), 223–262 (2002)
Aït-Sahalia, Y.: Closed-form likelihood expansions for multivariate diffusions. Ann. Stat. 36(2), 906–937 (2008)
Beskos, A., Papaspiliopoulos, O., Roberts, G., Fearnhead, P.: Exact and efficient likelihood-based inference for discretely observed diffusion processes (with discussion). J. R. Stat. Soc. Ser. B 68, 333–382 (2006)
Bibby, M., Sorensen, M.: Martingale estimation functions for discretely observed diffusion processes. Bernoulli 1, 1–39 (1995)
Black, F., Scholes, M.: The pricing of options and corporate liabilities. J. Polit. Econ. 81, 637–654 (1973)
Borodin, A., Salminen, P.: Handbook of Brownian Motion—Facts and Formulae, 2nd edn. Birkhäuser-Verlag, Basel (2002)
Cameron, R.H.: A “Simpson’s Rule” for the numerical evaluation of Wiener integrals in function space. Duke Math. J. 18, 111–130 (1951)
Cameron, R.H.: A family of integrals serving to connect the Wiener and Feynman integrals. J. Math. Phys. 39(126–140), 1961 (1960)
Cameron, R.H., Martin, W.T.: Transformations of Wiener integrals under a general class of linear transformations. Trans. Am. Math. Soc. 58(2), 184–219 (1945)
Chorin, A.J.: Accurate evaluation of Wiener integrals. Math. Comput. 27(121), 1–15 (1973)
Cox, J., Ross, S.: The valuation of options for alternative stochastic processes. J. Financ. Econ. 3, 145–166 (1976)
Dieterich, W., Peschel, I., Schneider, W.R.: Diffusion in periodic potentials. Z. Phys. B Condens. Matter 27, 177–187 (1977)
Durham, G.B., Gallant, A.R.: Numerical techniques for simulated maximum likelihood estimation of stochastic differential equations. J. Bus. Econ. Stat. 20, 297–316 (2002)
Emanuel, D., MacBeth, J.: Modelling the persistence of conditional variances. J. Financ. Quant. Anal. XVII 4, 533–554 (1982)
Feller, W.: Two singular diffusion problems. Ann. Math. 54, 173–182 (1951)
Feynman, R.: Space-time approach to nonrelativistic quantum mechanics. Rev. Mod. Phys. 20, 367–387 (1948)
Gelfand, I., Yaglom, A.: Integration in functional spaces and its application in quantum physics. J. Math. Phys. 1(1), 48–69 (1960)
Girolami, M., Calderhead, B.: Riemann manifold Langevin and Hamiltonian Monte Carlo methods. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 73(2), 123–214 (2011)
Glasserman, P.: Monte Carlo Methods in Financial Engineering. Springer, New York (2003)
Holden, H., Øksendal, B., Ubøe, J., Zhang, T.: Stochastic Partial Differential Equations: A Modeling, White Noise Functional Approach, 2nd edn, ch. 4.3, p. 146. Springer, Berlin (2009)
Kac, M.: On distributions of certain Wiener functionals. Trans. Am. Math. Soc. 65, 1–13 (1949)
Kac, M.: Wiener and integration in function spaces. Bull. Am. Math. Soc. 72, 52–68 (1966)
Karatzas, I., Shreve, S.: Brownian Motion and Stochastic Calculus. Springer, Berlin (1991)
Kloeden, P., Platen, E.: Numerical Solution of Stochastic Differential Equations. Springer, Berlin (1999) (corrected third printing)
Kolmogoroff, A.: Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung. Math. Ann. 104, 415–458 (1931)
Kulikov, G., Kulikova, M.: Accurate numerical implementation of the continuous–discrete extended Kalman filter. IEEE Trans. Autom. Control (2013) (in print)
Lighthill, M.: Introduction to Fourier Analysis and Generalised Functions. Cambridge University Press, Cambridge (1958)
Liptser, R., Shiryayev, A.: Statistics of Random Processes, vol. I, ch. 7.7, 2nd edn. Springer, Berlin (2001)
Lorenz, E.: Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130 (1963)
Nelson, E.: Feynman integrals and the Schrödinger equation. J. Math. Phys. 5(3), 332–343 (1964)
Ozaki, T.: Non-linear time series models and dynamical systems. In: Hannan, E.J., Krishnaiah, P.R., Rao, M.M. (eds.) Handbook of Statistics, vol. 5, pp. 25–83. North Holland, Amsterdam (1985). http://dx.doi.org/10.1016/S0169-7161(85)05004-0
Risken, H.: The Fokker–Planck equation, 2nd edn, ch. 4.2.1–2, p. 69. Springer, Berlin (1989)
Singer, H.: Finanzmarktökonometrie. Zeitstetige Systeme und ihre Anwendung in Ökonometrie und empirischer Kapitalmarktforschung [Econometrics of Financial Markets, in German], ch. 9.6, pp. 248. Physica-Verlag, Heidelberg (1999)
Singer, H.: Parameter estimation of nonlinear stochastic differential equations: simulated maximum likelihood vs. extended Kalman filter and Itô-Taylor expansion. J. Comput. Graph. Stat. 11(4), 972–995 (2002)
Singer, H.: Simulated maximum likelihood in nonlinear continuous–discrete state space models: importance sampling by approximate smoothing. Comput. Stat. 18(1), 79–106 (2003)
Singer, H.: Continuous–discrete state-space modeling of panel data with nonlinear filter algorithms. AStA Adv. Stat. Anal. 95, 375–413 (2011)
Singer, H.: Maximum maximum likelihood estimation of continuous–discrete state-space models: Langevin path sampling vs. numerical integration. (2013) (preprint)
Steinbauer, A.: Quadrature formulas for the Wiener measure. J. Complex. 15(4), 476–498 (1999)
Wagner, W.: Unbiased Monte Carlo evaluation of certain functional integrals. J. Comput. Phys. 71, 21–33 (1987)
Yosida, K.: Functional Analysis, pp. 419–430. Springer, Berlin (1980)
Author information
Authors and Affiliations
Corresponding author
Appendices
Appendix A: Backward and forward equations
In the case of terminal condition \(h(x)=\delta (y-x)\) one can rewrite (3) in the form
\(H(y,\tau ):=L(y,\tau )+v(y,\tau )\), using the formula \(A(x)\delta (y-x)=A^{\dagger }(y)\delta (y-x)\) (Risken 1989, pp. 67–70). Equation (40) is called a fundamental solution of backward equation (1) (see Karatzas and Shreve 1991, p. 368). Differentiating w.r.t. time \(t\) one obtains the forward Kolmogorov equation (Fokker–Planck equationFootnote 5)
where
is the forward operator.
In the case of \(h\ne \delta \), differentiation of (3) w.r.t. the forward time \(t\) does not lead to a PDE, due to the time ordering. In the autonomous case, however, i.e. \(f(x,t)=f(x)\), etc., one gets the simpler solutions
This implies both the forward equation
and the backward equation
The expression
fulfils a forward equation \(C_{t}=H^{\dagger }C\) with arbitrary initial condition \(C(s,y,s)=h(y)\). Using the fundamental solution (40), one can write
Appendix B: Feynman–Kac formula
1.1 B.1 Proof and finite-dimensional approximation
The Feynman–Kac formula (5) can be proved using the Itô lemma for the ‘discounted’ function \(C^{*}=C(t,X(\tau ),\tau )\beta (s,\tau )\), where \(C(t,x,s)\) is a solution of Eq. 1 and \(\beta (s,\tau )=\exp [\int _{s}^{\tau } v(X(u),u) \,\mathrm{d}u]\) (see Karatzas and Shreve 1991, ch. 5.7). One obtains
The first term in (43) is null due to the backward equation. Integrating from \(s\) to \(t\) one obtains
Taking the conditional expectation \(E[ \cdot | X(s)=x]\) the result (5) is proved. The fundamental solution (40) is thus given as
Instead of (5) the finite-dimensional approximation (11) is used for the MC estimators. A variant of (11) is
where \(p(x_{j+1},\tau _{j} | x_{j},\tau _{j})\) are transition densities on the grid (time slices) \(\tau _{j}=s+j\delta \tau \). One can set \(t=\tau _{T}, s=\tau _{0}\). In the exponent, the simple function \(x(\tau ):=\sum _{j=0}^{T-1}x_{j+1}\chi _{[\tau _{j},\tau _{j+1})}(\tau )\) appears, where \(\chi _{A}(x)\) is the indicator function of the set \(A\). One obtains, using the backward equation \(p_{s}(x_{1},\tau _{1}|x,s)+L(x,s)p(x_{1},\tau _{1}|x,s)=0\),
Inserting the Taylor expansion \(v(x_{1},s)=v(x,s)+v_{x}(x,s)(x_{1}-x)+O(x_{1}-x)^{2}\) one finds \(C_{s}=-v(x,s)C-L(x,s)C+O(\delta \tau )\), thus (45) approximately fulfils the backward equation (1). Note that the derivation is also valid for \(h=h(x_{T},\ldots ,x_{1})\).
One can also insert the transition densities
into (11). By partial integration we obtain
where \(\beta (\varvec{x},\varvec{\tau })=\mathrm{e}^{\sum _{j=1}^{T}v(x_{j},\tau _{j})\delta \tau }\). Inserting \(\exp [L(x,\tau _{j})\delta \tau ]\exp [v(x,\tau _{j})\delta \tau ] =\exp \{[L(x,\tau _{j})+v(x,\tau _{j})]\delta \tau \}+O(\delta \tau ^{2}) =\exp \{H(x,\tau _j)\delta \tau \}+O(\delta \tau ^{2})\) for small \(\delta \tau \) (Lie–Trotter formula; cf. Nelson 1964, app. B) one obtains the product formula
This coincides with (3) in the limit \(\delta \tau \rightarrow 0\) using the Lie–Trotter formula again. In a Banach space setting, similar product formulas were obtained by Kato (cf. Yosida 1980, ch. XIV.4).
The product formula can be used for numerical evaluation, when the differential operator \(H\) is represented as a integral operator with subsequent matrix approximations.
Appendix C: Product and quotient rule for the backward operator
Appendix D: Sampling from a reference density
One could also use a Langevin equation \(\mathrm{d}\varvec{X}(u)=(\partial _{\varvec{x}}\log \, p_{2})(\varvec{X}) \,\mathrm{d}u+\sqrt{2} \mathrm{d}W(u)\) to sample from \(p_{2}\) (17) since \(\partial _{\varvec{x}}\log C(x_{0},\tau _{0})=0\), but then the random vectors \(\varvec{X}(u)=[X_{T}(u),\ldots , X_{1}(u)]\) are autocorrelated over \(u\) (cf. Girolami and Calderhead 2011). Moreover, the density \(p_{2}\) appearing in (15) contains the unknown quantity \(C(x,s)=C(x_{0},\tau _{0})\). Instead, one can use a reference solution \(C_{r}\) with drift \(f_{r}\) and diffusion function \(\Omega _{r}\), which approximates \(C\). Then, the (suboptimal) importance density is
and the estimator (15) takes the form
\(X_{n}\sim p_{2}\). For practical calculations, one can use the ratio of Euler densities \(q_{1}/q_{r}\). One thus obtains a Monte Carlo correction to the reference solution \(C_{r}\). However, the measures \(q_{1}, q_{r}\) may be singular w.r.t. each other, since \(\Omega _{1}\ne \Omega _{r}\) in general.
Using an analogous argument as is Sect. 5, but with unequal diffusion matrices, it is possible to draw independent samples from \(p_{2}\) (47) using the drift correction (scalar notation)
and setting \(f_{2}=f_{r}+\delta f\). Equation (49) can be found by differentiating (47) w.r.t. time \(s\) and inserting the backward equations for \(p_{r}\), \(C_{r}\) and \(p_{2}\).
If \(\Omega _{r}=\Omega _{2}=\Omega \) (absolutely continuous case) and \(f_{r}=f\), one recovers the simple drift correction (19).
Rights and permissions
About this article
Cite this article
Singer, H. Importance sampling for Kolmogorov backward equations. AStA Adv Stat Anal 98, 345–369 (2014). https://doi.org/10.1007/s10182-013-0223-z
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10182-013-0223-z