Optimal Trade Execution Under Endogenous Pressure to Liquidate: Theory and Numerical Solutions

# Optimal Trade Execution Under Endogenous Pressure to Liquidate: Theory and Numerical Solutions

## Abstract

We study optimal liquidation of a trading position (so-called block order or meta-order) in a market with a linear temporary price impact (Kyle, 1985). We endogenize the pressure to liquidate by introducing a downward drift in the unaffected asset price while simultaneously ruling out short sales. In this setting the liquidation time horizon becomes a stopping time determined endogenously, as part of the optimal strategy. We find that the optimal liquidation strategy is consistent with the square-root law which states that the average price impact per share is proportional to the square root of the size of the meta-order (Bershova and Rakhlin, 2013; Farmer et al., 2013; Donier et al., 2015; Tóth et al., 2016).

Mathematically, the Hamilton-Jacobi-Bellman equation of our optimization leads to a severely singular and numerically unstable ordinary differential equation initial value problem. We provide careful analysis of related singular mixed boundary value problems and devise a numerically stable computation strategy by re-introducing time dimension into an otherwise time-homogeneous task.

###### keywords:
optimal liquidation, price impact, square-root law, singular boundary value problem, stochastic optimal control
###### Msc:
[2010] 34A12, 49J15, 91G80
1\cortext

[ca]Corresponding author

## 1 Introduction

We study optimal liquidation of an infinitely divisible asset when the execution price is subject to adverse price impact in proportion to the amount of the asset sold per unit of time, in line with Kyle (1985). The optimal liquidation strategy trades off expediency against the adverse price impact caused by a precipitous sale. However, our focus on liquidation is not fundamental; mutatis mutandis one can replace optimal liquidation with optimal acquisition in what follows. The novelty in our approach is that we rule out short sales in a falling market. This seemingly small change has a profound impact on the economics and mathematics of the problem. How and why this happens is the subject of the ensuing analysis.

Modelling of optimal execution with market impact is relatively new in the literature, going back to Almgren and Chriss (2000), Bertsimas and Lo (1998) and Subramanian and Jarrow (2001). Classical models Almgren and Chriss (2000), Bertsimas and Lo (1998), envisage a world where the unaffected price of the asset is a martingale and hence there is no pressure to trade quickly for an agent with linear utility. In these circumstances the incentive to trade is given by fiat – it is assumed that there is a fixed time limit by which the entire position must be liquidated.

The literature finds that optimal liquidation gives rise to ‘implementation shortfall’ (Perold, 1988) defined as the gap between the initial market value of the inventory and the expected revenue of the liquidation strategy; the latter always being lower due to the price impact. The shortfall itself is formed of two components, one due to ‘permanent price impact’ and another caused by ‘temporary impact’. The former cannot be influenced by the trading strategy, while the latter determines the optimal strategy and can be made arbitrarily small by making the liquidation time horizon longer. In this sense having more time is unambiguously beneficial to the trader.

A second strand of literature, Brown et al. (2010), Chen et al. (2014), Chen et al. (2015), identifies the motive to liquidate with a change in market conditions whereby tighter margin requirements lead to lower permitted amount of leverage. The change in market conditions occurs at discrete time points, while the optimal liquidation (deleveraging) is implemented continuously in time. Here for reasons of tractability the unaffected price is assumed constant during liquidation, although one could in principle use the results from the first strand of literature to make the modelling of the deleveraging phase more realistic.

In this paper we focus on the liquidation phase. Specifically, we study a situation where the unaffected price may be falling on average, which is highly plausible in a market with contracting liquidity. One expects that with the asset price decreasing the implementation shortfall should be more severe than in the martingale case. Surprisingly, the current literature finds that far from exhibiting a shortfall the optimal liquidation strategy may in this case record an expected surplus, see Schied (2013). On closer inspection one observes that the surplus arises due to short sale of the asset with subsequent acquisition at deflated price near the end of the allotted time horizon.

While strategic short sales in a bear market are not entirely implausible we feel it is important to examine a situation where such short sales are ruled out. The simplest way to achieve this is to stop the trading once the entire position has been liquidated. In doing so we recover the classical outcome from the martingale case whereby the price impact invariably leads to implementation shortfall. However, in a falling market without short sales it is no longer true that the shortfall can be made arbitrarily small by extending the liquidation time horizon.

Introduction of a stopping time is a novel feature in the optimal liquidation literature with a perfectly divisible asset. Previously, optimal stopping has appeared only in the context of optimal liquidation of an indivisible asset, see Mamer (1986) and Henderson and Hobson (2013). Although stopping on liquidation automatically precludes short sales, it does leave open the possibility of further intermediate acquisition. Ex-post it turns out that intermediate acquisition is not optimal, see Proposition 7.1 and Theorem 7.2. We show that the presence of the stopping time dramatically changes mathematical properties of the Hamilton-Jacobi-Bellman equation and leads to a severely singular and numerically unstable initial value problem. Part of our research contribution is in providing a comprehensive theoretical and numerical analysis of this HJB equation and related singular boundary value problems.

The paper is organized as follows. In Section 2 we survey related literature and present our model. In Section 3 we discuss reduction of our HJB partial differential equation (PDE) to an ordinary differential equation (ODE). Section 4 offers a probabilistic and control-theoretic interpretation of this reduction. Section 5 describes the singularity of the initial value problem (IVP) for the ODE of Section 3, while Section 6 shows how to obtain uniqueness from a related boundary value problem (BVP). In Section 7 we characterize the optimal strategy and its value function by means of the BVP of Section 6. In Section 8 we introduce and theoretically analyze a related PDE BVP which leads to a stable numerical scheme and present numerical results. Section 9 concludes.

## 2 Our model and related literature

We take the point of view of a trader with inventory whose initial value is given. The modelling is based on the premise that there is some price process – often called the ‘unaffected price’ – with exogenously given dynamics that governs the evolution of the asset price in the absence of our trading. In our case the unaffected price is a geometric Brownian motion

 dS(t)=λS(t)dt+σS(t)dB(t), (2.1)

where is a Brownian motion in its natural filtration.

The inventory attracts interest rate which becomes a storage cost when . We assume that the inventory is sold off continuously at a (stochastic) rate so that represents the amount of inventory sold per unit of time. Consequently, the inventory dynamics read

 dZ(t)=(rZ(t)−v(t))dt. (2.2)

Let be the first time when the entire inventory is disposed of. For a given pair of initial values

 s=S(0),z=Z(0), (2.3)

the expected discounted revenue from the disposal of the asset is given by

 J(s,z,v)=Es,z[∫T(Z=0)0e−ρt(S(t)−ηv(t))v(t)dt], (2.4)

where is the ‘affected price’ of the asset. In our setting measures the strength of ‘temporary impact’ the selling speed has on the price. The discount factor captures the opportunity cost of not holding alternative assets. The entire model is based on Černý (1999).

The task is to find optimal liquidation strategy that maximizes

 V(s,z):=supv∈AJ(s,z,v). (2.5)

We say that is an admissible control, and write , if process is predictable,

 E[∫t0|v(s)|mds]<∞ for all t>0 and m=1,2,…, (2.6)

and

 E(∫T(Z=0)0e−ρs|v(t)(S(t)−ηv(t))|dt)<∞. (2.7)

The optimization in our model can be seen, for specific parameter choices, as a special case of Ankirchner and Kruse (2013), Forsyth et al. (2012) and Schied (2013), with the crucial difference that in our case the liquidation time horizon is endogenous. We make a standing assumption that the time discounting is stronger than the expected appreciation and the interest on the asset combined,

 ρ>λ+r. (2.8)

To conclude this section we wish to make several observations that justify the choice of our modelling framework. The extant literature contains a number of variations on the model presented above. The trading may be discrete, rather than continuous, the unaffected price may be specified differently and the optimization criterion may involve a utility function. In common, existing models assume is fixed and exogenously given.

The most commonly considered specification for the affected price reads

 ~S:=S−γ(Z(0)−Z−−12ΔZ)−η1v+η2ΔZ, (2.9)

where is the ‘permanent’ price impact2 while and respectively, are known as ‘temporary’ price impacts in the continuous-time and discrete-time literature, respectively. It is assumed either that there is a finite number of fixed dates where is allowed to jump (discrete-time models) or that changes continuously at a stochastic time rate (continuous-time models). In each case is taken to be a predictable semimartingale with left limit process and jumps . Models in this category include Ankirchner et al. (2016), Brown et al. (2010), Chen et al. (2014), Gatheral and Schied (2011), Schied (2013), Schied and Schöneborn (2009) Ting et al. (2007) in continuous time and Almgren and Chriss (2000), Bertsimas and Lo (1998) in discrete time. Other impact specifications can be found, for example, in Chen et al. (2015), Cheridito and Sepin (2014), Forsyth (2011), Lorenz and Almgren (2011), Subramanian and Jarrow (2001) and Ting et al. (2007).

The revenue from liquidation over a fixed time horizon is given by . When the unaffected asset price process is a martingale, integration by parts together with suitable boundedness of and boundary condition yields

 E[R(T)]=Z(0)S(0)−γ2Z(0)2−η1∫T0v2(t)dt−η2N∑i=1(ΔZ(ti)2.

This equality offers several important insights:

1. Permanent impact (as defined here) has no strategic influence and in the absence of temporary impact () any strategy is optimal. The expected implementation shortfall is strictly positive.

2. With temporary impact it is optimal to liquidate at a constant rate, regardless of the strength of the permanent impact. The additional implementation shortfall equals in continuous time and in discrete time, respectively.

These observations suggest that temporary impact is responsible for the majority of strategic interaction also in the drifting market and we conjecture that the optimal strategy will therefore not change dramatically when the permanent impact is included. This is not to say that the implementation shortfall would be unaffected by the presence of permanent impact. Given the complexity of analysis to follow and the likely marginal gains to our understanding from the presence of permanent impact on the optimal trading strategy we feel justified in leaving out the permanent impact from our analysis.

More recent studies, excellently summarized in Gatheral (2010), consider an intermediate form of impact where the execution price is given by the formula

 St−∫t0f(vu)G(t−u)du.

Kernel is called the resiliency of the market and the two extreme cases, permanent impact and temporary impact, correspond to being constant or being the Dirac delta function, respectively. In Gatheral (2010) a case is made for a combination of power impact function, , with power law resiliency , , the latter tending to a Dirac delta function as . We note that our setup corresponds to the limiting case and we leave the analysis of the general impact function with general resiliency in the setup of this paper to future research.

## 3 HJB equation and dimension reduction

The value function defined in (2.5) formally solves the Hamilton-Jacobi-Bellman partial differential equation

 supv{12s2σ2Vss+λsVs+(rz−v)Vz−ρV+v(s−ηv)}=0, s>0,z>0,

with formal optimal control

 v∗=s−Vz2η,

giving rise to a quasilinear second order PDE

 12s2σ2Vss+λsVs+rzVz−ρV+(s−Vz)24η=0, (3.1)

with an initial condition

 V(s,0)=0. (3.2)

The self-similarity

 V(s,z)=s2u(x)/(ησ2), x=ησ2z/s, (3.3)

reduces (3.1, 3.2) to an initial value problem (IVP) for an ordinary differential equation (ODE)

 x2u′′ =axu′+bu−(u′−1)2/2, x>0, (3.4) u(0) =0, (3.5)

where

 a:=2(λ−r+σ2)/σ2, b:=−2(2λ−ρ+σ2)/σ2. (3.6)

The self-similarity reduces a problem with 7 independent parameters and to a problem with just three parameters: and .

## 4 Probabilistic interpretation of self-similarity

We begin by restating the HJB equation (3.1) in its variational form,

 supv(0){drift0(e−ρtV(S,Z))+v(0)(S(0)−ηv(0))}=0.

Plug in the self-similarity form of the value function and rearrange to obtain

 supv(0){drift0(e−ρtS2S(0)2u(ησ2ZS))+ησ2v(0)S(0)(1−ηv(0)S(0))}=0.

The next steps involve i) changing measure to given by where and are restrictions of and to ; ii) defining a new state variable ; and iii) reparametrizing the control to which yields

 supg(0){ˆdrift0(e(2λ+σ2−ρ)tu(X))+σ2g(0)(1−g(0))}=0. (4.1)

 dX=(rX−σ2g)dt+X(−dSS+d[S,S]S2),

while from the Girsanov theorem we obtain , which implies

 dX=((r−λ−σ2)X−σ2g)dt+σ2Xd^B, (4.2)

where is a Brownian motion under and denotes the stochastic logarithm of . In the final step we perform a time change from to defining and . This yields the dynamics

 d^X=(r−λ−σ2σ2^X−g)dt+^Xd^W, (4.3)

while (4.1) changes to

 supg(0){ˆdrift0(exp(2λ+σ2−ρσ2t)u(^X))+g(0)(1−g(0))}=0. (4.4)

With (4.3) in hand the optimality condition (4.4) explicitly reads

 0 = 12x2u′′(x)+r−λ−σ2σ2xu′(x) Missing or unrecognized delimiter for \left

and the (formal) optimal control equals . It is furthermore clear that (4.4) itself is a HJB equation of an optimal control problem

 u(x)=supgˆEx=^X(0)[∫T(^X=0)0exp(−ρ−2λ−σ2σ2t)g(t)(1−g(t))dt], (4.5)

with -dynamics of given by (4.3).

Note that the time in the transformed problem (4.5) is measured in terms of cumulative variance of the log return of the unaffected price, that is in ‘variance years’. One variance year corresponds to the physical time it takes to make . With one variance year is therefore equal to 25 calendar years. The new state variable corresponds to the size of temporary price impact as a percentage of current price, assuming inventory is completely liquidated at a constant rate over one variance year.

## 5 Singular initial value problem IVP0

Hereafter we refer to the IVP (3.4, 3.5) as . Note if and only if our standing assumption (2.8) holds. It has been shown in Brunovský et al. (2013) that is highly degenerate at For has infinitely many solutions with identical asymptotics near given by the formal power series

 hn(x)=x−23√2(a+b)x3/2+n∑i=2kix1+i/2,% n∈N, (5.1)

where are obtained recursively from

 kn+1 = 13(n+3)k1[kn((n+2)(2a−n)+4b) −12n−1∑j=1(3+j)(n−j+3)kj+1kn−j+1].

The series itself has zero radius of convergence for

 6a+4b−3=:K1>0>K2:=6a+2b−9,

see (Quittner, 2015, Remark 2). Asymptotic expansion of derivatives of is obtained by formal differentiation of the series in (5.1), ibid Theorem 1. Whenever or the power series ends at the 3rd element and constitutes a genuine solution of . This solution, however, is just one from a continuum and does not represent the optimal value function.

The highly degenerate nature of does not stem from the singularity of the linear terms in the ODE, which is well known and rather innocuous in the context of the Black-Scholes model, but from the singularity of the non-linear term. Liang (2009) studies singular IVPs of the form where is continuous. Note that the linear part of our ODE, belongs to Liang’s category, but the non-linear term does not.

Liang, too, observes multiplicity of solutions, but this multiplicity is less pronounced than in our case. In Liang’s work and uniquely determine the first derivatives of the solution, where , and, for non-integer , the solution becomes unique once the coefficient by has been specified. Therefore, in Liang’s case all solutions differ asymptotically by a multiple of near .

In contrast, has a continuum of solutions that differ asymptotically by

 xαexp(−β/√x),

where

 α:=2−23b,β:=√8(a+b),

see (Quittner, 2015, Theorems 2-5). These solutions invariably share their power series asymptotics to an arbitrary order as . A uniqueness result relevant for the current paper can be summarized as follows:

###### Proposition 5.1.

Under the assumption there is a unique solution of denoted by satisfying ,

 0≤u∞(x)≤x for x>0. (5.2)

The solution further satisfies for all as well as for .

###### Proof.

See Proposition 5.1 in Brunovský et al. (2013). ∎

Proposition 5.2 reveals certain qualitative characteristics of solutions of which can be observed empirically whenever an unstable numerical scheme is employed.

###### Proposition 5.2.

Any solution of (3.4) on with falls into one and only one of the following categories:

i) is constant;

ii) is strictly concave on ;

iii) is strictly convex on ;

iv) there is such that is strictly concave on , strictly convex on and for all ;

v) there is such that is strictly convex on , strictly concave on and for all .

###### Proof.

The conclusions follow readily from Brunovský et al. (2013), Lemma 4.1, applied to the equation

 x2y′′=(1+(a−2)x−y))y′+(a+b)y, (5.3)

with , obtained by differentiation and re-arrangement of (3.4). ∎

## 6 Boundary value problem BVP[0,∞)

In the context of the present paper it turns out to be advantageous to view Proposition 5.1 as a solution to a certain boundary value problem (BVP). We write whenever the limit on the right-hand side exists and complement the Dirichlet-type boundary condition with a Neumann-type boundary condition

 u′(∞)=0. (6.1)

Hereafter we refer to the mixed boundary value problem (3.4, 3.5, 6.1) as . It is seen below that the right-hand boundary condition (6.1) uniquely determines the solution found in Proposition 5.1.

###### Proposition 6.1.

Under the assumption has a unique solution which additionally satisfies , , , , as well as .

###### Proof.

possesses at least one solution, namely the solution identified in Proposition 5.1. Below we will prove uniqueness by showing that any solution of must also satisfy . By Lemma 3.1 in Brunovský et al. (2013) any local solution of the satisfies . Consider now the alternatives in Proposition 5.2 with and . Since any solution of also solves it cannot fall into the constant alternative i). Similarly, it cannot fall into category iii) with since then implies . Alternatives iv) and v) also imply . Therefore only category ii) remains as a possible alternative. One thus obtains globally, therefore is decreasing and implies We have thus proved and on integrating one obtains . This shows uniqueness by Proposition 5.1. ∎

The paper Brunovský et al. (2013) left two questions open. The first is whether the value function generated by the solution of from Proposition 5.1 via equation (3.3) is indeed the value function of the optimization problem (2.5). The second question concerns numerical computation of the solution to . We address both questions in turn, the former in Section 5 and the latter in Section 6.

## 7 Optimality

In this section we establish the precise connection between the boundary value problem and the optimal control and value function for the liquidation problem (2.5). We begin by formulating a natural sufficient condition for admissibility and investigate under what circumstances it is admissible to pursue further acquisition of the asset to be liquidated, .

###### Proposition 7.1.

Under the assumption (2.8) any predictable control satisfying is admissible. If additionally

 ρ>λ++r+, (7.1)

where then any predictable control satisfying for some is also admissible.

###### Proof.

i) We have and since is a GBM this implies for any finite and any which proves (2.6).

ii) To prove (2.7) first note that implies

 Z(t)≤zert+Kert−1r. (7.2)

To show integrability of the value function we first obtain an estimate of the integrand

 |v(t)(S(t)−ηv(t))|≤(v(t)++v(t)−)(S(t)+ηv(t)−)≤(K+v(t)+)(S(t)+ηK),

which for any bounded stopping time yields

 E[∫τ0e−ρt|v(t)(S(t)−ηv(t))|dt]≤K∫τ0e−ρt(E[S(t)]+ηK)dt+E[∫τ0e−ρtv(t)+(S(t)+ηK)dt]≤K∫∞0e−ρt(seλt+ηK)dtC<∞+E[∫τ0e−ρtv(t)+(S(t)+ηK)dt].

Continue with the integral inside the expectation in the second term, letting and integrating by parts. In preparation note which together with (7.2) implies for any bounded stopping time

 0 ≤ W(τ)=∫τ0v(t)+dt=∫τ0v(t)−dt+∫τ0rZ(t)dt+z−Z(τ) (7.3) ≤ Kτ+∫τ0r(zert+Kert−1r)dt+z=:g(τ).

Integration by parts yields

 ∫τ0e−ρtv(t)+(S(t)+ηK)dt = e−ρτW(τ)(S(τ)+ηK) +ρ∫τ0e−ρtW(t)(S(t)+ηK)dt−∫τ0e−ρtW(t)dS(t).

We continue with the second term on the right-hand side. Let then with from (7.3) which implies that is a (square-integrable) martingale. Hence for any bounded stopping time

 E[∫τ0e−ρtW(t)dS(t)]=E[∫τ0e−ρtλW(t)S(t)dt].

Pulling everything together

 E[∫τ0e−ρt|v(t)(S(t)−ηv(t))|dt]≤C+E[∫τ0e−ρtv(t)+(S(t)+ηK)dt].

The right-hand side is bounded for under the standing assumption (2.8). This is also true for if additionally and The last three inequalities together with the standing assumption (2.8) are equivalent to (7.1). Letting increase to we have by monotone convergence

 E[∫T0e−ρt|v(t)(S(t)−ηv(t))|dt]<∞.

The next theorem characterizes the optimal liquidation strategy and the corresponding value function. The inequality confirms the initial intuition that without short sales the implementation shortfall must be positive. We note that due to we have i.e. it is not optimal to buy more of the liquidated asset, even when (for ) strategies that involve further purchases are admissible.

###### Theorem 7.2.

Assume (2.8). Let be the unique solution of , with given by (3.6). Then the function is the value function of the optimization (2.5) and

 v∗(t):=12η(S(t)−Vz(S(t),Z∗(t)))=S(t)2η(1−u′∞(ησ2Z∗(t)S(t)))≥0 (7.4)

is the optimal control among all admissible controls defined in equations (2.6,2.7).

###### Proof.

To prove the theorem we apply the ‘Verification’ Theorem IV.5.1 of Fleming and Soner (2006). To this end, we have to check the following:

1. is and satisfies

 |V(s,z)|≤K(1+|(s,z)|m)

for some ;

2. for all admissible controls, where and ;

3. For any finite time

 limt→∞e−ρtE(s,z)[It≤T(Z∗=0)V(S(t),Z∗(t))]=0,

being the solution of

 dS(t) =λS(t)dt+σS(t)dB(t), dZ∗(t) =(rZ∗(t)−S(t)2η(1−u′(ησ2Z∗(t)S(t))))dt.

The regularity properties as well as the estimates of (i) and (ii) are immediate consequences of the properties of , which in particular imply

 0≤V(s,z)≤sz. (7.5)

The estimate (7.2) gives which in combination with inequality (7.5) and standing assumption (2.8) yields

 0 ≤ e−ρtE(s,z)[It≤T(Z∗=0)V(S(t),Z∗(t))] ≤ e−ρtzertE(s,z)[S(t)]=sze(r+λ−ρ)t↘0.

This proves item (iii). ∎

Observe that the optimal control deviates from the myopic strategy of maximizing the integrand of the objective function . In addition to the instantaneous impact on the execution price the current liquidation rate also affects future levels of the inventory . In (7.4) the optimal strategy at time differs from  by the amount , which is the marginal value of the optimal revenue with respect to the size of the remaining inventory. It follows that taking proper account of the role of future inventory level reduces the selling rate. By Proposition 5.1, is positive and decreasing to zero and so is in and therefore for large values of the selling rate is very close to the myopic strategy. For small values of the optimal rate of trading is non-linear, roughly proportional to as can be seen from the asymptotic expansion (5.1) and the formula for the optimal trading rate (7.4).

We remark that the classical martingale case with and fixed time horizon yields constant optimal liquidation speed . The resulting price impact per share, for fixed , is proportional to which is not consistent with broad empirical evidence that indicates power dependence roughly proportional to .

When estimating price impact empirically, an assumption has to be made about the rate of trading. In Almgren et al. (2005) this rate is assumed to be constant and the temporary impact of individual trades is estimated proportional to which yields per-share temporary price impact proportional to . Here, in contrast, the temporary impact is linear, proportional to , but the optimal rate of trading is non-linear, roughly proportional to for small values. ‘Small’ must be understood in context; we find that asymptotics is perfectly compatible with meta-orders whose optimal execution lasts several days, see Section 8.4.

We can also make qualitative conclusions about the optimized implementation shortfall by studying the asymptoptic expansion (5.1) whereby we find that for small the per-share price impact equals

 I(S(0),Z(0))=S(0)Z(0)−V(S(0),Z(0))S(0)Z(0)=43√η(ρ−λ−r)Z(0)/S(0)+O(Z(0)3/2),

which means that the price impact is proportional to the square root of the total trade size. There is a strong empirical evidence to support the square root law for meta-orders, see Bershova and Rakhlin (2013), Farmer et al. (2013), Donier et al. (2015) and Tóth et al. (2016) and references therein.

## 8 Computation of the solution

To make BVP amenable to numerical treatment we first truncate the spatial interval to with and solve the ODE (3.4) with mixed boundary conditions and . We refer to the truncated boundary value problem as BVP. In section 8.1 we prove that the solution of BVP is unique and that it converges pointwise upwards to the desired solution as .

Numerical solutions of BVPs for ordinary differential equations with singular coefficients have a well established literature, see for example Jamet (1969), Weinmüller (1984), Weinmüller (1986), and Auzinger et al. (1999) who consider BVPs with ODE of the form

 u′′=x−1A(x)u′+x−2B(x)u+F(x,u,u′), (8.1)

where and are continuous at and one of the boundaries is . Numerical solution of (8.1) can be computed by means of the Matlab function bvp5c after transformation , see Weinmüller (1986), equation (2.1a).

However, as we have mentioned already in the connection with IVP, our problem BVP is substantially more singular. This is not due to the singularity in the linear terms of ODE (3.4), which in fact can be accommodated in the ansatz (8.1), but because the non-linear part is not continuous in at zero. Attempts to compute the solution of BVP by some kind of shooting fail – both at and the trajectories blow up. Algorithm bvp5c is able to produce, with careful tuning of input parameters, a stable solution of BVP for not too close to zero. However, the quality of this solution near zero is poor, as can be seen in panel (b) of Figure 3.

To bypass the troublesome singularity at zero we introduce a time dimension into BVP in a strategy akin to the value function iteration method known from financial economics. This approach is also common in linear-quadratic optimal control problems where, however, it is not motivated by the presence of singularities, see Anderson and Moore (1989, Section 3.1).

We consider a parabolic PDE that corresponds to a finite horizon version of the time-homogeneous optimization (2.5). We formulate suitable boundary conditions on a finite spatial interval to obtain a parabolic problem BVP and show that its solution converges monotonically to the solution of BVP as . This is done in section 8.2. Unfortunately, BVP does not correspond to an optimal control problem due to the choice of boundary conditions.

In section 8.3 we formulate a finite difference scheme to solve BVP numerically. This scheme is well behaved with respect to the singularity at and produces a reliable approximation to , which for large enough is arbitrarily close to the desired solution .

### 8.1 Problem BVP[0,l]

###### Theorem 8.1.

Let . For given has a unique solution such that for all . The solution is strictly increasing, concave and satisfies for , and for , where is the unique solution of .

###### Proof.

Step 1) For any such that the function , resp. is a lower (resp. upper) solution of BVP in the sense of Definition II.1.1 in De Coster and Habets (2006), which crucially allows for the Neumann boundary condition at . Therefore by Theorem II.1.3 ibid the solution of the mixed boundary value problem BVP satisfies

 0≤uε(x)≤x for every ε>0. (8.2)

From here the proof proceeds as in Proposition 2.2 of Brunovský et al. (2013). From Bernstein’s condition Bernstein (1904) (see also Section I.4.3 of De Coster and Habets (2006) for related Nagumo condition) fixing we obtain a uniform (in ) a-priori bound on the derivative on . Together with (8.2) this yields via (3.4) an a-priori bound on on which means (as well as ) are equicontinuous on which in turn implies equicontinuity of via (3.4). One can thus extract a convergent subsequence of which convergences with its first two derivatives to some function on with and such that solves (3.4).

Step 2) By Brunovský et al. (2013), Lemma 3.1, . This, together with the conditions and excludes all alternatives of Proposition 5.2 except for ii). Therefore any solution of BVP must be concave and increasing on .

Step 3) To prove uniqueness of the solution assume that and are two solutions of BVP. Then solves

 x2p′′=axp′+bp−p′(u′−1)−12(p′)2, (8.3)

on which on differentiation yields

 x2p′′′=((a−2)x+1−u′−p′)p′′+(a+b−u′′)p′. (8.4)

Applying Lemma 4.1 of Brunovský et al. (2013) to (8.4) with , and , one obtains that obeys the same alternatives as in Proposition 5.2.

By construction we have , therefore alternatives (ii)-(v) of Proposition 5.2 are excluded and must be constant and thus necessarily equal to zero. Thus BVP has a unique solution which we denote by .

Step 4) Now we prove that the solutions grow with . Take and let , . Consider on which satisfies (8.3), (8.4) and therefore obeys the alternatives of Proposition 5.2.. As before we have . Since while we also have Hence in Proposition 5.2 (iii) is the only possible alternative, is strictly convex on and therefore on which implies and on .

Step 5) It remains to be proved that for , converges pointwise to the solution of BVP. Step 2) implies and by step 4)  is increasing in therefore for fixed the limit is well defined. Likewise and is increasing in hence we have a well-defined limit . Picking arbitrary and in we rewrite (3.4) in integral form

 uL(x) =uL(x0)+∫xx0u′L(ξ)dξ, (8.5) u′L(x) =u′L(x0)+∫xx0f(ξ,uL(ξ),u′L(ξ))dξ (8.6)

with

 f(x,u,v)=avx+bux2−12(v−1)2