Tutorials › Stochastic Differential Equations in Pharmacokinetics › Lamperti Coordinates and Boundary-Aware Simulation

Positivity-Preserving Models and Simulation · Tutorial 50 of 100

Lamperti Coordinates and Boundary-Aware Simulation

You will learn to transform a positive state-dependent diffusion to unit diffusion, then use the transformed domain and boundary behavior to assess whether a simulation step is valid.

Advanced 10 min read

What You'll Learn

  • Derive the scalar Lamperti transform and its transformed drift using Itô’s formula.
  • Identify how zero and other state-space boundaries map into transformed coordinates.
  • Relate scale functions and boundary-hitting probabilities before and after transformation.
  • Work through geometric Brownian, square-root, and CEV diffusion examples.
  • Diagnose why an unconstrained Gaussian step or arbitrary clipping can misrepresent boundary behavior.

Why Change Coordinates?

For a positive pharmacokinetic state, random fluctuations may grow or shrink with the current amount or concentration. A scalar diffusion model on \(X_t>0\) can be written \(dX_t=b(X_t)\,dt+\sigma(X_t)\,dW_t\), where \(W_t\) is standard Brownian motion. When \(\sigma\) depends on \(X_t\), a fixed-size increment of Brownian noise produces different-sized state changes at different concentrations. The Lamperti transform changes coordinates so that the transformed diffusion coefficient is constant. This can simplify analysis and simulation, but it does not remove the original model’s boundary behavior.

The key is to keep track of the transformed state space. A boundary at zero might map to a finite endpoint, to minus infinity, or to a coordinate where the transformed drift is singular. The transformed equation is valid in the interior; it does not, by itself, specify what to do when a numerical step reaches or crosses an endpoint.

Definition and Transformed Dynamics

Definition: Let \(X_t\) take values in an open interval \(I=(r_-,r_+)\), and suppose \(\sigma(x)>0\) on \(I\). Choose a reference point \(x_*\in I\) and define \(h(x)=\int_{x_*}^{x} \frac{1}{\sigma(z)}\,dz\). The Lamperti coordinate is \(Y_t=h(X_t)\). The function \(h\) is strictly increasing, so it maps \(I\) onto an interval \(J=(\ell_-,\ell_+)\), whose endpoints are the one-sided limits of \(h\) at the endpoints of \(I\).

Assume \(b\) and \(\sigma\) are sufficiently smooth on \(I\), with \(\sigma>0\), and that the SDE has a continuous solution up to its lifetime or first approach to the boundary. Applying Itô’s formula in the interior gives

$$ dY_t=\left(\frac{b(X_t)}{\sigma(X_t)}-\frac{1}{2}\sigma'(X_t)\right)dt+dW_t. $$

Since \(X_t=h^{-1}(Y_t)\), the transformed drift as a function of \(y\in J\) is \(\mu_Y(y)=b(x)/\sigma(x)-\sigma'(x)/2\), evaluated at \(x=h^{-1}(y)\). In particular, the transformed diffusion coefficient is one, but the drift can become singular near an endpoint.

Theorem: Under the conditions above, up to the lifetime or boundary approach of the solution, \(Y_t=h(X_t)\) satisfies the unit-diffusion SDE with drift \(\mu_Y\) displayed above. If \(h\) maps an endpoint \(r_-\) to \(\ell_-\), then continuous paths approaching \(r_-\) correspond under \(h\) to paths approaching \(\ell_-\). The analogous statement holds at the upper endpoint.

For the dynamics, \(h'(x)=1/\sigma(x)\) and \(h''(x)=-\sigma'(x)/\sigma(x)^2\). Itô’s formula therefore gives

$$ \begin{aligned} d h(X_t) &=h'(X_t)b(X_t)\,dt+h'(X_t)\sigma(X_t)\,dW_t +\frac{1}{2}h''(X_t)\sigma(X_t)^2\,dt\\ &=\left(\frac{b(X_t)}{\sigma(X_t)}-\frac{1}{2}\sigma'(X_t)\right)dt+dW_t. \end{aligned} $$

For the endpoint statement, \(h\) is continuous and strictly increasing on \(I\), and its limit at \(r_-\) is \(\ell_-\). Thus whenever a continuous path \(X_t\) approaches \(r_-\) from within \(I\), its image \(h(X_t)\) approaches \(\ell_-\); the converse follows from continuity of \(h^{-1}\) on \(J\) and its endpoint limit. This correspondence concerns paths in the continuous model; a discrete numerical step can cross an endpoint without the continuous process ever doing so. \(\square\)

Worked Coordinate Changes

Worked Example: Geometric Brownian Motion

Let \(dX_t=aX_t\,dt+\gamma X_t\,dW_t\), with \(X_t>0\), \(a=0.4\), and \(\gamma=0.3\). Here \(\sigma(x)=\gamma x\), so choose \(h(x)=\log(x)/\gamma\). The transformed coordinate is \(Y_t=\log(X_t)/0.3\), and its drift is

\(\frac{b(x)}{\sigma(x)}-\frac12\sigma'(x)=\frac{0.4x}{0.3x}-\frac12(0.3)=1.3333-0.15=1.1833\), to four decimal places. Thus \(dY_t=1.1833\,dt+dW_t\). If \(X_0=4\), then \(Y_0=\log(4)/0.3\approx4.6210\). The inverse map is \(X_t=\exp(0.3Y_t)\), which is positive for every finite \(Y_t\). In this example zero maps to \(-\infty\), not to a finite endpoint that a finite Gaussian increment can cross.

As a check, transforming back gives \(\log(X_t)=\log(4)+(a-\gamma^2/2)t+\gamma W_t\). Since \(a-\gamma^2/2=0.4-0.045=0.355\), the usual geometric Brownian solution agrees with the transformed equation.

Worked Example: A Square-Root Diffusion

Consider \(dX_t=k(\theta-X_t)\,dt+\xi\sqrt{X_t}\,dW_t\) with \(k=0.8\), \(\theta=5\), and \(\xi=0.6\). For \(x>0\), \(h(x)=2\sqrt{x}/\xi\), and \(h'(x)=1/(\xi\sqrt{x})\). The transformed drift is

\(\frac{0.8(5-x)}{0.6\sqrt{x}}-\frac12\frac{0.6}{2\sqrt{x}}=\frac{0.8(5-x)}{0.6\sqrt{x}}-\frac{0.15}{\sqrt{x}}\).

Writing \(y=2\sqrt{x}/0.6\), so \(x=(0.3y)^2\), gives \(\mu_Y(y)=\frac{2k\theta/\xi^2-1/2}{y}-\frac{k}{2}y\). Here \(2k\theta/\xi^2=8/0.36\approx22.2222\), so \(dY_t=(21.7222/Y_t-0.4Y_t)\,dt+dW_t\), for \(Y_t>0\).

At \(X_t=4\), the transformed value is \(y=2\sqrt{4}/0.6\approx6.6667\). The displayed drift there is \(21.7222/6.6667-0.4(6.6667)\approx0.5917\). Direct substitution in the untransformed formula gives \(0.8(5-4)/(0.6\sqrt4)-0.15/\sqrt4\approx0.5917\), confirming the calculation. Zero maps to the finite endpoint \(y=0\). For these parameters \(2k\theta=8>\xi^2=0.36\), so the standard Feller condition makes zero unattainable from a positive initial state. That fact does not license a numerical step to use a negative \(y\): the transformed drift is defined only for \(y>0\) and is singular at zero.

Worked Example: A CEV-Type Diffusion and a Finite Boundary

Take \(dX_t=(2-X_t)\,dt+0.5X_t^{3/4}\,dW_t\) on \(x>0\). Then \(h(x)=\int_1^x 2z^{-3/4}\,dz=8(x^{1/4}-1)\). The image of zero is the finite endpoint \(-8\), and \(x=(1+y/8)^4\) for \(y>-8\).

The transformed drift, first written in \(x\), is \(\frac{2-x}{0.5x^{3/4}}-\frac12(0.5)(3/4)x^{-1/4}=2(2-x)x^{-3/4}-0.1875x^{-1/4}\). At \(x=1\), \(y=0\), and the drift is \(2-0.1875=1.8125\). This is an interior calculation: the transformed equation has unit diffusion only while \(y>-8\). A Gaussian increment that puts \(y\) below \(-8\) has left the transformed state space. Applying the inverse fourth-power expression there would not repair the step; it would apply the inverse map outside its domain.

Scale and Boundary Hitting

The coordinate change also preserves one-dimensional hitting probabilities when the same physical boundaries are represented in the new coordinate. To state this precisely, suppose \(b,\sigma\) are continuous on an interval, \(\sigma>0\), and the diffusion is regular there. A scale density for \(X\) is \(s_X'(x)=\exp(-\int^x 2b(z)/\sigma(z)^2\,dz)\); its integral \(s_X\) is a scale function, defined up to a positive affine change. For endpoints \(a<x<c\) in the interval, the probability of hitting \(c\) before \(a\) is \((s_X(x)-s_X(a))/(s_X(c)-s_X(a))\).

Result: Under a Lamperti coordinate change, a scale function for the transformed diffusion can be chosen as \(s_Y(h(x))=s_X(x)\). Consequently, the probability of hitting one of two interior endpoints before the other is unchanged by the transformation.

To verify the result, the transformed drift is \(\mu_Y=b/\sigma-\sigma'/2\), and \(dy=dx/\sigma(x)\). A scale density for a unit-diffusion process is \(\exp(-\int 2\mu_Y(y)\,dy)\). Changing variables to \(x\), its exponent has differential

$$ 2\mu_Y(y)\,dy =\left(\frac{2b(x)}{\sigma(x)^2}-\frac{\sigma'(x)}{\sigma(x)}\right)dx. $$

Therefore the transformed scale density, expressed as a function of \(x\), is \(\sigma(x)\exp(-\int^x 2b(z)/\sigma(z)^2\,dz)\). Multiplying by \(dy=dx/\sigma(x)\) cancels the factor \(\sigma(x)\), so \(ds_Y=ds_X\) under the coordinate map. Integrating gives \(s_Y(h(x))=s_X(x)\), up to the permitted affine choice. Substitution into the hitting-probability formula proves the claim. \(\square\)

This result is useful because a coordinate change cannot turn an attainable physical boundary into a harmless modeling choice. It changes the description of the boundary, not the underlying continuous-path probability of reaching it. For boundaries at the edge of the state space, classification requires the relevant scale and speed integrals and the model’s boundary assumptions; merely seeing a finite or infinite transformed endpoint is not a classification.

Boundary-Aware Simulation: Diagnostics and Limits

A transformed SDE with unit diffusion may look like an unconstrained process with additive Gaussian noise. That description is incomplete whenever its state space is a proper interval \(J\). A numerical increment can leave \(J\), even when the continuous process does not reach the boundary. If the boundary is attainable, the continuous model’s specified behavior there matters; if it is unattainable, a crossing may indicate a coarse or unsuitable numerical step. In neither case does the Lamperti transform itself prescribe reflection, absorption, or clipping.

Before simulating, check that the transform is one-to-one on the intended state space, calculate both endpoint limits, and identify where the transformed drift is defined. Then inspect simulated paths for boundary crossings, unusually large transformed drifts, and sensitivity to reducing the step size. A run with no observed crossing is not proof that the underlying model or numerical method respects the boundary: rare crossings can be missed in a finite simulation.

ApproachWhat it offersBoundary caution
Lamperti coordinatesUnit diffusion in a scalar transformed SDE; useful for exposing drift and endpoint geometry.The transformed interval may be bounded, and the drift may be singular near its ends.
Direct state-coordinate simulationKeeps the model expressed in the original physical state.A numerical update can still leave the positive state space; positivity is not guaranteed by the coordinate choice.
Boundary-specific numerical methodsCan incorporate a model’s known boundary behavior.Must match the actual boundary classification and model assumptions; an arbitrary constraint is not a substitute.

The scalar transformation assumes a nonvanishing diffusion coefficient on the interior. It does not automatically provide a single coordinate transformation for a multidimensional system with correlated process noise. Nor does it supply a discrete-time transition density or guarantee that an approximate simulation has the correct boundary behavior. Those are separate modeling and numerical questions.

Common Pitfall

Common mistake: “The transformed diffusion coefficient is one, so I can take a Gaussian step and clip any value outside the allowed interval back to the boundary.” The unit diffusion describes the continuous stochastic differential equation only in its transformed domain. Clipping changes the simulated transition law and may create artificial boundary mass or an unintended boundary behavior. Instead, identify the original model’s boundary behavior, retain the correct transformed domain, and use a simulation method justified for that model; do not infer a boundary rule from the coordinate transform alone.

Check Your Understanding

Try each problem, then check the brief solutions.

  1. For \(dX_t=\alpha X_t\,dt+\gamma X_t\,dW_t\), with \(X_t>0\) and \(\gamma>0\), find a Lamperti coordinate and transformed drift.
  2. For \(dX_t=b(X_t)\,dt+\eta\,dW_t\), with constant \(\eta>0\), what does the Lamperti transformation do to the state coordinate and diffusion coefficient?
  3. If \(h(x)\) approaches a finite value as \(x\) approaches zero, does that fact alone show that the original diffusion hits zero?
  4. Why is clipping a simulated transformed value to a finite endpoint not automatically valid?

Solutions. (1) Choose \(h(x)=\log(x)/\gamma\). Its transformed drift is \(\alpha/\gamma-\gamma/2\), so \(dY_t=(\alpha/\gamma-\gamma/2)\,dt+dW_t\). (2) One may take \(h(x)=x/\eta\), up to an additive constant; the transformed diffusion coefficient is one and the drift is \(b(\eta y)/\eta\). (3) No. A finite transformed endpoint describes the coordinate map; boundary attainability depends on the diffusion’s dynamics, not endpoint finiteness alone. (4) Clipping changes the simulated law and imposes an implicit boundary behavior that may not match the continuous model.

Clinical Biostats