2. Formulation of problem

Throughout the paper, we assume that $(\Omega,\mathcal F, \mathbb P)$ be an atomless probability space. With a slight abuse of notation, we denote by $\mathcal M$ both the set of Borel probability measures on $\mathbb R$ and the set of distribution functions of real random variables. For $\Pi \in \mathcal M$ and $x \in \mathbb R$ , we have $\Pi(x)=\Pi((\!-\!\infty , x])$ . Let $\mathcal M ^p,\ p \in [1,\infty)$ , be the subset of $\mathcal M$ whose elements have finite p-th order moment. We write $X \sim \Pi$ if random variable X has distribution $\Pi$ , and $X \overset{d}= Y$ if X and Y have the same distribution. For $\Pi \in \mathcal M^2$ , we denote by $\mu(\Pi)$ and $\sigma ^2 (\Pi)$ the mean and variance of $\Pi$ , respectively.

2.1 Exploratory wealth process

In the financial market under study, we posit the inherent incompleteness of the market as in Basak and Chabakauri (Reference Basak and Chabakauri2010) and Dai et al. (Reference Dai, Jin, Kou and Xu2021), allowing for continuous trading in both a risk-free asset and a risky asset within the finite time horizon [0, T]. The price process $S_0(t)$ governing the risk-free asset is given by

\begin{align*} \mathrm{d} S_0(t) = rS_0(t)\mathrm{d} t,\end{align*}

where $r\gt0$ represents the constant risk-free interest rate. The price process S(t) of the risky asset is given by

(2.1) \begin{align} \dfrac{\mathrm{d} S(t)}{S(t)}=a(t,Y(t))\mathrm{d} t+b(t,Y(t))\mathrm{d} W(t).\end{align}

Here, W(t) is a Brownian motion defined on the filtered probability space $(\Omega, \mathcal F, \{\mathcal F_t\}_{0\leqslant t\leqslant T}, \mathbb P)$ adhering to standard conditions. Additionally, Y(t) is a diffusion process satisfying

\begin{align*} \mathrm{d} Y(t)=m(t,Y(t))\mathrm{d} t+v(t,Y(t))[\rho dW(t)+\sqrt{1-\rho ^2}d\widetilde W(t)],\end{align*}

where $\rho \in [\!-\!1,1]$ , and the Brownian motion $\widetilde W(t)$ is defined on $(\Omega, \mathcal F, \{\mathcal F_t\}_{0\leqslant t\leqslant T}, \mathbb P)$ independent of W(t). In the market context, $a(\cdot,\cdot)$ , $b(\cdot,\cdot)$ , $m(\cdot,\cdot)$ , and $v(\cdot,\cdot)$ are deterministic functions of t and y, with all randomness entering through Y(t). However, these functions themselves are unknown. Particularly, we are in the case of the complete market with constant parameters, where $a(t, y) \equiv a$ and $b(t, y) \equiv b$ . The process Y(t) is modeled as a diffusion representing macroeconomic or systemic risk factors. The parameter $\rho$ quantifies the correlation between the state process Y(t) and the stock price process S(t), thus characterizing the degree to which market-wide risk factors influence asset price dynamics.

Remark 1. Market incompleteness arises because the stochastic evolution of Y(t) affects the dynamics of the risky asset, yet cannot be fully hedged through trading in S(t) and $S_0(t)$ . This feature aligns with empirical observations, where financial markets exhibit unspanned stochastic volatility, unhedgeable economic risks, and other uncertainties. The market model described above exhibits substantial generality, encompassing numerous well-known models, such as the time-varying Gaussian mean return model and the stochastic volatility model; see Example 1 below. For a more detailed discussion on market incompleteness and the nonreplicability of contingent claims in such settings, we refer to Sections 1.4 and 5.4 of Föllmer and Schied (Reference Föllmer and Schied2011).

Example 1.

1. (i) If the stock price S(t) and the market price of risk Y(t) are governed by $\frac{\mathrm{d} S(t)}{S(t)} =(r+\sigma Y(t)) \mathrm{d} t+\sigma \mathrm{d} W(t)$ and $\mathrm{d} Y(t) =\iota \left(Y-Y(t)\right) \mathrm{d} t+v \mathrm{d}[\rho dW(t)+\sqrt{1-\rho ^2}d\widetilde W(t)],$ where $r, \sigma, \iota, v$ , Y are all positive constants and $\rho \in [\!-\!1,1]$ . This model, termed the time-varying Gaussian mean return model, effectively captures the intricate interplay between the dynamics of stock prices and market risk over time.

2. (ii) If the stock price S(t) and a state variable Y(t) follow $\frac{\mathrm{d} S(t)}{S(t)}=(r+\sigma Y(t)^{\frac{1+\alpha}{2 \alpha}}) \mathrm{d} t+Y(t)^{\frac{1}{2 \alpha}} \mathrm{d} W(t)$ and $\mathrm{d} Y(t)=\iota (Y-Y(t)) d t+v \sqrt{Y(t)} \mathrm{d} [\rho dW(t)+\sqrt{1-\rho ^2}d\widetilde W(t)],$ where $\alpha \neq 0$ is the constant elasticity of the market price of risk $\sigma Y(t)^{\frac{1}{2}}, \sigma \in \mathbb{R}, \iota\gt0, v\gt0$ , $Y \in \mathbb{R}$ and $\rho \in [\!-\!1,1]$ are all constants, then it is a stochastic volatility model. This formulation characterizes a stochastic volatility model, acknowledging the nonlinear relationship between the state variable and market risk elasticity.

Researchers, including Merton (Reference Merton1980), Kim and Omberg (Reference Kim and Omberg1996), Liu (Reference Liu2001), Wachter (Reference Wachter2002), Basak and Chabakauri (Reference Basak and Chabakauri2010), Dai et al. (Reference Dai, Jin, Kou and Xu2021), and Dai et al. (Reference Dai, Dong and Jia2023), have extensively explored dynamic portfolio choice problems within the framework of these two market settings and their specific variations.

In what follows, we consider a nonzero-sum game involving two competing agents or companies, referred as Agent 1 and Agent 2 for simplicity. Both agents have access to the financial market. For $i \in \{1,2\}$ , let $u_i(t)$ represent the discounted amount of Agent i invested in the risky asset at time t, and the rest of the wealth is invested in the risk-free asset. Define $\boldsymbol{u}_i=\{u_i(t),0\leqslant t\leqslant T\}$ and

(2.2) \begin{equation}\theta(t,Y(t))=\frac{a(t,Y(t))-r}{b(t,Y(t))}.\end{equation}

Then the dynamic of discounted wealth process of Agent i under strategy $u_i$ is given as

(2.3) \begin{align} \mathrm{d} X_i^{u_i}(t)&=u_i(t)b(t,Y(t))[\theta(t,Y(t))\mathrm{d} t+\mathrm{d} W(t)].\end{align}

Based on Wang et al. (Reference Wang, Zariphopoulou and Zhou2020a), we extend the control process (2.3) to a distributional control process and define the exploratory discounted wealth process for agents. Denote by $\Pi_i(t) \in \mathcal M^2$ the probability distribution function of the control $u_i$ at time t, where $\mathcal M^2$ represents the set of distribution functions with finite second-order moment on $\mathbb R$ . Let $\boldsymbol{\Pi}_i=\{\Pi_i(t),0\leqslant t\leqslant T\}$ and we write $\boldsymbol{\Pi}_i \in \mathcal M^2$ for simplicity. In contrast to the approach presented in Wang et al. (Reference Wang, Zariphopoulou and Zhou2020a), we additionally consider the correlation between the discounted wealth process and the market state, as well as the correlation between the two agents.

Now, we attempt to derive the exploratory version of the wealth process $X_i^{u_i}$ associated with randomized policy $\Pi_i$ . Such a structure can also be found in Dai et al. (Reference Dai, Dong and Jia2023), Wang et al. (Reference Wang, Zariphopoulou and Zhou2020a), Wang and Zhou (Reference Wang and Zhou2020). Let $W^n(t)$ and $\widetilde W^n(t),n=1,2,...,N,$ represent N paths independently sampled from W(t) and $\widetilde W(t)$ , respectively. Moreover, let $X_i^n(t)$ be the copies of the discounted wealth process of Agent i under strategy $u_i^n$ sampled from $\Pi_i$ . Then, for $n=1,2,...,N$ and $i \in \{1,2\}$ , the increments of $Y^n(t)$ and the corresponding $X_i^n(t)$ can be written as

(2.4) \begin{align}\begin{aligned} \Delta Y^n(t)&\equiv Y^n(t+\Delta t)-Y^n(t)\\&\approx m(t,Y^n(t))\Delta t+v(t,Y^n(t))[\rho \Delta W^n(t)+\sqrt{1-\rho ^2}\Delta \widetilde W^n(t)],\end{aligned}\end{align}

and

(2.5) \begin{align}\begin{aligned} \Delta X_i^n(t)&\equiv X_i^n(t+\Delta t)-X_i^n(t)\\&\approx u_i^n(t)b(t,Y^n(t))[\theta(t,Y^n(t))\Delta t+\Delta W^n(t)].\end{aligned}\end{align}

We denote the exploratory discounted wealth process of Agent i by $X_i^{\Pi_i}(t)$ . Consequently, $X_i^n(t)$ can be viewed as an independent sample from $X_i^{\Pi_i}(t)$ . By the law of large numbers and using (2.4) and (2.5), we observe that, as $N \to \infty$ ,

(2.6) \begin{align}\begin{split} &\dfrac{1}{N}\sum\limits_{n=1}^N \Delta X_i^n(t)\approx \dfrac{1}{N}\sum\limits_{n=1}^N u_i^n(t)b(t,Y^n(t))[\theta(t,Y^n(t))\Delta t+\Delta W^n(t)]\\ &\overset{a.s.}\longrightarrow \mathbb E\left[b(t,Y(t))\theta(t,Y(t))\int_U u\mathrm{d}\Pi_i(t,u)\right]\Delta t,\end{split}\end{align}

(2.7) \begin{align}\begin{split} &\dfrac{1}{N}\sum\limits_{n=1}^N (\Delta X_i^n(t))^2\approx \dfrac{1}{N}\sum\limits_{n=1}^N (u_i^n(t)b(t,Y^n(t)))^2\Delta t \\ &\overset{a.s.}\longrightarrow \mathbb E\left[b^2(t,Y(t))\int_U u^2\mathrm{d}\Pi_i(t,u)\right]\Delta t,\end{split}\end{align}

and

(2.8) \begin{align}\begin{split} &\dfrac{1}{N}\sum\limits_{n=1}^N \Delta X_i^n(t)\Delta Y^n(t)\approx \dfrac{1}{N}\sum\limits_{n=1}^N u_i^n(t)b(t,Y^n(t))\rho v(t,Y^n(t))\Delta t\\ &\overset{a.s.}\longrightarrow \mathbb E\left[\rho b(t,Y(t))v(t,Y(t))\int_U u\mathrm{d}\Pi_i(t,u)\right]\Delta t.\end{split}\end{align}

It is well known from the law of large numbers that

$$ \dfrac{1}{N}\sum\limits_{n=1}^N \Delta X_i^n(t)\overset{a.s.}\longrightarrow \mathbb E[\Delta X_i^{\Pi_i}(t)], \dfrac{1}{N}\sum\limits_{n=1}^N (\Delta X_i^n(t))^2\overset{a.s.}\longrightarrow \mathbb E[(\Delta X_i^{\Pi_i}(t))^2],$$

and

\begin{align*} \dfrac{1}{N}\sum\limits_{n=1}^N \Delta X_i^n(t)\Delta Y^n(t)\overset{a.s.}\longrightarrow \mathbb E[\Delta X_i^{\Pi_i}(t)\Delta Y(t)].\end{align*}

These together with (2.6)–(2.8), motivate our confidence in the dynamic of $\Delta X_i^{\Pi_i}(t)$ , which can be expressed as

(2.9) \begin{align}\begin{split} \mathrm{d} X_i^{\Pi_i}(t)&=b(t,Y(t))\theta(t,Y(t))\mu_i(t)\mathrm{d} t+b(t,Y(t))[\mu_i(t)\mathrm{d} W(t)+\sigma_i(t)\mathrm{d} \overline W_i(t)],\end{split}\end{align}

where

$$\mu_i(t)=\int_U u\mathrm{d}\Pi_i(t,u), \;\;\sigma^2_i(t)=\int_U u^2\mathrm{d}\Pi_i(t,u)-\mu^2_i(t)$$

with $\overline W_i$ being the Brownian motion independent of W(t) and $\widetilde W(t)$ . The remaining consideration involves the correlation between $\overline W_1(t)$ and $\overline W_2(t)$ . It is observed that, as $N\to \infty$ ,

\begin{align*} &\dfrac{1}{N}\sum\limits_{n=1}^N \Delta X_1^n(t)\Delta X_2^n(t)\approx \dfrac{1}{N}\sum\limits_{n=1}^N u_1^n(t)u_2^n(t)b^2(t,Y^n(t))\\ &\overset{a.s.}\longrightarrow \mathbb E\left[\mu_1(t)\mu_2(t)b^2(t,Y(t))\right]\Delta t,\end{align*}

leading to $\langle \overline W_1,\overline W_2 \rangle =0$ . By Lévy’s theorem, we conclude that $\overline W_1$ is independent of $\overline W_2$ .

Remark 2. Indeed, the above construction follows from the framework of relaxed stochastic control. In classical control theory, the strategy $u(\omega,t)$ is a stochastic process that takes a deterministic value for fixed $\omega$ and t. When we randomize the strategy $u(\omega,t)$ to $\Pi(\omega,t)$ , then for each fixed $\omega$ and t, $\Pi(\omega,t)$ is a distribution on U, which we write as $\Pi_t(\omega,\mathrm d u)$ . For any t, we can view $\Pi_t(\omega,B)$ with $B \in \mathcal B(U)$ , as a probability kernel with source $(\Omega,\mathcal F_t)$ and target $(U,\mathcal B(U))$ . By standard measure-theoretic results, $\Pi_t(\omega,\cdot)$ together with the probability $\mathbb P$ on $(\Omega,\mathcal F_t)$ uniquely induces a measure $\mathbb Q$ on $(\Omega\times U,\mathcal F_t\times \mathcal B(U))$ satisfying $ \mathbb Q(\mathrm d\omega,\mathrm d u)=\Pi_t(\omega,\mathrm d u)\mathbb P(\mathrm d\omega).$ The sampling described above is thus carried out on the probability space $(\Omega\times U,\mathcal F_t\times \mathcal B(U),\mathbb Q)$ . Therefore, by the generalized Fubini theorem, for any measurable function f on $(\Omega\times U,\mathcal F_t\times \mathcal B(U),\mathbb Q)$ , we have

\begin{align*} \mathbb E^{\mathbb Q}[f]=\int_{\Omega\times U}f\mathrm d \mathbb Q=\int_{\Omega}\int_U f(\omega,u)\Pi_t(\omega,\mathrm d u)\mathrm d \mathbb P=\mathbb E\left[\int_U f(\omega,u)\Pi_t(\omega,\mathrm d u)\right].\end{align*}

For $i,j\in\{1,2\}$ and $j\neq i$ , assume that Agent i takes into account not only his own wealth but also the wealth gap between himself and Agent j at the terminal time T. Given a strategy $\Pi_j\in \mathcal M^2$ employed by Agent j, Agent i will choose a strategy $\Pi_i$ to maximize the following objective

(2.10) \begin{align}\begin{split} &\;\mathbb E_t[(1-k_i)X_i^{\Pi_i}(T)+k_i(X_i^{\Pi_i}(T)-X_j^{\Pi_j}(T))]-\dfrac{\gamma_i}{2}\mathrm{Var}_t[(1-k_i)X_i^{\Pi_i}(T)+k_i(X_i^{\Pi_i}(T)-X_j^{\Pi_j}(T))]\\ &=\mathbb E_t[X_i^{\Pi_i}(T)-k_iX_j^{\Pi_j}(T)]-\dfrac{\gamma_i}{2}\mathrm{Var}_t[X_i^{\Pi_i}(T)-k_iX_j^{\Pi_j}(T)],\end{split}\end{align}

where $\mathbb E_t(\!\cdot\!)$ and $\mathrm{Var}_t(\!\cdot\!)$ denote the conditional expectation and variance with given $X_1^{\Pi_1}(t)$ , $X_2^{\Pi_2}(t)$ and Y(t), respectively, $\gamma_i\gt0$ represents the risk-aversion coefficient for Agent i, and $k_i\in (0,1)$ measures the sensitivity of Agent i to the performance of Agent j. Notably, a larger $k_i$ implies that Agent i emphasis on the relative performance against his opponent (Agent j), thereby intensifying the competitiveness of the game. Similar models addressing relative performance concerns, which do not incorporate RL, have been studied in the literature. For related discussions, see, for example, Bensoussan et al. (Reference Bensoussan, Siu, Yam and Yang2014), as well as the relevant references cited in our introduction.

Remark 3. Unlike the objective of maximizing the insurer’s utility in a nonzero-sum game (see, e.g., Browne, Reference Browne2000 and Bensoussan et al. Reference Bensoussan, Siu, Yam and Yang2014), we consider a mean-variance objective as in Wang et al. (Reference Wang, Zhang, Jin and Qian2019) and (Reference Wang, Zhang, Jin and Qian2021). On one hand, the agents seek to maximize the expected wealth gap; on the other hand, they prioritize stability and predictability over excessive volatility. By minimizing the variance of the wealth gap, the agent reduces uncertainty in their relative wealth position, aligning with standard assumptions of risk aversion.

Let $\hat X_i^{\Pi_i,\Pi_j}(t)=X_i^{\Pi_i}(t)-k_iX_j^{\Pi_j}(t)$ be the wealth difference of the two agents. It is obvious from (2.9) that $\hat X_i^{\Pi_i,\Pi_j}(t)$ follows the dynamic

(2.11) \begin{align}\begin{split} \mathrm{d} \hat X_i^{\Pi_i,\Pi_j}(t)&=b(t,Y(t))[\theta(t,Y(t))(\mu_i(t)-k_i\mu_j(t))\mathrm{d} t\\ &\quad +(\mu_i(t)-k_i\mu_j(t))\mathrm{d} W(t)+\sigma_i(t)\mathrm{d}\overline W_i(t)-k_i\sigma_j(t)\mathrm{d}\overline W_j(t)].\end{split}\end{align}

2.2 Objective function

We employ the Choquet regularizer $\Phi_h$ to quantify randomness. Given a concave function $h\;:\;[0,1]\to \mathbb{R}$ of bounded variation with $h(0)=h(1)=0$ and $\Pi\in \mathcal M$ , the Choquet regularizer $\Phi_h$ on $\mathcal M$ is defined as

(2.12) \begin{equation}\Phi_h(\Pi)= \int_\mathbb{R} h\circ \Pi([x,\infty))\mathrm{d} x.\end{equation}

We denote the set of $h\;:\;[0,1]\to \mathbb{R}$ by $\mathcal H$ . In fact, (2.12) is a signed Choquet integral characterized by Wang et al. (Reference Wang, Wei and Willmot2020c) via comonotonic additivity, which essentially builds on the seminal works of Schmeidler (Reference Schmeidler1989) and Yaari (Reference Yaari1987).

As stated in Lemma 2.2 of Han et al. (Reference Han, Wang and Zhou2023), the regularizer $\Phi_h$ is rigorously defined and serves as an important metric for quantifying the degree of randomness or exploration within the context of RL. Specifically, the concavity of h ensures that $\Phi_h$ is also concave. This means that $\Phi_h(\lambda \Pi_1 + (1-\lambda) \Pi_2 ) \geqslant \lambda \Phi_h( \Pi_1) + (1-\lambda) \Phi_h(\Pi_2 )$ for all $\Pi_1,\Pi_2\in \mathcal M$ and $\lambda\in [0,1],$ which intuitively implies that the linear combination of two distributions is more random. Moreover, the condition $h(0)=h(1)=0$ means for any $c\in\mathbb{R}$ , $\Phi_h(\delta_c)=0$ , where $\delta_c$ is the Dirac mass at c. This indicates that degenerate distributions do not have any randomness measured by $\Phi_h$ . Additionally, the agents have the flexibility to opt for different regularizers, contingent upon their preferences, as reflected by the distortion function h. For more detailed discussions about the properties associated with $\Phi_h$ , we refer to Han et al. (Reference Han, Wang and Zhou2023) and Guo et al. (Reference Guo, Han and Wang2025).

It is useful to note that $\Phi_h$ admits a quantile representation, see Lemma 1 of Wang et al. (Reference Wang, Wang and Wei2020b). For a distribution $\Pi\in\mathcal{M}$ , let its left-quantile for $p\in(0,1]$ be defined as

$$Q_\Pi(p)=\inf \left\{x\in \mathbb{R}\;:\; \Pi(x) \geqslant p\right\}, $$

then we have

(2.13) \begin{equation}\Phi_h(\Pi) =\int_0^1 Q_\Pi(1-p) \mathrm{d} h(p) \end{equation}

if h is left-continuous.

For any fixed $\Pi_j \in \mathcal M^2$ , we incorporate the Choquet regularizer $\Phi_{h_i}$ for Agent i, along with the exploration weight function $\lambda_i(t)$ , into the mean-variance criterion (2.10). This regularizer plays a role analogous to differential entropy in entropy-regularized reinforcement learning (e.g., Wang et al., Reference Wang, Zariphopoulou and Zhou2020a; Wang and Zhou, Reference Wang and Zhou2020), serving to prevent the learned policy from collapsing to a deterministic solution. Thus, each agent aims to achieve an exploratory mean-variance problem within the framework of RL. This yields the corresponding objective function

(2.14) \begin{align} J_i(t,\hat x_i,y;\;\boldsymbol{\Pi}_i,\boldsymbol{\Pi}_j)\;:\!=\;\mathbb E_t\left[\int_t^T\lambda_i(s)\Phi_{h_i}(\Pi_i(s))ds+\hat X_i^{\Pi_i,\Pi_j}(T)\right]-\dfrac{\gamma_i}{2}\mathrm{Var}_t[\hat X_i^{\Pi_i,\Pi_j}(T)],\end{align}

where $\hat X_i^{\Pi_i,\Pi_j}(t)=\hat x_i$ and $\ Y(t)=y$ with t representing the initial time.

Remark 4. We note that a larger $\lambda_i(t)$ promotes increased exploration, as it results in a higher weight on $\Phi_{h_i}(\Pi_i(t))$ . When $\lambda_i(t) \equiv \lambda_i$ , the exploration weight remains constant over time. It is often more realistic to set $\lambda_i(t)$ to decrease over time; for instance, $\lambda(t)$ may follow a power-decaying pattern, expressed as $\lambda_i(t) = \lambda_0 (T+\lambda)^{\lambda_0}/(t+\lambda)^{\lambda_0+1} $ with $\lambda_0, \lambda \gt 0$ . Alternatively, $\lambda_i(t)$ may decay exponentially with $\lambda_i(t) = \lambda_0 e^{\lambda_0(T-t)}$ with $\lambda_0 \gt 0$ . For further discussions on selecting $\lambda_i(t)$ , we refer to Section 3.4 of Jiang et al. (Reference Jiang, Saunders and Weng2022).

For $i,j\in\{1,2\}$ and $i\neq j$ , denote by $\mathcal A(\boldsymbol{\Pi}_j)$ the set of all admissible feedback control of $\boldsymbol{\Pi}_i$ . Given $\boldsymbol{\Pi}_j \in \mathcal M^2$ , $\boldsymbol{\Pi}_i$ is said to be in $\mathcal A(\boldsymbol{\Pi}_j)$ if the following conditions hold:

1. (i) For each $s\in[t,T]$ , $\Pi_i(s)\in \mathcal M^2(\mathbb R)$ ;

2. (ii) There exists a deterministic mapping $\xi_i\;:\;[t,T]\times\mathbb R \times\mathbb R\to \mathcal M^2(\mathbb R)$ , such that $\Pi_i(s)=\xi_i(s,\hat X_i^{\Pi_i,\Pi_j}(s),Y(s))$ ;

3. (iii) For any $A \in \mathcal B(\mathbb R)$ , $\{\int_A\Pi_i(s,u)\mathrm{d} u, t\leqslant s \leqslant T\} $ is $\mathcal F_s$ –progressively measurable;

4. (iv) $\mathbb E_t\int_t^T [|b(s,Y(s))\theta(s,Y(s))(\mu_i(s)-k_i\mu_j(s))|+b^2(s,Y(s))(\mu_i^2(s)+\sigma^2_i(s))]ds\lt \infty$ ;

5. (v) $\mathbb E_t\int_t^T|\lambda_i(s)\Phi_{h_i}(\Pi_i(s))|ds \lt\infty $ .

Next, we define the profile of the game, which encapsulates the comprehensive strategic behavior of both agents in the game.

Definition 1. A strategy pair $(\boldsymbol{\Pi}_1,\boldsymbol{\Pi}_2)$ is called the profile of the game if it satisfies the following conditions: given the strategy $\boldsymbol{\Pi}_2$ , $\boldsymbol{\Pi}_1$ is an admissible feedback control, and conversely, given the strategy $\boldsymbol{\Pi}_1$ , $\boldsymbol{\Pi}_2$ is also an admissible feedback control.

As mentioned in the introduction, to seek the time-consistent equilibrium strategy, we formulate the time-inconsistent dynamic optimization problem into a noncooperative game theoretic framework proposed by Björk and Murgoci (Reference Björk and Murgoci2014) and Björk et al. (Reference Björk, Khapko and Murgoci2017).

Definition 2. Let $(\boldsymbol{\Pi}_1,\boldsymbol{\Pi}_2)$ be a profile of the game. For $i,j\in\{1,2\}$ and $i\neq j$ , $\boldsymbol{\Pi}_i$ is said to be equilibrium response of Agent i if for a fixed $\Delta$ and for any initial state $(t,\hat x_i,y)$ and an arbitrary $\eta_i\in\mathcal M^2$ , $\boldsymbol{\Pi}_i^{\eta_i,\Delta}$ defined by

\begin{align*} \Pi_i^{\eta_i,\Delta}(s)=\left\{ \begin{aligned} &\eta_i,&t\leqslant s\leqslant t+\Delta,\\ &\Pi_i(s),&t+\Delta\leqslant s\leqslant T, \end{aligned} \right. \end{align*}

satisfying

\begin{align*} \liminf\limits_{\Delta\rightarrow 0^+}\dfrac{J_i(t,\hat x_i,y;\;\boldsymbol{\Pi}_i,\boldsymbol{\Pi}_j)-J_i(t,\hat x_i,y;\;\boldsymbol{\Pi}_i^{\eta_i,\Delta},\boldsymbol{\Pi}_j)}{\Delta}\geqslant 0. \end{align*}

The equilibrium response $\boldsymbol{\Pi}_i$ of Agent i can be viewed as a mapping of $\boldsymbol{\Pi}_j$ , and thus can be written as $\boldsymbol{\Pi}_i=\kappa_i(\boldsymbol{\Pi}_j)$ at this point. Furthermore, the equilibrium response value function of Agent i is defined as

\begin{align*} \widetilde V_i(t,\hat x_i,y;\;\boldsymbol{\Pi}_j)\;:\!=\;J_i(t,\hat x_i,y;\;\boldsymbol{\Pi}_i,\boldsymbol{\Pi}_j). \end{align*}

3. Time-consistent Nash equilibrium

In this section, we present the Nash equilibrium for a general incomplete market. It is important to emphasize that the uniqueness of the equilibrium policy remains an open question in the context of time-inconsistent optimization problems, as discussed by Ekeland and Pirvu (Reference Ekeland and Pirvu2008). Therefore, within the framework of game theory, we focus on a specific Nash equilibrium defined below.

Definition 3. A profile $(\boldsymbol{\Pi}_1^*,\boldsymbol{\Pi}_2^*)$ is called the time-consistent Nash equilibrium of the game if for $i,j\in\{1,2\}$ and $i\neq j$ , $\boldsymbol{\Pi}_i^*$ is the equilibrium response of Agent i, that is $\boldsymbol{\Pi}_1^*=\kappa_1(\boldsymbol{\Pi}_2^*)$ and $\boldsymbol{\Pi}_2^*=\kappa_2(\boldsymbol{\Pi}_1^*)$ . Furthermore, the equilibrium value function of Agent i is given by

\begin{align*} V_i(t,\hat x_i,y)\;:\!=\;\widetilde V_i(t,\hat x_i,y;\;\boldsymbol{\Pi}^*_j)=J_i(t,\hat x_i,y;\;\boldsymbol{\Pi}_i^*,\boldsymbol{\Pi}_j^*). \end{align*}

3.1 Verification theorem

Analogous to the work of Björk et al. (Reference Björk, Khapko and Murgoci2017) on optimization problems involving a broad class of objective functionals, we present the following verification theorem.

Let $\mathcal{D}=[0,T]\times\mathbb{R}^2$ and then for any $\varphi\in C^{1,2,2}(\mathcal{D})$ , we can denote the infinitesimal generator of $(\hat X^{\Pi_i,\Pi_j}(t),Y(t))$ by

(3.1) \begin{equation}\begin{aligned} \mathcal L^{\Pi_i,\Pi_j}\varphi(t, x, y)=&\displaystyle\frac{\partial\varphi}{\partial t}+b(t,y)\theta(t,y)(\mu_i(t)-k_i\mu_j(t)) \frac{\partial\varphi}{\partial x} +\frac{1}{2}b^2(t,y)\left((\mu_i(t)-k_i\mu_j(t))^2\right.\\&\left.+\sigma^2_i(t)+k^2_i\sigma^2_j(t)\right)\frac{\partial^2\varphi}{\partial x^2} + m(t,y)\frac{\partial\varphi}{\partial y}+\frac{1}{2}v^2(t,y) \frac{\partial^2\varphi}{\partial y^2}\\&+ \rho v(t,y) b(t,y)(\mu_i(t)-k_i \mu_j(t))\frac{\partial^2\varphi}{\partial x \partial y }.\end{aligned}\end{equation}

Let $\hat X_i^{u_i,u_j}(t)=X_i^{u_i}(t)-k_iX_j^{u_j}(t)$ and denote the infinitesimal generator of $(\hat X_i^{u_i,u_j}(t),Y(t))$ by $\mathcal L^{u_i,u_j}$ . It can be directly verified that $\mathcal L^{\Pi_i,\Pi_j}\varphi(t, x, y)=\int_{\mathbb R}\int_{\mathbb R}\mathcal L^{u_i,u_j}\varphi(t, x, y)\mathrm{d}\Pi_i(u_i)\mathrm{d}\Pi_j(u_j)$ .

Theorem 1 (Verification theorem). For $i,j\in\{1,2\}$ and $i\neq j$ , fix $\boldsymbol{\Pi}_j$ and suppose that functions $V_i(t,x,y)\in C^{1,2,2}(\mathcal{D})$ , $g_i(t,x,y)\in C^{1,2,2}(\mathcal{D})$ and strategy $\boldsymbol{\Pi}_i$ satisfy the following properties:

1. (i) $V_i$ and $g_i$ solve the extended HJB system

   (3.2) \begin{align} &\sup\limits_{\Pi_i\in \mathcal M^2}\left\{\mathcal L^{\Pi_i,\Pi_j}V_i(t,\hat x_i,y)-\dfrac{\gamma_i}{2}\mathcal L^{\Pi_i,\Pi_j}g_i^2(t,\hat x_i,y)+\gamma_ig_i\mathcal L^{\Pi_i,\Pi_j}g_i(t,\hat x_i,y)+\lambda_i(t)\Phi_{h_i}(\Pi_i) \right\}=0, \end{align}
   and
   (3.3) \begin{align} \mathcal L^{\Pi_i^*,\Pi_j}g_i(t,\hat x_i,y)=0 \end{align}
   with
   (3.4) \begin{align} V_i(T,\hat x_i,y)=\hat x_i,\ g_i(T,\hat x_i,y)=\hat x_i, \end{align}
   where $\mathcal L^{\Pi_i,\Pi_j}$ is the infinitesimal generator given by (3.1).

2. (ii) $\Pi_i^*$ realizes the supremum in (3.2) and $\boldsymbol{\Pi}^*_i=\{\Pi^*_i(t),0\leqslant t\leqslant T\}$ is admissible.

Then $\boldsymbol{\Pi}^*_i$ is the equilibrium response of Agent i. Furthermore, $V_i(t,\hat x_i,y)=J_i(t,\hat x_i,y;\;\boldsymbol{\Pi}^*_i,\boldsymbol{\Pi}_j)$ is the equilibrium response value function of Agent i and $g_i(t,\hat x_i,y)=\mathbb E_t[\hat X_i^{\Pi^*_i,\Pi_j}(T)]$ .

3.2 Solution to the general case

While the preservation of the uniqueness of equilibrium response remains uncertain, the aforementioned theorem provides a constructive framework for identifying a specific Nash equilibrium. We first concentrate on the equilibrium response of Agent 1, and analogous results for Agent 2 can be obtained using the same method. It is assumed that the equilibrium value function $V_1$ and the function $g_1$ can be precisely expressed as

(3.5) \begin{align} V_1(t,x,y)=x+D_1(t,y),\;\;\;g_1(t,x,y)=x+d_1(t,y).\end{align}

It is obvious that $D_1(T,y)=d_1(T,y)=0$ . Using (3.1) and substituting (3.5) into (3.2), we streamline (3.2) to

(3.6) \begin{align}\begin{split} \sup\limits_{\Pi_1\in \mathcal M^2}&\left\{ \dfrac{\partial D_1(t,y)}{\partial t}+b(t,y)\theta(t,y)(\mu_1(t)-k_1\mu_2(t))+m(t,y)\dfrac{\partial D_1(t,y)}{\partial y}\right.\\ &-\dfrac{\gamma_1}{2}b^2(t,y)[(\mu_1^2(t)+\sigma_1^2(t))+k_1^2(\mu_2^2(t)+\sigma_2^2(t))-2k_1\mu_1(t)\mu_2(t)]\\ &-\dfrac{\gamma_1}{2}v^2(t,y)\left(\dfrac{\partial d_1(t,y)}{\partial y} \right)^2+\dfrac{1}{2}v^2(t,y)\dfrac{\partial ^2 D_1(t,y)}{\partial y^2}\\ &\left.-\gamma_1\rho v(t,y)b(t,y)(\mu_1(t)-k_1\mu_2(t))\dfrac{\partial d_1(t,y)}{\partial y}+\lambda_1(t)\Phi_{h_1}(\Pi_1)\right\}=0.\end{split}\end{align}

We can see from (3.6) that the supremum only depends on the mean and variance of $\Pi_1$ except $\Phi_{h_1}(\Pi_1)$ . We proceed with our analysis relying on the crucial Lemma 1.

Lemma 1 (Theorem 3.1 of Liu et al. Reference Liu, Cai, Lemieux and Wang2020). If h is continuous and not constantly zero, then a maximizer $\Pi^*$ to the optimization problem

(3.7) \begin{align}\max_{\Pi \in \mathcal M^2}\Phi_h (\Pi) \mbox{ subject to } \mu (\Pi) =m \;\mbox{and}\; \sigma^2 (\Pi) = s^2\end{align}

has the following quantile function

(3.8) \begin{equation}Q_{\Pi^*}(p) = m + s\frac{ h'(1-p) }{ ||h'||_2}, \;\; \mbox{ a.e. }p\in (0,1),\end{equation}

and the maximum value of (3.7) is $\Phi_h(\Pi^*)= s||h'||_2$ .

Let $\Lambda_1(\Pi_1)$ denote the expression enclosed in braces in (3.6). We then have the following proposition.

Proposition 1. Let $h_1\in \mathcal H$ be a continuous function. For any strategy $\boldsymbol{\Pi}_1=\{\Pi_1(t),0\leqslant t\leqslant T\}\in \mathcal A(\boldsymbol{\Pi}_2)$ with mean process $\{\mu_1(t)\}_{0\leqslant t\leqslant T}$ and variance process $\{\sigma_1(t)^2\}_{0\leqslant t\leqslant T}$ , there exists a strategy $\boldsymbol{\Pi}^*_1=\{\Pi^*_1(t),0\leqslant t\leqslant T\}\in \mathcal A(\boldsymbol{\Pi}_2)$ defined by

\begin{equation*}Q_{\Pi^*_1(t)}(p) = \mu_1(t) + \sigma_1(t)\frac{ h_1'(1-p) }{ ||h_1'||_2}, \;\; \;\;\mbox{ a.e.}\;p\in (0,1),\;\;0\leqslant t\leqslant T,\end{equation*}

which shares the same mean and variance processes as $\boldsymbol{\Pi}_1$ and satisfies $\Lambda_1(\Pi_1(t)) \leqslant \Lambda_1(\Pi^*_1(t))$ .

Proof. From (3.6), $\Lambda_1(\Pi_1)$ depends on $\Pi_1$ only through $\mu_1(\Pi_1)$ , $\sigma^2_1(\Pi_1)$ and $\Phi_{h_1}(\Pi_1)$ . Since $\boldsymbol{\Pi}_1^*$ has the same mean and variance processes as $\boldsymbol{\Pi}_1$ , the problem reduces to finding $\boldsymbol{\Pi}_1^*$ that maximizes

(3.9) \begin{align} \max\limits_{\Pi_1(t)\in \mathcal{M^2}}\Phi_{h_1}(\Pi_1(t))\ \ \text{subject to }\mu_1(\Pi_1(t))=\mu_1(t),\ \sigma_1(\Pi_1(t))^2=\sigma_1(t)^2 \;\;\text{for all }t\in[0,T].\end{align}

By Lemma 1, the maximizer $\Pi_1^*(t)$ of this problem has quantile function

$$Q_{\Pi_1^*(t)}(p)=\mu_1(t)+\sigma_1(t)\dfrac{h_1'(1-p)}{\Vert h_1'\Vert_2},$$

and satisfies $\Phi_{h_1}(\Pi_1^*(t))=\sigma_1(t)\Vert h_1'\Vert_2.$ This completes the proof.

To find the equilibrium response of Agent 1, we need to solve problem (3.6). Let $\mu_1^*(t)$ and $\sigma_1^*(t)$ represent the mean and standard deviation of equilibrium response of Agent 1. As in the proof of Proposition 1, with given $\mu_1(t)$ and $\sigma_1(t)$ , $\Lambda_1(\Pi_1(t))$ only depends on $\Phi_{h_1}(\Pi_1(t))$ . Thus, we can rewrite (3.6) as

\begin{align*} \max\limits_{\Pi_1\in \mathcal{M^2}}\Lambda_1(\Pi_1)=\max\limits_{m\in \mathbb R,s\gt0}\max\limits_{\substack{\Pi_1\in \mathcal{M^2}\\ \mu_1(\Pi_1)=m,\sigma_1(\Pi_1)^2=s^2}}\Lambda_1(\Pi_1).\end{align*}

The inner maximization problem is equivalent to (3.7) or (3.9), and its maximizer $\Pi_1$ satisfies $\Phi_{h_1}(\Pi_1)= s\Vert h_1'\Vert_2$ . Substituting it back into $\Lambda(\Pi_1)$ , we obtain

\begin{align*} (\mu_1^*(t),\sigma_1^*(t))=&\mathop{\textrm{arg max}}\limits_{m\in \mathbb R,s\gt0}\left\{b(t,y)\theta(t,y)(m-k_1\mu_2(t))\right.\\ &-\dfrac{\gamma_1}{2}b^2(t,y)[(m^2+s^2)+k_1^2(\mu_2^2(t)+\sigma_2^2(t))-2k_1m\mu_2(t)]\\ &\left.-\gamma_1\rho v(t,y)b(t,y)(m-k_1\mu_2(t))\dfrac{\partial d_1(t,y)}{\partial y}+\lambda_1(t)s\Vert h_1'\Vert_2\right\}.\end{align*}

By the first-order condition, we deduce that

\begin{align*} b(t,y)\theta(t,y)-\gamma_1b^2(t,y)\mu_1^*(t)+k_1\gamma_1\mu_2(t)b^2(t,y)-\gamma_1\rho v(t,y)b(t,y)\dfrac{\partial d_1(t,y)}{\partial y}=0,\end{align*}

and

\begin{align*} -\gamma_1b^2(t,y)\sigma_1^*(t)+\lambda_1(t)\Vert h_1'\Vert_2=0.\end{align*}

So the mean and standard deviation of equilibrium response of Agent 1 are

(3.10) \begin{align} \mu_1^*(t)&=\dfrac{\theta(t,y)}{\gamma_1b(t,y)}+k_1\mu_2(t)-\dfrac{\rho v(t,y)}{b(t,y)}\dfrac{\partial d_1(t,y)}{\partial y},\end{align}

and

(3.11) \begin{align} \sigma_1^*(t)&=\dfrac{\lambda_1(t)\Vert h_1'\Vert_2}{\gamma_1b^2(t,y)}.\end{align}

By substituting (3.10) and (3.11) back into (3.3) and (3.6), we then get that $d_1(t,y)$ and $D_1(t,y)$ , respectively, satisfy

\begin{align*} \dfrac{\partial d_1(t,y)}{\partial t}+(m(t,y)-\rho v(t,y)\theta(t,y))\dfrac{\partial d_1(t,y)}{\partial y}+\dfrac{1}{2}v^2(t,y)\dfrac{\partial^2 d_1(t,y)}{\partial y^2}+\dfrac{\theta^2(t,y)}{\gamma_1}=0,\end{align*}

and

\begin{align*} &\dfrac{\partial D_1(t,y)}{\partial t}+m(t,y)\dfrac{\partial D_1(t,y)}{\partial y}+\dfrac{1}{2}v^2(t,y)\dfrac{\partial^2 D_1}{\partial y^2}(t,y)\\ &-\dfrac{\gamma_1}{2}(1-\rho^2)v^2(t,y)\left(\dfrac{\partial d_1}{\partial y}(t,y)\right)^2 -\rho v(t,y)\theta(t,y)\dfrac{\partial d_1}{\partial y}\\ &+\dfrac{\theta^2(t,y)}{2\gamma_1}-\dfrac{\gamma_1k_1^2}{2}b^2(t,y)\sigma_2^2(t)+\dfrac{\lambda_1^2(t)\Vert h_1'\Vert_2^2}{2\gamma_1b^2(t,y)}=0.\end{align*}

Repeating the procedure outlined above can yield results for Agent 2. So if $(\boldsymbol{\Pi}_1^*,\boldsymbol{\Pi}_2^*)$ constitutes a Nash equilibrium, we have

(3.12) \begin{align} \left\{ \begin{aligned} \mu_1^*(t)-k_1\mu_2^*(t)&=\dfrac{\theta(t,y)}{\gamma_1b(t,y)}-\dfrac{\rho v(t,y)}{b(t,y)}\dfrac{\partial d_1(t,y)}{\partial y},\\ \mu_2^*(t)-k_2\mu_1^*(t)&=\dfrac{\theta(t,y)}{\gamma_2b(t,y)}-\dfrac{\rho v(t,y)}{b(t,y)}\dfrac{\partial d_2(t,y)}{\partial y}. \end{aligned} \right.\end{align}

The solution to the system of equations can be directly calculated as

(3.13) \begin{align} \left\{ \begin{aligned} \mu_1^*(t)&=\dfrac{1}{1-k_1k_2}\left[\dfrac{\theta(t,y)}{b(t,y)}\left(\dfrac{1}{\gamma_1}+\dfrac{k_1}{\gamma_2}\right)-\dfrac{\rho v(t,y)}{b(t,y)}\left(\dfrac{\partial d_1(t,y)}{\partial y}+\dfrac{k_1\partial d_2(t,y)}{\partial y}\right)\right],\\ \mu_2^*(t)&=\dfrac{1}{1-k_1k_2}\left[\dfrac{\theta(t,y)}{b(t,y)}\left(\dfrac{1}{\gamma_2}+\dfrac{k_2}{\gamma_1}\right)-\dfrac{\rho v(t,y)}{b(t,y)}\left(\dfrac{\partial d_2(t,y)}{\partial y}+\dfrac{k_2\partial d_1(t,y)}{\partial y}\right)\right]. \end{aligned} \right.\end{align}

Summarizing the above, we get following theorem.

Theorem 2. For $i,j\in\{1,2\}$ and $i\neq j$ , let $d_i(t,y)$ and $D_i(t,y)$ be the solutions of

(3.14) \begin{align} \dfrac{\partial d_i(t,y)}{\partial t}+(m(t,y)-\rho v(t,y)\theta(t,y))\dfrac{\partial d_i(t,y)}{\partial y}+\dfrac{1}{2}v^2(t,y)\dfrac{\partial^2 d_i}{\partial y^2}+\dfrac{\theta^2(t,y)}{\gamma_i}=0, \end{align}

and

(3.15) \begin{align} \begin{split} &\dfrac{\partial D_i(t,y)}{\partial t}+m(t,y)\dfrac{\partial D_i(t,y)}{\partial y}+\dfrac{1}{2}v^2(t,y)\dfrac{\partial^2 D_i}{\partial y^2}(t,y)\\ &-\dfrac{\gamma_i}{2}(1-\rho^2)v^2(t,y)\left(\dfrac{\partial d_i}{\partial y}(t,y)\right)^2 -\rho v(t,y)\theta(t,y)\dfrac{\partial d_i}{\partial y}(t,y)\\ &+\dfrac{\theta^2(t,y)}{2\gamma_i}-\dfrac{\gamma_ik_i^2}{2}b^2(t,y)\sigma_j^2(t)+\dfrac{\lambda_i^2(t)\Vert h_i'\Vert_2^2}{2\gamma_ib^2(t,y)}=0, \end{split} \end{align}

with terminal condition $D_i(T,y)=d_i(T,y)=0$ . Then $(\boldsymbol{\Pi}_1^*,\boldsymbol{\Pi}_2^*)$ with quantile functions

(3.16) \begin{align} \begin{split} Q_{\Pi_i^*(t)}(p)=&\dfrac{1}{1-k_1k_2}\left[\dfrac{\theta(t,y)}{b(t,y)}\left(\dfrac{1}{\gamma_i}+\dfrac{k_i}{\gamma_j}\right)-\dfrac{\rho v(t,y)}{b(t,y)}\left(\dfrac{\partial d_i(t,y)}{\partial y}+\dfrac{k_i\partial d_j(t,y)}{\partial y}\right)\right]\\ &+\dfrac{\lambda_i(t)}{\gamma_ib^2(t,y)}h_i'(1-p) \end{split} \end{align}

is a Nash equilibrium with $p\in (0,1)$ , and $V_i(t,x,y)=x+D_i(t,y)$ is the equilibrium value function of Agent i.

From (3.16), it is evident that the equilibrium distribution of Agent i is uniquely determined by his own Choquet regularizer, $h'_i$ , and remains independent of his opponent’s regularizer, $h_j$ . Furthermore, (3.11) and (3.13) show that while the mean of Agent i’s distribution depends on both his own parameters and those of his opponent, the variance is solely determined by his own parameters, specifically $\lambda_i$ , $h_i$ , and $\gamma_i$ . These insights align with intuitive expectations in the context of RL. Although an opponent’s risk tolerance, sensitivity, or strategic decisions can influence the expected outcomes of the decision-making process, the degree of exploration, as reflected by variance, is solely a function of the agent’s intrinsic characteristics. Additionally, (3.11) shows that larger value of $\lambda_i$ indicates a stronger emphasis on exploration, leading to more dispersed exploration around the current position of Agent i. In contrast, an increase in the risk aversion parameter $\gamma_i$ reflects a more cautious approach, leading to reduced variance in the exploratory strategy.

3.3 Solution to Gauss mean return model

In this subsection, we examine the Gaussian mean return model as a special case of the state process Y(t) shown in Example 1, that is,

(3.17) \begin{align} a(t,y)=r+\sigma y,\ b(t,y)=\sigma,\ m(t,y)=\iota(Y-y),\ v(t,y)=v,\end{align}

where r, $\sigma$ , $\iota$ , v, and Y are positive constants. Thus, by (2.2), we have $\theta(t,y)=y$ . We formulate the following proposition as a direct consequence of Theorem 2.

Proposition 2. For the Gauss mean return model, $p\in (0,1)$ , $i,j\in\{1,2\}$ and $i\neq j$ , profile $(\boldsymbol{\Pi}_1^*,\boldsymbol{\Pi}_2^*)$ with quantile functions

(3.18) \begin{align} \begin{split} Q_{\Pi_i^*(t)}(p)=&\dfrac{1}{1-k_1k_2}\left[\dfrac{y}{\sigma}\left(\dfrac{1}{\gamma_i}+\dfrac{k_i}{\gamma_j}\right)-\dfrac{\rho v}{\sigma}\left((a_2^i(t)+k_ia_2^j(t))y+(a_1^i(t)+k_ia_1^j(t))\right)\right]\\ &+\dfrac{\lambda_i(t)}{\gamma_i\sigma^2}h_i'(1-p), \end{split} \end{align}

is a Nash equilibrium. Moreover, the corresponding equilibrium value function $V_i$ and $g_i$ have the following form

(3.19) \begin{equation} V_i(t,\hat x_i,y)=\hat x_i+\dfrac{1}{2}b_2^i(t)y^2+b_1^i(t)y+b_0^i(t), \end{equation}

and

(3.20) \begin{equation} g_i(t,\hat x_i,y)=\hat x_i+\dfrac{1}{2}a_2^i(t)y^2+a_1^i(t)y+a_0^i(t), \end{equation}

where $a^i_n(t)$ , $b^i_n(t)$ , $n=0,1,2$ , are continuously differentiable functions defined as

(3.21) \begin{align} \left\{ \begin{aligned} a_0^i(t)&=\dfrac{\iota^2 Y^2}{\gamma_i(\iota +\rho v)^2}\left(T-t+\frac{1-e^{-2(\iota+\rho v)(T-t)}}{2(\iota+\rho v)}-\frac{2(1-e^{-(\iota+\rho v)(T-t)})}{(\iota+\rho v)}\right)\\&\;\;\;\;+\dfrac{v^2}{2\gamma_i(\iota +\rho v)}\left(T-t-\frac{1-e^{-2(\iota+\rho v)(T-t)}}{2(\iota+\rho v)}\right), \\ a_1^i(t)&=\dfrac{\iota Y}{\gamma_i(\iota +\rho v)^2}[1-e^{-(\iota+\rho v)(T-t)}]^2,\\ a_2^i(t)&=\dfrac{1}{\gamma_i(\iota +\rho v)}[1-e^{-2(\iota+\rho v)(T-t)}], \end{aligned} \right. \end{align}

and

(3.22) \begin{align} \left\{ \begin{aligned} {b_2^i}'(t)=&\;2\iota b_2^i(t)+\gamma_i(1-\rho^2)v^2a_2^i(t)^2+2\rho va_2^i(t) -\dfrac{1}{\gamma_i},\\ {b_1^i}'(t)=&\;\iota b_1^i(t)-\iota Yb_2^i(t)+\gamma_i(1-\rho^2)v^2a_2^i(t)a_1^i(t)+\rho va_1^i(t),\\ {b_0^i}'(t)=&-\iota Yb_1^i(t)-\dfrac{v^2}{2}b_2^i(t)+\dfrac{\gamma_i(1-\rho^2)v^2}{2}a_1^i(t)^2\\ &+\dfrac{\gamma_ik_i^2\sigma^2}{2}\sigma_j(t)^2-\dfrac{\lambda_i^2(t)\Vert h_i'\Vert_2^2}{2\gamma_i\sigma^2}, \end{aligned} \right. \end{align}

with $b_0^i(T)=b_1^i(T)= b_2^i(T)=0$ .

Proof. For the Gauss mean return model, (3.14) can be simplified as

(3.23) \begin{align} \dfrac{\partial d_i(t,y)}{\partial t}+[\iota(Y-y)-\rho v y]\dfrac{\partial d_i(t,y)}{\partial y}+\dfrac{v^2}{2}\dfrac{\partial^2d_i(t,y)}{\partial y^2}+\dfrac{y^2}{\gamma_i}=0. \end{align}

By letting $d_i(t,y)=\dfrac{1}{2}a_2^i(t)y^2+a_1^i(t)y+a_0^i(t)$ and substituting it into (3.23), we obtain

(3.24) \begin{align} \left\{ \begin{aligned} {a_2^i}'(t)&=2a_2^i(t)(\iota +\rho v)-\dfrac{2}{\gamma_i},& a_2^i(T)=0,\\ {a_1^i}'(t)&=a_1^i(t)(\iota +\rho v)-a_2^i(t)\iota Y,& a_1^i(T)=0,\\ {a_0^i}'(t)&=-a_1^i(t)\iota Y-\dfrac{v^2}{2}a_2^i(t),& a_0^i(T)=0. \end{aligned} \right. \end{align}

It can be shown that (3.21) is the solution to (3.24). By substituting $d_i$ into (3.16), we derive (3.18). Consequently, $(\boldsymbol{\Pi}^*_1,\boldsymbol{\Pi}^*_2)$ is indeed a Nash equilibrium. Similarly, by simplifying (3.15), we can get $D_i(t,y)=\dfrac{1}{2}b_2^i(t)y^2+b_1^i(t)y+b_0^i(t)$ with $b^i_n(t)$ , $n=0,1,2$ , given by (3.22).

By setting $y=(a-r)/\sigma,$ $\iota=0$ and $v=0$ in (3.17), the price dynamics (2.1) reduces to the classical Black-Scholes model. In this case, the market becomes complete and the corresponding results are straightforward, as stated below.

Corollary 1. In the Black-Scholes model, for $p\in (0,1)$ , $i,j\in\{1,2\}$ and $i\neq j$ , profile $(\boldsymbol{\Pi}_1^*,\boldsymbol{\Pi}_2^*)$ with quantile functions

(3.25) \begin{align} Q_{\Pi_i^*(t)}(p)=&\dfrac{1}{1-k_1k_2}\left[\dfrac{a-r}{\sigma^2}\left(\dfrac{1}{\gamma_i}+\dfrac{k_i}{\gamma_j}\right)\right] +\dfrac{\lambda_i(t)}{\gamma_i\sigma^2}h_i'(1-p) \end{align}

is a Nash equilibrium.

In the Black-Scholes model, the influence of $k_1$ , $k_2$ , $\gamma_1$ , and $\gamma_2$ on the equilibrium strategies is clearly reflected in (3.25), aligning with the properties of Gauss mean return model discussed in Section 6.

4. Policy iteration

In this section, we employ the policy iteration method to find equilibrium strategies in two steps. For $i, j \in \{1, 2\}$ and $i \neq j$ , we first fix $\boldsymbol{\Pi}_j$ and estimate the associated value function $V_i^{\boldsymbol{\Pi}_i}$ given a policy $\boldsymbol{\Pi}_i$ . Then we update the previous policy $\boldsymbol{\Pi}_i$ to a new one $\tilde{\boldsymbol{\Pi}}_i$ based on the obtained value function $V_i^{\boldsymbol{\Pi}_i}$ . Despite the learning process not leading to a monotone iteration algorithm due to the “optimality” is in the sense of equilibrium, we demonstrate that the iterative process converges uniformly to the desired equilibrium policy.

Assuming $\boldsymbol{\Pi}_j$ is fixed, and letting $\boldsymbol{\Pi}_i$ be an admissible strategy for $i, j \in \{1,2\}$ with $i \neq j$ , we denote the value function under $\boldsymbol{\Pi}_i$ as $V_i^{\boldsymbol{\Pi}_i}(t,\hat x_i,y)=J_i(t,\hat x_i,y;\boldsymbol{\Pi}_i,\boldsymbol{\Pi}_j)$ . Similarly, we define $g_i^{\boldsymbol{\Pi}_i}(t,\hat x_i,y)=\mathbb E_t[\hat X_i^{\Pi_i,\Pi_j}(T)]$ . According to Björk et al. (Reference Björk, Khapko and Murgoci2017), the functions $V_i^{\boldsymbol{\Pi}_i}$ and $g_i^{\boldsymbol{\Pi}_i}$ satisfy the following equations

(4.1) \begin{align} &\mathcal L^{\Pi_i,\Pi_j}V^{\boldsymbol{\Pi}_i}_i-\dfrac{\gamma_i}{2}\mathcal L^{\Pi_i,\Pi_j}{g^{\boldsymbol{\Pi}_i}_i}^2+\gamma_ig^{\boldsymbol{\Pi}_i}_i\mathcal L^{\Pi_i,\Pi_j}g^{\boldsymbol{\Pi}_i}_i+\lambda_i(t)\Phi_{h_i}(\Pi_i)=0,\end{align}

and

(4.2) \begin{align} \mathcal L^{\Pi_i,\Pi_j}g^{\boldsymbol{\Pi}_i}_i(t,\hat x_i,y)=0,\end{align}

with

(4.3) \begin{align} V^{\boldsymbol{\Pi}_i}_i(T,\hat x_i,y)=\hat x_i,\;\;\ g^{\boldsymbol{\Pi}_i}_i(T,\hat x_i,y)=\hat x_i.\end{align}

Theorem 3. For $p\in (0,1)$ , $i,j\in\{1,2\}$ , and $i\neq j$ , with $\boldsymbol{\Pi}_j$ fixed, let $\boldsymbol{\Pi}_i^0$ be the initial policy of Agent i with quantile function given by

(4.4) \begin{align} Q_{\Pi_i^0(t)}(p)=\dfrac{y}{\gamma_i\sigma}+k_i\mu_j(t)-\dfrac{\rho v}{\sigma}(a_2^{i0}(t)y+a_1^{i0}(t))+\psi^0 h_i'(1-p).\end{align}

Choose one policy

$$\Pi_i(t)\in \mathop{\textrm{arg max}}\limits_{\Pi_i\in\mathcal M^2}\{\mathcal L^{\Pi_i,\Pi_j}V^{\boldsymbol{\Pi}_i^n}_i-\dfrac{\gamma_i}{2}\mathcal L^{\Pi_i,\Pi_j}{g^{\boldsymbol{\Pi}_i^n}_i}^2+\gamma_ig^{\boldsymbol{\Pi}_i^n}_i\mathcal L^{\Pi_i,\Pi_j}g^{\boldsymbol{\Pi}_i^n}_i+\lambda_i(t)\Phi_{h_i}(\Pi_i)\},$$

and denote this policy as $\boldsymbol{\Pi}_i^{n+1}$ , $n=0,1,2,...$ Then the following statements holds.

1. (i) The sequence of updated policies $\boldsymbol{\Pi}_i^{n}$ for $n\geqslant 1$ has the quantile function

   (4.5) \begin{align} Q_{\Pi_i^n(t)}(p)=\dfrac{y}{\gamma_i\sigma}+k_i\mu_j(t)-\dfrac{\rho v}{\sigma}(a_2^{in}(t)y+a_1^{in}(t))+\dfrac{\lambda_i(t)}{\gamma_i\sigma^2} h_i'(1-p), \end{align}
   where $a_1^{in}$ and $a_2^{in}$ satisfy
   (4.6) \begin{align} \left\{ \begin{aligned} {a_2^{in}}'(t)&=2\iota a_2^{in}(t)+2\rho va_2^{in-1}(t) -\dfrac{2}{\gamma_i},&a_2^{in}(T)=0,\\ {a_1^{in}}'(t)&=\iota a_1^{in}(t) +\rho va_1^{in-1}(t)-a_2^{in}(t)\iota Y,& a_1^{in}(T)=0. \end{aligned} \right. \end{align}

2. (ii) As $n\to \infty$ , $a_1^{in}(t)$ and $a_2^{in}$ uniformly converge to $a_1^i$ and $a_2^i$ in (3.21), respectively.

Theorem 3 shows that the iteration does not change the form of the policy (see (4.5)), and thus, it suffices to parameterize the iterative policy through two deterministic functions $(a_1^{(in)}(t), a_2^{(in)}(t))$ . In particular, when the initial policy chosen as the form of the equilibrium policy in Proposition 2, our algorithm is guaranteed to uniformly converge to the equilibrium policy.

The next result guarantees the convergence of policies as the two agents iterate simultaneously.

Proof. We just need to prove that the mean and variance of $\Pi_i$ converge to the mean and variance of $\Pi_i^*$ . The convergence of variance is obvious according to the proof of Theorem 3(i) and (4.10). Let $\mu_i^n(t)$ be the mean of $\Pi_i^n(t)$ . Based on (3.10), we have

\begin{align*} \begin{bmatrix} \mu_1^{n+1}(t)\\ \mu_2^{n+1}(t) \end{bmatrix}= \begin{bmatrix} 0\;\;\; & k_1 \\ k_2\;\;\; & 0 \end{bmatrix} \begin{bmatrix} \mu_1^{n}(t)\\ \mu_2^{n}(t) \end{bmatrix}+ \begin{bmatrix} \dfrac{y}{\gamma_1\sigma}-\dfrac{\rho v}{\sigma}(a_2^1(t)y+a_1^1(t))\\[4mm] \dfrac{y}{\gamma_2\sigma}-\dfrac{\rho v}{\sigma}(a_2^2(t)y+a_1^2(t)) \end{bmatrix}. \end{align*}

Consider the normed space $\mathbb R^2$ with $\Vert\cdot\Vert$ defined as $\Vert \vec{x}\Vert=\max\{x_1,x_2\}$ for $\vec{x}=[x_1,x_2]'\in \mathbb R^2$ . It is well known that $(\mathbb R^2,\Vert\cdot\Vert)$ is a Banach space. Define

\begin{align*} f(\vec{x})=\begin{bmatrix} 0\;\;\; & k_1 \\ k_2\;\;\; & 0 \end{bmatrix} \vec{x}+ \begin{bmatrix} \dfrac{y}{\gamma_1\sigma}-\dfrac{\rho v}{\sigma}(a_2^1(t)y+a_1^1(t))\\[4mm] \dfrac{y}{\gamma_2\sigma}-\dfrac{\rho v}{\sigma}(a_2^2(t)y+a_1^2(t)) \end{bmatrix}. \end{align*}

We have $\Vert f(\vec{x})-f(\vec{y})\Vert\leqslant \max\{k_1,k_2\}\Vert \vec{x}-\vec{y} \Vert$ , establishing f as a contraction mapping with a unique fixed point. By (3.12), the fixed point is precisely $[\mu_1^*(t),\mu_2^*(t)]'$ . Thus, the mean of $\Pi_i^n(t)$ also converges. Moreover, let M ^′ be an uniformly upper bound of $|\mu_i^0(t)-\mu_i^*(t)|$ . Then, we have

\begin{align*} \left\Vert\begin{bmatrix} \mu_1^{n}(t)\\ \mu_2^{n}(t) \end{bmatrix}- \begin{bmatrix} \mu_1^{*}(t)\\ \mu_2^{*}(t) \end{bmatrix}\right\Vert=& \left\Vert\begin{bmatrix} 0 & k_1 \\ k_2 & 0 \end{bmatrix} \begin{bmatrix} \mu_1^{n-1}(t)-\mu_1^*(t)\\ \mu_2^{n-1}(t)-\mu_2^*(t) \end{bmatrix}\right\Vert\\ =&\left\Vert\begin{bmatrix} 0 & k_1 \\ k_2 & 0 \end{bmatrix}^n \begin{bmatrix} \mu_1^{0}(t)-\mu_1^*(t)\\ \mu_2^{0}(t)-\mu_2^*(t) \end{bmatrix}\right\Vert\\ \leqslant & M'(\max\{k_1,k_2\})^n. \end{align*}

Thus, for $p\in (0,1)$ , $i,j\in\{1,2\}$ and $i\neq j$ , $Q_{\Pi_i^n(t)}(p)$ converges uniformly to $Q_{\Pi_i^*(t)}(p)$ of (3.18) as $n\to\infty$ .

5. RL algorithm design

In this section, we devise an algorithm to learn the Nash equilibrium. As mentioned in the introduction, game scenarios involve multiple agents in the environment, requiring the utilization of multi-agent reinforcement learning algorithms, which inherently introduce greater complexity. Specifically, in single-agent reinforcement learning, a basic assumption is the stability of the environment, wherein the transition probability and reward function remain constant. However, when other intelligent agents are introduced into the environment, this assumption no longer holds true. In a multi-agent context, any change in one agent’s strategy can significantly impact other agents, leading to dynamic evolution of the environment with their strategies. Moreover, as the number of agents increases, the complexity of training also escalates. Fortunately, in our model, each agent will only affect the mean of the other agent’s strategy. According to (3.12), the difference $\mu_i(t)-k_i\mu_j(t)$ for $i,j\in\{1,2\}$ and $i\neq j$ remains fixed, regardless of how the agents adjust their strategies. Thus, once the difference is known, we can directly apply the learning procedure outlined in Theorem 4. This allows us to reduce the game problem into two independent optimization problems. We then apply the method proposed in Section 4 to learn the difference. Below, we briefly introduce this method.

Assume that the risk-free interest rate r, the risk-aversion coefficient $\gamma_1, \gamma_2$ , the sensitivity coefficient $k_1, k_2$ , and the exploration weight $\lambda_1, \lambda_2$ are known. The agents have no prior knowledge of S(t) and Y(t), but can observe the pair (S(t),Y(t)) at each time t. In the continuous time setting, we discretize the interval [0, T] into N subintervals with equal length $\Delta t=t_{k+1}-t_k, k=0,1,...,N-1$ . Based on (2.1) and (2.3), we get

\begin{align*} \mathrm{d} X_i^{u_i}(t)=u_i(t)\dfrac{\mathrm{d} e^{-rt}S(t)}{e^{-rt}S(t)}.\end{align*}

Therefore, when the Agent i follows strategy $\Pi_i(t_k)$ at time $t_k$ and samples action $u_i(t_k)$ from $\Pi_i(t_k)$ , the discounted wealth at $t_{k+1}$ is

\begin{align*} X_i^{\Pi_i}(t_{k+1})\approx X_i^{\Pi_i}(t_k)+u_i(t_k)\dfrac{e^{-rt_{k+1}}S(t_{k+1})-e^{-rt_k}S(t_k)}{e^{-rt_k}S(t_k)}.\end{align*}

Assume that the policy $\boldsymbol{\Pi}_j$ of Agent j is fixed, and that both $\Pi_i(t)$ and $\Pi_j(t)$ have density $\pi_i(t)$ and $\pi_j(t)$ . The algorithm for learning the difference $\mu_i(t)-k_i\mu_j(t)$ is based on the standard idea of policy evaluation followed by policy update. We do not adopt the approach of generalized policy iteration, in which policy evaluation and policy update interact with each other. There are three main reasons for this. First, Theorem 3 guarantees a convergence rate that is sufficiently fast, making the generalized policy iteration unnecessary. Second, unlike classical reinforcement learning, our equilibrium strategy involves two value functions, $V_i$ and $g_i$ , which complicates the policy evaluation step, particularly in the context of generalized policy iteration. Third, as previously discussed, our policy update does not always yield an improved policy, which further limits the applicability of generalized policy iteration.

We represent the policy $\boldsymbol{\Pi}_i$ of Agent i using its quantile function. By Theorem 3, the quantile function can be defined as

(5.1) \begin{align}Q_{\Pi_i^{\Psi}(t)}(p)=k_i\mu_j(t)+\psi_0 \dfrac{y}{\gamma_i}-\psi_0\psi_1(a_2^{i0}(t)y+a_1^{i0}(t))+\dfrac{\lambda_i}{\gamma_i \psi_2^2} h_i'(1-p),\end{align}

where $\Psi = (\psi_0,\psi_1,\psi_2)\in\mathbb R^3$ denotes the set of parameters to be learned, and $a_1^{i0}$ and $a_2^{i0}$ are functions that will be updated according to Theorem 3(i). The parameters $\Psi$ can be initialized either randomly or using predetermined constants, and $a_1^{i0}$ and $a_2^{i0}$ can be initialized as zero functions.

For policy evaluation procedure, based on Proposition 2, we parameterize $V_i$ and $g_i$ in the form of (3.19) and (3.20), respectively, as follows

(5.2) \begin{align} \begin{split} V_i^{\Theta}(t,\hat x_i,y)&=\hat x_i+\dfrac{1}{2}p(\theta_i^{V,2},T-t)y^2+p(\theta_i^{V,1},T-t)y+p(\theta_i^{V,0},T-t),\\ g_i^{\Theta}(t,\hat x_i,y)&=\hat x_i+\dfrac{1}{2}p(\theta_i^{g,2},T-t)y^2+p(\theta_i^{g,1},T-t)y+p(\theta_i^{g,0},T-t), \end{split}\end{align}

where $p(\theta,t)$ is a suitable parametric function with coefficient vector $\theta\in \mathbb R^d$ . The full parameter set is denoted by $\Theta=(\theta_i^{V,0},\theta_i^{V,1},\theta_i^{V,2},\theta_i^{g,0},\theta_i^{g,1},\theta_i^{g,2})\in\mathbb R^{6d}$ . A typical choice for $p(\theta,t)$ is a linear combination of the first d terms of a basis expansion, such as a truncated Taylor or Fourier series, where $\theta$ denotes the corresponding coefficients. Since $g^{\boldsymbol{\Pi}_i}$ satisfies $\mathcal L^{\Pi_i,\Pi_j}g^{\boldsymbol{\Pi}_i}_i(t,\hat x_i,y)=0$ with terminal condition $g^{\boldsymbol{\Pi}_i}_i(T,\hat x_i,y)=\hat x_i$ for any $\boldsymbol{\Pi}_i$ , it can be interpreted via the Feynman-Kac formula as the value function of a time-consistent optimal problem. This allows it to be evaluated using continuous-time reinforcement learning. Theorem 3 in Jia and Zhou (Reference Jia and Zhou2022a) and Theorem 4 in Jia and Zhou (Reference Jia and Zhou2022b) show that $g^{\boldsymbol{\Pi}_i}$ can be estimated by minimizing the martingale loss function

\begin{align*} \mathrm{ML}_{g}(\Theta)\;:\!=\;\dfrac{1}{2}\mathbb{E}\left[\sum\limits_{k=0}^{N-1} \left(\hat X_i(T) - g_i^{\Theta}(t_k,\hat X_i(t_k),Y(t_k))\right)^2\Delta t\right],\end{align*}

which corresponds to the continuous-time analogue of the Monte Carlo policy evaluation with function approximation (e.g., Sutton and Barto, Reference Sutton and Barto2018). Then we can use the stochastic gradient descent method to minimize $\mathrm{ML}_{g}(\Theta)$ , which is a standard method in reinforcement learning. The gradient of $\mathrm{ML}_{g}(\Theta)$ with respect to $\Theta$ can be computed as

(5.3) \begin{align} \nabla_g \Theta=-\sum\limits_{k=0}^{N-1} \left(\hat X_i(T) - g_i^{\Theta}(t_k,\hat X_i(t_k),Y(t_k))\right)\dfrac{\partial g_i^{\Theta}}{\partial \Theta}\Delta t,\end{align}

and the parameters $\Theta$ are updated according to

\begin{align*} \Theta_{n+1}=\Theta_n-\alpha_g\nabla_g \Theta.\end{align*}

Once $g_i$ is evaluated, we can proceed to evaluate $V_i$ , which satisfies

(5.4) \begin{align} \mathcal L^{\Pi_i,\Pi_j}V_i(t,\hat x_i,y) + R(t,\hat x_i,y)= 0,\;\;\;V_i(T,\hat x_i,y)=\hat x_i,\end{align}

where the residual term is given by

$$R(t,\hat x_i,y)=-\dfrac{\gamma_i}{2}\mathcal L^{\Pi_i,\Pi_j}g_i^2+\gamma_i g_i\mathcal L^{\Pi_i,\Pi_j}g_i+\lambda_i(t)\Phi_{h_i}(\Pi_i). $$

The structure of (5.4), together with the Feynman-Kac formula, implies that $V_i$ can also be interpreted as the value function of a time-consistent optimal control problem. Consequently, $V_i$ can be evaluated analogously to $g_i$ , by minimizing a martingale loss function

\begin{align*} \mathrm{ML}_{V}(\Theta)\;:\!=\;\dfrac{1}{2}\mathbb{E}\left[\sum\limits_{k=0}^{N-1} \left(\hat X_i(T) - V_i^{\Theta}(t_k,\hat X_i(t_k),Y(t_k))+\sum\limits_{l=k}^{N-1}R(t_l,\hat X_i(t_l),Y(t_l))\Delta t\right)^2\Delta t\right],\end{align*}

and the gradient of $\mathrm{ML}_{V}(\Theta)$ with respect to $\Theta$ can be computed as

(5.5) \begin{align} \nabla_V \Theta=-\sum\limits_{k=0}^{N-1} \left(\hat X_i(T) - V_i^{\Theta}(t_k,\hat X_i(t_k),Y(t_k))+\sum\limits_{l=k}^{N-1}R(t_l,\hat X_i(t_l),Y(t_l))\Delta t\right)\dfrac{\partial V_i^{\Theta}}{\partial \Theta}\Delta t.\end{align}

Then $\Theta$ can be updated as

\begin{align*} \Theta_{n+1}=\Theta_n-\alpha_V\nabla_V \Theta.\end{align*}

By applying Itô’s formula, for any $\varphi\in C^{1,2,2}$ and strategy $\Pi_i$ of Agent i, we have

\begin{align*}\begin{split} &\mathbb E_t[\varphi(t+\Delta t,\hat X_i^{\Pi_i,\Pi_j}(t+\Delta t),Y(t+\Delta t))]-\varphi(t,\hat X_i^{\Pi_i,\Pi_j}(t),Y(t))\\ & =\mathbb E_t\int_t^{t+\Delta t}\mathcal L^{\Pi_i,\Pi_j}\varphi(s,\hat X_i^{\Pi_i,\Pi_j}(s),Y(s))\mathrm{d} s,\end{split}\end{align*}

and thus

(5.6) \begin{align} \mathcal L^{\Pi_i,\Pi_j}\varphi(t,\hat X_i^{\Pi_i,\Pi_j}(t),Y(t))\approx \dfrac{\varphi(t+\Delta t,\hat X_i^{\Pi_i,\Pi_j}(t+\Delta t),Y(t+\Delta t))-\varphi(t,\hat X_i^{\Pi_i,\Pi_j}(t),Y(t))}{\Delta t}.\end{align}

We can use (5.6) to approximate $R(t,\hat x_i,y)$ by replacing $\varphi$ with $V_i$ , $g_i$ and $g_i^2$ .

For policy update procedure, we want to maximize the following function

(5.7) \begin{align} L_i(\Psi;t,\hat x,y)\;:\!=\;\mathcal L^{\Pi^{\Psi}_i(t),\Pi_j}V^{\Theta}_i-\dfrac{\gamma_i}{2}\mathcal L^{\Pi^{\Psi}_i(t),\Pi_j}{g^{\Theta}_i}^2+\gamma_ig^{\Theta}_i\mathcal L^{\Pi^{\Psi}_i(t),\Pi_j}g^{\Theta}_i+\lambda_i\Phi_{h_i}(\Pi^{\Psi}_i(t))\end{align}

for all possible $t,\hat x,y$ . By the proof of Theorem 3, we only need to maximize $L_i(\Psi;\;t,\hat x,y)$ at initial state $(t_0,\hat x_0,y_0)$ . First, we can update policy (5.1) to

(5.8) \begin{align} Q_{\Pi_i^{\Psi}(t)}(p)=k_i\mu_j(t)+\psi_0 \dfrac{y}{\gamma_i}-\psi_0\psi_1(a_2^{i1}(t)y+a_1^{i1}(t))+\dfrac{\lambda_i(t)}{\gamma_i \psi_2^2} h_i'(1-p),\end{align}

where $a_2^{i1}(t)y+a_1^{i1}(t)$ comes from the derivative of evaluated $g^{\Theta}$ with respect to y. Then we only need to maximize $L_i(\Psi;\;t_0,\hat x_0,y_0)$ with respect to $\Psi$ . The maximization of $L_i(\Psi;\;t_0,\hat x_0,y_0)$ can also be achieved by stochastic gradient ascent method. Since

\begin{align*} \mathcal{L}^{\Pi_i,\Pi_j}\varphi(t,\hat x,y)=\int_R\int_R\mathcal{L}^{u_i,u_j}\varphi(t,\hat x,y)\pi_j(u_j)\pi_i(u_i) \mathrm{d} u_j\mathrm{d} u_i,\end{align*}

then we have

\begin{align*} \dfrac{dL_i}{d\Psi}=&\int_R\int_R\left[\left(\mathcal{L}^{u_i,u_j}V^{\Theta}_i-\dfrac{\gamma_i}{2}\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\gamma_ig^{\Theta}_i\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\lambda_i\Phi_{h_i}(\Pi^{\Psi}_i)\right)\pi_j(u_j)\dfrac{\partial \pi_i^{\Psi}(u_i)}{\partial \Psi}\right.\\ &+\left.\lambda_i\dfrac{\partial \Phi_{h_i}(\Pi_i^{\Psi})}{\partial \Psi}\pi_j(u_j)\pi_i(u_i)\right]\mathrm{d} u_j\mathrm{d} u_i.\\ =&\int_R\int_R\left[\left(\mathcal{L}^{u_i,u_j}V^{\Theta}_i-\dfrac{\gamma_i}{2}\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\gamma_ig^{\Theta}_i\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\lambda_i\Phi_{h_i}(\Pi^{\Psi}_i)\right)\pi_j(u_j)\pi_i(u_i)\dfrac{\frac{\partial \pi_i^{\Psi}(u_i)}{\partial \Psi}}{\pi_i(u_i)}\right.\\ &+\left.\lambda_i\dfrac{\partial \Phi_{h_i}(\Pi_i^{\Psi})}{\partial \Psi}\pi_j(u_j)\pi_i(u_i)\right]\mathrm{d} u_j\mathrm{d} u_i.\\ =&\int_R\int_R\left[\dfrac{\partial \log \pi_i^{\Psi}(u_i)}{\partial \Psi}\left(\mathcal{L}^{u_i,u_j}V^{\Theta}_i-\dfrac{\gamma_i}{2}\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\gamma_ig^{\Theta}_i\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\lambda_i\Phi_{h_i}(\Pi^{\Psi}_i)\right)\right.\\ &+\left.\lambda_i\dfrac{\partial \Phi_{h_i}(\Pi_i^{\Psi})}{\partial \Psi}\right]\pi_j(u_j)\pi_i(u_i)\mathrm{d} u_j\mathrm{d} u_i.\end{align*}

Thus, the gradient of $L_i(\Psi)$ with respect to $\Psi$ can be computed as

(5.9) \begin{align} \nabla\Psi=\dfrac{\partial \log \pi_i^{\Psi}(u_i)}{\partial \Psi}\left(\mathcal{L}^{u_i,u_j}V^{\Theta}_i-\dfrac{\gamma_i}{2}\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\gamma_ig^{\Theta}_i\mathcal{L}^{u_i,u_j}g^{\Theta}_i+\lambda_i\Phi_{h_i}(\Pi^{\Psi}_i)\right) +\lambda_i\dfrac{\partial \Phi_{h_i}(\Pi_i^{\Psi})}{\partial \Psi},\end{align}

and the infinitesimal generator $\mathcal{L}^{u_i,u_j}$ can be approximated by samples. Note that if the distribution $\Pi^{\Psi}$ is supported on $[S^{\Psi}_{\min},S^{\Psi}_{\max}]$ that depends on the parameter $\Psi$ , then the gradient computation must also account for the derivatives of the interval endpoints, $S^{\Psi}_{\min}$ and $S^{\Psi}_{\max}$ , with respect to $\Psi$ . Then we can update $\Psi$ as

\begin{align*} \Psi_{n+1}=\Psi_n+\alpha_{\Psi}\nabla\Psi.\end{align*}

By repeating policy evaluation and policy update, we can obtain the final policy. The complete procedure is summarized in Algorithm 1.

Algorithm 1

6. Numerical results

Given the variety of Choquet regularizers available, it is possible to select different regularizers for each agent. In Equation (2.13), h^′(x) represents the “probability weight” assigned to x when calculating the (nonlinear) Choquet expectation (see, e.g., Gilboa and Schmeidler, Reference Gilboa and Schmeidler1989 and Quiggin, Reference Quiggin1982). Consequently, the choice of the distortion function h can directly influence the agent’s attitude toward risk. As shown by Han et al. (Reference Han, Wang and Zhou2023), Choquet regularizers can generate several widely used exploratory samplers, such as the $\varepsilon$ -greedy strategy, exponential, uniform, and Gaussian. Below, we assume that the agents adopt different Choquet regularizers, resulting in their optimal exploration distributions being normal and exponential, respectively.

Assume that Agent 1 applies the the Choquet regularizer

(6.1) \begin{equation} \Phi_{h_1} (\Pi_1)=\int_0^1 Q_{\Pi_1}(p) z (p) \mathrm{d} p,\end{equation}

where z is the quantile function of a standard normal distribution, yielding $h_1(p)=\int_0^p z(1-s)\mathrm{d} s$ with $p\in [0,1]$ . Further, Agent 2 uses the Choquet regularizer

(6.2) \begin{equation} \Phi_{h_2} (\Pi_2) = -\int_0^1 Q_{\Pi_2}(1-p)(\log(p)+1) \mathrm{d} p.\end{equation}

It is known as the cumulative residual entropy (e.g., Rao et al., Reference Rao, Chen, Vemuri and Wang2004 and Hu and Chen, Reference Hu and Chen2020) and $h_2(p)= -p\log p$ with $p\in [0,1]$ . For Gauss mean return model, let

$$\mu_i^*(t)= \dfrac{1}{1-k_1k_2}\left[\dfrac{y}{\sigma}\left(\dfrac{1}{\gamma_i}+\dfrac{k_i}{\gamma_j}\right)-\dfrac{\rho v}{\sigma}\left((a_2^i(t)+k_ia_2^j(t))y+(a_1^i(t)+k_ia_1^j(t))\right)\right],$$

where $a^i_n(t)$ , $n=1,2$ , are given by (3.21). Based on (3.16), the equilibrium distribution $\Pi_1^*$ is a normal distribution given as

$${\Pi}_1^{*} (\!\cdot\!;\;t)= {\mathcal N}\left(\mu_1^*(t), \dfrac{\lambda_1^2(t)}{\gamma_1^2\sigma^4}\right),$$

and the equilibrium distribution $\Pi_2^*$ is an exponential distribution given as

$$ \Pi_2^*(u;\;t) = 1-\exp\left[\dfrac{\gamma_2 \sigma^2}{\lambda_2(t)}(\mu_2(t)-u)-1 \right].$$

We first investigate the influence of parameters $k_i$ , $\gamma_i$ and t on the equilibrium strategies of both agents. We assume that $\lambda_1(t)=\lambda_2(t) = \lambda_0 e^{\lambda_0(T-t)}$ . Unless otherwise specified, the parameters in (3.16) are set as in Table 1.Footnote ¹

Table 1.

Parameter values used in the model.

Figure 1.

The effects of t, $k_1$ , $k_2$ , $\gamma_1$ , and $\gamma_2$ on the Nash equilibrium.

In Figure 1, we set $T=20$ to make the differences between the plotted curves visually distinguishable. The first and second rows display the density functions for Agent 1 at $t = 0.1$ and $t = 18$ , respectively, while the third and fourth rows correspond to Agent 2’s density functions. For clarity, we focus on the characteristics of the parameters for Agent 1, as Agent 2’s performance with respect to these parameters is similar. The key observations are as follows.

1. (i) As $k_1$ increases, Agent 1 tends to adopt riskier strategies, leading to a higher mean investment in risky assets. This suggests that greater sensitivity to the opponent’s performance enhances Agent 1’s motivation to outperform.

2. (ii) As Agent 1’s risk aversion parameter $\gamma_1$ rises, the mean of the equilibrium distribution decreases. This indicates that higher levels of risk aversion prompt Agent 1 to adopt more cautious strategies, reducing their expected investment in risky assets.

3. (iii) As Agent 2 becomes more sensitive to Agent 1’s behavior (i.e., as $k_2$ increases), Agent 1 tends to adopt riskier strategies to increase the likelihood of achieving higher returns. This suggests that Agent 2’s increased sensitivity further motivates Agent 1 to excel. Economically, when relative wealth is taken into account, a larger $k_2$ leads Agent 2 to hold more risky assets, which in turn compels Agent 1 to increase its holdings of risky assets as well in order to maintain the current position and avoid falling behind.

4. (iv) As Agent 2 becomes more risk-seeking (i.e., as $\gamma_2$ decreases), Agent 1 adopts more aggressive strategies, investing a higher mean in risky assets. This adjustment is necessary to prevent Agent 1 from losing market share or competitive advantages if they fail to align their strategy with the increased risks.

5. (v) Comparing the equilibrium strategies between the first and second rows, the main difference lies in time t. Since a time-decaying temperature parameter $ \lambda_1(t)$ is employed, it means that as time progresses, the weight assigned to exploration decreases as time progresses. Consequently, the variance of the equilibrium strategies also decreases, as reflected in (3.11).

It is noteworthy that an agent’s own parameters consistently have a greater impact on their strategy than those of their opponent. This outcome is expected. For instance, when the opponent places greater emphasis on the wealth gap, it is the opponent who adjusts their behavior by increasing their own investment. The optimal response of the original decision-maker remains relatively stable, as the opponent’s concern for relative performance predominantly affects their own risk-taking behavior, rather than prompting significant changes in the other agent’s strategy.

Furthermore, Figure 1 shows that the mean of the equilibrium strategies decreases over time t. This trend can be attributed to the negative value chosen for the parameter $\rho$ . In the following, we investigate how the mean strategy evolves with respect to both $\rho$ and t.

Figure 2 shows that as the correlation $\rho$ increases, the mean of the strategies for both Agent 1 and Agent 2 decrease. This is because positive correlation between the Brownian shocks driving the asset price and the factor process Y(t) increases uncertainty in the asset dynamics, increasing effective risk exposure and making the investor more cautious. However, when risks are negatively correlated ( $\rho \in [\!-\!1,0)$ ), adverse movements in Y(t) tend to offset those in the asset price through the diffusion term, effectively providing a natural hedging effect that reduces risk. As a result, the investor is willing to allocate more on average under negative correlation.

Figure 2.

The effects of t and $\rho$ on the mean of the Nash equilibrium.

Additionally, for negatively correlated risks, we find that the mean investment decreases as time t increases. This is because as time advances toward the terminal time T, the opportunity to exploit this hedging effect through dynamic rebalancing diminishes. The shorter the remaining horizon, the less effective the negative correlation becomes at mitigating risk over time, leading the investor to reduce the position. Conversely, when risks are positively correlated, the mean investment increases with time t. At early stages, the cumulative risk from positively correlated shocks is higher, prompting caution. However, as t approaches T, the remaining exposure horizon shortens, reducing the impact of correlated shocks on total risk and allowing the investor to increase the allocation. Specifically, when $\rho = 0$ , the shocks driving the factor process Y(t) and the asset price are independent. In this case, the mean investment is time-invariant, resulting in identical values at different t, as can be seen directly from (3.18).

Next, we conduct numerical experiments with simulated data to demonstrate our Algorithm 1 for Gauss mean return model and Black-Scholes model. We first emphasize that there are some factors that affect the accuracy of the algorithm. First of all, Theorems 3 and 4 give theoretical convergence, which provides support for algorithm design, but the theoretical results depend on the parameters of the market model and these parameters are actually unknown. Secondly, as mentioned earlier, a very important property in classical reinforcement learning is that the strategy obtained after each update is better than before, but this property no longer holds in time-inconsistent problem. These increases the impact of errors on the convergence in each algorithm iteration, which may affect the convergence to the true equilibrium strategy. Finally, due to the use of maximizing $L_i(\Psi;\;t_0,\hat X_i(t_0),Y(t_0))$ instead of maximizing $L_i(\Psi;\;t,\hat x,y)$ for all possible $t,\hat x,y$ , each iteration highly relies on the current sample, and different samples may cause the algorithm to converge to different strategies. Meanwhile, similar to classical reinforcement learning, due to the time horizon is finite, each sample only has $N+1$ time points, which is far less than the number required to make the strategy maximizing $L_i(\Psi;\;t_0,\hat X_i(t_0),Y(t_0))$ and maximizing $L_i(\Psi;\;t,\hat x,y)$ for all possible $t,\hat x,y$ close enough. In summary, due to various reasons, algorithms for time-inconsistency problems rely more on model settings, especially parameter selection, than the classical reinforcement learning.

We use the stock process parameters detailed in Table 1. Additionally, the other parameter settings for the algorithm are presented in Table 2. For Gauss mean return model, the value function is parameterized as in (5.2) with $p(\theta,t)$ chosen as

\begin{align*} p(\theta,t)=\theta_1t^2+\theta_0t.\end{align*}

The policy is parameterized as (5.1) with suitable initial values selected based on the problem context. From Theorem 3, we know that the optimal values of $\psi_0$ and $\psi_2$ are reciprocals of each other, and only once policy update is needed to get the optimal value of $\psi_2$ . Therefore, before starting training, we can perform a pretraining to get optimal $\psi_2$ to reduce the amount of subsequent training. We examine the behavior of the mean value under Nash equilibrium when Agent 1 adopts a normal distribution and Agent 2 follows an exponential distribution. The corresponding Choquet regularizers are specified in Equations (6.1) and (6.2). Figure 3 presents the mean of the discounted value invested in the risky asset under Nash equilibrium, evaluated at the market state $Y(t) = 0.273$ , while Figure 4 illustrates the learning trajectory of $\psi_2$ , which captures the evolution of the variance. By comparing the mean and variance of learned policy with the true policy in Figures 3 and 4, we observed that our experimental results closely approximate the theoretical values. This underscores the effectiveness of our approach in approximating Nash equilibrium solutions. Owing to the asymmetry and heavier tail of the exponential distribution compared to the normal distribution, Figure 4 shows that the trajectory under exponential-distribution-based exploration exhibits greater volatility, though it still converges to the true value. From this standpoint, the normal distribution may appear to yield more stable learning dynamics. However, this interpretation should be viewed with caution, as the regularizers used here capture agents’ subjective preferences rather than purely statistical features of the distributions.

Table 2.

Parameter settings for the algorithm.

Figure 3.

The mean value of Nash equilibrium.

Figure 4.

Convergence of the learned $\psi_2$ .

For Black-Scholes model, policy evaluation and policy updates can still be implemented using the methods mentioned above. The difference is that the quantile function of the policy $\boldsymbol{\Pi}_i$ is given by

\begin{align*} Q_{\Pi_i^{\Psi}(t)}(p)=k_i\mu_j(t)+\psi_0 \dfrac{y}{\gamma_i}+\dfrac{\lambda_i}{\gamma_i \psi_1^2} h_i'(1-p),\end{align*}

and the value function is parameterized as

\begin{align*} V_i^{\Theta}(t,\hat x_i)=\hat x_i+\theta_i^{V}(T-t)\ \text{ and }\ g_i^{\Theta}(t,\hat x_i)=\hat x_i+\theta_i^{g}(T-t).\end{align*}

From Corollary 1, we know that the optimal strategy does not change over time and the optimal value of $1/\psi_0$ and $\psi_1$ is $\sigma$ . Thus we are actually learning the optimal $\sigma$ . We plot the density function of the equilibrium policy in Figure 5 when $Y(t)=0.273$ and $\sigma=0.15$ . We can see that the output of our algorithm closely approximates the theoretical value. Note that the “true” value or policy in Figures 3–5 is computed under the assumption of full knowledge of all model parameters shown in Tables 1 and 2.

Figure 5.

Learned policy for Black-Scholes model.