Dynamic intertemporal utility optimization by means of Riccati transformation of Hamilton-Jacobi Bellman equation

TL;DR

This paper introduces a Riccati transformation method to solve a complex Hamilton-Jacobi-Bellman equation in stochastic portfolio optimization, enabling more accurate computation of optimal strategies considering intertemporal utility.

Contribution

The paper presents a novel Riccati method for transforming nonlinear HJB equations into a quasi-linear form, facilitating numerical solutions in portfolio optimization.

Findings

The numerical scheme achieves second-order convergence.

Optimal strategies are significantly affected by intertemporal utility considerations.

Application to German DAX 30 data demonstrates practical relevance.

Abstract

In this paper we investigate a dynamic stochastic portfolio optimization problem involving both the expected terminal utility and intertemporal utility maximization. We solve the problem by means of a solution to a fully nonlinear evolutionary Hamilton-Jacobi-Bellman (HJB) equation. We propose the so-called Riccati method for transformation of the fully nonlinear HJB equation into a quasi-linear parabolic equation with non-local terms involving the intertemporal utility function. As a numerical method we propose a semi-implicit scheme in time based on a finite volume approximation in the spatial variable. By analyzing an explicit traveling wave solution we show that the numerical method is of the second experimental order of convergence. As a practical application we compute optimal strategies for a portfolio investment problem motivated by market financial data of German DAX 30 Index and show the effect of considering intertemporal utility on optimal portfolio selection.

Figures

Tables

Table 1. Table 1. The L ∞ ( ( 0 , T ) : L 2 ( x L , x R ) ) L_{\infty}((0,T):L_{2}(x_{L},x_{R})) and L ∞ ( ( 0 , T ) : L ∞ ( x L , x R ) ) L_{\infty}((0,T):L_{\infty}(x_{L},x_{R})) norm of the error of the numerical solution with the spatial step h ℎ h and time step k = h 2 𝑘 superscript ℎ 2 k=h^{2} and the exact traveling wave solution. Corresponding experimental orders of convergence.

$h$	$L_{\infty }((0,T):L_2)$-err	$EO{C}_{k=h^2}$	$L_{\infty }((0,T):L_{\infty }))$-err	$EO{C}_{k=h^2}$
0.05	1.1886e-01	–	5.8577e-02	–
0.025	3.2102e-02	1.8885	1.5919e-02	1.8796
0.0125	8.1969e-03	1.9695	4.0718e-03	1.9670
0.01	5.2598e-03	1.9882	2.6133e-03	1.9874
0.005	1.3196e-03	1.9949	0.6558e-03	1.9945

Equations

$$d{x}_t=\mu (x_t,t,{\mathbf{\theta }}_t)dt+\sigma (x_t,t,{\mathbf{\theta }}_t)d{W}_t,$$

$$d{x}_t^i={\mu }^{i}dt+\sum _{k=1}^n{\sigma }^{ki}d{W}_t^k$$

$$d{x}_t^{\mathbf{\theta }}={\mathbf{\mu }}^T\mathbf{\theta }dt+\sigma (\mathbf{\theta })d{W}_t$$

$$\sigma (x,t,\mathbf{\theta }{)}^2={\mathbf{\theta }}^T\mathbf{\Sigma }\mathbf{\theta },$$

$$\mu (x,t,\mathbf{\theta })={\mathbf{\mu }}^T\mathbf{\theta }-\frac{1}{2}\sigma (x,t,\mathbf{\theta }{)}^2+\epsilon e^{-x}+r,$$

$$d{y}_t^{\tilde{\mathbf{\theta }}}=(\epsilon +(r+\mu (\tilde{\mathbf{\theta }}))y_t^{\tilde{\mathbf{\theta }}})dt+\sigma (\tilde{\mathbf{\theta }})y_t^{\tilde{\mathbf{\theta }}}d{W}_t,$$

$$\sigma (x,t,\mathbf{\theta }{)}^2=\underset{\mathbf{\Sigma }\in \mathcal{K}}{\max }{\mathbf{\theta }}^T\mathbf{\Sigma }\mathbf{\theta },$$

$$\mu (x,t,\mathbf{\theta })=\underset{\mathbf{\mu }\in \mathcal{E}}{\min }{\mathbf{\mu }}^T\mathbf{\theta }-\frac{1}{2}\sigma (x,t,\mathbf{\theta }{)}^2+\epsilon e^{-x}+r,$$

$$\max_{{\bm{\theta}}|_{[0,T)}}\mathbb{E}\left[u(x_{T}^{\bm{\theta}})\,+\int_{0}^{T}c(x_{s}^{\bm{\theta}},s)ds\big{|}\,x_{0}^{\bm{\theta}}=x_{0}\right],$$

$$V(x,t):=\underset{\mathbf{\theta }{|}_{[t,T)}}{\sup }\mathbb{E}\left[u(x_T^{\mathbf{\theta }})+{\int }_t^{T}c(x_s^{\mathbf{\theta }},s)ds|x_t^{\mathbf{\theta }}=x\right]$$

$${\partial }_{t}V+\underset{\mathbf{\theta }\in \Delta }{\max }\left(\mu (x,t,\mathbf{\theta }){\partial }_{x}V+\frac{1}{2}\sigma (x,t,\mathbf{\theta }{)}^2{\partial }_x^{2}V\right)+c(x,t)=0,$$

$$V(x,T)=u(x),$$

$$u(x)=-e^{-ax},$$

$$a(x)=-\frac{u^{\prime \prime }(x)}{u^{\prime }(x)}forx\in \mathbb{R}.$$

$$c(x,t)=-\kappa e^{-dx-\varrho (T-t)},$$

$$\phi (x,\tau )=-\frac{{\partial }_x^{2}V(x,t)}{{\partial }_{x}V(x,t)},where\tau =T-t.$$

$${\partial }_{t}V-\alpha (\cdot ,\phi ){\partial }_{x}V+c=0,V(\cdot ,T)=u(\cdot ),$$

$$\alpha (x,\tau ,\phi )=\underset{\mathbf{\theta }\in \Delta }{\min }\left(-\mu (x,t,\mathbf{\theta })+\frac{\phi }{2}\sigma (x,t,\mathbf{\theta }{)}^2\right),\tau =T-t.$$

$$\alpha(\varphi)=\left\{\begin{array}[]{ll}A\varphi-\frac{B}{\varphi}+C,&\hbox{if}\ \varphi>\varphi_{*},\\ E\varphi+D,&\hbox{if}\ \varphi\leq\varphi_{*},\end{array}\right.$$

$${\partial }_x\alpha (x,\tau ,\phi )={\alpha }_x^{\prime }(x,\tau ,\phi )+{\alpha }_{\phi }^{\prime }(x,\tau ,\phi ){\partial }_x\phi ,$$

$$-{\partial }_{\tau }\phi +{\partial }_x\left({\partial }_x\alpha (\cdot ,\phi )-\alpha (\cdot ,\phi )\phi \right)=\frac{1}{b(T-\tau )}{\partial }_x\left(e^{{\int }_{x_{\ast }}^x\phi (\xi ,\tau )d\xi }{\partial }_{x}c\right),$$

$$\phi (x,0)=-u^{\prime \prime }(x)/u^{\prime }(x),(x,\tau )\in \mathbb{R}\times (0,T),$$

$$V(x,t)=a(t)+b(t){\int }_{x_{\ast }}^x{e}^{-{\int }_{x_{\ast }}^{\xi }\phi (\eta ,\tau )d\eta }d\xi ,t=T-\tau ,$$

$$\frac{d}{dt}a(t)$$

$$\frac{d}{dt}b(t)$$

$$-{\partial }_{\tau }\phi =-\frac{{\partial }_x^2{\partial }_{t}V}{{\partial }_{x}V}+\frac{{\partial }_x^{2}V{\partial }_x{\partial }_{t}V}{({\partial }_{x}V{)}^2}=-\frac{{\partial }_x^2{\partial }_{t}V}{{\partial }_{x}V}-\phi \frac{{\partial }_x{\partial }_{t}V}{{\partial }_{x}V},$$

$${\partial }_x^{2}V=-\phi {\partial }_{x}V,and{\partial }_x^{3}V=-{\partial }_x(\phi {\partial }_{x}V)=({\phi }^2-{\partial }_x\phi ){\partial }_{x}V,$$

$$-{\partial }_{\tau }\phi$$

$$V(x,T)$$

$$\frac{d}{dt}\left(b(t)e^{-{\int }_t^T\omega (\eta )d\eta }\right)=-{\partial }_{x}c(x_{\ast },t)e^{-{\int }_t^T\omega (\eta )d\eta }\le 0.$$

$$b(t)\ge b(T)e^{{\int }_t^T\omega (\eta )d\eta }=u^{\prime }(x_{\ast })e^{{\int }_t^T\omega (\eta )d\eta }>0,$$

$${\int }_{x_{\ast }}^x{\partial }_{\xi }V\left({\partial }_{\xi }\alpha -\alpha \phi \right)d\xi$$

Full text

AMS-MOS Classification: 35K55, 34E05, 70H20, 91B70, 90C15, 91B16

Keywords: Dynamic stochastic portfolio optimization, Dynamic utility, Hamilton-Jacobi-Bellman equation, Riccati transformation, finite volume scheme.

1. Introduction

In this paper, we investigate the impact of presence of a nontrivial intertemporal utility function on stochastic optimal portfolio selection problem. The problem can be formulated in terms of the expected terminal and intertemporal utility maximization problem, in which the underlying stochastic process is controlled by a vector of time-dependent weights of assets entering a financial portfolio.

To solve the expected terminal and intertemporal utility maximization problem, we follow a methodology based on solving a fully nonlinear parabolic Hamilton-Jacobi-Bellman equation for the intermediate value function of the corresponding optimal control problem. A similar problem was investigated by Federico, Gassiat and Gozzi in [7], where they studied a problem of terminal and intertemporal utility maximization in an investment-consumption portfolio setting and the current utility being dependent also on the wealth process. They studied properties of solutions to a dual control problem.

The novelty of our paper is generalization of the transformation method proposed and analyzed by Abe and Ishimura [1], Ishimura and Ševčovič [10] and Kilianová and Ševčovič [12, 13] for the case of a nontrivial intertemporal utility function. The transformation is also referred to as the Riccati transformation as it involves the ratio between the second and the first derivatives of the value function. The transformed function can be viewed as the absolute risk aversion coefficient of an investor. Secondly, we generalize the underlying stochastic process to more general processes with arbitrary drift and volatility functions. Such a general setting can include, in particular, processes arising in the so-called worst-case portfolio optimization studied recently by Kilianová and Trnovská in [14]. In contrast to the problem involving the terminal utility maximization only (cf. [1],[10], [12], [13]), the resulting transformed equation is a non-local quasi-linear parabolic equation containing non-local terms involving the intertemporal utility function. The non-local parabolic equation can be further transformed into a coupled system of two quasi-linear local parabolic equations. We analyze these governing equations and show how their solutions are related to solving the original HJB equation. As a tool for solving the associated terminal and intertemporal utility maximization problem, we generalize the numerical method proposed by Kilianová and Ševčovič [12, 13] for the case when a non-local term appears in the quasi-linear parabolic equation. We furthermore derive a-priori lower and upper bounds of a solution which are given in terms of the risk aversion coefficients of the terminal and intertemporal utility functions. The main advantage of the Riccati transformation method is twofold. First, the transformed function has a practical representation and meaning as an intertemporal risk aversion of an investor and it is a globally bounded function even in the case when the utility function is unbounded. Moreover, there are natural boundary conditions for a solution defined on a truncated numerical domain. Secondly, the nonlinearity appearing as a diffusion function in the transformed equation can be computed in a fast and efficient way using modern tools of conic convex programming.

As a practical application we compute optimal strategies for a portfolio investment problem motivated by the market financial data from German DAX 30 Index. We compare the optimal portfolio selection strategies for the case of absence of intertemporal utility maximization to the case when a non-trivial intertemporal utility is considered. We illustrate the effect of an intertemporal utility function on the optimal portfolio selection.

The paper is organized as follows. In the next section we introduce and discuss basic model assumptions made on the underlying stochastic process. The process of logarithmized portfolio wealth $x_{t}$ at time $t$ is controlled by a vector of weights ${\bm{\theta}}_{t}$ belonging to a compact convex subset of $\mathbb{R}^{n}$. A dynamic stochastic optimization problem with intertemporal utility is formulated in Section 3. Following the Bellman optimality principle, we present a fully nonlinear backward parabolic equation for the intermediate value function satisfying a given terminal condition. In Section 4 we present the so-called Riccati transformation of the value function, leading us from the fully nonlinear HJB equation to a single quasi-linear parabolic equation in the divergent form containing a non-local term. We furthermore analyze qualitative properties of an auxiliary value function arising from a parametric convex programming problem. Existence of a classical Hölder smooth solution and its a-priori bounds are also derived in this section. Section 5 is devoted to a numerical approximation scheme for solving the transformed quasi-linear parabolic function. The scheme is based on the finite volume approximation method involving dual finite volumes. We compare the numerical scheme to the fixed policy iteration method for solving HJB equations as investigated by e.g. by Huang et al. [8], Reisinger and Witte [28]. Finally, in Section 6 we test the accuracy of the proposed numerical method on an explicit traveling wave example and compute the experimental order of convergence. We show that the experimental order of convergence is approximately two which indicates the second order of convergence of the numerical method. Subsequently, we apply the proposed method to optimal portfolio selection problem and present corresponding results with and without intertemporal utility maximization.

2. Underlying stochastic process with control

Throughout the paper we will assume that the underlying stochastic process $\{x_{t}\}$ satisfies the stochastic differential equation (SDE)

[TABLE]

where the control process $\{{\bm{\theta}}_{t}\}$ is adapted to the process $\{x_{t}\}$, $\{W_{t}\}$ is the standard one-dimensional Wiener process and functions $(x,t,{\bm{\theta}})\mapsto\mu(x,t,{\bm{\theta}})$ and $(x,t,{\bm{\theta}})\mapsto\sigma(x,t,{\bm{\theta}})^{2}$ are $C^{1,1}$ smooth in $x,t$ and ${\bm{\theta}}$ variable, i.e. their first derivatives are Lipschitz continuous functions.

Remark 1**.**

The motivation for studying the SDE (1) as an underlying stochastic process controlled by ${\bm{\theta}}$ arises from the stochastic dynamic optimal portfolio management. Let $x_{t}^{i}=\ln y_{t}^{i}$ denote the logarithm of the asset value $y_{t}^{i}$ entering a portfolio consisting of $n$ assets with vector of weights ${\bm{\theta}}$. Then $dx_{t}^{i}=dy_{t}^{i}/y_{t}^{i}$ is the return on the asset $i$. Suppose that the process of each such a return is driven by

[TABLE]

where $W_{t}^{j}$ is a one-dimensional Wiener process such that the increments $dW_{t}^{j}$ and $dW_{t}^{i}$ are independent for $j\not=i$. The mean return on the increment of the portfolio $x^{\bm{\theta}}=\sum_{i=1}^{n}\theta^{i}x^{i}$ with the vector of weights ${\bm{\theta}}$ is ${\bm{\mu}}^{T}{\bm{\theta}}\,dt$ and its variance is equal to $\sum_{i,j,k=1}^{n}\theta^{i}\sigma^{ki}\sigma^{kj}\theta^{j}\,dt$.

Following Merton [23, 24] we can describe the stochastic process $x_{t}^{\bm{\theta}}$ by the following one-dimensional SDE of the form (1):

[TABLE]

where $W_{t}$ is the one-dimensional Wienner process, $\sigma({\bm{\theta}})^{2}={\bm{\theta}}^{T}{\bm{\Sigma}}{\bm{\theta}}$ and ${\bm{\Sigma}}$ is the covariance matrix, ${\bm{\Sigma}}_{ij}=\sum_{k=1}^{n}\sigma^{ki}\sigma^{kj}$.

Example 1**.**

As an example of the stochastic process (1), one can consider a portfolio optimization problem with regular cash inflow (e.g. pension planning). In this example, the volatility function is given by

[TABLE]

where ${\bm{\Sigma}}$ is a positive definite covariance matrix of asset returns. The drift function is given by

[TABLE]

where ${\bm{\mu}}$ is the vector of mean returns of assets, $\varepsilon$ is an inflow ($\varepsilon>0$)/outflow ($\varepsilon<0$) to/from the portfolio, $r\geq 0$ is an interest rate of a risk-free bond. The stochastic process $\{x^{\bm{\theta}}_{t}\}$ controlled by $\{{\bm{\theta}}_{t}\}$ is a logarithmic transformation of the stochastic process $\{y_{t}^{\tilde{\bm{\theta}}}\}_{t\geq 0}$ driven by the stochastic differential equation

[TABLE]

where $\tilde{\bm{\theta}}(y,t)={\bm{\theta}}(x,t)$ with $x=\ln y$ (cf. Kilianová and Ševčovič [12]).

Another example stems from the so-called worst-case portfolio optimization problem investigated by Kilianová and Trnovská [14]. The volatility function is given by

[TABLE]

where $\mathcal{K}$ is an uncertainty convex set of positive definite covariance matrices. Typically, only a part of the covariance matrix is exactly prescribed while other entries are not precisely determined. For instance, if only the diagonal $d$ is known, we have ${\mathcal{K}}=\{{\bm{\Sigma}}\succ 0,\ diag({\bm{\Sigma}})=d\}$. The drift function is given by

[TABLE]

where $\mathcal{E}$ is a given uncertainty convex set of mean returns.

3. Dynamic stochastic optimization problem with intertemporal utility function

Our goal is to extend the model of terminal utility maximization studied previously in Kilianová and Ševčovič [12] by including an intertemporal utility function. Maximization of dynamic utility has been investigated in vast literature in the past by a number of methods. In this paper, we assume that the investor has a certain utility $c$ from intertemporal wealth but a different utility $u$ from terminal wealth. We assume the overall utility to be time-additive. Then we can formulate the problem of dynamic utility maximization as follows:

[TABLE]

(c.f. [7] where they included consumption as well). Here $\{x_{t}^{{\bm{\theta}}}\}$ is Itō’s stochastic process of the form (1) on a finite time horizon $[0,T]$, $u:\mathbb{R}\to\mathbb{R}$ is a given terminal utility function and $x_{0}$ a given initial state condition of $\{x_{t}^{\bm{\theta}}\}$ at $t=0$. The function ${\bm{\theta}}:\mathbb{R}\times[0,T)\to\mathbb{R}^{n}$ maps $(x,t)\mapsto{\bm{\theta}}(x,t)$ and it represents an unknown control function governing the underlying stochastic process $\{x_{t}^{\bm{\theta}}\}$. The function $c:\mathbb{R}\times[0,T)\to\mathbb{R}$ is the intertemporal utility function. In what follows, we will assume $c$ is a $C^{2}$ smooth function and it is non-decreasing in the $x$ variable. Clearly, one can add a time discounting factor to both utility functions in (5). We assume that the control parameter ${\bm{\theta}}$ belongs to a closed convex subset $\Delta$ of the compact simplex $\mathcal{S}^{n}=\{{\bm{\theta}}\in\mathbb{R}^{n}\ |\ {\bm{\theta}}\geq\mathbf{0},\mathbf{1}^{T}{\bm{\theta}}=1\}\subset\mathbb{R}^{n}$, where $\mathbf{1}=(1,\cdots,1)^{T}\in\mathbb{R}^{n}$.

If we introduce the value function

[TABLE]

then $V(x,T):=u(x)$. Following Bertsekas [5], the value function $V=V(x,t)$ satisfies the fully nonlinear Hamilton-Jacobi-Bellman (HJB) parabolic equation:

[TABLE]

for $(x,t)\in\mathbb{R}\times[0,T)$; see also [12] and Kossaczký, Ehrhardt, Günther [17, 18].

As an example of the terminal utility function, one can consider e.g. a constant absolute risk aversion (CARA) function:

[TABLE]

with constant absolute risk aversion $a\equiv a(x)>0$, where

[TABLE]

We note that a CARA utility function in the variable $x$ corresponds to a CRRA (constant relative risk aversion) utility function $u(y)=-y^{-a}$ in the variable $y=e^{x}$. In practical applications, $y$ can stand for a portfolio value and $x=\ln y$ its log-transform, for which (1) holds.

Another choice for the utility function $u$ can be a decreasing/increasing absolute risk aversion (DARA/IARA) function with $a(x)$ decreasing/increasing in the $x$ variable. Typically, the intertemporal utility function $c$ is a non-decreasing concave discounted function, i.e.

[TABLE]

where $\kappa,d\geq 0$ and $\varrho$ is a discounting factor. We note that including a discount factor $e^{-rT}$ into the terminal utility function $u$ does not play any role in the solution, as one can transfer this constant into the coefficient $\kappa$ of the intertemporal utility function $c$ simply by multiplying (5) by $e^{rT}$.

4. The Riccati transformation of the HJB equation with intertemporal utility function

Following the papers by Abe and Ishimura [1], Ishimura and Ševčovič [10] and Kilianová and Ševčovič [12], the Riccati transformation of the value function $V$ can be introduced as follows:

[TABLE]

Suppose for a moment that the value function $V(x,t)$ is increasing in the $x$-variable. This is a natural assumption in the case when the terminal utility function $u(x)$ is increasing itself. Then the HJB equation (7) can be rewritten as follows:

[TABLE]

where $\alpha(x,\tau,\varphi)$ is the value function of the following parametric optimization problem:

[TABLE]

The proof of the following result is a straightforward generalization of [12, Theorem 4.1] for a more general drift and volatility functions depending on $x$ and $t$ and therefore it is omitted.

Theorem 1**.**

Assume that the functions $(x,t,{\bm{\theta}})\mapsto\mu(x,t,{\bm{\theta}})$ and $(x,t,{\bm{\theta}})\mapsto\sigma(x,t,{\bm{\theta}})^{2}$ are $C^{1,1}$ smooth in $x,t$ and ${\bm{\theta}}$ variables, and such that the objective function $f(x,t,\varphi,{\bm{\theta}}):=-\mu(x,t,{\bm{\theta}})+\frac{\varphi}{2}\sigma(x,t,{\bm{\theta}})^{2}$ is strictly convex in the variable ${\bm{\theta}}\in\Delta$ for any $\varphi\in(\varphi_{min},\infty)$ where $\Delta\subset\mathbb{R}^{n}$ is a compact convex set. Then the value function $\alpha$ is $C^{1,1}$ smooth for $x\in\mathbb{R},\tau\in[0,T],\varphi\in(\varphi_{min},\infty)$. Moreover, $\varphi\mapsto\alpha(\cdot,\varphi)$ is a strictly increasing function. For the derivative of $\alpha$ we have $\alpha^{\prime}_{\varphi}(x,\tau,\varphi)=(1/2)\sigma(x,T-\tau,\hat{\bm{\theta}}(x,\tau,\varphi))^{2}$ where $\hat{\bm{\theta}}(x,\tau,\varphi)\in\Delta$ is the argument of the minimum of $\alpha(x,\tau,\varphi)$ with respect to ${\bm{\theta}}$.

Example 2**.**

If we consider an example of the decision set $\Delta=\{{\bm{\theta}}\in\mathbb{R}^{2},\,\theta_{1},\theta_{2}\geq 0,\theta_{1}+\theta_{2}=1\},n=2,\mu={\bm{\mu}}^{T}{\bm{\theta}},\sigma^{2}={\bm{\theta}}^{T}{\bm{\Sigma}}{\bm{\theta}}$, then the value function $\alpha$ has the form:

[TABLE]

where constants $A>0,B>0,C,D,E>0,\varphi_{*}>0$ depend on the mean return vector ${\bm{\mu}}$ and covariance matrix ${\bm{\Sigma}}$ and are such that $\alpha$ is $C^{1,1}$ continuous function having one point of discontinuity of the second derivative $\alpha^{\prime\prime}$ at $\varphi_{*}$. The minimizer $\hat{\bm{\theta}}=\hat{\bm{\theta}}(\varphi)$ increases the number of positive weights when $\varphi$ passes through $\varphi_{*}$. For $n>2$, the number of discontinuities of $\alpha_{\varphi}^{\prime\prime}$ increases (cf. [12]).

In what follows, we shall denote by $\partial_{x}\alpha$ the total differential of the function $\alpha(x,\tau,\varphi)$ where $\varphi=\varphi(x,\tau)$, that is

[TABLE]

where $\alpha^{\prime}_{x}$ and $\alpha^{\prime}_{\varphi}$ are partial derivatives of $\alpha$ with respect to variables $x$ and $\varphi$, respectively.

The relationship between the transformed function $\varphi$ and the value function $V$ is given by the following theorem.

Theorem 2**.**

Assume that the utility function $u(x)$ and the intertemporal utility function $c(x,t)$ are $C^{2}$ smooth functions and such that $u$ is increasing and $c$ is non-decreasing in the $x$ variable. Then an increasing function $V(x,t)$ in the $x$ variable is a solution to the Hamilton-Jacobi-Bellman equation (7) if and only if the transformed function $\varphi(x,\tau)=-\partial_{x}^{2}V(x,t)/\partial_{x}V(x,t),\,\ t=T-\tau$, is a solution to the quasi-linear parabolic non-local PDE:

[TABLE]

and

[TABLE]

where the functions $a(t)$ and $b(t)$ are solutions to the system of ODEs:

[TABLE]

Here $x_{*}\in\mathbb{R}$ is a fixed real number, $\gamma(t):=\alpha(x_{*},\tau,\varphi(x_{*},\tau))$, and $\omega(t):=\partial_{x}\alpha(x_{*},\tau,\varphi(x_{*},\tau))-\alpha(x_{*},\tau,\varphi(x_{*},\tau))\varphi(x_{*},\tau)$ where $\tau=T-t$.

P r o o f. Let $V$ be a solution to the HJB equation (7) satisfying the terminal condition $V(x,T)=u(x)$ and such that $\partial_{x}V(x,t)>0$ for each $(x,t)\in\mathbb{R}\times[0,T)$. Thus $V$ solves (10), i.e. $\partial_{t}V=\alpha(x,\tau,\varphi)\,\partial_{x}V-c$ where $\varphi=-\partial^{2}_{x}V/\partial_{x}V$. Therefore, $V$ is given by (14) with $a(t)=V(x_{*},t)$ and $b(t)=\partial_{x}V(x_{*},t)$.

Since

[TABLE]

[TABLE]

it follows from the equation $\partial_{t}V-\alpha(\cdot,\varphi)\partial_{x}V+c=0$ that $\varphi$ satisfies:

[TABLE]

It means that the function $\varphi$ is a solution to the Cauchy problem (12)–(13). By differentiating (10) with respect to $x$ we obtain $\partial_{t}\partial_{x}V=\partial_{x}(\alpha\partial_{x}V)-\partial_{x}c=\partial_{x}\alpha\,\partial_{x}V+\alpha\partial^{2}_{x}V-\partial_{x}c=(\partial_{x}\alpha-\alpha\,\varphi)\partial_{x}V-\partial_{x}c$. Taking $x=x_{*}$ we conclude $\partial_{t}\partial_{x}V(x_{*},t)=\omega(t)\partial_{x}V(x_{*},t)-\partial_{x}c(x_{*},t)$. As $\partial_{x}V(x_{*},T)=U^{\prime}(x_{*})$ we obtain $b(t)=\partial_{x}V(x_{*},t)$ is the solution to the ODE (16). Furthermore, as $\partial_{t}V(x_{*},t)=\alpha(x_{*},\tau,\varphi(x_{*},\tau))\partial_{x}V(x_{*},t)-c(x_{*},t)=\gamma(t)b(t)-c(x_{*},t)$, $\tau=T-t$, and $V(x_{*},T)=u(x_{*})$ we conclude $a(t)=V(x_{*},t)$ solves (15), as claimed.

On the other hand, suppose that a function $\varphi$ solves (12)–(13) and functions $a,b$ solve (15)–(16). Then $V(x,t)$ given by (14) satisfies $-\partial^{2}_{x}V(x,t)/\partial_{x}V(x,t)=\varphi(x,\tau)$, $\tau=T-t$, and

[TABLE]

The function $b(t)$ satisfying (16) is a positive function. Indeed, as $c$ is non-decreasing at $x_{*}$ we have

[TABLE]

Integrating the above inequality over $(t,T)$ we obtain $b(T)-b(t)e^{-\int_{t}^{T}\omega(\eta)d\eta}\leq 0$ and so

[TABLE]

for any $t\in[0,T]$. Furthermore, as $\partial_{x}V(x,t)=b(t)e^{-\int_{x_{*}}^{x}\varphi(\xi,\tau)d\xi}>0$, the function $V(x,t)$ is increasing in the $x$ variable.

Note that for any $\xi$ we have

[TABLE]

Moreover, as $\varphi$ solves (12), we have

[TABLE]

Differentiating (14) with respect to $t$ we obtain

[TABLE]

which means that $V(x,t)$ solves equation (10). Since $\partial_{x}V>0$, the function $V$ solves the HJB equation (7), as claimed. $\diamondsuit$

Notice that the system of parabolic-ordinary differential equations (12)–(16) can also be rewritten as a system of two quasi-linear parabolic equations. Indeed, let us denote

[TABLE]

Then, by (14) we have $\psi(x,\tau)=1/\partial_{x}V(x,t)$. With regard to (10) we obtain

[TABLE]

As $\partial_{x}\psi=\varphi\,\psi$, we obtain $\partial_{x}\varphi\,\psi=\partial_{x}^{2}\psi-\varphi\,\partial_{x}\psi$. Since $\partial_{x}\alpha=\alpha^{\prime}_{\varphi}\,\partial_{x}\varphi+\alpha^{\prime}_{x}$, we conclude that the function $\psi$ satisfies the following parabolic equation:

[TABLE]

The terminal condition $\psi(x,T)$ can be deduced from the terminal condition $V(x,T)=u(x)$. That is, $\psi$ satisfies the initial condition $\psi(x,0)=1/u^{\prime}(x)$. In summary, we have shown the following theorem:

Theorem 3**.**

Assume that the terminal utility function $u(x)$ and the intertemporal function $c(x,t)$ are $C^{2}$ smooth functions and such that $u$ is increasing and $c$ is non-decreasing in the $x$ variable. Then an increasing function $V(x,t)$ in the $x$ variable is a solution to the Hamilton-Jacobi-Bellman equation (7) if and only if the pair $(\varphi,\psi)$ of transformed functions $\varphi(x,\tau)=-\partial_{x}^{2}V(x,t)/\partial_{x}V(x,t)$ and $\psi(x,\tau)=1/\partial_{x}V(x,t),t=T-\tau$, is a solution to the system of quasi-linear parabolic PDEs:

[TABLE]

and the value function $V(x,t)$ is given by (14).

At this point, we would like to emphasize the advantage of the suggested approach. By defining $\alpha$ as in (11) and subsequently setting up the PDEs in (12) or (18)–(19), one can compute the function $\alpha$ beforehand and then plug it into the corresponding PDEs as a known function. In this way, we do not have to deal with the maximization operator from (7) in each $x$ and $t$ separately which significantly simplifies the computation process.

Next we derive a-priori bounds on a solution $\varphi(x,\tau)$ of (12). We will use parabolic comparison principle (cf. [27]). To this end, we need to restrict the form of the value function $\alpha$ and utility functions $u,c$ by the following assumptions:

1. A1.

The value function $\alpha(x,\tau,\varphi)$ is separable in the following sense:

[TABLE]

where $\tilde{\alpha}$ is a $C^{1,1}$ smooth strictly increasing function with a bounded and Lipschitz continuous derivative for $\varphi\in(\varphi_{min},\infty)$ and $\alpha_{0}$ is a $C^{2}$ smooth function in $x\in\mathbb{R}$ and $\tau\in[0,T]$ variables. 2. A2.

There exist constants $\underline{\varphi},\overline{\varphi}\in\mathbb{R}$ such that $\varphi_{min}\leq\underline{\varphi}\leq 0\leq\overline{\varphi}$, and the functions $\alpha_{0}$ and $c$ satisfy the estimates:

[TABLE]

[TABLE]

for any $x\in\mathbb{R}$ and $\tau\in[0,T]$.

Example 3**.**

If $\alpha(x,\tau,\varphi)=\tilde{\alpha}(\varphi)-\varepsilon(\tau)e^{-x}-r(\tau)$ with $\varepsilon\geq 0$ is the value function introduced in Section 2 then $\tilde{\alpha}$ is defined on $[-1,\infty)$ and $\alpha_{0}(x,\tau)=-\varepsilon(\tau)e^{-x}-r(\tau)$ satisfies the assumption (A2) with $\underline{\varphi}=-1$ and any $\overline{\varphi}\geq 0$.

The intertemporal utility function $c$ of the form $c(x,t)=-\kappa e^{-dx-\varrho(T-t)}$ with $\kappa,d\geq 0$ satisfies (A2) if $\overline{\varphi}\geq d$.

Theorem 4**.**

Assume that the utility function $u(x)$ is a $C^{2}$ smooth strictly increasing function for $x\in\mathbb{R}$. Assume that the value function $\alpha$ and intertemporal utility function $c$ satisfy Assumption (A) with constants $\underline{\varphi}\leq 0\leq\overline{\varphi}$.

If the utility function $u$ satisfies the inequalities $\underline{\varphi}\leq\varphi(x,0)=-u^{\prime\prime}(x)/u^{\prime}(x)\leq\overline{\varphi}$ for any $x\in\mathbb{R}$, then, for the bounded solution $\varphi$ to (12) we have a-priori estimate: $\underline{\varphi}\leq\varphi(x,\tau)\leq\overline{\varphi}$ for any $\tau\in[0,T)$ and $x\in\mathbb{R}$.

P r o o f. Let $\psi(x,\tau)$ be a $C^{2}$ smooth nonnegative function, $\psi(x,\tau)\geq 0$. Let us define the parabolic operator:

[TABLE]

Since $\alpha(x,\tau,\varphi)=\tilde{\alpha}(\varphi)+\alpha_{0}(x,\tau)$, we have ${\mathcal{L}}(\tilde{\varphi})=\partial_{x}^{2}\alpha_{0}-\tilde{\varphi}\partial_{x}\alpha_{0}-\tilde{\varphi}\,\psi\,\partial_{x}c$ for any constant function $\tilde{\varphi}\in\mathbb{R}$. Thus, for constant functions $\underline{\varphi},\overline{\varphi}$ and nondecreasing function $c$, and $\alpha_{0}$ satisfying assumption (A2) we have

[TABLE]

Now let $\varphi$ be a solution to (12). Then ${\mathcal{L}}(\varphi)=\psi\partial^{2}_{x}c$ where $\psi(x,\tau)=\frac{e^{\int_{x_{*}}^{x}\varphi(\xi,\tau)d\xi}}{b(T-\tau)}>0$. Note the fact that $\partial_{x}\psi=\varphi\,\psi$. Hence the bounded solution $\varphi$ satisfies the following inequalities: ${\mathcal{L}}(\underline{\varphi})\geq{\mathcal{L}}(\varphi)\geq{\mathcal{L}}(\overline{\varphi})$.

If the initial condition satisfies the inequalities $\underline{\varphi}\leq\varphi(x,0)=-u^{\prime\prime}(x)/u^{\prime}(x)\leq\overline{\varphi}$ for any $x\in\mathbb{R}$ then, applying the parabolic comparison principle (cf. [27]), we obtain that the bounded solution $\varphi$ to (12) satisfies the inequalities $\underline{\varphi}\leq\varphi(x,\tau)\leq\overline{\varphi}$ for any $\tau\in[0,T]$ and $x\in\mathbb{R}$, as claimed. $\diamondsuit$

By $H^{k+\lambda}(\Omega)$, $0<\lambda<1$, we denote the Banach space consisting of all uniformly continuous functions $\varphi$ on $\bar{\Omega}=[x_{L},x_{R}]$ whose $k$-th derivative is uniformly $\lambda$-Hölder continuous, i.e. the Hölder semi-norm $\langle\varphi\rangle^{(\lambda)}_{k}=\sup_{x,y\in\Omega,x\not=y}|\partial^{k}_{x}\varphi(x)-\partial^{k}_{x}\varphi(y)|/|x-y|^{\lambda}$ is finite. Let $Q_{T}=\Omega\times(0,T)$ be a bounded cylinder. Following Ladyzhenskaya et al. [20] we introduce the parabolic Hölder space $H^{2k+\lambda,k+\lambda/2}(Q_{T})$ consisting of all continuous functions $\varphi:Q_{T}\to\mathbb{R}$ such that functions $\partial^{k}_{\tau}\varphi,\partial^{2k}_{x}\varphi$ are $\lambda$-Hölder continuous in the $x$-variable and $\lambda/2$-Hölder continuous in the $\tau$-variable.

Theorem 5**.**

Let $\Omega=(x_{L},x_{R})$ be a bounded interval. Assume $\alpha(x,\tau,\varphi)$ is $C^{2}$ smooth in the $t\in[0,T]$ and $x\in\Omega$ variables, $C^{1,1}$ smooth in the $\varphi$ variable, and such that $0<\alpha^{\prime}_{-}\leq\alpha^{\prime}_{\varphi}(x,\tau,\varphi)\leq\alpha^{\prime}_{+}<\infty$ for any $x\in\Omega,\tau\in[0,T]$, and $\varphi\geq\varphi_{min}$. Assume $c\in H^{2+\lambda,1+\lambda/2}(\Omega_{T})$ for some $0<\lambda<1/2$, and $c(x,t)$ is a non-decreasing function in the $x$ variable.

If the initial conditions $\varphi(\cdot,0),\psi(\cdot,0)\in H^{2+\lambda}(\Omega)$, then there exists a classical solution $(\varphi,\psi)$ to the system of quasi-linear parabolic equations (18)–(20) satisfying prescribed Dirichlet boundary conditions at $x_{L},x_{R}$. Moreover, $\eta,\psi\in H^{2+\lambda,1+\lambda/2}(\Omega_{T})$ where $\eta(x,\tau)=\alpha(x,\tau,\varphi(x,\tau))$. The function $\tau\mapsto\partial_{\tau}\varphi(x,\tau)$ is $\lambda/2$-Hölder continuous for all $x\in\mathbb{R}$ whereas $x\mapsto\partial_{x}\varphi(x,\tau)$ is Lipschitz continuous for all $\tau\in[0,T]$.

Proof.

The proof is analogous to that of [12, Theorem 5.3] where we proved existence of Hölder smooth solutions in the case when $c\equiv 0$. The methodology of the proof is based on the Schauder type of estimates (c.f. [20]).

Since the diffusion function $\alpha$ need not be $C^{2+\lambda}$ smooth in the $\varphi$ variable we first rewrite the system (18)–(20) using the auxiliary function $\eta=\alpha(\cdot,\varphi)$. Then $\varphi=\beta(\cdot,\eta)$ where $\beta$ is the inverse function to the strictly increasing function $\alpha$, i.e. $\varphi=\beta(\cdot,\alpha(\cdot,\varphi))$. Moreover, $\beta^{\prime}_{\eta}=1/\alpha^{\prime}_{\varphi}$ and $\partial_{\tau}\varphi=\beta^{\prime}_{\eta}\partial_{\tau}\eta+\beta^{\prime}_{\tau}$. Then system (18)–(20) can be rewritten in the form:

[TABLE]

where $\Phi=(\eta,\psi)$, $\delta(\cdot,\eta)=\alpha^{\prime}_{\varphi}(\cdot,\beta(\cdot,\eta))=1/\beta^{\prime}(\cdot,\eta)$, and

[TABLE]

Note that there are constants $\alpha^{\prime}_{\pm}$ such that $0<\alpha^{\prime}_{-}\leq\delta\leq\alpha^{\prime}_{+}<\infty$. Applying a $C^{2+\lambda}$ regularization of the function $\eta\mapsto\delta(\cdot,\eta)$ and following the proof of [12, Theorem 5.3] and the result on existence of classical solutions to the regularized equation (c.f. [20, Ch. V, pp. 495-496]) we conclude existence of a weak solution $\Phi\in W^{2,1}_{2}(Q_{T})$ of (21) satisfying the prescribed Dirichlet boundary and initial conditions. Recall that the parabolic Sobolev space $W^{2,1}_{2}(Q_{T})$ consists of all square integrable functions $\Phi\in L_{2}(Q_{T})$ such that $\partial_{x}\Phi,\partial^{2}_{x}\Phi,\partial_{\tau}\Phi\in L_{2}(Q_{T})$. The space $W^{2,1}_{2}(Q_{T})$ is continuously embedded into the Hölder space $H^{\lambda,\lambda/2}(Q_{T})$ for $0<\lambda<1/2$ (c.f. [20]). The rest of the proof is based on a simple bootstrap argument. A weak solution $\Phi\in W^{2,1}_{2}(Q_{T})$ is a solution to the linear equation:

[TABLE]

where the diffusion coefficient $\tilde{\delta}(\cdot)=\delta(.,\eta(\cdot))$, and $2\times 2$ matrices $\tilde{B}_{1},\tilde{B}_{0}$ belong to $H^{\lambda,\lambda/2}(Q_{T})$ because $\eta,\psi\in H^{\lambda,\lambda/2}(Q_{T})$ and $\delta$ is a Lipschitz continuous function. According to [20, Theorem 12.2, Chapter III] we have $(\eta,\psi)\equiv\Phi\in H^{2+\lambda,1+\lambda/2}(Q_{T})$ where $\eta=\alpha(\cdot,\phi)$. The proof of theorem now follows. ∎

5. Numerical approximation scheme

A numerical scheme that we propose for solving quasi-linear parabolic equation (12) is based on a semi-implicit in time approximation method. Spatial discretization is based on a finite volume approximation scheme (cf. LeVeque [21]) combined with a nonlinear equation iterative solver method proposed by Mikula and Kútik in [19]. This methodology for solving the transformed HJB equation was proposed and analyzed in [12] and [13] for the case of absence of an intertemporal utility function, i.e. $c=0$. In such a case, analysis of the experimental order of convergence suggested the second order of convergence with respect to the spatial discretization step (see [19], [12]).

Equation (12) belongs to a wide class of quasi-linear parabolic equations of the general form:

[TABLE]

satisfying the initial condition $\varphi(x,0)=-u^{\prime\prime}(x)/u^{\prime}(x)$ where $x\in\mathbb{R},\tau\in(0,T)$. Here $A(x,\tau,\varphi)=\alpha(x,\tau,\varphi),\quad B(x,\tau,\varphi)=-\alpha(x,\tau,\varphi)\varphi$, and

[TABLE]

5.1. A semi implicit time-space discretization of the transformed non-local parabolic equation

Since the original spatial domain for the $x$ variable is unbounded, we first truncate it into a bounded computational domain $[x_{L},x_{R}]$ and we use uniform spatial discretization mesh points $x_{i}=x_{L}+ih$ for $i=0,\cdots,n+1$ where $h=(x_{R}-x_{L})/(n+1)$. Thus $x_{0}=x_{L}$ and $x_{n+1}=x_{R}$. Following the idea of dual finite volumes (cf. [21]), the inner mesh points $x_{i}$, $i=1,\cdots,n$, are the centers of dual finite volumes $(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}})$. In what follows, the dual volumes will be denoted by $(x_{i-},x_{i+})$, i.e. $x_{i\pm}=x_{i\pm\frac{1}{2}}$. Clearly, $h=x_{i+}-x_{i-}$. The time discretization levels are set to $\tau^{j}=jk,j=0,\cdots,m$, where $k=T/m$ and $m$ is the number of time discretization steps. If we integrate equation (12) over dual finite volumes, apply the midpoint rule on the left-hand side integral and approximate the time derivative by the Euler forward finite difference, we arrive at the following system of equations:

[TABLE]

where

[TABLE]

and the integral $\int_{x_{i-}}^{x_{i+}}C{\hbox{d}}x$ over the interval $(x_{i-},x_{i+})$ is approximated by means of the mid-point rule integration, that is

[TABLE]

Let us denote

[TABLE]

and approximate the derivatives $\partial_{x}\varphi$ at dual mesh points $x_{i\pm}$ by the central differences:

[TABLE]

Let us fix a point $x_{*}=x_{i_{*}}$ for some spatial index $i_{*}$. As for the integral $\int_{x_{*}}^{x}\varphi(\eta,\tau)d\eta$ appearing in the non-local term $J_{i}$ at time layer $j$ (denoted as $J_{i}^{j}$) we use the trapezoidal integration rule:

[TABLE]

where $\Phi^{j}_{i}=\frac{h}{2}(\varphi^{j}(x_{L})+2\varphi^{j}(x_{1})+\dots+2\varphi^{j}(x_{i-1})+\varphi^{j}(x_{i}))$, which can be efficiently calculated recursively as follows:

[TABLE]

Hence

[TABLE]

Here $b^{j}$ is a discrete explicit/implicit Euler approximation of the solution $b(T-\tau^{j})$ to the ODE (16): $-db/d\tau=\omega b-\partial_{x}c(x_{*},T-\tau)$, i.e.

[TABLE]

when treated explicitly, or

[TABLE]

when treated implicitly. Here $\omega^{j}=\omega(T-\tau^{j})=(\partial_{x}\alpha-\alpha\varphi)|_{x=x_{*}}$ can be approximated by

[TABLE]

To compute a solution at the new time layer $\tau^{j+1}$, we take the terms $D_{i\pm}^{j},E_{i\pm}^{j},F_{i\pm}^{j}$ from the previous time layer $\tau^{j}$ and the terms $\partial_{x}\varphi|_{i\pm}^{j+1}$ from the new layer $\tau^{j+1}$. Rearranging the new-layer terms to the left-hand side and the old-layer terms to the right-hand side, we obtain a tridiagonal system of linear algebraic equations:

[TABLE]

which can be efficiently and fastly solved by the Thomas algorithm.

We assume Neumann boundary conditions at the boundaries $x_{L},x_{R}$. More precisely, $\partial_{x}\varphi(x,\tau)=0\text{ at }x=x_{L},x_{R}$, for all $\tau\in(0,T]$. The boundary conditions can be deduced from the asymptotic behavior of equation (12) for $x\to\pm\infty$. After discretization, these boundary conditions attain the form:

[TABLE]

5.2. Comparison with policy iteration method for solving HJB equations

In this section we discuss comparison of the numerical approximation scheme (30) and the fixed policy iteration method for solving HJB equation investigated by Huang et al. [8] and Reisinger and Witte [28].

Denote $V^{j}=V(\cdot,T-\tau^{j}),c^{j}=c(\cdot,T-\tau^{j})$. Then the time implicit time discretization of the HJB equation (7) can be written as follows:

[TABLE]

for $j=1,\cdots,m$. That is,

[TABLE]

where

[TABLE]

The fixed policy iteration method consists of replacing ${\bm{\theta}}^{j}$ by ${\bm{\theta}}^{j-1}$ in (32) and solving a linear equation for $V^{j}$, i.e.

[TABLE]

Since

[TABLE]

the fixed policy iteration method for solving HJB equation (10) corresponds to the numerical solution of the transformed equation (12) by means of the semi-implicit scheme (30) in which $\alpha$ is approximated by its linearization at $\varphi^{j-1}$ from the previous time step $\tau^{j-1}$.

Remark 2**.**

The main advantage of our method is twofold. First, we work with a transformed function $\varphi$ representing risk aversion of the investor. The boundary conditions for a truncated domain can be set up in a natural way, e.g. homogeneous Neumann boundary conditions. For the original problem formulated in terms of the intertemporal value function $V$, one can expect unbounded exponential like solution and numerical problems when working small and large values of $V$ and treatment of boundary conditions for $V$.

Secondly, the advantage consists of the possibility of evaluation of the value function $\alpha$ in a fast and efficient way instead of computation of ${\bm{\theta}}^{j}$ in (33). This might be useful when treating problems involving convex conic programming, e.g. worst-case portfolio selection problem, for which one can use efficient tools for solving convex conic programming optimization problems (cf. [14]).

6. Computational results and conclusions

6.1. Numerical benchmark to a traveling wave solution

Suppose that the value function $\alpha$ depends only on the $\varphi$ variable (e.g. we set $\varepsilon=0$, $r=0$ in (2)–(3)) and the intertemporal utility function is given as follows:

[TABLE]

Note that $c(x,T-\tau)=W(x-v\tau)$. Here $v\in\mathbb{R}$ is a given constant traveling wave speed. Then the function $V(x,t)=u(x-v(T-t))$ satisfies the equation:

[TABLE]

i.e. $V(x,t)$ is a solution to the HJB equation (10) and $V(x,T)=u(x)$. Hence the function

[TABLE]

is a traveling wave solution to (12) satisfying the initial condition $\varphi(x,0)=-u^{\prime\prime}(x)/u^{\prime}(x)$. The explicit solution of the form (34) can be used to test our numerical approximation scheme. As a testing example one can consider utility and value functions of the form:

[TABLE]

Then $u$ represents a convex-concave utility function with variable absolute risk aversion $a(x)$ given by

[TABLE]

If we set $c(x,t)=W(x-v(T-t))$ then $\varphi(x,\tau)=a(x-v\tau)=-u^{\prime\prime}(x-v\tau)/u^{\prime}(x-v\tau)$ is a solution to (12) satisfying the initial condition $\varphi(x,0)=a(x)=2x/(1+x^{2})$. Consequently, $V(x,t)=u(x-v(T-t))$ is the traveling wave solution to the HJB equation (10). The traveling wave solution $\varphi$ is depicted in Fig. 1 (left) for times $\tau^{j}=jT/10,j=0,\cdots,10$, where $T=1$ and $v=5$. As for a numerical solution, we considered the truncated computational domain $[x_{L},x_{R}]=[-20,20]$ and Dirichlet boundary conditions $\varphi(x_{L},\tau)=a(x_{L}-v\tau),\ \varphi(x_{R},\tau)=a(x_{R}-v\tau)$, for all $\tau>0$, which coincide with exact values of the explicit solution.

Let $\varphi_{expl}$ be the explicit traveling wave solution and $\varphi_{num}$ be the numerical solution constructed by means of our approximation scheme presented in Section 5. The $L_{2}$ and $L_{\infty}$ discrete norms are defined as follows:

[TABLE]

and the error between solutions as

[TABLE]

We consider the following relation between spatial and time discretization steps: $k=h^{2}$. Supposing $error(h)=O(h^{\delta})$, estimation of the order parameter $\delta$ can be obtained by means of the experimental order of convergence (or convergence ratio). It can be defined with respect to the norm of the space $L_{\infty}((0,T):L_{2}(x_{L},x_{R}))$ as follows:

[TABLE]

where $h_{1}>h_{2}>\cdots>h_{J}$. The results of computation of EOCs are summarized in Table 1.

The numerical results indicate the second order of convergence of the proposed numerical method. This is in accordance with the earlier result of the authors in [12] where we showed the same order of experimental convergence for the special case when $c=0$.

An example of the difference of explicit and numerical solution $\varphi_{expl}-\varphi_{num}$ is depicted in Fig. 1 (right). We can observe that the error is largest where the function $\varphi$ is steep.

6.2. Dynamic portfolio optimization example

Now we illustrate the solution of the proposed scheme on an example of dynamic portfolio optimization. Following Kilianová and Ševčovič [12, 13], we consider a stochastic dynamic portfolio optimization problem for a portfolio consisting of 30 stocks forming the German DAX30 stock index from August 2010 to April 2012. We chose the same data set as in [12, 13] for the purpose of comparison. As for the drift and volatility function, we will assume their form:

[TABLE]

where ${\bm{\Sigma}}$ is a positive definite covariance matrix. The function $\alpha(x,\tau,\varphi)$ can be rewritten as follows: $\alpha(x,\tau,\varphi)=\tilde{\alpha}(\varphi)-\varepsilon e^{-x}$, where $\tilde{\alpha}$ is the value function of the parametric quadratic optimization problem

[TABLE]

A graphical example of the function $\tilde{\alpha}$ in which ${\bm{\mu}}$ and ${\bm{\Sigma}}$ were obtained from the DAX30 data set is depicted in Figure 2. We can observe jumps in the graph of the second derivative of $\alpha$. Indeed, according to Theorem 1, the function $\alpha$ is $C^{1,1}$ continuous only. Furthermore, jumps in $\alpha^{\prime\prime}$ correspond to the points $\varphi$ where the set of indices $\{i:{\bm{\theta}}_{i}>0\}$ with positive weights is enlarged by a new index (cf. [12, 13]).

As for the utility functions, we use

[TABLE]

for the terminal utility and

[TABLE]

for the intertemporal utility.

Utility functions parameters used in our example are: $a=9$, $\kappa=1$, $d\in\{0,8,11\}$, $\varrho=0$. Parameters corresponding to model data are: $\varepsilon=1$. Parameters of the numerical scheme are: $h=0.01$, $k=0.5h^{2}$, $x_{L}=-4$, $x_{R}=8$, $x_{*}=x_{200}=-2.01,i_{*}=200$. Note that the solution does not depend on $x_{*}$ so one can choose arbitrary $x_{*}$. Nevertheless, a suitable choice of $x_{*}$ is very important in order to stabilize numerical computation. The free parameter $x_{*}$ enters integral term in (12) as well as the ODE for $b$, i.e. (16). We calculate the function $\alpha$ for $\varphi\in(-1,15)$ with a fine division step $h_{\varphi}=0.05$. The investment period is $T=1$.

Figure 3 presents numerical results for $d=0$ (the case of a trivial intertemporal function $c\equiv 0$), $d=8$ and $d=11$. The main difference we can observe is that while for the problem without intertemporal utility function we obtain a solution $\varphi(x,\tau)$ which is increasing on the interval $[x_{L},x_{R}]$, in problems with a nontrivial intertemporal utility function the solution $\varphi(x,\tau)$ turns out to be non-monotone in $x$. It eventually becomes increasing in the $x$ variable when $\tau$ is approaching the maturity $T$. Furthermore, the range of values of $\varphi$ is a smaller interval when compared to the case without intertemporal utility function. This has a practical consequence: as $\varphi$ has a small variation in the $x$ variable for $d\approx a$, so does the vector of optimal weights ${\bm{\theta}}$ (see Figure 3, right column). Notice that in the case when $\varepsilon=0,c\equiv 0$ there is a constant solution $\varphi(x,\tau)\equiv\varphi(x,0)=a$ to (12) corresponding to the so-called Merton solution to the optimal portfolio selection problem (cf. [24]). Note that the solution $\varphi$ satisfies a-priori estimates derived in Theorem 4.

In summary, there is a non-trivial effect on optimal portfolio selection when considering a non-trival intertemporal utility function $c(x,t)$ which has a similar behavior as the terminal utility function $u(x)$. We furthermore showed that the optimal solution $\varphi(x,\tau)$ to the transformed HJB equation (12) need not be monotonically increasing. In terms of the optimal portfolio selection vector ${\bm{\theta}}$, the optimal weight $\ theta_{i}$ for some of the stocks entering the active set can attain local minimum with respect to the $x$ variable. Such a behavior cannot be observed in models without intertemporal utility and non-trivial portfolio inflow $\varepsilon>0$ investigated in the recent papers [12, 13].

As far as numerical aspects of the Riccati transformation method are concerned, we showed that using this transformation enables to solve a quasi-linear parabolic equation by means of modern numerical methods based on finite volume approximation in a more efficient way when compared to traditional numerical methods based on the fixed policy iteration method or other explicit numerical approximation approaches the original HJB equation.

Acknowledgements

The authors were supported by VEGA 1/0062/18 and DAAD ENANEFA-2018 grants.