This paper introduces a novel stochastic control framework to enhance the capabilities of automated investment managers, or robo-advisors, by accurately inferring clients' investment preferences from past activities. Our approach leverages a continuous-time model that incorporates utility functions and a generic discounting scheme of a time-varying rate, tailored to each client's risk tolerance, valuation of daily consumption, and significant life goals.  We address the resulting time inconsistency issue through state augmentation and the establishment of the dynamic programming principle and the verification theorem. Additionally, we provide sufficient conditions for the identifiability of client investment preferences. To complement our theoretical developments, we propose a learning algorithm based on maximum likelihood estimation within a discrete-time Markov Decision Process framework, augmented with entropy regularization. We prove that the log-likelihood function is locally concave, facilitating the fast convergence of our proposed algorithm. Practical effectiveness and efficiency are showcased through two numerical examples, including Merton's problem and an investment problem with unhedgeable risks.

Our proposed framework not only advances financial technology by improving personalized investment advice but also contributes broadly to other fields such as healthcare, economics, and artificial intelligence, where understanding individual preferences is crucial.

§ INTRODUCTION

Automated investment managers, commonly known as robo-advisors, have emerged as a modern alternative to traditional financial advisors in recent years [Capponi et al., 2022, D’Acunto et al., 2019, Rossi and Utkus, 2020]. The effectiveness and viability of robo-advisors depend significantly on their ability to provide customized financial guidance tailored to the unique needs of each client. To provide impactful personalized advice, two critical steps must be undertaken:  first, accurately estimate the client’s investment preferences, and second, formulate investment recommendations that align with these preferences. This paper focuses on the first step, involving a detailed analysis of the client's investment preferences.

More often than not, it is difficult for the automated investment manager to have full access to clients' investment preferences. Therefore, it is worth exploring whether it is possible to infer relevant information by observing the clients' past investment activities. On the other hand, inferring a client's investment preferences is typically challenging, as it involves several complex aspects that vary from individual to individual. For example, clients may have short-term or long-term investment objectives [Jin and Yu Zhou, 2008]. Additionally, they might exhibit varying utility functions [Musiela and Zariphopoulou, 2006, Zariphopoulou, 2001], reflecting distinct risk tolerance related to profit-and-loss (PnL) outcomes and valuation of daily consumption. Furthermore, individuals often demonstrate diverse time preferences in terms of the trade-off between immediate and deferred outcomes [Bjork and Murgoci, 2010]. Finally, clients may have specific life goals [Capponi and Zhang, 2024], such as saving for their children's education or building a retirement nest egg, rather than focusing solely on generating the highest possible portfolio return or beating the market.

The inference of preferences in sequential decision-making is a critical component not only for financial investments but also in other fields, leveraging insights into individual behaviors to optimize decisions and predict outcomes. In economics, utility functions are inferred to model consumer behavior, guiding businesses in product development and pricing strategies [Derbaix and Abeele, 1985, Shin and Yu, 2021].  Healthcare professionals use inferred utility functions to evaluate patient preferences regarding different treatment options, which is essential for effective healthcare management and policy-making [Chewning et al., 2012, Richesson and Vehik, 2010]. Additionally, in artificial intelligence, particularly in areas like reinforcement learning (RL) and game theory, inferring utility functions helps in designing algorithms that can predict and mimic human decision-making processes, enhancing the interaction between humans and machines [Boularias et al., 2011, Christiano et al., 2017].

Our framework, results, and contributions.
We propose a novel stochastic control framework in continuous time that incorporates all the aforementioned investment preferences. This framework includes two utility functions that allow the client to define their risk tolerance related to the PnL outcomes and their valuation of daily consumption. Additionally, it allows for a generic discounting scheme under a time-varying rate, enabling the clients to balance immediate and deferred outcomes. This time-varying discounting scheme further incorporates specific life goals by assigning greater importance to times of significant expenditures, such as college tuition fees for children. Lastly, we address control problems on both finite-time and infinite-time horizons to accommodate clients' preferred investment duration. The control problem is time inconsistent under the generic time-varying discounting scheme. We address this issue by state augmentation to account for the cumulative discount rate. We study the well-definedness of the augmented control framework by establishing the regularity of the value function, the dynamic programming principle (DPP), and the verification theorem (see Propositions <ref>, <ref>, <ref> for finite-time horizon and Propositions <ref>, <ref>, <ref> for infinite-time horizon). In addition, we identify sufficient conditions for identifying both the utility functions and the discounting scheme by solely observing the optimal policies provided by the client (see Theorem <ref> for finite-time horizon and Theorem <ref> for infinite-time horizon).

To complement the above theoretical framework, we propose an inference procedure based on maximum likelihood estimation. To demonstrate the effectiveness of this procedure, along with tractable theoretical guarantees, we focus on a specific case within the discrete-time Markov Decision Process (MDP), featuring Shannon's entropy regularization over an infinite-time horizon. The discrete-time MDP is especially relevant in the context of statistical inference and machine learning. The entropy term encourages full exploration of the state-action space and simultaneously introduces smoothness into the analysis [Haarnoja et al., 2017]. We employ a parametric framework where the client uses an exponential discounting scheme, parameterized by $\bar{\rho}$, and a utility function parameterized by $\bar{\theta}$. Both sets of parameters are unknown to the automated investment manager. Mathematically, we show that the true preference parameter $(\bar{\rho},\bar{\theta})$ is a stationary point of the log-likelihood function and the log-likelihood function is locally concave near $(\bar{\rho},\bar{\theta})$; see Proposition <ref> and Theorem <ref>. This landscape property facilitates the design of a gradient-based algorithm to update the inferred preference parameter. We demonstrate the promising performance of our algorithm through two examples–Merton's problem and an investment problem under unhedgeable risks.

Considering the wide-ranging applications and the versatility of our proposed framework, we use the term “inference agent” instead of “automated investment manager” to describe the individual who interacts with the clients and infers their preferences.

Related literature and comparisons to our results. Our developments are associated with several lines of literature as follows.

Utility inference.
Back in 1964, Kalman [Kalman, 1964] asked the question of whether it is possible to recover the quadratic cost by observing an optimal linear policy; a similar question was also considered by Boyd et al., 1994. In fact, economists have long been interested in such questions within the context of determining utility functions from observations  such as Samuelson, 1948 and Richter, 1966. For instance, Keeney and Raiffa, 1976 studied the proper rank of actions based on some deterministic evaluations under a static setting. Sargent, 1978 later extended this question into a dynamic setting where the actions were specified as labor demand and evaluations as wages. Dybvig and Rogers, 1997 paid special attention to the recoverability or identifiability of utility and showed that Von Neumann-Morgenstern preferences over terminal consumption can be inferred from wealth process of a discrete-time, binomial model or continuous-time Gaussian model.

Cox et al., 2014 studied the utility inference problem for the optimal consumption and allocation of wealth in continuous time by observing the actions of the client. The authors observed that there are infinitely many utility functions generating a given consumption pattern when the dynamic is deterministic and the consumption and investment strategies have to satisfy a consistency condition in the stochastic setting.
El Karoui and Mrad, 2021 took a “forward-looking” perspective of the connection between the observable process \(\{\mathscr{X}_t(x)>0:\mathscr{X}_0(x)=x>0\}\) (i.e., the characteristic process) and the corresponding utility process \(\{U(t,z):z>0,U(0,z)=u(z)\}\) (i.e., the dynamic utility); see the concept forward utility proposed by Musiela and Zariphopoulou, 2007. Different than the backward-looking perspective where the connection between the observable and utility is governed by some Markovian decision-making rule, the authors interpreted such a connection as the martingale property of the process \(\{U(t,\mathscr{X}_t(x)\}\), since Markov property no longer existed under the forward-looking viewpoint. To fully explore the concavity of utilities, the authors introduced an adjoint process of \(\mathscr X\), \(\{Y_t(y):Y_0(y)=y\}\), representing the decreasing marginal utility \(\{U_z(t,z)\}\) so that \(\{Y_t(u_z(x))=U_z(t,\mathscr X_t(x))\}\). Given the initial utility \(u\), the observable process \(\mathscr X\) and its adjoint process \(Y\), the authors fully characterized the martingale dynamic utility and its dual form via the Itô-Ventzel formula and showed that they are solutions to some Hamilton-Jacobi-Bellman-type stochastic partial differential equations; this set of analytical tools was introduced in [Nicole and Mohamed, 2013] and [El Karoui et al., 2019]. In [El Karoui et al., 2024], the authors also extended the result of [El Karoui and Mrad, 2021] to allow an exogenous default time \(\tau\).

In recent years, utility inference has been integrated with machine learning to embrace the potential of the big data era (and the progress is summarized in the next paragraph). In addition, inference problems in sequential decision-making for modern applications are more complex than inferring solely the utility function. Other preferences such as time preferences and specific investment goals should also be included, leading to the main formulation of our paper.

Theory of inverse optimal control. Inverse optimal control aims at inferring the underlying reward function that motivates the observed behavior of a rational agent in a sequential decision-making framework; within the context of MDP, inverse optimal control is also known as inverse reinforcement learning (IRL). In this area,  Ng et al., 2000 considered a particular setting that the true reward function is some linear combination of several action-free basis functions and that the true reward function maximally distinguishes the observed policy from the rest. They reformulated this question into a constrained linear programming problem eventually leading to a well-defined solution. In [Abbeel and Ng, 2004], the reward was assumed to be a linear combination of several features that best distinguish the demonstrated policy from other policies. The key assumption in both works is that the true reward function should maximize the margin between observations and the other policies. It also played a central role in the model of the well-known GAIL (generative adversarial imitation learning) algorithm [Ho and Ermon, 2016]. Other than the “maximum margin” setting, another commonly adopted setting in IRL is to assume that an observed randomized policy should maximize the causal entropy of an underlying regularized MDP. For instance, Ziebart et al., 2008 studied the maximum entropy IRL based on known features. They assumed that the reward is a linear function of such features. Ziebart, 2010
extended this approach to a selected set of non-linear rewards; see also [Levine et al., 2011] and [Boularias et al., 2011] for similar settings. Wulfmeier et al., 2015 followed this approach but with rewards represented by neural networks. Finn et al., 2016 combined the idea of adversarial training and IRL. They trained a discriminator to recover the reward function. Reddy et al., 2019 proposed a soft Q imitation learning algorithm to imitate the expert's policy by learning her Q function. Garg et al., 2021 proposed an algorithm to learn the soft Q function which implicitly represents both the reward function and the policy. Zeng et al., 2022 adopted the maximum likelihood estimator and showed that their algorithm converges to a stationary point under a finite-time guarantee.  Back to our preference inference problem, since it is to infer the utility functions and the time preferences of the client simultaneously, these existing IRL algorithms are not directly applicable. Such a multi-facet inference problem motivates our main algorithm. Furthermore,
we are able to provide a loss landscape analysis that facilitates fast convergence of our proposed algorithm; see Proposition <ref> and Theorem <ref>.

Identifiability issues in IRL. In 1998, Russell [Russell, 1998] pointed out the ill-posedness of inverse optimal control or IRL problems under a generic setting. Both the “maximum margin” and the “maximum entropy” settings mentioned above are reasonable assumptions to ameliorate this ill-posedness. Nonetheless, without prior access to the underlying true reward function, it is difficult to verify either one of them. To guarantee identifiability in IRL, alternative and more verifiable conditions and assumptions are required. Under an entropy regularized MDP setting, Cao et al., 2021 pointed out two possible remedies for the identifiability issue. One way is to provide additional observations of the same agent (i.e., keeping the underlying reward function the same) under different environments; see also a repeated IRL setting proposed in [Amin and Singh, 2016] and [Amin et al., 2017]. It was shown in [Cao et al., 2021] that under proper technical conditions on the transition kernels, observations from two distinct environments would suffice. Another approach is to provide additional structural assumptions on the MDP environment or the family of candidate reward functions based on prior domain knowledge; see also the identification of an action-free reward in [Fu et al., 2018]. Both Cao et al., 2021 and Kim et al., 2021 provided sufficient structural conditions for the MDP environment that guarantee identifiability.

However, as pointed out by Schlaginhaufen and Kamgarpour, 2023, the identifiability may no longer hold without the entropy regularization. In addition, the majority of these previous studies rely on the full disclosure of the MDP environment, including the transition kernel, time horizon, and the rate of an exponential discounting scheme. Though Dong and Wang, 2024 provided a mathematical formulation and an algorithm for the partial information setting, it remains to be explored whether identifiability of both the unknown MDP information and the true reward function is viable. In this paper, we establish such identifiability for our preference inference problem, which is also one of the major theoretical contributions; see Theorems <ref> and <ref>.

Time inconsistency in stochastic control. Unlike assuming an exponential discounting scheme for the client, a general discounting scheme may lead to a time-inconsistent policy. In economics, one of the earliest studies on the inconsistency in dynamic utility maximization is [Strotz, 1955], where the optimality of the problem derived today is different from that of tomorrow due to some non-exponential discounting mechanism. Later Pollak, 1968 proposed a game-theoretic consistent planning approach for the discrete-time problem, where the game is among decision makers at different time steps and the optimal decision path is considered to be the Nash equilibrium. There has been a line of works following this consistent planning approach under both discrete- and continuous-time settings; see, for instance, [Björk and Murgoci, 2014, Björk et al., 2017, Ekeland and Lazrak, 2010, Hu et al., 2012, Hu et al., 2017, Yong, 2012], and more recently, [Dai et al., 2023, Hernández and Possamaï, 2023]. Apart from this game-theoretic approach, Karnam et al., 2017 introduced the idea of “dynamic utility” to a family of time-inconsistent optimization problems over a finite-time horizon. By modeling the utility as the solution to a backward stochastic differential equation (BSDE), the DPP could be revived. For an infinite-time horizon setting which is suitable to model a long-run investment planning problem though, this BSDE approach can no longer be applied. Hence we propose a different way to revive DPP; see Propositions <ref> and <ref> in Section <ref>.

Robo-advising. Robo-advising has emerged over the last two decades as an alternative to traditional human financial advising, addressing limitations such as the human advisors' limited knowledge and high service fees [Capponi et al., 2022, D’Acunto and Rossi, 2021, D’Acunto et al., 2019]. Here, we mainly review some papers that explore the machine learning and inference aspects of this subject. The first RL algorithm for a robo-advisor was proposed by Alsabah et al., 2021, where the authors designed an exploration-exploitation algorithm to learn a constant risk appetite parameter and then applied a follow-the-leader type of algorithm to invest. Wang and Yu, 2021 introduced a framework consisting of two agents: the first, an inverse portfolio optimization agent, infers a risk preference parameter and the corresponding expected return; the second aggregates the learned information to formulate a new multi-period portfolio optimization problem solved by deep learning. To transcend the rather single-facet inference settings above, the theoretical framework and the numerical procedure in our paper are designed to capture the multiple investment needs of a client.

§ CONTINUOUS-TIME FRAMEWORK

In this section, we study a continuous-time framework of the joint consumption-allocation problem of an investing client. Her wealth consists of a risk-free asset and a risky asset. What distinguishes this framework from the classical ones is that the client holds a general preference of time, that is, the discounting scheme is not necessarily exponential. This could possibly lead to time-inconsistent decision-making. First, for the optimal control problem, we analyze the time-inconsistent dynamical decision-making problem for such a client, assuming the client's utility functions of consumption and wealth as well as her time preference are fully disclosed. The optimal decision relies on reviving a suitable DPP under this framework. Then, for the inverse optimal control problem, we establish an identifiability result for both the utility functions and the time preference of the client, assuming instead her optimal joint consumption-allocation plan is disclosed. Such an identifiability result provides inspirations for the algorithm proposed in Section <ref>.

§.§ Finite-time Horizon

We first focus on a finite-time horizon setting, with a decision horizon $T$, to address scenarios where the client has a short-term investment plan.

Market dynamics and client's wealth.  Let $(\Omega, \mathcal{F}, \mathbb{F} = (\mathcal{F}_t)_{t\ge 0},\mathbb{P})$ be a filtered probability space, supporting a one-dimensional $\mathbb{F}$-Brownian motion $W$. Assume there is a bond and a stock in the investment universe. The price of the bond follows

\begin{eqnarray}\label{eq:bond_price}
{\rm d} S_t^0 = r {\rm d}  t,
\end{eqnarray}

and the price of the stock follows

\begin{eqnarray}\label{eq:asset_price}
{\rm d }S_t = S_t(\mu {\rm d}t + \sigma {\rm d}W_t).
\end{eqnarray}

Assume the client choose an allocation process \(\balpha=\{\alpha_t\}_{t\in[0,T]}\)  and a consumption process \(\bc=\{c_t\}_{t\in[0,T]}\) with \(c_t\geq0\).  Namely, the client allocates $\alpha_t$ proportion of wealth to the stock and $1-\alpha_t$ proportion of wealth to the bond at time $t$. In addition,  the client is also making consumption $c_t$ to achieve certain satisfaction in life.

Fixing a sufficiently large constant \(M\in\R^+\) and introducing a compact space \(\mathcal{K}=[-M,M]\times[0,M]\), define

\begin{equation}
\label{eq:adm-set}
\begin{aligned}
\calA\coloneqq&\left\{(\balpha,\bc)\{(\alpha_t,c_t)\}_{t\geq0}\,\biggl|\,(\alpha_t,c_t)\in\mathcal{K},\, (\alpha_t,c_t)\in\bar\calF_t\coloneqq\sigma\left(\sigma(\beta_s,0\leq s\leq t)\times\mathcal{F}^W_t\right),\right.\\
&\hspace{100pt}{}\E\left[|\ctrldX_t|^2\right]<\infty,\,\forall t\geq0\biggl\}
\end{aligned}
\end{equation}

as the admissible set of all possible joint consumption-allocation processes. Hence the wealth process follows:

\begin{eqnarray}
\end{eqnarray}

Client's preference.
In the finite-time horizon, the preference of the client can be characterized by a pair of utility functions and a discount scheme.
More specifically, consider utility functions $U_1, U_2$ that belong to the following class

\begin{equation}
\label{eq:gen-util}
\begin{aligned}
\mathcal U\coloneqq\biggl\{U:\R\to[-\infty,+\infty)\biggl|\,\,&U\text{ is strictly positive, increasing and concave on }(0,+\infty),\\
&\text{there exists a sufficiently large constant }C\in\R^+\text{ such that}\\
&|U(x)|\leq C\left(1+x^2\right)\text{ for all }x\in(0,\infty),
\\&
U\in\calC^2((0,+\infty)),\,U(x)=U(0)\text{ for }x\leq 0\biggl\}.
\end{aligned}
\end{equation}

Here $U_1$ quantifies the the client's evaluation regarding the consumption whereas $U_2$ quantifies her evaluation regarding the terminal wealth at the end of the investment plan.

General discounting scheme. We are particularly interested in a client that is subject to a general discounting scheme \(\beta=\{\beta_t=\beta(t)\}_{t\geq0}\), where

* \(\beta_t\in[0,1]\) for all \(t\in[0,T]\); and

* there exists \(\dot{\beta}:[0,\infty)\to\R\) such that \(\dot{\beta}\) is bounded and integrable on \([0,t]\) with \(\beta_t=\int_0^t\dot{\beta}_sds+\beta_0\) for any \(t>0\).

Such a discounting scheme \(\{\beta_t\}_{t\geq0}\) reflects a generic time preference of the client. A time-varying discounting rate could account for different levels of appreciation for the immediate outcome and the delayed fulfillment. It could also provide the flexibility of assigning greater importance to times of significant expenditures, such as college tuition for children and down-payment of a house.

Then for any \((t,x,z)\in[0,T]\times\R\times[0,1]\), define the total reward as

\begin{equation}
\label{eq:acc-util}
\end{equation}

subject to the wealth process (<ref>) and

\begin{align}
\dd \beta_t=\dot{\beta}_t \dd t.\label{eq:gen-discnt}
\end{align}

For any \((x,z)\in\R\times[0,1]\), define the value function as follows,

\begin{equation}
\label{eq:gen-value}
V(t,x,z)=\sup_{(\balpha,\bc)\in\calA}J(t,x,z,\balpha,\bc),\,t\in[0,t);\ \ V(T,x,z)=zU_2(x).
\end{equation}

subject to (<ref>) and (<ref>).

In this section, we also aim to recover the DPP to the above time-inconsistent utility optimization problem (<ref>), where the time-inconsistency is particularly due to the general discounting scheme. We take a different approach than the BSDE characterization of dynamic utility in [Karnam et al., 2017]; instead, we extend the state space to incorporate the discounting process (similar to [Bäuerle and Rieder, 2014]) and then re-establish DPP accordingly.

§.§.§ Preliminary Analysis

First, we establish the well-definedness of the control problem (<ref>)-(<ref>) and introduce some analytical properties associated with it.

Assume that \(U_1,U_2\in\mathcal{U}\). Moreover, assume that \(U_1(0)=0\) and \(U_2(0)=-\infty\). For any \((t,x,z)\in[0,T]\times(0,\infty)\times[0,1]\), if the policy $\balpha^*,\bc^*$ satisfies that \(J(t,x,z,\balpha^*,\bc^*)=V(t,x,z)\), then it holds almost surely that

\begin{eqnarray}
\ctrldX[\balpha^*,\bc^*]_s\in(0,\infty)\quad \text{for all}\quad s\in[t,T],
\end{eqnarray}

where \(\ctrldX[\balpha^*,\bc^*]\) solves (<ref>) on \([t,T]\) given \((\balpha,\bc)=(\balpha^*,\bc^*)\) and \(\ctrldX_t=x\).

For any \((\balpha,\bc)\in\calA\), we have \(\ctrldX_s\leq \ctrldX[\balpha,\pmb{0}]_s\) for all \(s\in[t,T]\) almost surely. Notice that \(\ctrldX[\balpha,\pmb{0}]_{t'}=\ctrldX[\balpha,\pmb{0}]_t\exp\left\{\int_t^{t'}\alpha_l(\mu-r)+r-\frac{\sigma^2\alpha_l^2}{2}dl+\int_t^{t'}\sigma\alpha_ldW_l\right\}\) for \({t'}\geq t\). If \(\ctrldX_{t'}\leq0\), then
\[\ctrldX_T\leq\ctrldX[\balpha,\pmb{0}]_T\leq0,\]
and hence \(J(t,x,z,\balpha,\bc)=-\infty\). On the other hand,

\[J(t,x,z,\balpha,\pmb{0})=\beta_T\E{U_2\left(x\exp\left\{\int_t^{T}\alpha_l(\mu-r)+r-\frac{\sigma^2\alpha_l^2}{2}dl+\int_t^{T}\sigma\alpha_ldW_l\right\}\right)}>0.\]

Then if \(J(t,x,z,\balpha^*,\bc^*)=V(t,x,z)\), then \(\ctrldX[\balpha^*,\bc^*]_s\in(0,\infty)\) for all \(s\in[t,T]\) almost surely.

Assume that \(U_1,U_2\in\mathcal{U}\). Moreover, assume that \(U_1(0)=0\) and \(U_2(0)=-\infty\). Then it holds that the value function \(V:[0,T]\times\R\times[0,1]\to[-\infty,+\infty)\) defined in (<ref>) is strictly concave and strictly increasing in \(x\in(0,\infty)\) given any \((t,z)\in[0,T]\times(0,1]\).

Fix any \((t,x,z)\in[0,T]\times(0,\infty)\times(0,1]\).

* Strictly concave and positive. Take \(y\in(0,\infty)\setminus\{x\}\) and \(\lambda\in(0,1)\). Define \(u=\lambda x+(1-\lambda)y\) and \(u\in(0,\infty)\). Take any \((\balpha^x,\bc^x),(\balpha^y,\bc^y)\in\calA\) and define \((\balpha^u,\bc^u)=\lambda(\balpha^x,\bc^x)+(1-\lambda)(\balpha^y,\bc^y)\). Then it immediately follows that \((\balpha^u,\bc^u)\in \calA\). Let \(\ctrldX[\balpha^u,\bc^u]\) (resp. \(\ctrldX[\balpha^x,\bc^x]\) or \(\ctrldX[\balpha^y,\bc^y]\)) be the solution to the SDE (<ref>) over \([t,T]\) given \((\balpha,\bc)=(\balpha^u,\bc^u)\) (resp. \((\balpha,\bc)=(\balpha^x,\bc^x)\) or \((\balpha,\bc)=(\balpha^y,\bc^y)\)) and \(\ctrldX_t=u\) (resp. \(\ctrldX_t=x\) or \(\ctrldX_t=y\)). Then we have
\[\ctrldX[\balpha^u,\bc^u]_s = \lambda\ctrldX[\balpha^x,\bc^x]_s+(1-\lambda)\ctrldX[\balpha^y,\bc^y]_s,\ \ s\in[t,T].\]

By Lemma <ref>, we can assume that both \(\ctrldX[\balpha^x,\bc^x]_s\) and \(\ctrldX[\balpha^y,\bc^y]_s\) are strictly positive for \(s\in[t,T]\) almost surely. Since \(U_1,U_2\in\mathcal{U}\), then

\[J(t,u,z,\balpha^u,\bc^u)>\lambda J(t,x,z,\balpha^x,\bc^x)+(1-\lambda)J(t,y,z,\balpha^y,\bc^y)>0.\]

Taking the supremum over both \((\balpha^x,\bc^x)\) and \((\balpha^y,\bc^y)\),

\[V(t,u,z)>\lambda V(t,x,z)+(1-\lambda)V(t,y,z)>0.\]

* Strictly increasing. Fix any \(\Delta x>0\) take any \((\balpha,\bc)\in\calA\) such that \(\ctrldX_s>0\) for \(s\in[t,T]\) almost surely. Let \(\widehat{\ctrldX}\) be the solution to (<ref>) given \(\ctrldX_t=x+\Delta x\). Then,

\[\Delta_T :=\widehat{\ctrldX_T}-\ctrldX_T=\Delta x\exp\left\{\int_t^{T}\alpha_l(\mu-r)+r-\frac{\sigma^2\alpha_l^2}{2}dl+\int_t^{T}\sigma\alpha_ldW_l\right\}>0\text{ a.s.},\]

and therefore

\[J(t,x+\Delta x,z,\balpha,\bc)-J(t,x,z,\balpha,\bc)=\beta_T\E\left[U_2\left(\widehat{\ctrldX_T}\right)-U_2(\ctrldX_T)\right]>0.\]

Hence, \(V(t,x+\Delta x,z)>V(t,x,z)\).

Having established some preliminary properties of the value function, we first show a necessary condition for the value function (<ref>).

Take the same assumptions on \(U_1,U_2\) as in Lemma <ref>. For any \((t,x,z)\in[0,T)\times\R\times[0,1]\) and \(\tau\in\mathbb T_t\) where \(\mathbb T_t\) denotes all \(\{\bar\calF_t\}_{t\geq0}\)-adapted stopping times \(\tau\) such that \(\tau\in[t,T]\) a.s., the value function \(V\) defined in (<ref>) satisfies

\begin{equation}
\label{eq:gen-dpp}
\tag{DPP}
V(t,x,z)=\sup_{(\balpha,\bc)\in\calA}\E\left[\int_t^{\tau}\beta_sU_1(c_s)ds+V(\tau, \ctrldX_{\tau},\beta_{\tau})\biggl|\ctrldX_t=x,\beta_t=z\right],
\end{equation}

with \(V(T,x,z)=zU_2(x)\).

Fix any \((t,x,z)\in[0,T)\times\R\times[0,1]\), \((\balpha,\bc)\in\calA\) and \(\tau\in\mathbb T_t\). We have

\begin{align*}
J(t,x,z,\balpha,\bc)&= \E\left[\int_t^{\tau}\beta_sU_1(c_s)ds+\int_{\tau}^{T}\beta_sU_1(c_s)ds+\beta_TU_2(\ctrldX_T)\biggl|\ctrldX_t=x,\beta_t=z\right]\\
&\hspace{20pt}+\E\left[\E\left[J(\tau,X,Z,\balpha,\bc)\biggl|\tau, X=\ctrldX_{\tau},Z=\beta_{\tau}\right]\biggl|\ctrldX_t=x,\beta_t=z\right].
\end{align*}

By the definition given by (<ref>), for any \(\epsilon>0\) and \(\Delta t\in[0,T-t]\), there exists \((\balpha^{\epsilon,t+\Delta t,x,z},\bc^{\epsilon,t+\Delta t,x,z})\in\calA\) such that

\begin{equation}\label{eq:gen-upper}
J(t+\Delta t,x,z,\balpha^{\epsilon,t+\Delta t,x,z},\bc^{\epsilon,t+\Delta t,x,z})>\sup_{(\balpha,\bc)\in\calA}J(t+\Delta t,x,z,\balpha,\bc)-\epsilon=V(t+\Delta t, x,z)-\epsilon,
\end{equation}

\begin{equation}\label{eq:gen-lower}
J(t+\Delta t, x,z,\balpha,\bc)\leq J(t+\Delta t,x,z,\balpha^{\epsilon,t+\Delta t,x,z},\bc^{\epsilon,t+\Delta t,x,z})\leq V(t+\Delta t, x,z).
\end{equation}

Then consider \((\bar\balpha,\bar\bc)=\{(\bar\alpha_s,\bar c_s)\}_{s\in[t,T]}\) such that

\begin{equation*}
(\bar\alpha_s,\bar{c}_s)=(\alpha_s,c_s)\mathbbm{1}\{\tau>s\}+ \left(\alpha^{\epsilon,\tau,\ctrldX_{\tau},\beta_{\tau}}_s,c^{\epsilon,\tau,\ctrldX_{\tau},\beta_{\tau}}_s\right)\mathbbm{1}\{\tau\leq s\},
\end{equation*}

\begin{equation*}
\beta_{\tau}=z+\int_t^{\tau}\dot{\beta}_sds,\ \
\ctrldX_{\tau}=x+\int_t^{\tau}d\ctrldX_s,
\end{equation*}

according to (<ref>) and (<ref>). Notice that \((\bar\balpha,\bar\bc)\in\calA\). By (<ref>) and (<ref>), for any \((\balpha,\bc)\in\calA\), we have

\[\begin{aligned}
V(t,x,z)&\geq J(t,x,z,\bar\balpha,\bar\bc)>\E\left[\int_t^{\tau}\beta_sU_1(c_s)ds\biggl|\ctrldX_t=x,\beta_t=z\right]\\
&\hspace{50pt}+\E\left[V(\tau, \ctrldX_{\tau},\beta_{\tau})\biggl|\ctrldX_t=x,\beta_t=z\right]-\epsilon
\end{aligned}\]

for any \(\epsilon>0\), and

\[J(t,x,z,\balpha,\bc)\leq \E\left[\int_t^{\tau}\beta_sU_1(c_s)ds+V(\tau, \ctrldX_{\tau},\beta_{\tau})\biggl|\ctrldX_t=x,\beta_t=z\right].\]

It follows that for any \((t,x,z)\in[0,T)\times\R\times[0,1]\) and \(\tau\in\mathbb T_t\),

\begin{equation*}
V(t,x,z)=\sup_{(\balpha,\bc)\in\calA}\E\left[\int_t^{\tau}\beta_sU_1(c_s)ds+V(\tau, \ctrldX_{\tau},\beta_{\tau})\biggl|\ctrldX_t=x,\beta_t=z\right].
\end{equation*}

For any \(\alpha\in\R\) and \(c\in\R^+\), define the following operator

\[\infgen \phi(t,x,z)=\left\{\left[\alpha(\mu-r)x\right]-c\right\}\partial_x\phi(t,x,z)+\frac{\sigma^2\alpha^2}{2}x^2\partial_x^2\phi(t,x,z),\]

for any test function \(\phi\in\calC^{\infty}_b([0,T)\times\R\times\R^+)\bigcap\calC_b^0([0,T)\times\R\times\R^+)\). Following the DPP under a generic discounting scheme (<ref>), we have the following Hamilton-Jacobi-Bellman (HJB) equation

\begin{equation}
\label{eq:gen-hjb}
\tag{HJB}\left\{
\begin{aligned}
&\partial_tV(t,x,z)+\dot{\beta}_t\partial_zV(t,x,z)+rx\partial_xV(t,x,z)+\sup_{(\alpha,c)\in\mathcal K}\left\{zU_1(c)+\infgen V(t,x,z)\right\}=0,\ \ t\in[0,T);\\
\end{aligned}\right.
\end{equation}

The next result provides sufficient conditions for the value function in (<ref>) regarding classical solutions to (<ref>).

Take the same assumptions on \(U_1,U_2\) as in Lemma <ref>. Let \(w:[0,T]\times\R\times[0,1]\to\R\) be a function such that

\[w\in\calC^{1,2,1}\left([0,T]\times\R\times[0,1]\right)\bigcap\calC^0\left([0,T]\times\R\times[0,1]\right),\]

and there exists a constant \(C>0\) with

\[w(t,x,z)\leq C(1+|x|^2),\ \ \forall (t,x,z)\in[0,T]\times\R\times[0,1].\]

* Assume that for any \((\alpha,c)\in\mathcal K\),

\[
\left\{
\begin{aligned}
&\partial_tw(t,x,z)+\dot{\beta}_t\partial_zw(t,x,z)+rx\partial_xw(t,x,z)+zU_1(c)+\infgen w(t,x,z)\leq0,\\
&\hspace{180pt}\forall (t,x,z)\in[0,T)\times\R\times[0,1];\\
&w(T,x,z)\geq zU_2(x), \ \  \forall (x,z)\in\R\times[0,1].
\end{aligned}\right.
\]

Then \(w\geq V\) on \([0,T]\times\R\times[0,1]\).

* Assume further that there exists \(\hat{\alpha}:[0,T]\times\R\times[0,1]\to[-M,M]\) and \(\hat{c}:[0,T]\times\R\times[0,1]\to[0,M]\) such that

\[
\left\{
\begin{aligned}
&\partial_tw(t,x,z)+\dot{\beta}_t\partial_zw(t,x,z)+rx\partial_xw(t,x,z)+zU_1(\hat c(t,x,z))+\infgen[\hat{\alpha}(t,x,z),\hat{c}(t,x,z)]w(t,x,z)=0,\\
&\hspace{180pt}\forall (t,x,z)\in[0,T)\times\R\times[0,1];\\
&w(T,x,z)=zU_2(x), \ \  \forall (x,z)\in\R\times[0,1],
\end{aligned}\right.
\]

also, with \(\beta_t=\beta_0+\int_0^t\dot{\beta}_sds\in[0,1]\) for all \(t\in[0,T]\), the following SDE,

\[dX_t=
\left\{X_t\left[\hat{\alpha}(t,X_t,\beta_t)(\mu-r)+r\right]-\hat{c}(t,X_t,\beta_t)\right\}dt
\]

admits a unique solution \(\ctrldX[\hat{\balpha},\hat{\bc}]\) given \(X_0=x\) for any \(x\in\R\), and

\[\left(\hat{\balpha}=\left\{\hat{\alpha}_t\right\}_{t\in[0,T]}=\left\{\hat{\alpha}(t,\ctrldX[\hat{\balpha},\hat{\bc}]_t,\beta_t)\right\}_{t\in[0,T]}, \hat{\bc}=\left\{\hat{c}_t\right\}_{t\in[0,T]}=\left\{\hat{c}(t,\ctrldX[\hat{\balpha},\hat{\bc}]_t,\beta_t)\right\}_{t\in[0,T]}\right)\in\calA.\]

Then \(w= V\) on \([0,T]\times\R\times[0,1]\), with \((\hat{\balpha},\hat{\bc})\) being an optimal joint allocation-consumption process.

The assumptions on \(U_1\) and \(U_2\) guarantee a quadratic growth rate in \(x\). Consider arbitrary \((t,x,z)\in[0,T)\times\R\times[0,1]\) and \((\balpha,\bc)\in\calA\).

* Define

\[\tau_n\coloneqq\inf\left\{s\geq t\biggl|\int_t^s|\partial_xw(u,\ctrldX_u,\beta_u)|^2du\geq n\right\},\ \ \forall n\in\mathbb N^+.\]

Then we have \(\lim_{n\uparrow\infty}\tau_n\overset{a.s.}{=}\infty\) and the stopped process \(\left\{\int_t^{s\wedge\tau_n}\partial_xw(u,\ctrldX_u,\beta_u)dW_u\right\}_{s\in[t,T]}\) is a martingale for all \(n\in\mathbb N^+\). The for any \(s\in[t,T]\), by Itô's formula, we have

\begin{align*}
&.{}+r\ctrldX_u\partial_xw(u,\ctrldX_u,\beta_u)+\infgen[\alpha_u,c_u] w(u,\ctrldX_u,\beta_u)\Big\}\,du\\
\end{align*}

Therefore, taking expectations on both sides we have

\begin{align*}
&\E\left[w(s\wedge\tau_n,\ctrldX_{s\wedge\tau_n},\beta_{s\wedge\tau_n})\biggl|\ctrldX_t=x,\beta_t=z\right]=w(t,x,z)+\mathbb E\left[\int_t^{s\wedge\tau_n}\partial_tw(u,\ctrldX_u,\beta_u)\right.\\
&\hspace{20pt}\left.{}+\dot{\beta}_u\partial_zw(u,\ctrldX_u,\beta_u)+r\ctrldX_u\partial_xw(u,\ctrldX_u,\beta_u)+\infgen[\alpha_u,c_u] w(u,\ctrldX_u,\beta_u)du\biggl|\ctrldX_t=x,\beta_t=z\right]\\
&\hspace{20pt}\leq w(t,x,z)-\E\left[\int_t^{s\wedge\tau_n}\beta_uU_1(c_u)du\biggl|\ctrldX_t=x,\beta_t=z\right],
\end{align*}

where the well-posedness of \(\E\left[\int_t^{s\wedge\tau_n}\beta_uU_1(c_u)du\biggl|\ctrldX_t=x,\beta_t=z\right]\) is guaranteed by the quadratic growth rate condition on \(U_1\) and the fact that \((\balpha,\bc)\in\calA\). The quadratic growth rate assumption on \(w\) together with \((\balpha,\bc)\in\calA\) allows us to apply dominated convergence theorem and get

\[\begin{aligned}
\leq&w(t,x,z)-\E\left[\int_t^T\beta_uU_1(c_u)du\biggl|\ctrldX_t=x,\beta_t=z\right]\\
\implies&w(t,x,z)\geq\E\left[\int_t^T\beta_uU_1(c_u)du+\beta_TU_2(\ctrldX_T)\biggl|\ctrldX_t=x,\beta_t=z\right]=J(t,x,z,\balpha,\bc).\end{aligned}\]

Hence, \(w(t,x,z)\geq V(t,x,z)\) by taking the supreme of \((\balpha,\bc)\) over \(\calA\).

* Applying a similar localization-and-Itô argument as in the previous part, we have that for any \(s\in[t,T]\)

\begin{align*}
&\E\left[w(s,\ctrldX[\hat{\balpha},\hat{\bc}]_{s},\beta_{s})\biggl|X_t=x,\beta_t=z\right]=w(t,x,z)+\mathbb E\left[\int_t^{s}\partial_tw(u,\ctrldX[\hat{\balpha},\hat{\bc}]_u,\beta_u)\right.\\
&\hspace{20pt}\left.{}+\dot{\beta}_u\partial_zw(u,\ctrldX[\hat{\balpha},\hat{\bc}]_u,\beta_u)+r\ctrldX[\hat{\balpha},\hat{\bc}]_u\partial_xw(u,\ctrldX[\hat{\balpha},\hat{\bc}]_u,\beta_u)+\infgen[\hat{\alpha}_u,\hat c_u] w(u,\ctrldX[\hat{\balpha},\hat{\bc}]_u,\beta_u)du\biggl|\ctrldX[\hat{\balpha},\hat{\bc}]_t=x,\beta_t=z\right]\\
&\hspace{20pt}= w(t,x,z)-\E\left[\int_t^{s}\beta_uU_1(\hat{c}_u)du\biggl|\ctrldX[\hat{\balpha},\hat{\bc}]_t=x,\beta_t=z\right].
\end{align*}

In particular, when \(s=T\),

\[\begin{aligned}
\implies&w(t,x,z)=\E\left[\int_t^T\beta_uU_1(\hat{c}_u)du+\beta_TU_2(\ctrldX[\hat{\balpha},\hat{\bc}]_T)\biggl|\ctrldX[\hat{\balpha},\hat{\bc}]_t=x,\beta_t=z\right]=J(t,x,z,\hat{\balpha},\hat{\bc}).\end{aligned}\]

Then we have \(w(t,x,z)=J(t,x,z,\hat{\balpha},\hat{\bc})\leq V(t,x,z)\). Combined with the result from the previous part, we have \(w=V\) on \([0,T]\times\R\times\R^+\), with \((\hat{\balpha},\hat{\bc})\in\calA\) being a corresponding optimal joint allocation-consumption process.

Without assuming the existence of a classical solution to (<ref>), we could instead consider its viscosity solution.

Denote \(\calD=[0,T)\times(0,\infty)\times[0,1]\).

* A lower semi-continuous function \(\underline{v}:\calD\to\R\) is a viscosity subsolution to (<ref>) if for any \((t_0,x_0,z_0)\in\calD\) and any \(\phi\in\calC^{1,2,1}(\calD)\) such that

\[\min_{(t,x,z)\in B(t_0,x_0,z_0)}(\phi-\underline{v})(t,x,z)=(\phi-\underline{v})(t_0,x_0,z_0)=0\]

for some neighborhood \(B(t_0,x_0,z_0)\subset\calD\),

\begin{equation}
\label{eq:gen-subsoln}
-\partial_t\phi(t_0,x_0,z_0)-\dot{\beta}_t\partial_z\phi(t_0,x_0,z_0)-rx_0\partial_x\phi(t_0,x_0,z_0)-\sup_{(\alpha,c)\in\mathcal K}\Big\{z_0U_1(c)+\mathcal{L}^{\alpha,c}\phi(t_0,x_0,z_0)\Big\}\leq0.
\end{equation}

* An upper semi-continuous function \(\overline{v}:\calD\to\R\) is a viscosity supersolution to (<ref>) if for any \((t_0,x_0,z_0)\in\calD\) and any \(\psi\in\calC^{1,2,1}(\calD)\) such that

\[\max_{(t,x,z)\in B(t_0,x_0,z_0)}(\psi-\underline{v})(t,x,z)=(\phi-\underline{v})(t_0,x_0,z_0)=0\]

for some neighborhood \(B(t_0,x_0,z_0)\subset\calD\),

\begin{equation}
\label{eq:gen-supersoln}
-\partial_t\psi(t_0,x_0,z_0)-\dot{\beta}_t\partial_z\psi(t_0,x_0,z_0)-rx_0\partial_x\psi(t_0,x_0,z_0)-\sup_{(\alpha,c)\in\mathcal K}\Big\{z_0U_1(c)+\mathcal{L}^{\alpha,c}\psi(t_0,x_0,z_0)\Big\}\geq0.
\end{equation}

* A continuous function \(v:\calD\to\R\) is a viscosity solution to (<ref>) if it is both a viscosity subsolution and a viscosity supersolution to (<ref>).

Take the same assumptions on \(U_1,U_2\) as in Lemma <ref>. The value function \(V\) in (<ref>) is the unique viscosity solution to (<ref>) over the any  \(\bar{\mathcal D}=[0,T]\times \mathcal D_1\times \mathcal D_2\subset\calD\) with \(\mathcal D_i\) compact, \(i=1,2\).

First, notice that under the assumptions on admissible control specified in (<ref>) as well as those on the utility functions specified in (<ref>), the continuity of the value function \(V\) in (<ref>) over the domain \(\calD\) can be established following the classical results of [Linos, 1983, Lions, 1982]. Therefore, \(V\) is bounded and uniformly continuous on \(\bar{\calD}\). Combining Proposition <ref> and similar arguments of Itô's formula in its proof, the viscosity solution property in Definition <ref> can be established. The uniqueness result follows a classical comparison principal <cit.>.

Given \(V\) being a \(\mathcal{C}^{1,2,1}\left([0,T)\times(0,\infty)\times[0,1]\right)\bigcap\calC\left([0,T]\times(0,\infty)\times[0,1]\right)\), define the Hamiltonian as

\[H(t,x,z,\alpha,c,p,q)\coloneqq zU_1(c)+\left\{x\left[\alpha\mu+(1-\alpha)r\right]-c\right\}p+\frac{\sigma^2\alpha^2}{2}x^2q.\]

Then, we have that for any \((t,x,z)\in[0,T)\times(0,\infty)\in[0,1]\),

\[
(\alpha_t^*,c_t^*)=(\alpha^*(t,x,z),c^*(t,x,z))=\argmax_{(\alpha,c)\in\mathcal K}H\Big(t,x,z,\alpha,c,\partial_xV(t,x,z),\partial_x^2V(t,x,z)\Big),
\]

\begin{align}
\alpha^*_t&=-\frac{\partial_xV(t,x,z)}{\partial_x^2V(t,x,z)}\cdot\frac{(\mu-r)}{\sigma^2x}\wedge M,\label{eq:gen-opt-alpha}\\
\end{align}

note that by Lemma <ref>, \(\alpha_t^*>0\).

§.§.§ The Inverse Problem: Identifiability of the Utility Functions

In this section, we focus on the “inverse” problem with respect to the optimal asset allocation-consumption scenario and study the “identifiability” of the utility functions as well as the discounting scheme out of the optimal investment policies. More specifically, we assume the client provides her decision policies (in the sense of the allocation-consumption processes) to the inference agent.

Following practical protocols, we assume that the inference agent does not know the discounting scheme \(\dot{\beta}\) nor the utility functions \(U_1\) and \(U_2\). Nevertheless, the inference agent tries to infer these characteristic functions based on available information, namely the joint allocation-consumption process (i.e., control policy) provided by the client.

To start, let

\begin{equation}
\label{eq:gen-inv-opt-ctrl}
({\bar\alpha},{\bar c}):[0,T]\times(0,\infty)\times[0,1]\to\mathcal{K}
\end{equation}

be some allocation and consumption policies of a client such that

\[\left(\bar{\balpha}=\left\{\bar\alpha_t\right\}_{t\in[0,T]}=\left\{\bar\alpha(t,\ctrldX[\bar{\balpha},\bar{\bc}]_t,\beta_t)\right\}_{t\in[0,T]},\bar{\bc}=\left\{\bar c_t\right\}_{t\in[0,T]}=\left\{\bar c(t,\ctrldX[\bar{\balpha},\bar{\bc}]_t,\beta_t)\right\}_{t\in[0,T]}\right)\in\calA,\]

where \(\left(\left\{\beta_t\right\}_{t\in[0,T]},\left\{\ctrldX[\bar{\balpha},\bar{\bc}]_t\right\}_{t\in[0,T]}\right)\) solves

\begin{align*}
&d\ctrldX[\bar{\balpha},\bar{\bc}]_t=\left\{\ctrldX[\bar{\balpha},\bar{\bc}]_t\left[\bar\alpha\left(t,\ctrldX[\bar{\balpha},\bar{\bc}]_t,\beta_t\right)(\mu-r)+r\right]-\bar c\left(t,\ctrldX[\bar{\balpha},\bar{\bc}]_t,\beta_t\right)\right\}dt+\sigma \bar\alpha\left(t,\ctrldX[\bar{\balpha},\bar{\bc}]_t,\beta_t\right)\ctrldX[\bar{\balpha},\bar{\bc}]_tdW_t,
\end{align*}

for any given \((\beta_0,\ctrldX[\bar{\balpha},\bar{\bc}]_0)\in[0,1]\times(0,\infty)\), and

\[\left(\bar{\balpha},\bar{\bc}\right)\in\argmax_{(\balpha,\bc)\in\calA}J(t,x,z,\balpha,\bc),\ \ \forall(t,x,z)\in[0,T]\times(0,\infty)\times[0,1]\]

subject to (<ref>) and (<ref>) . Write

\begin{equation}
\label{eq:inv-value}
\bar{V}(t,x,z)=J(t,x,z,\bar\balpha,\bar\bc),\ \ \forall (t,x,z)\in[0,T]\times(0,\infty)\times[0,1].
\end{equation}

Note that the inference agent has full access to (<ref>).

Assume that

* \(\bar\alpha,\bar c\in\calC^{1,1,1}([0,T)\times(0,\infty)\times[0,1])\bigcap\calC^{0}([0,T]\times(0,\infty)\times[0,1])\);

* \(\bar\alpha(t,x,z)\in(0,M)\) for all \((t,x,z)\in[0,T]\times(0,\infty)\times[0,1]\);

* for any \((t,z)\in[0,T]\times[0,1]\),

\[\bar c(t,x,z)<x,\ \  \forall x>0;\]

* both \(\bar\alpha(T,\cdot,\cdot)\) and \(\bar c(T,\cdot,\cdot)\) are “\(z\)-free”, denoted by

\[\bar\alpha(T,x,z)\equiv\bar\alpha_T(x),\ \ \bar c(T,x,z)\equiv\bar c_T(x),\ \ \forall (x,z)\in(0,\infty)\times[0,1];\]

* for \(x\in(0,\infty)\), \(\bar\alpha_T(x)>0\), \(\bar c_T\) is invertible, and the following difference for any \((t,z)\in[0,T)\times(0,1]\),

\[\Delta(t,x,z)\coloneqq\int_x^1\frac{dy}{y\bar\alpha(t,y,z)}-\int_{\bar c_T^{-1}(\bar c(t,x,z))}^1\frac{dy}{y\bar\alpha_T(y)}\]

depends only on \((t,z)\), namely \(\Delta(t,x,z)\equiv\Delta(t,z)\).

Then both the discounting scheme characterized by \(\dot{\beta}\) and the utility functions \(U_i\in\mathcal U\) for \(i=1,2\), with \(U_1(0)=0\) and \(U_2(0)=-\infty\), are identifiable up to an affine transform.

This result is also consistent with the finding in [Cao et al., 2021] that the identifiability of the unknown utility function in an inverse optimal control problem is equivalent to the identifiability of the corresponding value function under the observed optimal policy. Assumptions <ref> and <ref> enunciate the precise dependency of value function \(\bar V\) in (<ref>) and the observed policy \((\bar\alpha,\bar c)\).

First, by Assumption <ref>, for all  \((t,x,z)\in[0,T]\times(0,\infty)\times[0,1]\), (<ref>) is equivalent to

\begin{equation}
\label{eq:gen-hjb-pde}
\left\{\begin{aligned}
&\hspace{100pt}-\frac{(r-\mu)^2}{2\sigma^2}\frac{\left[\partial_xV(t,x,z)\right]^2}{\partial_x^2V(t,x,z)}=zU^*_1\left(\frac{\partial_xV(t,x,z)}{z}\right),\ \ t\in[0,T);\\
\end{aligned}
\right.
\end{equation}

where \(U_1^*\) is the Legendre transform of the concave utility function \(U_1:[0,\infty)\to\R\),

\[U_1^*(\kappa)\coloneqq \inf_{c\in[0,M]}\left\{\kappa c-U_1(c)\right\},\ \ \forall\kappa\in\R.\]

Now, we construct the value function \(\bar V\) in (<ref>) from (<ref>). By (<ref>) and Assumption <ref>,

\[
\bar\alpha_T(x)=-\frac{U_2'(x)}{U_2''(x)}\frac{\mu-r}{\sigma^2x}\implies U_2(x)=k_1\int_1^x\exp\left\{\int_y^1\frac{\mu-r}{\sigma^2u\bar\alpha_T(u)}du\right\}dy+k_2,
\]

for some \(k_1,k_2>0\); in particular,

\[U_2'(x)=k_1\exp\left\{\int_x^1\frac{\mu-r}{y\bar\alpha_T(y)}dy\right\}.\]

By (<ref>), for any \(x>0\),

\[U_1'(\bar c_T(x))=U_2'(x)\implies U_1'(x)=U_2'(\bar c_T^{-1}(x))=k_1\exp\left\{\int_{\bar c_T^{-1}(x)}^1\frac{\mu-r}{\sigma^2y\bar\alpha_T(y)}dy\right\} \]

\[U_1(x)=k_1\int_1^x\exp\left\{\int_{\bar c_T^{-1}(y)}^1\frac{\mu-r}{\sigma^2u\bar\alpha_T(u)}du\right\}dy+k_3\]

for some \(k_3>0\).

For any \((t,z)\in[0,T)\times(0,1]\) and \(x>0\), by (<ref>)–(<ref>),

\[\begin{cases}
\bar\alpha(t,x,z)=-\frac{\partial_x\bar V(t,x,z)}{\partial_x[\partial_x\bar V(t,x,z)]}\frac{\mu-r}{\sigma^2x},\\
\partial_x\bar V(t,x,z)=zU_1'(\bar c(t,x,z));
\end{cases}\implies\partial_x \bar V(t,x,z)=K_1(t,z)\exp\left\{\int_x^1\frac{\mu-r}{\sigma^2y\bar\alpha(t,y,z)}dy\right\},\]

\[K_1(t,z)=k_1z\exp\{-\frac{\mu-r}{\sigma^2}\Delta(t,z)\}.\]

Rewrite (<ref>) as

\[
\partial_t\bar V(t,x,z)+\dot{\beta}_t\partial_z\bar V(t,x,z)=-\left\{\left[r+\frac{\bar\alpha(t,x,z)(\mu-r)}{2}\right]x-\bar c(t,x,z)\right\}\partial_x\bar V(t,x,z)-zU_1(\bar c(t,x,z)).
\]

Differentiating with respect to \(x\) on both sides, we have

\[\dot{\beta}_t=-\frac{\partial_t[\partial_x\bar V(t,x,z)]+\partial_x\left\{\left\{\left[r+\frac{\bar\alpha(t,x,z)(\mu-r)}{2}\right]x-\bar c(t,x,z)\right\}\partial_x\bar V(t,x,z)+zU_1(\bar c(t,x,z))\right\}}{\partial_z[\partial_x\bar V(t,x,z)]}.\]

We conclude the analysis on finite-time horizon by discussing a special case with an explicit solution.

Set $\beta_0=1$, $\dot\beta_t = 0$ $(0\leq t \leq T)$, $U_1(c) = 0$ and a CRRA (power) utility $U_2(x) = \frac{x^\theta}{\theta}$ with $0<\theta<1$.  Also set the constant $M$ such that $M > \frac{\mu-r}{\sigma^2}$.
In this case, we face a classic control problem that is time-consistent. Hence state augmentation is not necessary.  In addition, both the optimal control and the inverse problem have explicit representations. The goal here is to identify the parameter $\theta$ from the client.

Consequently, define the value function:

\begin{eqnarray}
V(t,x) = \sup_{\alpha \in \mathcal{A}}\mathbb{E}\Big[ U_2(X_T)\,\left|\,X_t = x\Big].\right.
\end{eqnarray}

The value function satisfies the following HJB equation:

\begin{eqnarray}
-\partial_t V - \sup_{\alpha \in \mathcal{A}} \Big[\mathcal{L}^\alpha V(t,x)\Big] = 0,
\end{eqnarray}

with boundary condition $V(T,x) = U_2(x)= \frac{x^\theta}{\theta}$,
where the generator is defined as $\mathcal{L}^\alpha V(t,x) = x\Big( \alpha \mu + (1-\alpha)r\Big)\partial_{x}V+ \frac{1}{2} x^2 \alpha^2\sigma^2 \partial_{x}^2 V$.
The optimal policy follows:

\begin{eqnarray}
\bar\alpha(t,x) = -M\vee \left(- \frac{(\mu-r)\partial_x V}{x \sigma^2 \partial_x^2 V}\right)\wedge M.
\end{eqnarray}

For now assume that $-M\leq \left(- \frac{(\mu-r)\partial_x V}{x \sigma^2 \partial_x^2 V}\right)\leq M$ (to be checked later). Then plugging it back into the HJB equation, we have

\begin{eqnarray}
-\frac{\partial V}{\partial t} = x r \partial_{x}V - \frac{1}{2} \frac{(\mu-r)(\partial_{x}V)^2}{\sigma^2 \partial^2_{x}V}
\end{eqnarray}

with boundary condition $V(T,x) =  \frac{x^\theta}{\theta}$.
We take the ansatz $V(t,x) = \phi(t)\frac{x^\theta}{\theta}$. Hence $\phi(t)$ satisfies:

\begin{eqnarray}
\phi'(t) = \rho \phi(t) = 0, \quad \phi(T)=1,
\end{eqnarray}

where $\rho = p\times \sup_{\alpha\in[-M,M]}\Big[ \alpha(\mu-r) + r -\frac{1}{2}a^2(1-p)\sigma^2\Big]$. Hence $\bar\alpha = \frac{
\mu-r}{\sigma^2(1-\theta)}\in[-M,M]$. In this case, it is obvious that condition $-M\leq \left(- \frac{(\mu-r)\partial_x V}{x \sigma^2 \partial_x^2 V}\right)\leq M$ is satisfied. Therefore we can recover the preference parameter by using $1-\frac{\mu-r}{\bar \alpha \sigma^2}$.

§.§ Infinite-time Horizon

Now we shift our focus to an infinite-time horizon setting that accommodates a long-run investment planning scenario.

Recall that the investing client is holding a general discounting scheme \(\beta=\{\beta_t\}_{t\geq0}\) where

* \(\beta_t\in[0,1]\) for all \(t\in[0,\infty)\) such that \(\lim_{t\to\infty}\beta_t=0\); and

* there exists \(\dot{\beta}:[0,\infty)\to\R\) such that \(\dot{\beta}\) is integrable on \([0,t]\) with \(\beta_t=\int_0^t\dot{\beta}_sds+\beta_0\) for any \(t>0\).

For any \((t,x,z)\in[0,\infty)\times\R\times[0,1]\), define the total reward function as

\begin{equation}
\label{eq:acc-util-inf}
\end{equation}

subject to  (<ref>) and (<ref>), under a given allocation process \(\balpha=\{\alpha_t\}_{t\geq0}\) and a given consumption process \(\bc=\{c_t\}_{t\geq0}\) with \(c_t\geq0\). For any \((t,x,z)\in[0,\infty)\times\R\times[0,1]\), define the value function as follows,

\begin{equation}
\label{eq:gen-value-inf}
V_\infty(t,x,z)=\sup_{(\balpha,\bc)\in\calA}J_\infty(t,x,z,\balpha,\bc),\,t\in[0,t);\ \ \lim_{t\to\infty}V_\infty(t,x,z)=0,
\end{equation}

subject to (<ref>) and (<ref>).

It is easy to show that the value function \(V_\infty\) in (<ref>) will have similar results as specified in Section <ref> therefore here we state these results without proofs.
First, we have the necessary condition for \(V_\infty\).

For any \((t,x,z)\in[0,\infty)\times\R\times[0,1]\) and \(\tau\in\bar{\mathbb T}_t\) where \(\bar{\mathbb T}_t\) denotes all \(\{\bar\calF_t\}_{t\geq0}\)-adapted stopping times \(\tau\) such that \(\tau\in[t,\infty)\) a.s.. Then the value function \(V_\infty\) defined in (<ref>) satisfies

\begin{equation}
\label{eq:gen-dpp-inf}
\tag{DPP'}
V_\infty(t,x,z)=\sup_{(\balpha,\bc)\in\calA}\E\left[\int_t^{\tau}\beta_sU_1(c_s)ds+V_\infty(\tau, \ctrldX_{\tau},\beta_{\tau})\,\biggl|\,\ctrldX_t=x,\beta_t=z\right].
\end{equation}

The corresponding HJB equation is given by

\begin{equation}
\label{eq:gen-hjb-inf}
\tag{HJB'}\left\{
\begin{aligned}
&\partial_tV_\infty(t,x,z)+\dot{\beta}_t\partial_zV_\infty(t,x,z)+rx\partial_xV_\infty(t,x,z)+\sup_{\alpha\in\R,c\geq0}\Big\{zU_1(c)+\infgen V_\infty(t,x,z)\Big\}=0,\ \ t\in[0,\infty);\\
\end{aligned}\right.
\end{equation}

Likewise, combining Proposition <ref> and Itô's formula, we have the following result.

If the value function \(V_\infty\) in (<ref>) is jointly continuous on \(\calD'=[0,\infty)\times(0,\infty)\times[0,1]\), then it is a viscosity solution to (<ref>) over the domain \(\calD'\).

With a classical solution to (<ref>), we have the following verification theorem serving as sufficient conditions for \(V_\infty\).

Suppose that \(U_1:[0,\infty)\to\R^+\in\mathcal U\) continuous at \(0\). Let \(w:[0,\infty)\times\R\times[0,1]\to\R\) be a function such that

\[w\in\calC^{1,2,1}\left([0,\infty)\times\R\times[0,1]\right),\]

and there exists a constant \(C>0\) with

\[w(t,x,z)\leq C(1+|x|^2),\ \ \forall (t,x,z)\in[0,\infty)\times\R\times[0,1].\]

* Assume that for any \((\alpha,c)\in\mathcal K\),

\[
\left\{
\begin{aligned}
&\partial_tw(t,x,z)+\dot{\beta}_t\partial_zw(t,x,z)+rx\partial_xw(t,x,z)+zU_1(c)+\infgen w(t,x,z)\leq0,\\
&\hspace{180pt}\forall (t,x,z)\in[0,\infty)\times\R\times[0,1];\\
&\lim_{t\to\infty}w(t,x,z)=\infty, \ \  \forall (x,z)\in\R\times[0,1].
\end{aligned}\right.
\]

Then \(w\geq V_\infty\) on \([0,\infty)\times\R\times[0,1]\).

* Assume further that there exists \(\hat{\alpha}:[0,\infty)\times\R\times[0,1]\to[-M,M]\) and \(\hat{c}:[0,\infty)\times\R\times[0,1]\to[0,M]\) such that

\[
\left\{
\begin{aligned}
&\partial_tw(t,x,z)+\dot{\beta}_t\partial_zw(t,x,z)+rx\partial_xw(t,x,z)+zU_1(\hat c(t,x,z))+\infgen[\hat{\alpha}(t,x,z),\hat{c}(t,x,z)]w(t,x,z)=0,\\
&\hspace{180pt}\forall (t,x,z)\in[0,T)\times\R\times[0,1];\\
&\lim_{t\to\infty}w(t,x,z)=\infty, \ \  \forall (x,z)\in\R\times[0,1],
\end{aligned}\right.
\]

also, with \(\beta_t=\beta_0+\int_0^t\dot{\beta}_sds\in[0,1]\) for all \(t\geq0\), the following SDE,

\[dX_t=
\left\{X_t\left[\hat{\alpha}(t,X_t,\beta_t)(\mu-r)+r\right]-\hat{c}(t,X_t,\beta_t)\right\}dt
\]

admits a unique solution \(\ctrldX[\hat{\balpha},\hat{\bc}]\) given \(X_0=x\) for any \(x\in\R\), and

\[\left(\hat{\balpha}=\left\{\hat{\alpha}_t\right\}_{t\geq0}=\left\{\hat{\alpha}(t,\ctrldX[\hat{\balpha},\hat{\bc}]_t,\beta_t)\right\}_{t\geq0}, \hat{\bc}=\left\{\hat{c}_t\right\}_{t\geq0}=\left\{\hat{c}(t,\ctrldX[\hat{\balpha},\hat{\bc}]_t,\beta_t)\right\}_{t\geq0}\right)\in\calA.\]

Then \(w= V_\infty\) on \([0,\infty)\times\R\times[0,1]\), with \((\hat{\balpha},\hat{\bc})\) being an optimal joint allocation-consumption process.

Given that \(V_\infty\in\calC^{1,2,1}(\calD')\) being a classical solution to (<ref>), then the optimal policy is given by

\begin{align}
\bar\alpha^*(t,x,z)&=-M\vee-\frac{\partial_xV_\infty(t,x,z)}{\partial_x^2V_\infty(t,x,z)}\cdot\frac{(\mu-r)}{\sigma^2x}\wedge M,\label{eq:gen-opt-alpha-inf}\\
\bar c^*(t,x,z)&=\argmax_{c\in[0,M]}\Big\{-\partial_xV_\infty(t,x,z)c+zU_1(c)\Big\}.\label{eq:gen-opt-c-inf}
\end{align}

Accordingly, for the inverse problem, we also have the following identifiability result.

Assume that

* \(\bar\alpha,\bar c\in\calC^{1,1,1}([0,\infty)\times(0,\infty)\times[0,1])\bigcap\calC^{0}([0,\infty)\times(0,\infty)\times[0,1])\);

* \(\bar\alpha(t,x,z)\in(-M,M)\) for all \((t,x,z)\in[0,T]\times(0,\infty)\times[0,1]\);

* for any \((t,z)\in[0,\infty)\times[0,1]\),

\[\bar c(t,x,z)<x,\ \  \forall x>0;\]

* \(\exists(t_0,z_0)\in[0,\infty)\times(0,1]\) such that  \(\bar c_0(\cdot)\coloneqq \bar c(t_0,\cdot,z_0)\) is invertible, and the following difference for any \((t,x,z)\in[0,T)\times(0,\infty)\times(0,1]\),

\[\Delta(t,x,z)\coloneqq\int_x^1\frac{dy}{y\bar\alpha(t,y,z)}-\int_{\bar c_0^{-1}(\bar c(t,x,z))}^1\frac{dy}{y\bar\alpha_T(y)}\]

depends only on \((t,z)\), namely, \(\Delta(t,x,z)\equiv\Delta(t,z)\).

Then both the discounting scheme characterized by \(\dot{\beta}\) and the utility function \(U_1\in\mathcal U\) with \(U_1(0)=0\) are identifiable up to an affine transform.

§.§.§ Special case: exponential discount rate

When $\beta(t) = \exp(-\rho t)$ for some $\rho>0$. The problem becomes time-homogeneous and time-consistent, hence the value function can be defined as

\begin{eqnarray}
V(x)= \sup_{(\alpha,c)\in \mathcal{A}}   \mathbb{E}\Bigg[\int_0^\infty e^{-\rho t} U_1(c_t)\dd t\left|X_0=x\Bigg].\right.
\end{eqnarray}

The stationary HJB Equation follows:

\begin{eqnarray}\label{eq:HJB3}
\rho V = \max_{c,\alpha}\Bigg\{U_1(c) + \Big((\alpha(\mu-r) + r )x -c\Big)  V_x + \frac{1}{2} x^2\sigma^2 \alpha^2 V_{xx} \Bigg\},
\end{eqnarray}

with $V(0) =0$.
First order condition leads to:

\begin{eqnarray}
\bar \alpha(x) = - \frac{\mu - r}{x\sigma^2} \frac{ V_x}{V_{xx}},\quad \bar c (x) = \left(U_1'\right)^{-1}(V_x),
\end{eqnarray}

or equivalently (for the second equation), $U_1'(\bar c) = V_x$.
Plug in  $\bar a(x)$ and $\bar c(x)$ back into the HJB  equation (<ref>):

\begin{eqnarray}\label{eq:HJB4}
{\rho} V &=& U_1(c^*) + \Big( r x -c^* + \frac{1}{2}(\mu-r)a^*(x)\Big)  V_x.
\end{eqnarray}

If we denote $W(x) := \int_0^xU'_1(c^*(y)){\rm d}y$ and $G(x) = U_1(c^*(x))$, then (<ref>) becomes the following first-order nonlinear equation:

\begin{eqnarray}
\label{eq:HJB6}
\rho W(x) &=& G(x) + \Big( r x -c^*(x) + \frac{1}{2}(\mu-r)a^*(x)\Big)W{'}(x),
\\
G'(x) &=& W'(x)\,(c^*)'(x)\nonumber
\end{eqnarray}

with initial conditions $W(0)=0$ and $G(0)=0$.
Assuming $\bar \alpha$ is differentiable and taking derivative on both side of the first equation in (<ref>):

\begin{eqnarray}\label{eq:HJB7}
0 &=& \Big( r \rho + \frac{1}{2}(\mu-r)\bar\alpha'(x)\Big)H(x)+ \Big( r x -\bar c(x) + \frac{1}{2}(\mu-r)\bar \alpha(x)\Big)H'(x),
\end{eqnarray}

where $H(x) = W{'}(x) = U_1'(\bar c(x))$.

Note that

\begin{eqnarray}
-   \frac{\Big( r \rho + \frac{1}{2}(\mu-r)(\bar a)'(x)\Big)}{\Big( r x -\bar c(x) + \frac{1}{2}(\mu-r)\bar a(x)\Big)} = \frac{H'(x)}{H(x)}
\end{eqnarray}

and $ -  \frac{x\sigma^2 a^*(x)} {\mu - r}=  \frac{ V_x}{V_{xx}}$. Therefore $\Bar{\rho}$ must satisfies

\begin{eqnarray}
\frac{\Big( r \rho + \frac{1}{2}(\mu-r)(a^*)'(x)\Big)}{\Big( r x -c^*(x) + \frac{1}{2}(\mu-r)a^*(x)\Big)} =\frac{x\sigma^2 a^*(x)} {\mu - r}.
\end{eqnarray}

Then we plug $\Bar{\rho}$ into (<ref>) to solve for $H(x)$. This recovers $U_1'$.

§ DISCRETE-TIME MDP WITH ENTROPY REGULARIZATION

The continuous-time framework in Section <ref> emphasizes the well-definedness of the mathematical framework when the client is subject to a generic discounting scheme with a time-varying rate, and it outlines conditions necessary for ensuring identifiability for both the utility functions and the discounting scheme. Building on these insights, this section explores a practical scenario focusing on the inference procedure. We adopt a parametric framework in which the client utilizes an exponential discounting scheme, parameterized by $\bar{\rho}$, alongside a utility function parameterized by $\bar{\theta} \in \mathbb{R}^d$. The client's preference parameter is summarized as $(\Bar{\rho},\Bar{\theta}
)$, which is unknown to the inference agent. The inference agent employs a maximum likelihood estimation method to infer the parameters $(\bar{\rho}, \bar{\theta})$. This analysis is conducted within a discrete-time MDP setting under the regularization of a Shannon entropy type. It encourages the client to fully explore the state-action space and introduces smoothness to the analysis at the same time; see [Haarnoja et al., 2017], for instance.

Mathematically, let us consider the entropy regularized MDP with state space $\calS$ and action space $\calA$, which can be finite or infinite. The state process follows $s_{t+1}\sim P(\cdot|s_t,a_t)$, with $P:\mathcal{S}\times \mathcal{A}\rightarrow\mathcal{P}(\mathcal{A})$ the transition kernel that maps from the joint state-action space to the distribution over the state space. After taking action $a$ at state $s$, we assume the client receives a deterministic reward $R(s,a)\in[0,1]$. Throughout the remainder of this paper, we will use the notation $\SumInt$ to denote the summation or integration over the action space, emphasizing that our framework accommodates both finite and infinite action spaces.

Under a generic preference parameter $(\theta,\rho)$, consider the entropy regularized objective:

\begin{align}
Q_{\rho,\theta}^*(s,a)=\max_\pi\mathbb{E}^\pi\left.\left[\sum_{t=1}^\infty \rho^ t \Big(U_\theta\big(R(s_t,a_t)\big)+\mathcal{H}\big(\pi(\cdot|s_t)\big)\Big)\right|s_0=s,a_0=a,a_t\sim \pi(s_t)\right],
\end{align}

where $\calH\big(\pi(\cdot|s)\big):= -\SumInt_{a\in\calA}\pi(a|s)\log\big(\pi(a|s)\big)$ is the Shannon's entropy. The optimal policy is then given by:

\begin{align}
\label{equ:optimal_policy}
\pi_{\rho,\theta}(a|s) =\frac{e^{Q_{\rho,\theta}^*(s,a)}}{\SumInt_{a\in\calA} e^{Q_{\rho,\theta}^*(s,a')}},
\end{align}

and the soft Bellman equation holds:

\begin{align}
\label{equ:soft_Bellman_equation}
V_{\rho,\theta}^*(s) = \log\Big(\SumInt_{a\in\calA} e^{Q_{\rho,\theta}^*(s,a)} \Big).
\end{align}

§.§ Maximum Likelihood Estimation

With a trajectory $\tau=\{(s_t,a_t)\}_{t=0}^\infty$ following the client's policy $\pi_{\Bar{\rho},\Bar{\theta}}$, we adopt a maximum likelihood estimation method to infer the client's preference parameter $(\Bar{\rho},\Bar{\theta})$, which is unknown to the inference agent. Specifically, the discounted likelihood of a trajectory $\tau=\{(s_t,a_t)\}_{t=0}^\infty$ following the client's policy $\pi_{\Bar{\rho},\Bar{\theta}}$ is defined as

\begin{eqnarray}
&&\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\left[\log \left(\prod_{t=0}^\infty (P(s_{t+1}|s_t,a_t) \pi_{\rho,\theta}(a_t|s_t))^{\gamma^t}\right)\right]\nonumber\\
&&\qquad=\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t\log \pi_{\rho,\theta}(a_t|s_t) \Big] + \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t\log  P(s_{t+1}|s_t,a_t)  \Big],\label{eq:likelihood_obj}
\end{eqnarray}

where the notation $\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}$ represents that the trajectory $\tau$ is sampled from applying policy $\pi_{\Bar{\rho},\Bar{\theta}}$ and $\gamma$ is a discount factor specified by the inference agent, which is potentially different from $\Bar\rho$.

Note that in our case $\gamma\neq \Bar{\rho}$ because $\Bar{\rho}$ is the client's discount factor and is unknown to the inference agent.
This distinguishes us from the usual IRL literature, where $\Bar{\rho} $ is always assumed to be known [Bloem and Bambos, 2014]. For example, Zeng et al., 2022 studied the IRL problem using a maximum likelihood estimator by setting  $\gamma = \rho$ in (<ref>) and showed that their algorithm
converges to a stationary point with a finite-time guarantee. Note that this stationary point may not be the ground-truth solution.

The maximum likelihood inference problem can be written as:

\begin{align}
\max_{(\rho,\theta)\in\Theta}\quad & \mathcal{L}(\rho,\theta):=
\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \log \pi_{\rho,\theta}(a_t|s_t)\Big],
\end{align}

where $\pi_{\rho,\theta}$ is the optimal policy under the preference parameter $(\rho,\theta)$ defined in (<ref>). Here for simplicity we set $\Theta := (0,1)\times \mathbb{R}^{d}$. The maximum likelihood problem is to find a preference parameter $(\Bar{\rho},\Bar{\theta})$ that generates the client's trajectory with the highest likelihood.

The goal is to investigate the landscape of the log-likelihood function $\mathcal{L}(\rho,\theta)$ with respect to $(\rho,\theta)$ and understand the possibility of recovering $(\Bar{\rho},\Bar{\theta})$, which is also referred to as the inverse problem. To proceed, we first show that $(\Bar{\rho},\Bar{\theta})$ is a stationary point of the likelihood function (see Proposition <ref>) and then show that likelihood function is concave near $(\Bar{\rho},\Bar{\theta})$ (see Theorem <ref>). Interestingly, the results of the landscape analysis is independent of the choice of $\gamma$, making our proposed method robust and practical.

It holds that

\begin{eqnarray}\label{eq:first-order}
\nabla_\theta \mathcal{L}(\rho,\theta)
\Big|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})} = 0,\quad \nabla_\rho \mathcal{L}(\rho,\theta)\Big|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})} = 0.
\end{eqnarray}

Proposition <ref> suggests that the gradient of the likelihood function equals zero at the client's preference parameter value $(\Bar{\rho},\Bar{\theta})$, and hence $(\Bar{\rho},\Bar{\theta})$ is a stationary point of the likelihood function $\mathcal{L}(\rho,\theta)$.

Our proof can be divided into three steps. We first provide some useful formulas regarding the first-order and second-order derivatives of $Q,V$ with respect to $(\rho,\theta)$. With such formulas, we next derive the derivatives of the log-likelihood function. Finally, we show that (<ref>) holds.

Step 1. To begin with, for any $(\rho,\theta)\in \Theta$, for any $(s_t,a_t)\in\calS\times\calA$,

\begin{align}
\nabla_\theta Q_{\rho,\theta}(s_t,a_t)
&= \nabla_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[\nabla_\theta V_{\rho,\theta}(s_{t+1}) \Big]\nonumber\\
&= \nabla_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[\nabla_\theta \log({\SumInt_a e^{Q_{\rho,\theta}(s_{t+1},a)}}  ) \Big]\label{equ:likelihood_gradient_0.1}\\
&= \nabla_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[
\SumInt_a \pi_{\rho,\theta}(a|s_{t+1})  \nabla_\theta  Q_{\rho,\theta}(s_{t+1},a) \Big]\label{equ:likelihood_gradient_0.2}\\
&= \nabla_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t),a_{t+1}\sim \pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[ \nabla_\theta  Q_{\rho,\theta}(s_{t+1},a_{t+1}) \Big],\label{equ:likelihood_gradient_1}
\end{align}

where (<ref>) holds by the soft Bellman equation and (<ref>) holds  because $\pi_{\rho,\theta}$ is the optimal policy.
Applying (<ref>) recursively yields:

\begin{eqnarray}\label{eq:recursion_result}
\nabla_\theta Q_{\rho,\theta}(s_t,a_t) = \E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta U_\theta(R(s_{k}',a_{k}'))\Big|s_t,a_t  \Big],
\end{eqnarray}

where $\tau'=\{s'_k,a'_k\}_{k=0}^{\infty}$ denotes a trajectory following $\pi_{\rho,\theta}$. Similarly,

\begin{align}
\nabla_\rho Q_{\rho,\theta}(s_t,a_t)
&= \nabla_\rho U_\theta(R(s_t,a_t))+\nabla_\rho (\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ V_{\rho,\theta}(s_{t+1}) \Big])\label{equ:likelihood_gradient_0.25}\\
&=  \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ V_{\rho,\theta}(s_{t+1}) \Big] + \rho\E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\rho \log({\SumInt_a e^{Q_{\rho,\theta}(s_{t+1},a)}}  )  \Big]\label{equ:likelihood_gradient_0.3}\\
&=  \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ V_{\rho,\theta}(s_{t+1}) \Big]+ \rho\E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[
\SumInt_a \pi_{\rho,\theta}(a|s_{t+1})  \nabla_\rho  Q_{\rho,\theta}(s_{t+1},a) \Big]\label{equ:likelihood_gradient_0.4}\\
&=  \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ V_{\rho,\theta}(s_{t+1}) \Big]+ \rho\E_{s_{t+1}\sim P(\cdot|s_t,a_t),\ a_{t+1}\sim \pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[ \nabla_\rho  Q_{\rho,\theta}(s_{t+1},a_{t+1}) \Big],\label{equ:likelihood_gradient_0}
\end{align}

where (<ref>) holds by the Bellman equation and (<ref>) holds because $\pi_{\rho,\theta}$ is the optimal policy.  Applying (<ref>) recursively yields:

\begin{eqnarray}\label{eq:recursion_result2}
\nabla_\rho Q_{\rho,\theta}(s_t,a_t) = \E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t+1}^\infty \rho^{k-t-1}V_{\rho,\theta}(s'_{k})\Big|s_t,a_t  \Big].
\end{eqnarray}

Furthermore, we have for any $s_t\in\calS$,

\begin{align}
&\nabla_\theta V_{\rho,\theta}(s_t) = \nabla_\theta (\log \SumInt_a e^{Q_{\rho,\theta}(s_t,a)})\nonumber\\
&= \E_{a_t\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta {Q_{\rho,\theta}(s_t,a_t)}\Big]= \E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta U_\theta(R(s_{k}',a_{k}'))\Big|s_t \Big],\label{equ:likelihood_gradient_3.5}
\end{align}

where the last equation holds by (<ref>). In addition,

\begin{align}
&\nabla_\rho V_{\rho,\theta}(s_t)
= \nabla_\rho (\log \SumInt_a e^{Q_{\rho,\theta}(s_t,a)})\nonumber\\
&= \E_{a_t\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\rho {Q_{\rho,\theta}(s_t,a_t)}\Big]= \E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t+1}^\infty \rho^{k-t-1}V_{\rho,\theta}(s_{k}')\Big|s_t \Big],\label{equ:likelihood_gradient_3.6}
\end{align}

where the last equation holds by (<ref>).

In summary, it holds that for any $(\rho,\theta)\in\Theta,s_t\in\calS,a_t\in\calA$,

\begin{align}
\nabla_\theta Q_{\rho,\theta}(s_t,a_t) &= E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta U_\theta(R(s_{k}',a_{k}'))\Big|s_t,a_t  \Big],\label{equ:gradient_Q_theta}\\
\nabla_\rho Q_{\rho,\theta}(s_t,a_t) &= E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t+1}^\infty \rho^{k-t-1}V_{\rho,\theta}(s_{k}')\Big|s_t,a_t  \Big],\label{equ:gradient_Q_eho}\\
\nabla_\theta V_{\rho,\theta}(s_t) &=\E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta U_\theta(R(s_{k}',a_{k}'))\Big|s_t \Big],\label{equ:gradient_V_theta}\\
\nabla_\rho V_{\rho,\theta}(s_t) &= \E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t+1}^\infty \rho^{k-t-1}V_{\rho,\theta}(s_{k}')\Big|s_t \Big].\label{equ:gradient_V_rho}
\end{align}

Step 2. Next, we derive the gradients of the log-likelihood function. For any $(\rho,\theta)\in\Theta$,

\begin{align}
\mathcal{L}(\rho,\theta) &= \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \log \pi_{\rho,\theta}(a_t|s_t)\Big]
=\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \log \frac{e^{Q_{\rho,\theta}(s_t,a_t)}}{\SumInt_a e^{Q_{\rho,\theta}(s_t,a)}}\Big]\label{equ:likelihood_gradient_4.5}\\
&=\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \Big(Q_{\rho,\theta}(s_t,a_t) -V_{\rho,\theta}(s_t)\Big)\Big]\nonumber\\
&= \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t U_\theta(R(s_t,a_t))\Big]-\E_{s_0\sim\mu(\cdot)}\Big[V_{\rho,\theta}(s_0)\Big]+(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} V_{\rho,\theta}(s_{t})\Big],\label{equ:likelihood_gradient_1.5}
\end{align}

where $\mu$ is the distribution of the initial state $s_0$. (<ref>) holds by the optimality of the policy, and (<ref>) holds by the soft Bellman equation. Taking the gradient of (<ref>) with respect to $\theta$ gives

\begin{align}
\nabla_\theta \mathcal{L}(\rho,\theta)& =    \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \nabla_\theta U_\theta(R(s_t,a_t))\Big] - \E_{s_0\sim\mu(\cdot)}\Big[\nabla_\theta V_{\rho,\theta}(s_0)\Big]\nonumber\\
&\qquad+(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\theta V_{\rho,\theta}(s_{t})\Big].\label{equ:likelihood_gradient_3}
\end{align}

Combining (<ref>) with (<ref>) gives:

\begin{align}
\nabla_\theta \mathcal{L}(\rho,\theta) &= \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \nabla_\theta U_\theta(R(s_t,a_t)) \Big]-\E_{\tau\sim\pi_{\rho,\theta}}\Big[\sum_{t=0}^\infty {\rho}^t \nabla_\theta U_\theta(R(s_t,a_t)) \Big]\nonumber\\
&\qquad+ (\rho-\gamma) \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}
\E_{\tau'\sim\pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t} \nabla_\theta U_\theta(R(s'_{k},a'_{k}))\,\Big|\,s_t\Big]\Big].\label{equ:likelihood_gradient_6}
\end{align}

Similarly, taking the gradient of (<ref>) with respect to $\rho$ gives

\begin{align}
\nabla_\rho \mathcal{L}(\rho,\theta)
&= - \E_{s_0\sim\mu(\cdot)}\Big[\nabla_\rho V_{\rho,\theta}(s_0)\Big]
+ \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} V_{\rho,\theta}(s_t)\Big]\nonumber\\
&\qquad +(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\rho V_{\rho,\theta}(s_t)\Big].\label{equ:likelihood_gradient_5}
\end{align}

Combining (<ref>) with (<ref>) yields:

\begin{align}
\nabla_\rho \mathcal{L}(\rho,\theta)  &={\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} V_{\rho,\theta}(s_t)\Big]} -\mathbb{E}_{\tau\sim \pi_{\rho,\theta}}\Big[\sum_{t=1}^\infty \rho^{t-1}V_{\rho,\theta}(s_{t}) \Big]  \nonumber\\
&\qquad+(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \E_{\tau'\sim \pi_{\rho,\theta}}\Big[\sum_{k=t+1}^\infty \rho^{k-t-1}V_{\rho,\theta}(s'_{k})\,\Big|\, s_t  \Big]\Big].\label{equ:likelihood_gradient_7}
\end{align}

Step 3. Finally, when $(\rho,\theta) =(\Bar{\rho},\Bar{\theta})$, by (<ref>) we have

\begin{eqnarray}
\nabla_\theta \mathcal{L}(\rho,\theta)|_{(\rho,\theta) =(\Bar{\rho},\Bar{\theta})} &= \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \nabla_\theta U_\theta(R(s_t,a_t)) \Big]-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty {{\rho}^t} \nabla_\theta U_\theta(R(s_t,a_t)) \Big]\nonumber\\
&\qquad+ (\rho-\gamma) \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}
\sum_{k=t}^\infty \rho^{k-t} \nabla_\theta U_\theta(R(s_{k},a_{k}))\Big].\label{eq:grad0}
\end{eqnarray}

Note that for the last line of the above equation,

\begin{eqnarray}
&&\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}
\sum_{k=t}^\infty \rho^{k-t} \nabla_\theta U_\theta(R(s_{k},a_{k}))\Big]= \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \left(\frac{\gamma}{\rho}\right)^{t-1}
\sum_{k=t}^\infty \rho^{k-1} \nabla_\theta U_\theta(R(s_{k},a_{k}))\Big]\nonumber\\
&=&   \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\sum_{t=1}^{k} \left(\frac{\gamma}{\rho}\right)^{t-1}
\rho^{k-1} \nabla_\theta U_\theta(R(s_{k},a_{k}))\Big]\label{eq:inter1}\\
\rho^{k-1} \nabla_\theta U_\theta(R(s_{k},a_{k}))\Big] = \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho}
\nabla_\theta U_\theta(R(s_{k},a_{k}))\Big],\label{eq:inter2}
\end{eqnarray}

where (<ref>) holds by changing the order of summations. Plugging (<ref>) back into (<ref>), we have the desired result that $ \nabla_\theta \mathcal{L}(\rho,\theta)|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})} = 0$.

Similarly we have $ \nabla_\rho \mathcal{L}(\rho,\theta)|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})} = 0$.

We next show results on the Hessian matrix.

It holds that

\begin{eqnarray*}
\nabla^2_\theta \mathcal{L}(\rho,\theta){|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})}}
&=& - \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty\gamma^{t}
\mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a) \Big] \Big],\\
\nabla^2_\rho \mathcal{L}(\rho,\theta)  {|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})}}
&=& -  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big],\\
\nabla_\theta\nabla_\rho \mathcal{L}(\rho,\theta){|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})}}
&=& -\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty\gamma^t{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a),\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big],
\end{eqnarray*}

in which we define

\begin{align}
\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]
& :=\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\nabla_\theta Q_{\rho,\theta}(s_t,a)^\top\Big]\nonumber\\
&\qquad- \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]^\top,\label{equ:definition_of_variance}
\end{align}

\begin{align}
&{\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t+1},a),\nabla_\rho Q_{\rho,\theta}(s_{t+1},a)\Big]\nonumber\\
&:=\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]\nonumber
\\
&\qquad \qquad - \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]\in\mathbb{R}^d.\label{equ:definition_of_covariance}
\end{align}

In addition,

\begin{eqnarray*}
\mathcal H(\Bar{\rho},\Bar{\theta})&&:=
\begin{pmatrix}
\nabla^2_\theta \mathcal{L}(\rho,\theta) &  \nabla_\theta \nabla_\rho \mathcal{L}(\rho,\theta)\\
\nabla_\theta \nabla_\rho \mathcal{L}(\rho,\theta)^\top &  \nabla^2_\rho \mathcal{L}(\rho,\theta)
\end{pmatrix}{\Bigg|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})} }
\end{eqnarray*}

is negative semi-definite.

Theorem <ref> suggests that the log-liklihood function $\mathcal{L}(\rho,\theta)$ is concave near the client's preference parameter $(\Bar{\rho},\Bar{\theta})$. As mentioned earlier, an interesting finding is that the negative semi-definite property of the Hessian  does not rely on  the choice $\gamma$, making the likelihood estimation method particularly suitable for inference problems.

Our proof consists of two parts. We first derive formulas for the second-order derivatives of $Q,V$ with respect to $\theta$ and $\rho$. Then we calculate the second-order derivatives of the log-likelihood function and study its Hessian matrix when $(\rho,\theta) = (\Bar{\rho},\Bar{\theta})$.

Step 1. To begin with, for any $(\rho,\theta)\in\Theta$ by taking the derivative of (<ref>), we have:

\begin{align}
\nabla_\theta^2V_{\rho,\theta}(s_t)
&= \nabla_\theta \Big(\SumInt_a \pi_{\rho,\theta}(a|s_t) \nabla_\theta Q_{\rho,\theta}(s_t,a)^\top\Big)\nonumber\\
&=\SumInt_a \nabla_\theta \Big( e^{Q_{\rho,\theta}(s_t,a)-V_{\rho,\theta}(s_t) } \Big)\nabla_\theta Q_{\rho,\theta}(s_t,a)^\top
+\SumInt_a \pi_{\rho,\theta}(a|s_t) \nabla_\theta^2 Q_{\rho,\theta}(s_t,a) \label{equ:second_order_derivative_1}\\
&= \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\nabla_\theta Q_{\rho,\theta}(s_t,a)^\top\Big]\nonumber\\
&\qquad- \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]^\top
+\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta^2 Q_{\rho,\theta}(s_t,a)\Big]\label{equ:second_order_derivative_2}\\
&= \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]
+\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta^2 Q_{\rho,\theta}(s_t,a)\Big],\label{equ:second_order_derivative_2.5}
\end{align}

where the covariance matrix $\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]$ is defined in (<ref>).
In particular, (<ref>) holds because:

\begin{align*}
\pi_{\rho,\theta}(a|s_t) = \frac{e^{Q_{\rho,\theta}(s_t,a)}}{\SumInt_{a'}e^{Q_{\rho,\theta}(s,a')}} = e^{Q_{\rho,\theta}(s_t,a)-V_{\rho,\theta}(s_t) },
\end{align*}

(<ref>) holds by (<ref>).
In addition, we have:

\begin{align}
\nabla^2_\theta Q_{\rho,\theta}(s_t,a_t)
&= \nabla_\theta^2 U_\theta (R(s_t,a_t)) + \rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[\nabla_\theta^2 V_{\rho,\theta}(s_{t+1})\Big]\nonumber\\
&= \nabla_\theta^2 U_\theta (R(s_t,a_t))
+\rho\E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[\mathbb{V}_{a_{t+1}\sim\pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big]\Big]\nonumber\\
+\rho \E_{\tau\sim\pi_{\rho,\theta}}\Big[\nabla_\theta^2 Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big|s_t,a_t\Big],\label{equ:second_order_derivative_3}
\end{align}

where (<ref>) holds by (<ref>). Applying (<ref>) recursively yields:

\begin{align}
\nabla^2_\theta Q_{\rho,\theta}(s_t,a_t)
& = \E_{\tau'\sim\pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta^2 U_\theta (R(s_{k}',a_{k}'))\Big|s_t,a_t \Big]\nonumber\\
&\qquad+ \E_{\tau'\sim\pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t+1} \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k+1}')}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k+1}',a)\Big]\Big|s_t,a_t\Big].\label{equ:second_order_derivative_2.6}
\end{align}

Similarly, for any $\rho\in(0,1)$, $\theta\in\Theta$, by taking the gradient  of (<ref>) with respect to $\rho$, we have

\begin{align}
\nabla_\rho^2V_{\rho,\theta}(s_t)
&= \nabla_\rho \Big(\SumInt_a \pi_{\rho,\theta}(a|s_t) \nabla_\rho Q_{\rho,\theta}(s_t,a)\Big)\nonumber\\
&=\SumInt_a \nabla_\rho \Big( e^{Q_{\rho,\theta}(s_t,a)-V_{\rho,\theta}(s_t) } \Big)\nabla_\rho Q_{\rho,\theta}(s_t,a)
+\SumInt_a \pi_{\rho,\theta}(a|s_t) \nabla_\rho^2 Q_{\rho,\theta}(s_t,a) \nonumber\\
&= \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[(\nabla_\rho Q_{\rho,\theta}(s_t,a))^2\Big]
- \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]^2
+\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\rho^2 Q_{\rho,\theta}(s_t,a)\Big]\nonumber\\
&= \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]
+\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\rho^2 Q_{\rho,\theta}(s_t,a)\Big].
\label{equ:second_order_derivative_2.7}
\end{align}

Similarly, by taking the gradient of (<ref>),

\begin{align}
\nabla_\rho^2 Q_{\rho,\theta}(s_t,a_t)
&= \nabla_\rho\Big(\E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ V_{\rho,\theta}(s_{t+1})\Big]
+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\rho V_{\rho,\theta}(s_{t+1})\Big]\Big)\nonumber\\
&= 2\E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\rho V_{\rho,\theta}(s_{t+1})\Big]
+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\rho^2 V_{\rho,\theta}(s_{t+1})\Big]\nonumber\\
& = 2\E_{\tau\sim\pi_{\rho,\theta}}\Big[ \nabla_\rho Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big|s_t,a_t\Big]
+\rho \E_{s_{t+1}}\Big[\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t+1},a)\Big]\Big] \nonumber\\
&\qquad+\rho \E_{a_{t+1}\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\rho^2 Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big|s_t,a_t\Big].
\label{equ:second_order_derivative_6}
\end{align}

Applying (<ref>) recursively yields:

\begin{align}
\nabla_\rho^2 Q_{\rho,\theta}(s_t,a_t)
&= 2\E_{\tau'\sim\pi_{\rho,\theta}}\Big[ \sum_{k=t}^\infty \rho^{k-t} \nabla_\rho Q_{\rho,\theta}(s_{k+1}',a_{k+1}')\Big|s_t,a_t\Big] \nonumber\\
\Big[\sum_{k=t}^\infty \rho^{k-t+1}\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k+1}')}\Big[\nabla_\rho Q_{\rho,\theta}(s_{k+1}',a)\Big]\Big|s_t,a_t\Big].\label{equ:second_order_derivative_2.8}
\end{align}

Furthermore, for any $(\rho,\theta)\in\Theta$, by taking the gradient of (<ref>) with respect to $\theta$ , we have

\begin{align}
&\nabla_\theta\nabla_\rho V_{\rho,\theta}(s_t)
= \nabla_\theta \Big(\SumInt_a \pi_{\rho,\theta}(a|s_t) \nabla_\rho Q_{\rho,\theta}(s_t,a)\Big)\nonumber\\
&=\SumInt_a \nabla_\theta \Big( e^{Q_{\rho,\theta}(s_t,a)-V_{\rho,\theta}(s_t) } \Big)\nabla_\rho Q_{\rho,\theta}(s_t,a)
+\SumInt_a \pi_{\rho,\theta}(a|s_t) \nabla_\theta\nabla_\rho Q_{\rho,\theta}(s_t,a) \nonumber\\
&= \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta Q_{\rho,\theta}(s_t,a)\nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]
- \E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\theta Q_{\rho,\theta}(s_t,a)\Big]\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[ \nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]\nonumber\\
&\qquad+\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta\nabla_\rho Q_{\rho,\theta}(s_t,a)\Big]\nonumber\\
& = {\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t+1},a),\nabla_\rho Q_{\rho,\theta}(s_{t+1},a)\Big]\nonumber\\
&\qquad+\E_{a\sim\pi_{\rho,\theta}(\cdot|s_t)}\Big[\nabla_\theta\nabla_\rho Q_{\rho,\theta}(s_t,a)\Big],\label{equ:second_order_derivative_2.9}
\end{align}

where the “covariance” between $\nabla_\theta Q$ and $\nabla_\rho Q$ is defined in (<ref>).
Note that ${\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t+1},a),\nabla_\rho Q_{\rho,\theta}(s_{t+1},a)\Big] \in \mathbb{R}^d$, as we have $\theta\in \mathbb{R}^d$ and $\rho \in \mathbb{R}$.

Lastly, by taking the gradient of (<ref>), we obtain

\begin{align*}
&\nabla_\theta\nabla_\rho Q_{\rho,\theta}(s_t,a_t)\\
&= \nabla_\theta\Big(\E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ V_{\rho,\theta}(s_{t+1})\Big]
+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\rho V_{\rho,\theta}(s_{t+1})\Big]\Big)\\
&= \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\theta V_{\rho,\theta}(s_{t+1})\Big]
+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\theta\nabla_\rho V_{\rho,\theta}(s_{t+1})\Big]\\
& = \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ \nabla_\theta V_{\rho,\theta}(s_{t+1})\Big]
+ \rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[ {\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t+1})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t+1},a),\nabla_\rho Q_{\rho,\theta}(s_{t+1},a)\Big]\Big]\\
&\qquad+ \rho \E_{\tau\sim\pi_{\rho,\theta}}\Big[\nabla_\theta\nabla_\rho Q_{\rho,\theta}(s_{t+1},a)\Big|s_t,a_t\Big].
\end{align*}

Applying the last equation recursively yields:

\begin{align}
\nabla_\theta\nabla_\rho Q_{\rho,\theta}(s_t,a_t)
&=\E_{\tau'\sim\pi_{\rho,\theta}}\Big[ \sum_{k=t}^\infty \rho^{k-t}\nabla_\theta V_{\rho,\theta}(s_{k+1}')\Big|s_t,a_t\Big]\nonumber\\
&\qquad+ \E_{\tau'\sim\pi_{\rho,\theta}}\Big[ \sum_{k=t}^\infty\rho^{k-t+1}{\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k+1})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k+1}',a),\nabla_\rho Q_{\rho,\theta}(s_{k+1}',a)\Big]\Big|s_t,a_t\Big].\label{equ:second_order_derivative_2.95}
\end{align}

In summary, by combining (<ref>), (<ref>), (<ref>), (<ref>),(<ref>),
and (<ref>), we have the following formulas of the second-order gradients of the value function: for any $(\rho,\theta)\in\Theta$, for any $(s_t,a_t)\in\calS\times\calA$,

\begin{align}
\nabla^2_\theta V_{\rho,\theta}(s_t)
&= \E_{\tau'\sim\pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta^2 U_\theta (R(s_{k}',a_{k}')) \Big|s_t\Big]\nonumber\\
&\qquad+ \E_{\tau'\sim\pi_{\rho,\theta}}\Big[\sum_{k=t}^\infty \rho^{k-t} \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k}')}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k}',a)\Big]\Big|s_t\Big],\label{eq:v_secondorder_1}\\
\nabla^2_\rho V_{\rho,\theta}(s_t)
&= 2\E_{\tau'\sim\pi_{\rho,\theta}}\Big[ \sum_{k=t}^\infty \rho^{k-t} \nabla_\rho Q_{\rho,\theta}(s_{k+1}',a_{k+1})\Big|s_t\Big] \nonumber\\
\Big[\sum_{k=t}^\infty \rho^{k-t}\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k}')}\Big[\nabla_\rho Q_{\rho,\theta}(s_{k}',a)\Big]\Big|s_t\Big],\label{eq:v_secondorder_2}\\
\nabla_\theta\nabla_\rho V_{\rho,\theta}(s_t)
&= \E_{\tau'\sim\pi_{\rho,\theta}}\Big[ \sum_{k=t}^\infty \rho^{k-t}\nabla_\theta V_{\rho,\theta}(s_{k+1}')\Big|s_t\Big]\nonumber\\
&\qquad+ \E_{\tau'\sim\pi_{\rho,\theta}}\Big[ \sum_{k=t}^\infty\rho^{k-t}{\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k}')}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k}',a),\nabla_\rho Q_{\rho,\theta}(s_{k}',a)\Big]\Big|s_t\Big].\label{eq:v_secondorder_3}
\end{align}

Step 2. Next, we calculate the derivatives of the log-likelihood function. By straight-forward calculations using (<ref>) and combining with (<ref>), the second-order derivative of the log-likelihood function to $\theta$ satisfies:

\begin{align*}
\nabla^2_\theta \mathcal{L}(\rho,\theta) & =    \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \nabla^2_\theta U_\theta(R(s_t,a_t))\Big] - \E_{s_0\sim\mu(\cdot)}\Big[\nabla^2_\theta V_{\rho,\theta}(s_0)\Big]\nonumber\\
&\qquad+(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla^2_\theta V_{\rho,\theta}(s_{t})\Big]\\
& =  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \nabla^2_\theta U_\theta(R(s_t,a_t))\Big]
-\E_{\tau\sim\pi_{\rho,\theta}}\Big[\sum_{t=0}^\infty \rho^{t}\nabla_\theta^2 U_\theta (R(s_{t},a_{t})) \Big]\nonumber\\
&\qquad- \E_{\tau\sim\pi_{\rho,\theta}}\Big[\sum_{t=0}^\infty \rho^{t} \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a)\Big]\Big]\\
&\qquad+(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla^2_\theta V_{\rho,\theta}(s_{t})\Big].
\end{align*}

Note that when $(\rho,\theta)=(\Bar{\rho},\Bar{\theta})$,

\begin{eqnarray*}
\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla^2_\theta V_{\rho,\theta}(s_{t})\Big]
&=& \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \sum_{k=t}^\infty \rho^{k-t}\nabla_\theta^2 U_\theta (R(s_{k},a_{k})) \Big]\\
&&+  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}\sum_{k=t}^\infty \rho^{k-t} \mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{k})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k},a) \Big]\Big]\\
&=&  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho}
\nabla_\theta^2 U_\theta (R(s_{k},a_{k})) \Big]\\
&&+  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho}
\mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{k})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k},a) \Big]\Big],
\end{eqnarray*}

where the first equality holds by (<ref>) and the last equality holds by changing the order of summations. Plugging the above result to $\nabla^2_\theta \mathcal{L}(\rho,\theta)$, we obtain that,

\begin{eqnarray*}
\nabla^2_\theta \mathcal{L}(\rho,\theta) |_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})}
& =&  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \nabla^2_\theta U_\theta(R(s_t,a_t))\Big]
-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \rho^{t}\nabla_\theta^2 U_\theta (R(s_{t},a_{t})) \Big]\nonumber\\
&\qquad& - \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \rho^{t} \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a)\Big]\Big]\\
\nabla_\theta^2 U_\theta (R(s_{t},a_{t})) \Big] - \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\gamma^{t}
\nabla_\theta^2 U_\theta (R(s_{t},a_{t})) \Big]\\
\mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a) \Big]\Big]\\
&\qquad&- \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\gamma^{t}
\mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a) \Big] \Big]\\
&=& - \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{{t=0}}^\infty\gamma^{t}
\mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a) \Big] \Big].
\end{eqnarray*}

Similarly, using the result in (<ref>) and (<ref>),

\begin{eqnarray}
\nabla^2_\rho \mathcal{L}(\rho,\theta)
&=& -2\E_{\tau\sim\pi_{\rho,\theta}}\Big[ \sum_{t=0}^\infty \rho^{t} \nabla_\rho Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big] -\E_{\tau\sim\pi_{\rho,\theta}}
\Big[\sum_{t=0}^\infty \rho^{t}\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big]\nonumber\\
&&+2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\rho V_{\rho,\theta}(s_t)\Big]\nonumber\\
&& +(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla^2_\rho V_{\rho,\theta}(s_t)\Big].\label{equ:likelihood_hessian3}
\end{eqnarray}

Note that when $(\rho,\theta) = (\Bar{\rho},\Bar{\theta})$,

\begin{eqnarray*}
\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla^2_\rho V_{\rho,\theta}(s_t)\Big]
&=&2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}\sum_{k=t}^\infty \rho^{k-t} \nabla_\rho Q_{\rho,\theta}(s_{k+1},a_{k+1})\Big]\\
&& + \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \sum_{k=t}^\infty \rho^{k-t} \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{k},a)\Big]\Big]\\
&=&2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho} \nabla_\rho Q_{\rho,\theta}(s_{k+1},a_{k+1})\Big]\\
&& + \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho} \mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{k})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{k},a)\Big]\Big],
\end{eqnarray*}

where the last equality holds by changing the order of summations.
Plugging the above result to $\nabla^2_\rho \mathcal{L}(\rho,\theta)$, we have

\begin{eqnarray*}
&\nabla^2_\rho& \mathcal{L}(\rho,\theta)  |_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})}\\
&=& -2\E_{\tau\sim\pi_{\rho,\theta}}\Big[ \sum_{t=0}^\infty \rho^{t} \nabla_\rho Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big] -\E_{\tau\sim\pi_{\rho,\theta}}
\Big[\sum_{t=0}^\infty \rho^{t}\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big]\nonumber\\
&&+2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\rho V_{\rho,\theta}(s_t)\Big]\nonumber\\
&& + 2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \rho^{t} \nabla_\rho Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big]-2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t} \nabla_\rho Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big] \nonumber\\
&& + \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \rho^t \mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big] -  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^t \mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big] \\
&=& {2\, \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\rho V_{\rho,\theta}(s_t)\Big] - 2\, \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t} \nabla_\rho Q_{\rho,\theta}(s_{t+1},a_{t+1})\Big]}
\\
&&-  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^t \mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big] \\
&&{-2 \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\nabla_\rho Q_{\rho,\theta}(s_1,a_1)\Big]
-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\mathbb{V}_{a\sim\pi_{\rho,\theta}(\cdot|s_{0})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{0},a)\Big]\Big]}\\
&=& {-  \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \gamma^t \mathbb{V}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big],}
\end{eqnarray*}

where the last equality holds by the expressions of $\nabla_\rho V_{\rho,\theta}$ and $\nabla_\rho Q_{\rho,\theta}$ in (<ref>) and (<ref>).

\begin{eqnarray}
\nabla_\theta\nabla_\rho \mathcal{L}(\rho,\theta)& =&    - \E_{s_0\sim\mu(\cdot)}\Big[\nabla_\theta \nabla_\rho \, V_{\rho,\theta}(s_0)\Big]\nonumber+\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\theta V_{\rho,\theta}(s_{t})\Big]\\
&&\qquad +(\rho-\gamma)\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\theta\nabla_\rho V_{\rho,\theta}(s_{t})\Big].\label{eq:second_cross}
\end{eqnarray}

Note that when $(\rho,\theta)=(\Bar{\rho},\Bar{\theta})$,

\begin{eqnarray*} &&\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\theta\nabla_\rho V_{\rho,\theta}(s_{t})\Big]\\
&=&\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}\sum_{k=t}^\infty \rho^{k-t}\nabla_\theta V_{\rho,\theta}(s_{k+1})\Big] \\
&&+ \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1}\sum_{k=t}^\infty \rho^{k-t}{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{k})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k},a),\nabla_\rho Q_{\rho,\theta}(s_{k},a)\Big]\Big]\\
&=&\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho}\nabla_\theta V_{\rho,\theta}(s_{k+1})\Big] \\
&&+ \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{k=1}^\infty\frac{\gamma^{k}-\rho^{k}}{\gamma-\rho}{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{k})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{k},a),\nabla_\rho Q_{\rho,\theta}(s_{k},a)\Big]\Big].
\end{eqnarray*}

The last equality holds by changing the order of summations. Plugging the last equality back to (<ref>) and applying (<ref>) for $t=0$, we have

\begin{eqnarray}
&\nabla_\theta\nabla_\rho& \mathcal{L}(\rho,\theta)|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})}\nonumber\\
&=&   - \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty \rho^{t}\nabla_\theta V_{\rho,\theta}(s_{t+1})\Big]
-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty\rho^{t}{\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a),\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big]\nonumber\\
&&+\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\theta V_{\rho,\theta}(s_{t})\Big]\nonumber\\
&&+\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\rho^t\nabla_\theta V_{\rho,\theta}(s_{t+1})\Big]-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\gamma^t\nabla_\theta V_{\rho,\theta}(s_{t+1})\Big]\nonumber \\
&&+ \E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\rho^t{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a),\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big]\nonumber\\
&& -\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\gamma^t{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a),\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big]\nonumber\\
&=& -\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\gamma^t{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a),\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big] \nonumber\\
&&{-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[{\rm \bf Cov}_{a\sim\pi_{\rho,\theta}(\cdot|s_{0})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{0},a),\nabla_\rho Q_{\rho,\theta}(s_{0},a)\Big]\Big]}\label{equ:landscape_1}\\
&=& {-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=0}^\infty\gamma^t{\rm \bf Cov}_{a\sim\pi_{\Bar{\rho},\Bar{\theta}}(\cdot|s_{t})}\Big[\nabla_\theta Q_{\rho,\theta}(s_{t},a),\nabla_\rho Q_{\rho,\theta}(s_{t},a)\Big]\Big]},\nonumber
\end{eqnarray}

where (<ref>) holds by the fact that

\begin{align*}
&\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty \gamma^{t-1} \nabla_\theta V_{\rho,\theta}(s_{t})\Big]
-\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[\sum_{t=1}^\infty\gamma^t\nabla_\theta V_{\rho,\theta}(s_{t+1})\Big]
=\E_{\tau\sim\pi_{\Bar{\rho},\Bar{\theta}}}\Big[  \nabla_\theta V_{\rho,\theta}(s_{1})\Big].
\end{align*}

To summarize,

\begin{eqnarray*}
\mathcal H&&:=
\begin{pmatrix}
\nabla^2_\theta \mathcal{L}(\rho,\theta) &  \nabla_\theta \nabla_\rho \mathcal{L}(\rho,\theta)\\
\nabla_\theta \nabla_\rho \mathcal{L}(\rho,\theta)^\top &  \nabla^2_\rho \mathcal{L}(\rho,\theta)
\end{pmatrix}{\Bigg|_{(\rho,\theta) = (\Bar{\rho},\Bar{\theta})} }\\
\begin{pmatrix}
\nabla_\theta Q_{\rho,\theta}(s_{t},a)\\
\nabla_\rho Q_{\rho,\theta}(s_{t},a)
\end{pmatrix}\Big]\Big].
\end{eqnarray*}

Therefore, $\mathcal H$ is negative-semi definite by the definition of the covariance notation ${\mathbb{V}}$ in (<ref>).

§.§ Algorithm Design and Implementation

Motivated by the landscape analysis in Section <ref>, we design an algorithm that iteratively updates $\rho$ and $\theta$ to maximize the likelihood function; see Algorithm <ref>.
At each iteration $k$, the value function $V_{\rho^k,\theta^k}$ is first computed by the soft Q iteration (see e.g. [Reddy et al., 2019])  in lines <ref>-<ref>, and the parameters $\rho^k,\theta^k$ are then updated in line <ref> using the gradient computed in line <ref>.

Maximum likelihood update

[1]
Initialize $\rho^0,\theta^0$.
$k=1,2,\cdots, K$
Set $Q_{\rho^k,\theta^k}^{0}(s,a)=\frac{1}{1-\rho^k} $ for all $(s,a)\in\calS\times\calA$,
and $V_{\rho^k,\theta^k}^{0}(s)=\frac{1}{1-\rho^k} $ for all $s\in\calS$.
$i=1,2,\cdots, I$

\begin{align*}
Q_{\rho^k,\theta^k}^{i}(s,a) = U_{\theta^k}(R(s,a)) +\rho^k\, \E_{s'\sim P(\cdot|s,a)}\Big[V_{\rho^k,\theta^k}^{i}(s')\Big].
\end{align*}

Compute $V_{\rho^k,\theta^k}^{I}$ using the soft Bellman equation (<ref>).
With the value of $V^{I}_{\rho^k,\theta^k}$, compute $\nabla \mathcal{L}(\rho^k,\theta^k)$ using (<ref>) and (<ref>).

Update $(\rho^{k+1},\theta^{k+1}) = (\rho^{k},\theta^{k}) + \zeta^k \nabla \mathcal{L}(\rho^k,\theta^k)$.

Exponential case. Do we keep this experiment? We first test the performance using the exponential utility function $U_\theta(x) = 1 - e^{-\theta x}$ under the generic MDP setting. Recall that we denote $(\Bar{\theta},\Bar{\rho})$ as the preference parameter from the client and under which the optimal policy is denoted as $\pi_{\Bar{\rho},\Bar{\theta}}:=\pi_{\Bar{\theta},\Bar{\rho}}$. The goal is to recover $\Bar{\theta}$ and $\Bar{\rho}$ by using $\pi_{\Bar{\rho},\Bar{\theta}}$ provided by the client.

In the experiment, we set $\calS = 20$, $\calA = 10$, $\Bar{\rho} = 0.5$, $\Bar{\theta} = 3$, and $\gamma = 0.6$. Before showing the performance, we visualize the log-likelihood function and its gradient in Figure <ref>. One can see that the likelihood function is locally concave in $\theta$ and $\rho$ in a sufficiently large area around $(\Bar{\theta},\Bar{\rho})$, enabling us to find the true parameters under fast convergence by Algorithm <ref>.

Visualization of the log-likelihood function and its gradients (Left columns: visualization with respect to $\theta$ (under $\rho = \Bar{\rho}$) Right columns: visualization with respect to $\rho$ (under $\theta = \Bar{\theta}$)).

When implementing Algorithm <ref>, we initialize $\theta$ and $\rho$ randomly with $\theta^0$ sampled uniformly from $[0,1]$ and $\rho^0$ sampled uniformly from $[0.1,0.2]$. As shown in Figure <ref>, both $\theta$ and $\rho$ converge to the ground-truth value within 25 iterations.

The convergence result of Algorithm <ref>. The left plot is the value of $\theta$ at each iteration, and the right plot is for $\rho$.

Numerical example one: Merton's problem

We implement the discrete-time version of Merton's problem introduced in Section <ref>. The price of the bond follows $S_{t+1}^0=S_t^0+r\,\Delta$ and the price of the stock follows $S_{t+1}-S_t = S_t(\nu\Delta+\sigma \sqrt{\Delta} B_t)$, where $B_t$ are iid sampled from $\mathcal{N}(0,1)$.
Denote $(\alpha_t,c_t)\in \mathcal{A}:=[0,1]\times[0,2]$ as the pair of the consumption-allocation policy at time $t$, then the wealth process follows:

\begin{eqnarray}
X_{t+1}-X_t = \Big[ X_t (\alpha_t \nu + (1-
\alpha_t)r)-c_t\Big] \Delta + X_t\alpha_t \sigma \sqrt{\Delta}B_t.
\end{eqnarray}

The client provides a time-homogeneous policy $\pi_{\Bar{\rho},\Bar{\theta}}\in \mathcal{P}(\mathcal{A})$ to the inference agent, which solves

\begin{eqnarray}
\sup_{\pi}   \mathbb{E}\Big[\sum_{t=1}^\infty (\Bar{\rho})^t\Big( 1-\exp(-\Bar{\theta}\,c_t)\Big)+\mathcal{H}\big(\pi(\cdot|X_t)\big)\Big)\Big]
\end{eqnarray}

with $c_t = c(X_t)$ and $\alpha_t = \alpha(X_t)$. Here both $\Bar{\rho}$ and $\Bar{\theta}$ are unknown.

Visualization of the log-likelihood function and its gradients (Left columns: visualization with respect to $\theta$ (under $\rho = \Bar{\rho}$). Right columns: visualization with respect to $\rho$ (under $\theta = \Bar{\theta}$)).

In the experiment, we set
$\Bar{\rho} = 0.3$, $\Bar{\theta} = 2$, $\gamma = 0.6$, $r=1.05$, $\Delta=1$, $\nu=1.06$, and $\sigma=0.05$.  We discretize and truncate the state space of the wealth process as $\calS = \{0.13,  0.39,  \dots,  2.23,  2.5\}$, with evenly distanced values such that $|\calS| = 10$. In addition, we discretize the joint space of the allocation and consumption processes as $\calA = \{0.1, 0.11, \dots, 0.98, 1\}\times\{0    , 0.22,\dots , 1.77, 2\}$, with evenly distanced values such that $|\calA|=50$.

We visualize the log-likelihood function and its gradient in Figure <ref>. One can see that the likelihood function is locally concave in $\theta$ and $\rho$  around $(\Bar{\theta},\Bar{\rho})$ in a sufficiently large area, enabling us to find the true parameters by Algorithm <ref> under fast convergence rate.

When implementing Algorithm <ref>, we initialize the parameters randomly with $\theta^0$ sampled uniformly from $[0,1]$ and $\rho^0$ sampled uniformly from $[0.1,0.2]$. We set the learning rate as $\zeta^k = \frac{1000}{k}$ and the total steps of the soft Q update as  $I=100$. As shown in Figure <ref>, both $\theta$ and $\rho$ converge to the ground-truth value within 100 iterations.

Additionally, we analyze the behaviors of the client under different $\Bar{\rho}$ values. Figure <ref> suggests that the client opts for an overall higher consumption when $\Bar{\rho} = 0.1$ and an overall lower consumption when $\Bar{\rho} = 0.75$, indicating a bigger emphasis on deferred outcomes for the latter case.

The convergence result of Algorithm <ref>. The left plot shows the value of $\theta$ at each iteration, while the right plot displays the values for $\rho$.

Visualization of the client's consumption policy. The left plot illustrates consumption at various wealth levels under $\Bar{\rho} = 0.1$, while the right plot corresponds to $\Bar{\rho} = 0.75$.

Numerical example two: Investment under unhedgeable risk.
We consider a more complex investment problem, where the price of the primitive asset is modeled as a diffusion process whose coefficients evolve according to a correlated diffusive factor [Zariphopoulou, 2001]. The price of the bond follows the same dynamics as in Example One:
On the other hand, the price of the stock follows
$$S_{t+1}-S_t = S_t(\nu(Y_t,t)\Delta+\sigma(Y_t,t) \sqrt{\Delta} B^1_t),$$

with $Y_t$ the “stochastic factor model” and it is assumed to satisfy
$$Y_{t+1}-Y_t = b(Y_t,t)\Delta + d(Y_t,t)\sqrt{\Delta} B^1_t.$$
Here $B^1_t$  and  $B^2_t$ are iid sampled from $\mathcal{N}(0,1)$. We assume the correlation between $B_t^1$ and $B_t^2$ is $\eta\in(0,1)$.

Consider a problem with only investment and no consumption. Then the wealth process follows:

\begin{eqnarray}
X_{t+1}-X_t = \Big[ X_t (\alpha_t \,\nu(t,Y_t) + (1-
\alpha_t)r)\Big] \Delta + X_t\alpha_t \sigma (t,Y_t)\,\sqrt{\Delta}B^1_t,
\end{eqnarray}

under the investment strategy $\alpha_t\in \mathcal{A}=[0,1]$. The client provides a time-homogeneous policy $\pi_{\Bar{\rho},\Bar{\theta}}(x,y)\in \mathcal{P}(\mathcal{A})$ to the inference agent, which solves

\begin{eqnarray}
\label{equ:unhegeable}
\sup_{\pi}   \mathbb{E}\Big[\sum_{t=1}^\infty (\Bar{\rho})^t \frac{1}{\Bar{\theta}_1 \Bar{\theta}_2}(X_t)^{\Bar{\theta}_1} (Y_t)^{\Bar{\theta}_2} + \mathcal{H}\Big(\pi(\cdot|(X_t,Y_t))\Big)\Big]
\end{eqnarray}

for some  $\Bar{\rho}>0$ and $\Bar{\theta}_1,\Bar{\theta}_2\in(0,1)$ that  are unknown to the inference agent.

Visualization of the log-likelihood function and its gradients (Left columns: visualization with respect to $\theta_1$ (under $\theta_2 = \Bar{\theta}_2$ and $\rho = \Bar{\rho}$).
Middle columns: visualization with respect to $\theta_2$ (under $\theta_1 = \Bar{\theta}_1$ and $\rho = \Bar{\rho}$).
Right columns: visualization with respect to $\rho$ (under $\theta_1 = \Bar{\theta}_1$ and $\theta_2 = \Bar{\theta}_2$).)

In the experiment, we discretize and truncate the state space for the wealth process and the stochastic factor model as $\calS=\{0.1, 0.7, 1.3, 1.9, 2.5\}\times \{0.1 , 0.32, 0.55 , 0.77, 1\}$,  with evenly distanced values such that $|\calS|=25$. In addition, we discretize the action space of the allocation process as $\calA= \{0.1  , 0.32, 0.55 , 0.77, 1 \}$, with evenly distanced values such that  $|\calA|=5$.
We set $r = 1.05$, $\Delta = 1$,
$\Bar{\theta_1} = 3$, $\Bar{\theta}_2 = 2$, $\Bar{\rho} = 0.3$, and $\gamma = 0.6$. For the drift and diffusion terms, we set $b(y,t) = -0.6y+0.2$,
$d(y,t)=0.3 y + 0.3$, $\nu(t,y)=y$, and $\sigma(t,y)=0.5y+0.3$.

As shown in Figure <ref>,
we visualize the log-likelihood function and its gradient. One can see that the likelihood function is locally concave in $\theta$ and $\rho$ in an area around $(\Bar{\theta},\Bar{\rho})$.

When implementing Algorithm <ref>, we initialize the parameters randomly with $\theta_1^0$, $\theta_2^0$ sampled uniformly from $[1,2]$ and $\rho^0$ sampled uniformly from $[0.1,0.2]$. We set the learning rate as $\zeta^k = \frac{1000}{\sqrt{k}}$ and the total steps of the soft Q update as $I=100$. As shown in Figure <ref>, both $\theta_1,\theta_2$ and $\rho$ converge to the ground-truth values within 1500 iterations.

The convergence result of Algorithm <ref>. The left plot shows the value of $\theta_1$ at each iteration, the middle plot is for $\theta_2$, and the right plot is for $\rho$.

Furthermore, Figure <ref> illustrates the client's investment allocation policy $\alpha$ across various wealth levels (under fixed factor  value $1$), considering $\Bar{\rho} = 0.1$ and $\Bar{\rho} = 0.75$. The influence of $\Bar{\rho}$ on the investment decisions in this example is less pronounced compared to Merton's problem. This difference arises because, in Merton's problem, the client confronts a trade-off between higher consumption for instantaneous rewards and lower consumption for better future rewards. Conversely, the client addressing (<ref>) strives for a higher $X_t$ regardless of her $\Bar{\rho}$. Our algorithm consistently finds the optimal parameters, although the convergence speed here is slower compared to that for Merton's problem due to the above-mentioned reasons.

Visualization of the client's allocation policy (under fixed factor value $1$). The left plot illustrates her allocation at various wealth levels with $\Bar{\rho} = 0.1$, while the right plot is for $\Bar{\rho} = 0.75$.

§ STATEMENTS AND DECLARATIONS

Funding. R. X. is partially supported by the NSF CAREER award DMS-2339240 and a JP Morgan Faculty Research Award. H. C. is partially supported by the departmental startup fund.

Competing interests. The authors have no relevant financial or non-financial interests to disclose.

Author contributions. All authors have contributed equally to this study.

[Abbeel and Ng, 2004]
P. Abbeel and A. Y. Ng.
Apprenticeship learning via inverse reinforcement learning.
In Proceedings of the twenty-first international conference on
Machine learning, page 1, 2004.

[Alsabah et al., 2021]
H. Alsabah, A. Capponi, O. Ruiz Lacedelli, and M. Stern.
Robo-advising: Learning investors’ risk preferences via portfolio
Journal of Financial Econometrics, 190 (2):0
369–392, 2021.

[Amin and Singh, 2016]
K. Amin and S. Singh.
Towards resolving unidentifiability in inverse reinforcement
arXiv preprint arXiv:1601.06569, 2016.

[Amin et al., 2017]
K. Amin, N. Jiang, and S. Singh.
Repeated inverse reinforcement learning.
Advances in Neural Information Processing Systems,
30:0 1815–1824, 2017.

[Bäuerle and Rieder, 2014]
N. Bäuerle and U. Rieder.
More risk-sensitive markov decision processes.
Mathematics of Operations Research, 390 (1):0
105–120, 2014.

[Bjork and Murgoci, 2010]
T. Bjork and A. Murgoci.
A general theory of markovian time inconsistent stochastic control
Available at SSRN 1694759, 2010.

[Björk and Murgoci, 2014]
T. Björk and A. Murgoci.
A theory of markovian time-inconsistent stochastic control in
discrete time.
Finance and Stochastics, 18:0 545–592, 2014.

[Björk et al., 2017]
T. Björk, M. Khapko, and A. Murgoci.
On time-inconsistent stochastic control in continuous time.
Finance and Stochastics, 21:0 331–360, 2017.

[Bloem and Bambos, 2014]
M. Bloem and N. Bambos.
Infinite time horizon maximum causal entropy inverse reinforcement
In 53rd IEEE conference on decision and control, pages
4911–4916. IEEE, 2014.

[Boularias et al., 2011]
A. Boularias, J. Kober, and J. Peters.
Relative entropy inverse reinforcement learning.
In Proceedings of the Fourteenth International Conference on
Artificial Intelligence and Statistics, pages 182–189. JMLR Workshop and
Conference Proceedings, 2011.

[Boyd et al., 1994]
S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan.
Linear matrix inequalities in system and control theory.
SIAM, 1994.

[Cao et al., 2021]
H. Cao, S. Cohen, and L. Szpruch.
Identifiability in inverse reinforcement learning.
Advances in Neural Information Processing Systems,
34:0 12362–12373, 2021.

[Capponi and Zhang, 2024]
A. Capponi and Y. Zhang.
A continuous time framework for sequential goal-based wealth
Management Science, 2024.

[Capponi et al., 2022]
A. Capponi, S. Olafsson, and T. Zariphopoulou.
Personalized robo-advising: Enhancing investment through client
Management Science, 680 (4):0 2485–2512,

[Chewning et al., 2012]
B. Chewning, C. L. Bylund, B. Shah, N. K. Arora, J. A. Gueguen, and G. Makoul.
Patient preferences for shared decisions: a systematic review.
Patient education and counseling, 860 (1):0
9–18, 2012.

[Christiano et al., 2017]
P. F. Christiano, J. Leike, T. Brown, M. Martic, S. Legg, and D. Amodei.
Deep reinforcement learning from human preferences.
Advances in neural information processing systems, 30, 2017.

[Cox et al., 2014]
A. M. Cox, D. Hobson, and J. OBłÓJ.
Utility theory front to back—inferring utility from agents'choices.
International Journal of Theoretical and Applied Finance,
170 (03):0 1450018, 2014.

[Dai et al., 2023]
M. Dai, Y. Dong, and Y. Jia.
Learning equilibrium mean-variance strategy.
Mathematical Finance, 330 (4):0 1166–1212,

[Derbaix and Abeele, 1985]
C. Derbaix and P. V. Abeele.
Consumer inferences and consumer preferences. the status of cognition
and consciousness in consumer behavior theory.
International Journal of Research in Marketing, 20
(3):0 157–174, 1985.

[Dong and Wang, 2024]
C. Dong and Y. Wang.
Towards generalized inverse reinforcement learning.
arXiv preprint arXiv:2402.07246, 2024.

[Dybvig and Rogers, 1997]
P. H. Dybvig and L. C. G. Rogers.
Recovery of preferences from observed wealth in a single realization.
The Review of Financial Studies, 100 (1):0
151–174, 1997.
ISSN 08939454, 14657368.
URL <http://www.jstor.org/stable/2962259>.

[D’Acunto and Rossi, 2021]
F. D’Acunto and A. G. Rossi.
Springer, 2021.

[D’Acunto et al., 2019]
F. D’Acunto, N. Prabhala, and A. G. Rossi.
The promises and pitfalls of robo-advising.
The Review of Financial Studies, 320 (5):0
1983–2020, 2019.

[Ekeland and Lazrak, 2010]
I. Ekeland and A. Lazrak.
The golden rule when preferences are time inconsistent.
Mathematics and Financial Economics, 4:0 29–55, 2010.

[El Karoui and Mrad, 2021]
N. El Karoui and M. Mrad.
Recover dynamic utility from observable process: Application to the
economic equilibrium.
SIAM Journal on Financial Mathematics, 120
(1):0 189–225, 2021.

[El Karoui et al., 2019]
N. El Karoui, C. Hillairet, and M. Mrad.
Construction of an aggregate consistent utility, without pareto
optimality. application to long-term yield curve modeling.
In Frontiers in Stochastic Analysis–BSDEs, SPDEs and their
Applications: Edinburgh, July 2017 Selected, Revised and Extended
Contributions 8, pages 169–199. Springer, 2019.

[El Karoui et al., 2024]
N. El Karoui, C. Hillairet, and M. Mrad.
Bi-revealed utilities in a defaultable universe: A new point of view
on consumption.
Probability, Uncertainty and Quantitative Risk, 90
(1):0 13–34, 2024.

[Finn et al., 2016]
C. Finn, P. Christiano, P. Abbeel, and S. Levine.
A connection between generative adversarial networks, inverse
reinforcement learning, and energy-based models.
arXiv preprint arXiv:1611.03852, 2016.

[Fleming and Soner, 2006]
W. H. Fleming and H. M. Soner.
Controlled Markov processes and viscosity solutions,
volume 25.
Springer Science & Business Media, 2006.

[Fu et al., 2018]
J. Fu, K. Luo, and S. Levine.
Learning robust rewards with adverserial inverse reinforcement
In International Conference on Learning Representations, 2018.

[Garg et al., 2021]
D. Garg, S. Chakraborty, C. Cundy, J. Song, and S. Ermon.
Iq-learn: Inverse soft-q learning for imitation.
Advances in Neural Information Processing Systems,
34:0 4028–4039, 2021.

[Haarnoja et al., 2017]
T. Haarnoja, H. Tang, P. Abbeel, and S. Levine.
Reinforcement learning with deep energy-based policies.
In Proceedings of the 34th International Conference on Machine
Learning - Volume 70, ICML'17, page 1352–1361. JMLR.org, 2017.

[Hernández and Possamaï, 2023]
C. Hernández and D. Possamaï.
Me, myself and i: a general theory of non-markovian time-inconsistent
stochastic control for sophisticated agents.
The Annals of Applied Probability, 330 (2):0
1396–1458, 2023.

[Ho and Ermon, 2016]
J. Ho and S. Ermon.
Generative adversarial imitation learning.
In Advances in neural information processing systems, pages
4565–4573, 2016.

[Hu et al., 2012]
Y. Hu, H. Jin, and X. Y. Zhou.
Time-inconsistent stochastic linear–quadratic control.
SIAM journal on Control and Optimization, 500
(3):0 1548–1572, 2012.

[Hu et al., 2017]
Y. Hu, H. Jin, and X. Y. Zhou.
Time-inconsistent stochastic linear-quadratic control:
Characterization and uniqueness of equilibrium.
SIAM Journal on Control and Optimization, 550
(2):0 1261–1279, 2017.

[Jin and Yu Zhou, 2008]
H. Jin and X. Yu Zhou.
Behavioral portfolio selection in continuous time.
Mathematical Finance: An International Journal of Mathematics,
Statistics and Financial Economics, 180 (3):0 385–426,

[Kalman, 1964]
R. E. Kalman.
When Is a Linear Control System Optimal?
Journal of Basic Engineering, 860 (1):0
51–60, 03 1964.

[Karnam et al., 2017]
C. Karnam, J. Ma, and J. Zhang.
Dynamic approaches for some time-inconsistent optimization
The Annals of Applied Probability, 270 (6):0
3435 – 3477, 2017.
URL <https://doi.org/10.1214/17-AAP1284>.

[Keeney and Raiffa, 1976]
R. L. Keeney and H. Raiffa.
Decisions with Multiple Objectives: Preferences and Value
Wiley, 1976.

[Kim et al., 2021]
K. Kim, S. Garg, K. Shiragur, and S. Ermon.
Reward identification in inverse reinforcement learning.
In M. Meila and T. Zhang, editors, Proceedings of the 38th
International Conference on Machine Learning, volume 139 of
Proceedings of Machine Learning Research, pages 5496–5505. PMLR,
18–24 Jul 2021.
URL <http://proceedings.mlr.press/v139/kim21c.html>.

[Levine et al., 2011]
S. Levine, Z. Popovic, and V. Koltun.
Nonlinear inverse reinforcement learning with gaussian processes.
Advances in neural information processing systems,
24:0 19–27, 2011.

[Linos, 1983]
P. L. Linos.
Optimal control of diffustion processes and hamilton-jacobi-bellman
equations part i: the dynamic programming principle and application.
Communications in Partial Differential Equations, 80
(10):0 1101–1174, 1983.
URL <https://doi.org/10.1080/03605308308820297>.

[Lions, 1982]
P. L. Lions.
Optimal stochastic control of diffusion type processes and
hamilton-jacobi-bellman equations.
In W. H. Fleming and L. G. Gorostiza, editors, Advances in
Filtering and Optimal Stochastic Control, pages 199–215, Berlin,
Heidelberg, 1982. Springer Berlin Heidelberg.
ISBN 978-3-540-39517-1.

[Musiela and Zariphopoulou, 2006]
M. Musiela and T. Zariphopoulou.
Investments and forward utilities.
preprint, 2006.

[Musiela and Zariphopoulou, 2007]
M. Musiela and T. Zariphopoulou.
Investment and Valuation Under Backward and Forward Dynamic
Exponential Utilities in a Stochastic Factor Model, pages 303–334.
Birkhäuser Boston, Boston, MA, 2007.
ISBN 978-0-8176-4545-8.
URL <https://doi.org/10.1007/978-0-8176-4545-8_16>.

[Ng et al., 2000]
A. Y. Ng, S. Russell, et al.
Algorithms for inverse reinforcement learning.
In Icml, volume 1, page 2, 2000.

[Nicole and Mohamed, 2013]
E. K. Nicole and M. Mohamed.
An exact connection between two solvable sdes and a nonlinear utility
stochastic pde.
SIAM Journal on Financial Mathematics, 40
(1):0 697–736, 2013.
URL <https://doi.org/10.1137/10081143X>.

[Pollak, 1968]
R. A. Pollak.
Consistent Planning1.
The Review of Economic Studies, 350 (2):0
201–208, 04 1968.
ISSN 0034-6527.
URL <https://doi.org/10.2307/2296548>.

[Reddy et al., 2019]
S. Reddy, A. D. Dragan, and S. Levine.
Sqil: Imitation learning via reinforcement learning with sparse
arXiv preprint arXiv:1905.11108, 2019.

[Richesson and Vehik, 2010]
R. Richesson and K. Vehik.
Patient registries: utility, validity and inference.
Rare diseases epidemiology, pages 87–104, 2010.

[Richter, 1966]
M. K. Richter.
Revealed preference theory.
Econometrica, 340 (3):0 635–645, 1966.
ISSN 00129682, 14680262.
URL <http://www.jstor.org/stable/1909773>.

[Rossi and Utkus, 2020]
A. G. Rossi and S. P. Utkus.
Who benefits from robo-advising? evidence from machine learning.
Evidence from Machine Learning (March 10, 2020), 2020.

[Russell, 1998]
S. Russell.
Learning agents for uncertain environments.
In Proceedings of the Eleventh Annual Conference on
Computational Learning Theory, pages 101–103, 1998.

[Samuelson, 1948]
P. A. Samuelson.
Consumption theory in terms of revealed preference.
Economica, 150 (60):0 243–253, 1948.
ISSN 00130427, 14680335.
URL <http://www.jstor.org/stable/2549561>.

[Sargent, 1978]
T. J. Sargent.
Estimation of dynamic labor demand schedules under rational
Journal of Political Economy, 860 (6):0
1009–1044, 1978.

[Schlaginhaufen and Kamgarpour, 2023]
A. Schlaginhaufen and M. Kamgarpour.
Identifiability and generalizability in constrained inverse
reinforcement learning.
In International Conference on Machine Learning, pages
30224–30251. PMLR, 2023.

[Shin and Yu, 2021]
J. Shin and J. Yu.
Targeted advertising and consumer inference.
Marketing Science, 400 (5):0 900–922, 2021.

[Strotz, 1955]
R. H. Strotz.
Myopia and inconsistency in dynamic utility maximization.
The Review of Economic Studies, 230 (3):0
165–180, 1955.
ISSN 00346527, 1467937X.
URL <http://www.jstor.org/stable/2295722>.

[Wang and Yu, 2021]
H. Wang and S. Yu.
Robo-advising: Enhancing investment with inverse optimization and
deep reinforcement learning.
In 2021 20th IEEE international conference on machine learning
and applications (ICMLA), pages 365–372. IEEE, 2021.

[Wulfmeier et al., 2015]
M. Wulfmeier, P. Ondruska, and I. Posner.
Maximum entropy deep inverse reinforcement learning.
arXiv preprint arXiv:1507.04888, 2015.

[Yong, 2012]
J. Yong.
Time-inconsistent optimal control problems and the equilibrium hjb
Mathematical Control and Related Fields, 20
(3):0 271–329, 2012.
ISSN 2156-8472.

[Zariphopoulou, 2001]
T. Zariphopoulou.
A solution approach to valuation with unhedgeable risks.
Finance and stochastics, 5:0 61–82, 2001.

[Zeng et al., 2022]
S. Zeng, C. Li, A. Garcia, and M. Hong.
Maximum-likelihood inverse reinforcement learning with finite-time
Advances in Neural Information Processing Systems,
35:0 10122–10135, 2022.

[Ziebart, 2010]
B. D. Ziebart.
Modeling Purposeful Adaptive Behavior with the Principle of
Maximum Causal Entropy.
PhD thesis, Carnegie Mellon University, 2010.

[Ziebart et al., 2008]
B. D. Ziebart, A. L. Maas, J. A. Bagnell, A. K. Dey, et al.
Maximum entropy inverse reinforcement learning.
In Aaai, volume 8, pages 1433–1438. Chicago, IL, USA, 2008.

§.§ IRL References

* Maximum margin:
Ng et al., 2000.

* Bayesian inverse RL:
<cit.>, <cit.>, <cit.>.

* Maximum entropy:
Ziebart et al., 2008, <cit.>, Reddy et al., 2019, Finn et al., 2016, Zeng et al., 2022.

* Survey:

§.§ Bayesian Inverse Reinforcement Learning Algorithm Design

Given a $\texttt{MDP}(\calS,\calA,P,R)$, we consider the following problem

* The optimal Q function is given by

\begin{eqnarray}
Q^*(s,a,\rho,U)=\max_\pi\mathbb{E}^\pi\Big[\sum_{t=1}^\infty \rho^ t U(r_t(s_t,a^*_t))|s,a\Big].
\end{eqnarray}

* The agent takes a stochastic policy with a Boltzman-type distribution:

\begin{align}
\pi(a|s,\rho,U) = \frac{1}{Z}e^{\alpha Q^*(s,a,\rho,U)},
\end{align}

where $\alpha$ is a known behavioral constant of the agent and  $Z$ is the appropriate normalizing factor. This is equivalent to the setting of a MDP with an entropy regularization.

* Let $\tau=\{s_t,a_t\}_{t=0}^T$ denote the expected trajectory of the agent with $\rho$ and $U$. Bayesian inference is to calculate the posterior:

\begin{align}
\prob(\tau|\rho,U) &= \Pi_t \prob(s_t,a_t|\rho,U)\\
&= \frac{1}{Z}e^{\alpha \sum_t Q^*(s_t,a_t,\rho,U)}.
\end{align}

* The posterior is given by:

\begin{align}
\prob(\rho,U|\tau) &= \frac{\prob(\tau|\rho,U) \prob(\rho,U) }{\prob(\tau)}\\
&=\frac{1}{Z'}e^{\alpha \sum_t Q^*(s_t,a_t,\rho,U)}\prob(\rho,U),
\end{align}

where $\prob(\rho,U)$ is the prior distribution.

* Assume that $U$ is parameterized by $\theta\in\mathbb{R}$. Then the above posterior is:

\begin{align}
\prob(\rho,\theta|\tau)
&=\frac{1}{Z'}e^{\alpha \sum_t Q^*(s_t,a_t,\rho,\theta)}\prob(\rho,\theta) ,
\end{align}

and we can assume the prior $\prob(\rho,\theta) $ is a uniform distribution.

* The algorithm is to find the $\rho,\theta$ that maximizes the posterior. Let $\hat{\rho}$. $\hat{\theta}$ denotes the optimal solution with maximum posterior.

* It was shown that Bayesian inverse RL with a Laplacian prior is equivalent to the maximum margin IRL.

\begin{eqnarray*}
\nabla^2_\theta Q_{\rho,\theta}(s_t,a_t)
&=& \nabla^2_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[\nabla^2_\theta V_{\rho,\theta}(s_{t+1}) \Big]\\
&=& \nabla^2_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[\nabla^2_\theta \log({\SumInt_a e^{Q_{\rho,\theta}(s_{t+1},a)}}  ) \Big]\\
&=& \nabla^2_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\nabla_\theta\Big[
\SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}}{\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}}\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a) \Big]
\end{eqnarray*}

Note that

\begin{eqnarray}
&&  \nabla_\theta\Big[
\SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}}{\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}}\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a) \Big] \nonumber\\
&=& \SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}}{\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}} (\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a))^2 - \SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a')}{(\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')})^2}\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a)\nonumber\\
&&+\SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}}{\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}}\nabla^2_\theta  Q_{\rho,\theta}(s_{t+1},a)\nonumber\\
&=&  \frac{1}{2}  \frac{1}{(\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')})^2} \sum_{a,a'}e^{Q_{\rho,\theta}(s_{t+1},a)} e^{Q_{\rho,\theta}(s_{t+1},a')}\Big(\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a')-\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a)\Big)^2\nonumber\\
&&+\SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}}{\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}}\nabla^2_\theta  Q_{\rho,\theta}(s_{t+1},a)\nonumber\\
&\ge &\SumInt_a \frac{e^{Q_{\rho,\theta}(s_{t+1},a)}}{\sum_{a'} e^{Q_{\rho,\theta}(s_{t+1},a')}}\nabla^2_\theta  Q_{\rho,\theta}(s_{t+1},a).\label{equ:likelihood_hesian11}
\end{eqnarray}

Applying (<ref>) recursively yields:

\begin{align}
\nabla^2_\theta Q_{\rho,\theta}(s_t,a_t) \ge \E_{\tau\sim \pi_{\rho,\theta}}\Big[\sum_{t'=t}^\infty \rho^{t'-t}\nabla^2_\theta U_\theta(R(s_{t'},a_{t'}))\Big|s_t,a_t  \Big].
\end{align}

On the other hand,

\begin{eqnarray*}
\nabla^2_\theta Q_{\rho,\theta}(s_t,a_t)
&=& \nabla^2_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\nabla_\theta\Big[
\SumInt_a \pi_{\theta,\rho}(s_{t+1},a)\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a) \Big]\\
&=& \nabla^2_\theta U_\theta(R(s_t,a_t))+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[
\SumInt_a \nabla_\theta\pi_{\theta,\rho}(s_{t+1},a)\nabla_\theta  Q_{\rho,\theta}(s_{t+1},a) \Big]\\
&&+\rho \E_{s_{t+1}\sim P(\cdot|s_t,a_t)}\Big[
\SumInt_a \pi_{\theta,\rho}(s_{t+1},a)\nabla^2_\theta  Q_{\rho,\theta}(s_{t+1},a) \Big]
\end{eqnarray*}

$i=1,2,\cdots, I$

\begin{align*}
Q_{\rho^k,\theta^k}^{i+1}(s,a) =
+\rho^k \E\Big[V_{\rho^k,\theta^k}^{i+1}(s')\Big]
\end{align*}