{'fraction_non_alphanumeric': 0.057891267589487864, 'fraction_numerical': 0.01584845097870633, 'mean_word_length': 4.049782923299566, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 16, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal exercise of American options under time-dependent Ornstein-Uhlenbeck processes

We study the barrier that gives the optimal time to exercise an American option written on a time-dependent Ornstein-Uhlenbeck process, a diffusion often adopted by practitioners to model commodity prices and interest rates. By framing the optimal exercise of the American option as a problem of optimal stopping and relying on probabilistic arguments, we provide a non-linear Volterra-type integral equation characterizing the exercise boundary, develop a novel comparison argument to derive upper and lower bounds for such a boundary, and prove its differentiability and Lipschitz continuity in any closed interval that excludes the expiration date. We implement a Picard iteration algorithm to solve the Volterra integral equation and show illustrative examples that shed light on the boundary's dependence on the process's drift and volatility.


Introduction
The extra flexibility of being able to exercise at any time before maturity makes American options a popular financial instrument among traders, but it makes this type of vanilla options significantly harder to price when compared to their European counterparts. Indeed, while the arbitrage-free price of the European option written on a geometric Brownian motion was settled down in the notorious works of Merton (1973) and Black and Scholes (1973), the valuation of the American option had to go through a much longer process that took 40 years to complete, starting from the early works of Samuelson (1965) and McKean (1965), reaching important milestones in the works of Kim (1990), Jacka (1991), and Carr et al. (1992), and concluding in Peskir (2005). This pricing journey is registered in the surveys of Myneni (1992) and Barone-Adesi (2005).
The classic geometric Brownian motion is not a suitable model in many financial contexts. Examples are the so-called mean-reverting strategies, which arise when trading assets that tend to revert to an average price. Pair-trading strategies are the most popular of these scenarios. They take place when the spread between correlated assets is pulled towards a baseline level. In such a case, it is common to model the spread by means of an Ornstein-Uhlenbeck (OU) process to capture the mean-reverting behavior. See Leung and Li (2015) for the theory and empirical evidence, Gatev et al. (2006) for a practitioner perspective, and the handbook of Ehrman (2006) for a quick reference guide. Mean-reverting strategies also find room when trading options. Indeed, imagine that a trader holds two American options, one put and one call, written on two different stocks whose price at time t is X (1) t and X (2) t , with the same maturity date T and strike prices given by A 1 and A 2 . At any time t before maturity, the trader can exercise both options, hence making the profit X (1) . The trader can then carve out a new American put option written on the spread X t = X (1) t − X (2) t and with strike price given by A = A 2 − A 1 , and the spread can be modeled as an OU process. that holds an explicit and relatively simple dependency on the OU parameters, which allows deriving a clean intuition on how those affect the OSB shape. We comment on two algorithms and implement one of them to work out the solution of the integral equation. Our solution method, which is primarily probabilistic, yields important properties about the OSB, such as its differentiability and Lipschitz continuity. This order of smoothness of the OSB, which goes beyond its typical differentiability almost everywhere, is usually difficult to prove, especially for non-convex OSBs and non-differentiable gain functions. We also derive explicit upper and lower bounds for the OSB by relying on a comparison method that, to the best of our knowledge, has not been used within the framework of optimal stopping, and could potentially be used in generic processes with drifts that are monotone in the space component. The value function is characterized by means of the earlyexercise premium representation, which has a clear economic interpretation. As a matter of showing that the seemingly restrictive framework of continuously differentiable drifts and volatilities could be escaped by considering limit cases, we numerically explore the applicability of our work to approximate the OSB of diffusion bridges with drifts that are linear in the space component, like the Brownian bridge and the OU bridge. We also cover the American call option scenario by proving that its OSB is the reflection of the OSB of its put counterpart with respect to the strike price.
Recently, by using a different approach, Carr and Itkin (2021) also tackled the problem of pricing an American option with an OU process with differentiable time-dependent coefficients underneath. They performed time-space transformations to cast the partial differential operator of the associated free-boundary problem into that of the heat equation, which they solved by means of the heat potential method (the recent paper of Lipton and Kaushansky (2020) explains this method within the context of the American put option for a geometric Brownian motion). As a result, the OSB is proved to satisfy an integral equation of the Fredholm type, although the uniqueness of the solution was not addressed, nor numerical methods for finding approximate solutions were discussed. In addition, the dependence of the integral equation with respect to the process' coefficients does not seem to be evident, as neither is the interpretation of the pricing formula in financial terms. Remarkably, Carr and Itkin (2021) also allow for a time-dependent discount rate and repulsion drifts, which we do not account for in our setting.
The rest of the paper is structured as follows. Section 2 sets the problem and introduces basic notation. Section 3 gathers all the required properties to come up, in Section 4, with the pricing formula and free-boundary equation. Numerical experiments are displayed in Section 5, while concluding remarks are relegated to Section 6.

Formulation of the problem
Let X = {X s } s∈[0,T ] be a stochastic process satisfying the stochastic differential equation , where θ, α : [0, T ] → R are continuously differentiable functions and θ(s) > 0 for all s ∈ [0, T ]. We refer to θ and α as the slope and pulling level (functions), respectively. From (1), we get that the infinitesimal generator L of the process {(s, X s )} s∈[0,T ] is given by Above and henceforth, we use ∂ t and ∂ x to denote, respectively, the usual partial derivatives with respect to time and space of a function with arguments (t, x), and ∂ xx is used as a shorthand for Consider the finite-horizon, discounted Optimal Stopping Problem (OSP) where V is the value function, is the gain function for some A ∈ R that represents the strike price of the option, and λ ≥ 0 is the discounting rate. The supremum above is taken under all random times τ such that t + τ is a stopping time in {F s } s∈[0,T ] and E t,x represents the expectation under the probability measure P t,x defined as P t,x (·) := P (· | X t = x). For simplicity and in an abuse of notation, in the sequel we will refer to τ as a stopping time while keeping in mind that t + τ is the actual stopping time in the filtration {F s } s∈[0,T ] . The particular form of G implies that for all x, y ∈ R, inequalities which will be recurrently used throughout the paper. It is useful to keep track of the condition X t = x in a way that does not change the probability measure whenever t or x change. To do so, we denote the process X t, where the underlying Brownian motion {W s } s∈R + remains the same for all t ∈ [0, T ]. Theorem 38 from Chapter V in Protter (2003) guarantees that X t,x is well defined and is continuous with respect to x. Theorem 39 from Chapter V in Protter (2003) ensures that X t,x is continuously differentiable with respect to x and t, and combined with (1), it states that the processes ∂ t X t, , defined as the following P-a.s. limits and Mind the difference between the differential operators ∂ t and ∂ x and the processes ∂ t X t,x s and ∂ x X t,x s . From (7), it follows that the process ∂ x X t,x s is deterministic, i.e.
and that ∂ x X t,x s ∈ (0, 1] for all s ∈ [0, T − t]. Denote by D := {V = G} and C := D c = {V > G} the so-called stopping set and continuation set respectively. If for all (t, x) ∈ [0, T )×R, then we can guarantee that, under P t,x , the first hitting time of {X t+s } s∈[0,T −t] into D, denoted by τ * = τ * (t, x), is optimal in (3) (see, e.g., Karatzas and Shreve (1998, Appendix D) and Peskir and Shiryaev (2006, Section 2.2)), meaning that If there is another Optimal Stopping Time (OST) τ , then τ * ≤ τ P t,x -a.s. Solving the OSP (3) means providing tractable expressions for both the value function V and the OST τ * . The boundary of D (or C), denoted by ∂D (or ∂C), is called the Optimal Stopping Boundary (OSB). It turns out that both the value function V and the OSB ∂D are the solution of a Stefan problem with the infinitesimal generator L being the differential operator acting on V in C. For this reason, the OSB is also referred to as the free boundary (see Peskir and Shiryaev (2006) for more on the relation between OSPs and free-boundary problems). If the OSB can be depicted by the graph of a function b : Finally, it is convenient to recall the martingale and supermartingale properties of V , for all 0 ≤ s < T − t, as they will be used to prove results in Section 3.

Regularities of the boundary and the value function
In this section we derive regularity conditions about the OSB and the value function. These conditions allow obtaining a solution for the OSP (3), later addressed in Section 4, by using an extension of the Itô formula to derive a characterization of the OSB via a Volterra integral equation. We first point out that (9) holds true, as it is a consequence of Lemma 1 below. This allows us to prove the existence of the OST as well as its characterization in terms of the stopping set claimed in Section 2. We then shed light on the geometry of the stopping and continuation regions in the next proposition, which entails that the stopping set lies below the continuation set and that the boundary between them can be seen as the graph of a bounded function. We highlight the comparison argument used to derive the lower bound of the OSB. Besides its simplicity, it relies only on the monotonicity of the drift with respect to the space component, which makes it usable beyond the specifics of the OU process.
We prove now that x 1 ) ∈ D, then, for any x 2 < x 1 we have that, for τ * = τ * (t, x 2 ), where in the second inequality we used (5) and the fact that s ≥ 0 for all s according to (8)). We also used the mean value theorem alongside (8) to come up with the equality. We have proved that (t, x 2 ) belongs to D. Finally, being D closed, (t, b(t)) ∈ D for all t ∈ [0, T ), which guarantees that D has the claimed shape.
We use now the Itô-Tanaka formula to get that, for all (t, , and l A s (X) is the local time of the process X at A and up to time s, that is, We have already proved that b(t) ≤ γ(t) whenever γ(t) ≥ A, so we will assume that γ(t) < A. Fix (t, x) such that γ(t) < x < A and notice that, in such a case, Take a ball B ⊂ [0, T ]×(−∞, A) centered at (t, x) and sufficiently small such that F t (u, y) > 0 for all (u, y) ∈ B. In (14), replacing s by the first exit time from B (denoted by τ B c ), taking P t,x -expectation (cancels the martingale term), and noticing that X t+u never touches A for u ≤ τ B c under P t,x (i.e., the local-time term is null), we get Analogous arguments yield that (t, x) ∈ D for all x < min {A, γ(T )} and t sufficiently close to T . However, we are not including a proof of this assertion, but rather assuming that it holds and refer to Peskir and Uys (2005) for further details. Hence, we conclude that b(T ) := lim t→T b(t) = min{γ(T ), A}.
The lower boundedness of b is now addressed. Define m : where B x t = x + W t and {W t } t≥0 is the underlying standard Brownian motion in (1). The second inequality holds since the drift of the reflection of X t,x with respect to α is always higher than the drift of the reflection of B x with respect to α for all (t, x) ∈ [0, T ) × R, and therefore we can assume that the first process is greater than the last one pathwise P-a.s. (see Peng and Zhu, 2006, Corollary 3.1). The above inequalities, and the fact that G(x) is a decreasing function, guarantee that the value function (3) is lower than the value function associated with a reflected (with respect to m) Brownian motion, and hence the respective OSBs hold the reverse inequality, that is, the stopping region of the latter process is contained in the stopping region of the former. It is easy to show that the free-boundary for the discounted OSP with the gain function G and m-reflected Brownian motion is finite. Actually, one can obtain explicitly the (constant) OSB by directly solving the associated free-boundary problem. Therefore b is bounded from below.
The value function satisfies the regularity properties listed in the next proposition. The method used to get the Lipschitz continuity of V , based on properties (11) and (12), is a powerful technique in solving OSPs framed within Markovian processes. For example, the work of De Angelis and Stabile (2019) provides Lipschitz continuity of the value function in a similar fashion, but for smooth gain functions and time-homogeneous diffusions. For a neat summary of further developments on the regularity of the value function, see the work of De Angelis and Peskir (2020), which proves V ∈ C 1 ([0, T ] × R) in a fairly general setting, and references therein.
Proposition 2 (Regularity of V ). The value function V satisfies the following properties: where τ * = τ * (t, x). Additionally, for some positive constant L.
is optimal for the OSP with initial conditions (t, x) and suboptimal for any different initial condition, and due to inequality (5), it follows that, for δ > 0 small enough, where the last equality holds for a certain (random) time t δ ∈ (t − δ, t) due to the mean value theorem. In (18) and thereafter we use the notation By (39) below, the value of L x is bounded. Applying similar arguments as the ones used in (18), relying on the martingale and supermartingale properties of the value function, and noticing that where ρ * := τ * T − δ, X t,x T −t−δ and, in the second inequality, we added and subtracted the term . We now analyze the last term in (19) in detail. By the definition of V in (3), using the Itô-Tanaka formula, and acknowledging that where ). Then, arguing as in (19), it follows that with the first and second addend being lower bounded by δL x and −δ (L x + L), respectively. For 0 ≤ δ ≤ T − t and τ δ = τ * (t + δ, x) one gets the following by proceeding as in (18): So far, as L x is finite (see Lemma 1), we have proved that, for any x ∈ R, t → V (t, x) is Lipschitz continuous in [0, T ]. We will now prove that V is also Lipschitz continuous with respect to x ∈ R for all t ∈ [0, T ], which will complete the proof of (i).
Since G is decreasing and x), and argue as in (13) (1) is Lipschitz continuous (actually, it suffices to require local α-Hölder continuity) in [0, T ] × R, and the diffusion coefficient is constant, one can borrow a classic result from the theory of parabolic partial differential equations (see Friedman, 1964, Section 3, Theorem 9) to guarantee that, for an open rectangle R ⊂ C, the first initial-boundary value problem has a unique solution f ∈ C 1,2 (R). Fix (t, x) ∈ R and let τ R c be the first time X t,x s exits R. Therefore, we can use the Itô formula on f (X t,x s ) at s = τ R c , together with (23) and (24), to get the equality E t,x [V (X t+τ R c )] = f (t, x). Finally, notice that, due to the strong Markov property, , meaning that f = V (see Peskir and Shiryaev, 2006, Section 7.1, for further details). The fact that V is C 1,2 in D follows straightforwardly from V = G on D.
(iii) We already mentioned that x → V (t, x) is decreasing since G is decreasing as well while x → X t,x s is increasing. For (t, x) ∈ [0, T ) × R such that x = b(t), (15) follows after recalling that V is differentiable with respect to x in C and D, dividing by δ in (22), and taking δ → 0, while using the dominated convergence theorem. The same procedure used in (18) and (21) yields (16) and (17).
Smoothness of the free boundary is essential to prove the smooth-fit condition, in which the uniqueness of the value function V that solves the OSP is based. The works of De Angelis (2015) (2015) provides the continuity of the free boundary for one-dimensional processes with locally Lipschitz-continuous drift and volatility. The two-dimensional case (including time-space diffusions) is addressed by Peskir (2019) in fairly general settings, which proves the impossibility of first-type discontinuities of the OSB for Mayer-Lagrange OSPs, yet does not address second-type discontinuities. Unfortunately, this is not enough to prove the smooth-fit condition (piecewise monotonicity is sufficient; see Proposition 4). De Angelis and Stabile (2019) go further in proving the Lipschitz continuity of the free boundary in a higher dimensional framework, also for Mayer-Lagrange OSPs, although only tackling timehomogeneous processes and imposing somehow restrictive conditions. For example, the gain function is assumed to be C 3 in the spatial coordinate, which is not satisfied by our gain function G.
The next proposition proves the differentiability and Lipschitz continuity of the OSB by adapting the methodology used by De Angelis and Stabile (2019) to deal with the non-differentiable gain function G and the time-inhomogeneous OU process. Proof. Consider the function W (t, x) = V (t, x) − G(x) and the closed interval I = [ t,t ] ⊂ [0, T ). Proposition 1 guarantees that b is bounded from above and hence we can choose r > sup {b(t) : t ∈ I}. Hence, I × {r} ⊂ C, and then W (t, r) > 0 for all t ∈ I. Since W is continuous in C (see Proposition 2), there exists a constant a > 0 such that W (t, r) ≥ a for all t ∈ I. Therefore, for all δ such that 0 < δ ≤ a, the equation W (t, x) = δ has a solution in C for all t ∈ I. This solution is unique for each t, as ∂ x W > 0 in C (see (15)). Hence we can denote it by b δ (t), where b δ : I → (b(t), r]. Away from the boundary, W is regular enough to apply the implicit function theorem, which guarantees that b δ is differentiable and Notice that b δ is increasing in δ and therefore it converges pointwise to some limit function b 0 , which satisfies b 0 ≥ b in I as b δ > b for all δ. Since W (t, b δ (t)) = δ and W is continuous, it follows that W (t, b 0 (t)) = 0 after taking δ ↓ 0, which means that b 0 ≤ b in I and hence b 0 = b in I. For (t, x) ∈ C such that t ∈ I and x < r, consider the stopping times τ * = τ * (t, x) and Using (38) from Lemma 1, alongside the fact that τ * ≥ T − t, it readily follows after (16) that where L is a positive constant coming from the Lipschitz continuity of ∂ t µ and m is the function defined as Due to the tower property of conditional expectation and the strong Markov property, we have that Notice that, for b(t) < x < r, t + τ r , X t,x t+τr ∈ Γ t := ((t,t) × {r}) ∪ ({t} × (b(t), r]) on the set {τ r < τ * }. Then, we get the following upper bound on {τ r < τ * }: where we used the Cauchy-Schwartz inequality and (36) from Lemma 1. Plugging (28) into (27), recalling (26), and noticing that, for u ≤ τ * ∧ τ r , |X t,x u | ≤ max{|b|, r} with b := inf{b(t) : t ∈ [0, T ]}, we get that for some positive constant K (1) I . Using (15) and (37) from Lemma 1, and proceeding in the same way as for (27), we get the following lower bound for ∂ x W (t, x): Now, by putting together equations (25), (29), and (30), we obtain the following inequalities for a constant K (2) where in the last inequality we used that X t,x δ τr = r in the set {τ r ≤ τ δ , τ r <t − t}. Since I × {r} ⊂ C, there exists ε > 0 such that R ε := [ t,t + ε] × (r − ε, r + ε) ⊂ C. Notice that τ * (t, r) > τ ε (t, r) for all t ∈ I, with τ ε (t, r) = inf u ≥ 0 : t + u, X t,r u / ∈ R ε . Hence, Therefore, going back to (31), we have that, for all t ∈ I = [ t,t − ε] and ε > 0 small enough, To find a bound in the opposite direction that is also uniform with respect to δ and for all t ∈ I , we consider (17) and (26) and use the Markov inequality to get for some positive constant L. Hence, relying on the same arguments as the ones used to get (32), but with (33) instead of (29), we obtain b δ (t) ≤ K 1 + ε −1 + K I,ε + L(T −t) −1 .
We have proved that |b δ (t)| is bounded by a constant, uniformly in δ and for all t ∈ I . Then the Arzelà-Ascoly's theorem gives that b δ converges uniformly to b in I as δ → 0, which implies that b is differentiable and Lipschitz continuous in I . Since ε > 0 and I can be chosen arbitrarily, we conclude that b is differentiable anywhere on [0, T ).
The next proposition gives the smooth-fit condition for our OSP. Its proof relies on the work of Cox and Peskir (2015), which implicitly shows that the condition holds for piecewise monotone boundaries for general recurrent diffusions.
Besides, reasoning as in (22), we get that The supremum is finite (actually, it is lower than 1) since ∂ x µ < 0. Hence, we can apply the dominated convergence theorem to obtain with τ 0 := lim δ→0 τ δ (τ 0 is well-defined since the sequence τ δ decreases with respect to δ P-a.s.).
The smooth-fit condition arises by recalling that ∂ − x V (t, b(t)) = −1 since V = G on D.
We conclude this section with a technical lemma on the boundedness of several functionals of the underlying process X that have been used throughout the previous proofs.
Lemma 1 (Some bounds on functionals of X). The following inequalities hold for positive constants L, M (i) Moreover, s → M (i) s is an increasing function for i = 1, 3, and it is decreasing for i = 2.
To obtain (37), we use that µ is Lipschitz continuous and ∂ x µ < 0, together with representation (8) and (15), which lead to where L is the Lipschitz constant. Let us prove now (39). To do so, notice first that (6) implies Therefore, an application of Gronwall's inequality yields Hence, due to the Lipschitz continuity of µ and ∂ t µ, and using (36) alongside the Cauchy-Schwartz inequality, we have that which entails (38) and (39).

The pricing formula and the free-boundary equation
Our main result is collected in the following theorem. It gives the characterization of the OSB as the unique solution, up to regularity conditions, of a type-two Volterra integral equation, and the early-exercise premium representation for the value function.
Theorem 1 (Free-boundary equation and pricing formula).
The OSB related to the OSP (3) satisfies the integral equation where, for c 1 , c 2 , x 1 , x 2 ∈ R and t 1 , t 2 ∈ [0, T ] such that t 2 ≥ t 1 , where φ and Φ are the density and distribution functions of a standard normal, and ν and γ depend on the time-dependent OU parameters: Moreover, (40) has a unique solution among the class of continuous functions of bounded variation c : [0, T ] → R that satisfy c(t) < A for all t ∈ [0, T ). Also, the associated value function is given by Proof. Propositions 1-4 allow to apply an extension of the Itô formula (see, e.g., D'Auria et al., 2020, Lemma A2) to e −λs V (t + s, X t+s ). After setting s = T − t, taking P t,x -expectation, and shifting the integrating variable t units, it follows that where the martingale term vanishes after taking P t,x -expectation, and the local time term does not appear due to the smooth-fit condition. Recall that LV = λV on C and V = G on D. Therefore, (45) turns into the pricing formula The first term of the right-hand sum in (46) stands for the price of the European put option written on the same asset and expiring on the same date, while the second term is the early exercise premium, i.e., the cost of being able to exercise the option before the expiration time T .
By taking x ↑ b(t) in (46) we get the free-boundary equation The uniqueness of the solution of (47), up to the conditions stated in the theorem, follows well-known arguments described by (Peskir, 2005, Theorem 3.1).
We can provide a more tractable expression for both the pricing formula (46) and the freeboundary equation (47) by taking advantage of the Gaussianity of X. Indeed, it is easy to get (see, e.g., Gardiner, 2004, Section 4.4.9) that, under P t,x , X u has a normal distribution with mean and variance given by (42) and (43) for all u ∈ [t, T ]. This leads, after some algebra, to equations (40) and (44).
Further flexibility of the underlying process can be achieved by allowing a time-dependent volatility, as it is the case in the popular Hull-White model. We show in the next remark that such an extension requires no extra analysis: rather, it can be reduced to the constant-volatility model (1) by means of a change of time.
Remark 1 (Time-dependent volatility). Let X t (θ, α, σ, T ) = {X t } t∈[0,T ] be a stochastic process satisfying the stochastic differential equation where θ, α, σ : [0, T ] → R are continuously differentiable functions such that θ(t) > 0 and σ(t) > 0 for all t ∈ [0, T ], and T 0 σ 2 (t) dt < ∞. Define the change of timet =h(t) = t 0 σ 2 (u) du and the inverse transformation h =h −1 . Consider the time-changed process X = Xt t ∈[0, T ] defined by Xt = X h(t) (θ, α, σ) and T =h(T ). Then, where θ(t) = σ −2 (h(t))θ(h(t)) and α(t) = α(h(t)), that is, X is a time-dependent OU process with constant volatility. Hence, we can rely on equations (44) and (40) to compute the value function and OSB, V andb, associated to X, from where we can simply obtain the value function and OSB related to X(θ, α, σ, T ) as The optimal exercise of the put and call American option are fundamentally the same problem due to their symmetry with respect to the strike price, as shown in the next remark.
Remark 2 (Put-call parity). Let V p be the value function associated to an American put option written on the time-dependent OU processes X (p) = X (p) t t∈[0,T ] and with strike price A. Then, (p) , and V c is the value function of the American call option with strike price A and written on X (c) . Similarly, the optimal exercise strategy for the call option is to stop the process X (c) the first time it is below the function b c :

Numerical experiments
Non-linear Volterra-type integral equations are hard to solve, and the free-boundary equation (40) is not an exception. Mainly, there are two methods to solve this type of integral equations arising from OSPs, both relying on the contraction principle: (1) a backward induction approach that recursively builds the boundary by leveraging that b(t) depends only on {b(s)} s∈ [t,T ] ; and (2) the well-known method of Picard iterations that, given an initial candidate boundary b 0 , consequently computes subsequent boundaries until a convergence criterion is fulfilled. See Detemple (2005, Chapter 8) for more details on the method of backward induction and Pedersen and Peskir (2002) for implementations. The Picard iteration scheme has been used by, e.g., Detemple and Kitapbayev (2020) and De Angelis and Milazzo (2020). As far as we know, the convergence of both approaches has not been formally addressed, and authors tend to provide numerical evidence of the error decreasing as the number of iterations grows. We adopted the Picard method to display visual insight about the OSB's shape, as our numerical experiments suggested similar accuracy and faster computations compared to the backward induction approach (within the set of explored settings).
Take the partition 0 = t 0 < t 1 < · · · < t N = T for N ∈ N. We initialize the Picard scheme by starting with the constant boundary b (0) (t) = b(T ) for all t ∈ [0, T ]. Recall from Proposition 1 that b(T ) = min(A, γ(T ) = (θ(T ) + λ) −1 (θ(T )α(T ) + λA)) is known. Each iteration of the algorithm updates the boundary according to the following formula, which is a right Riemann sum version of the integral in (40) We stop the iterations of the fixed-point algorithm (49) at the first k = 1, 2, . . . such that the (squared) distance d k : 2 is less than δ = 10 −3 . Computational experiments suggested that a mesh that is denser at the terminal time T performs better and that it is preferable that the distances between consecutive nodes narrow smoothly. For that reason we used the logarithmically-spaced partition t i = log (1 + i(e − 1)/N ) with N = 200 (unless is otherwise specified).
For Figures 1-3, the images on top show computer drawings of the obtained OSB from the Picard iteration (49). The images below represent the sequences of errors of the algorithm, that is, the x-axis accounts for the iteration number k = 1, 2, . . . and the y-axis for d k . The decrease of these error sequences corroborates empirically the convergence of the algorithm.  We set T = 1, λ = 1, and A = 0. The smaller images below provide the errors log 10 (d k ) between consecutive boundaries for each iteration k = 1, 2, . . . of (49). Figure 1 provides insights into the shape of the OSB for different sets of slopes and pulling levels. We also allow the process to have non-constant volatility, as discussed in Remark 1. In each image, two of the three functions that define the dynamics in (48) are fixed, and the remaining one varies to make clear the marginal effect it has on the OSB's shape. We use functions Φ and φ to model abrupt changes in a smooth fashion. Figure 1(a) models a regime switch of the stock price, where the switching happens smoothly, deterministically, and depends only on the time variable. It also reflects the attracting behavior of the OSB towards α. Figure 1(b) exemplifies how this attraction strengthens as θ increases, and also depicts a sudden increase of the pulling strength, which could represent a sharp increase of confidence in a belief of the price evolution. Figure 1(c) introduces periods where the volatility spikes to unusual levels before returning to a baseline. It also shows how the OSB gets repelled from the pulling level as the volatility gets higher, which, in alignment with Remark 1, coincides with the effect of reducing the slope.
Recent interest in OSPs related to diffusion bridges whose drift is linear in the space component (see, e.g., D'Auria et al. (2020) and D'Auria et al. (2021)) led us to investigate the applicability of our OSB within this context. The difficulty of dealing with these processes is that their slopes θ explode as the time approaches the horizon, hence they do not fit our assumption of a bounded slope. However, we suggest in Figure 2 that this inconvenience could be circumvented in practice by arbitrarily approximating the explosion with non-exploding smooth drifts of time-dependent OU processes. We do so by relying on the OSB of a Brownian Bridge (BB) and an Ornstein-Uhlenbeck Bridge (OUB) in the non-discounted scenario, whose solutions are available for some particular settings.
The OSB for the BB with pinning point A and gain function x → (A − x) + takes the form A − Bσ √ T − t, for B ≈ 0.8399, a result that goes back to the work of Shepp (1969) and its connection to the American option by D' Auria et al. (2020). We use, in Figure 2(a), the Taylorexpansion approximation (at t = 0) of the BB slope (1 − t) −1 , namely, n i=0 t i , for different values of the order n indicated in the legend. The results in Figure 2(a) show that the application of the Picard iteration (49) on the time-dependent OU that order-n approximates a BB converges to the expected, known, OSB.  The OSB of the OUB is not known in closed-form, but rather as a characterization via the freeboundary equation provided in D' Auria et al. (2021), which can be also solved by means of a Picard algorithm. In that paper the authors treated the OUB case with the identity as the gain function. However, a simple translation with respect to A provides the OSB for x → (A − x) + as the gain function whenever the pinning point is lower than A. This is easy to note as both gain functions, x → A−x and x → (A−x) + coincide for all x ≤ A. Then, there is no fundamental difference in both OSPs (for λ = 0) whenever the boundary associated with the OUB lies below A. In Figure 2(b), to avoid the explosion of the OUB slope, t → a coth(a(1 − t)), we use instead the smooth function t → 1(t ≤ 1 − ε) a coth(a(1 − t)) + 1(t > 1 − ε) 1 − exp{−a 2 (t − 1 + ε)/ sinh 2 (aε)} + a coth(aε) , for a = 5, and different values of ε specified in the legend. Figure 2(b) shows that the application of (49) on the time-dependent OU that approximates an OUB with pinning point z = 2 converges to the (numerically-computed) OSB of the OUB.  Finally, Figure 3 illustrates how the algorithm's output seems to converge as the partition size increases. Observe the relatively good performance with just a few points in the partition. Additionally, we show how changing the discounting rate affects the OSB's shape: larger discounts decrease the separation of the OSB with respect to the strike price.
The repository https://github.com/aguazz/AmOpTDOU contains all the required R scripts to reproduce the numerical experiments.

Concluding remarks
If a time-dependent Ornstein-Uhlenbeck model fits well the price of the underlying asset, the strategy that maximizes the mean gain of a holder of an American put option on that asset is that of exercising once the price lies below the solution of the integral equation in Theorem 1. This integral equation can be solved numerically using the algorithm from Section 5. The fair price for the option is the one given at (44), which splits into the price of the European option written on the same asset plus an early-exercise premium.
The problem of the optimal exercise of the American option was treated in this paper as an optimal stopping problem. It is known that not having a monotonic optimal stopping boundary makes it hard to derive the smooth-fit condition, which is key to proving the uniqueness of the solution for the free-boundary equation. We overcame this handicap by proving its Lipschitz continuity instead. Also, it is worth highlighting the comparison method used to get the lower bound of the optimal boundary in Proposition 1, as it seems to be extensible to other settings.
The numerical experiments performed in Section 5 reveal a wide flexibility in the optimal boundary shape inherited from changing the form of the time-dependent parameters. Especially, Figure 2 suggests that the limiting differentiability assumption on the process' coefficients can be escaped to approximate the optimal boundary of a Brownian bridge and an Ornstein-Uhlenbeck bridge, and, potentially, other diffusion bridges whose drift is linear in the space component.