\section{Introduction}
Optimal transport (OT) is a method for quantifying the discrepancy between two probability distributions. The origins of OT can be traced back to Monge's mass transfer problem in the 18th century. However, it was not until the advent of linear programming in the mid 20th century that OT was revitalized by Kantorovich. In the 1990s, \cite{brenier1991polar} and \cite{gangbo1996geometry} provided new insights into OT through the lens of convex analysis and geometry. Nowadays, OT is widely recognized as an interdisciplinary tool with a wide range of applications in fields such as mathematics \citep{villani2009optimal}, machine learning \citep{arjovsky2017wasserstein,peyre2019computational}, economics \citep{galichon2016optimal}, statistics \citep{torous2021optimal}, biology \citep{schiebinger2019optimal}, operations research (OR), and management science (MS).

In particular, OT has a close connection with linear programming, which is a widely used OR tool. Additionally, if we regard the two marginals in OT as distributions of workers and jobs, OT becomes the optimal assignment problem in OR/MS \citep[Chapter 3]{galichon2016optimal}. Another major application of OT in OR/MS is robust optimization. Many decision-making tasks are based on a nominal probability measure, which may not be accurate. Wasserstein distances were used in distributionally robust optimization by considering alternative measures within a certain Wasserstein distance from the nominal distribution and selecting decision rules that perform well even in the worst-case scenario \citep{blanchet2019quantifying,blanchet2021dist,gao2016dist,mohajerin2018data,nguyen2022dist,kuhn2019tutorial,blanchet2021stat}.

While most literature on OT has focused on static data without a chronological structure, temporal data are ubiquitous in finance, statistics, and optimization. To address this, \cite{lassalle2013causal} introduced causal optimal transport (COT), which imposes a causality constraint on the transport plans between discrete stochastic processes. In a less formal sense, COT requires that when the past of one process, $X$, is given, the past of another process, $Y$, should be independent of the future of $X$ under the transport plan. COT has found applications in mathematical finance \citep{backhoff2020adapted}, video modeling \citep{xu2020cot}, mean-field games \citep{acciaio2021cournot}, and stochastic optimization \citep{pflug2012distance,acciaio2020causal}.

 However, when considering optimal decisions over time, a well-documented and classic human tendency is time inconsistency. There are empirical findings that human behavior may contradict the axiom of rationality and deviate from the previous optimal strategy at later stages. This is observed and studied in various areas such as non-exponential discounting \citep{laibson1997golden}, state-dependent preferences and mean-variance portfolio selection \citep{bjork14FS,bjork14MF,bjork17FS,han2021robust,dai2021dynamic,kovavcova2021time}, optimal stopping \citep{bayraktar2021equilibrium,pichler2022risk}, robust decision-making \citep{epstein2022optimal,arora2022}, prospect theory \citep{kahneman1979prospect}, probability distortion \citep{ma2021time}, and gambling behaviors \citep{barberis2012model}.

In this paper, we generalize the classic Monge-Kantorovich problem to a dynamic setting with time-inconsistent costs. Suppose transport plans are bi-causal in the sense that two processes are causal with respect to each other. The cost function is allowed to include time-inconsistent preferences such as nonlinear objectives and state dependence. Moreover, we notice that regularization with general $f$-divergence can also be a new source of time inconsistency. To address the self-contradictory behavior and the violation of dynamic programming principle, we adopt the subgame Nash equilibria solution as in \cite{bjork14FS,basak10}. This concept, given formally in Definition \ref{def:equitrans}, is well-defined and easy to understand in a discrete-time setting. Informally, the agent who is aware of time inconsistency can consider herself at different future times as different planners. A bi-causal equilibrium transport plan $\pi^*$ with successive regular kernels $\{ \bar{\pi}^*, \pi^*_1, ..., \pi^*_{T-1} \}$ should satisfy that, if all future selves at time $t+1, ..., T-1$ agree to use $\pi^*_{t+1}, ..., \pi^*_{T-1}$, then it is optimal for the planner at time $t$ to use $\pi^*_t$. Thus, there is no motivation for the planner to deviate from $\pi^*$ at any time $t$. The first main contribution of this paper is the characterization of the equilibrium transport under an extended dynamic programming principle, see Theorems \ref{thm:nlin_V} and \ref{thm:state_V}. With a strict quasi-convexity assumption, we also prove the existence and uniqueness of the equilibrium transport.

Our framework introduces a new methodology with a wide range of potential interdisciplinary applications. As an example, we consider the matching between employers and employees, quantified by firm net sales and wages. This problem has been studied extensively for the labor market of top executives \citep{gabaix2008has,demerjian2012quantifying,taylor2013ceo,ai2021unified}. A stylized empirical effect, known as Roberts' law in \cite{gabaix2008has}, states that larger firms pay more to top executives represented by chief executive officers (CEOs). This firm-wage matching is also called positive assortative matching (PAM), where firms ranked number $n$ in terms of sizes (net sales or market values) should pay qualified executives with wages ranked number $n$. However, empirical data do not always show a perfect correspondence between firm sizes and wages. This raises questions why PAM may fail and some high-ranked firms pay less to executives.

We conjecture that there are incentives for firms to maintain similar firm-wage matching over time, regardless of whether it aligns with PAM or not. These incentives can stem from regulatory requirements, momentum in financial performance, and status quo bias. Regulators require public companies to provide a year-over-year comparison of executive compensations, implying that regulators and stakeholders want to detect abnormal changes in wages. Additionally, wages are often tied to financial performance, such as stock prices, which may also exhibit momentum. Furthermore, people tend to keep the current situation stable and may be resistant to change, leading to a status quo bias in the firm-wage matching. Motivated by these factors, in Section \ref{sec:model}, we propose a novel job market model that incorporates a state-dependent term capturing the potential persistence or inertia of maintaining previous matching over time. This allows us to quantify the underlying incentives that contribute to the observed patterns in the labor market.

We call a job market more efficient if it is more consistent with the Roberts' law. Through our empirical analysis in Sections \ref{sec:exec} and \ref{sec:prof}, we find that a job market is less efficient when its inertia represented by state dependence is stronger. Specifically, over the last five years, the Consumer Discretionary, Real Estate, and Information Technology industries have the weakest correlations between net sales and wages, while they exhibit stronger inertia or state dependence. In addition to the well-studied executive job market, we also examine the academic job market which received less attention in the literature.

Thanks to the Transparent California project, we are able to analyze wage data from the University of California (UC) system. In Section \ref{sec:prof}, wages for faculty in Business, Economics, and Engineering departments, including full professors, associate professors, and assistant professors, are found to be almost perfectly consistent with U.S. News National University Rankings. There is no evidence of any inertia for faculty wages. However, we find that the median wages of postdocs from all departments are highly inconsistent with university rankings. The state dependence or inertia in the postdoc job market is even stronger than that of top executives in any industry. In contrast to the attention on executive compensation, the inefficiency and inertia in the postdoc job market are rarely discussed in the literature.

Our paper also offers practical techniques to enhance the robustness of the results. In Section \ref{sec:num_group}, we adopt a data-driven method to determine the optimal number of clusters in the discretization of the OT algorithm. Additionally, in Section \ref{sec:valid}, we employ a model validation procedure by testing the model with synthetic perfectly matched firm-wage paths. Besides, we utilize the bootstrapping technique for data augmentation. These procedures ensure that our findings are robust across multiple simulation runs and varying numbers of clusters, despite the inherent noise in the net sale and wage data. The Python code is publicly available on the GitHub repository: \url{https://github.com/hanbingyan/equitrans}.

The rest of the paper is organized as follows. We introduce the notation and bi-causal OT in Sections \ref{sec:notation} and \ref{sec:bic}. Section \ref{sec:ti} presents equilibrium transport and characterization with extended dynamic programming principle. Section \ref{sec:model} introduces a new job market model with inertia measured by state dependence. Sections \ref{sec:exec} and \ref{sec:prof} apply the model to executive and academic job markets. Technical proofs are deferred to Appendix \ref{sec:append}. Appendix \ref{sec:supp_figs} gives some supplementary Tables and Figures.

\subsection{Notation}\label{sec:notation}
Denote the finite number of periods as $T$. For each $t \in \{1, 2, ..., T-1, T\}$, suppose ${\mathcal X}_t$ is a closed (but not necessarily bounded) subset of $\mathbb{R}^d$. ${\mathcal X}_t$ is interpreted as the range of the process at time $t$. ${\mathcal X}_{1:T} := {\mathcal X}_1 \times \ldots \times {\mathcal X}_T$ is a closed subset of $\mathbb{R}^{T \times d}$. With some abuse of notation, we simply write ${\mathcal X}$ while we mean ${\mathcal X}_{1:T}$. With a fixed $p \in [1, \infty)$, we introduce the metric as $d_{\mathcal X}(x_{1:T}, x'_{1:T}) = \left[ \sum^T_{t=1} d_{{\mathcal X}_t}(x_t, x'_t)^p \right]^{1/p}$ for $x, x' \in {\mathcal X}$, where $d_{{\mathcal X}_t}(\cdot, \cdot)$ is a metric on ${\mathcal X}_t$.

Denote the set of all Borel probability measures on ${\mathcal X}$ as ${\mathcal P}({\mathcal X})$. The Wasserstein space of order $p$ is given by
\begin{equation*}
	{\mathcal P}_p({\mathcal X}) := \left\{ \mu \in {\mathcal P}({\mathcal X}) \Big| \int_{{\mathcal X}} d_{\mathcal X}(x_{1:T}, \bar{x}_{1:T})^p \mu(dx) < \infty \right\}
\end{equation*}
for some fixed $\bar{x}_{1:T} \in {\mathcal X}$.

The Wasserstein distance from the OT theory has been widely used as a metric on probability spaces. Consider two Borel probability measures $\mu, \mu' \in {\mathcal P}({\mathcal X})$. A coupling $\gamma \in {\mathcal P}({\mathcal X} \times {\mathcal X})$ is a Borel probability measure that admits $\mu$ and $\mu'$ as its marginals on ${\mathcal X}$. Denote $\Pi(\mu, \mu')$ as the set of all the couplings. The Wasserstein distance of order $p$ between $\mu$ and $\mu'$ is given by
\begin{equation}\label{eq:Wp}
	W_p(\mu, \mu') = \left( \inf_{\gamma \in \Pi(\mu, \mu')} \int_{{\mathcal X} \times {\mathcal X}} d_{\mathcal X}(x, x')^p \gamma (dx, dx') \right)^{1/p}.
\end{equation}
Let
\begin{equation*}
	C_p({\mathcal X}) := \left\{ f \in C({\mathcal X}) \Big|  \exists\, \text{ constant } C > 0, \; |f(x)| \leq C(1 + d_{\mathcal X}(x, \bar{x})^p)\right\}
\end{equation*}
be the set of continuous functions from ${\mathcal X}$ to $\mathbb{R}$ with a growth rate of order $p$ in $d_{\mathcal X}(x, \bar{x})$. Consider a sequence of probability measures $\{ \mu^k \}_{k \in {\mathbb N}}$ in ${\mathcal P}_p({\mathcal X})$ and $\mu$ as another probability measure in ${\mathcal P}_p({\mathcal X})$. Recall that $W_p(\mu^k, \mu) \rightarrow 0$ is equivalent to $\mu_k$ converges weakly in ${\mathcal P}_p({\mathcal X})$ to $\mu$; see \citet[Definition 6.7 and Theorem 6.8]{villani2009optimal}.

In the OT theory, we need to consider another closed set ${\mathcal Y} = {\mathcal Y}_{1:T} = {\mathcal Y}_1 \times ... \times {\mathcal Y}_T$. Similarly, we denote $d_{\mathcal Y}(y_{1:T}, y'_{1:T}) = \left[ \sum^T_{t=1} d_{{\mathcal Y}_t}(y_t, y'_t)^p \right]^{1/p}$ for $y, y' \in {\mathcal Y}$ and $d_{{\mathcal Y}_t}(\cdot, \cdot)$ is a metric on ${\mathcal Y}_t$. Under the metric $d_{\mathcal Y}$, we denote ${\mathcal P}_p({\mathcal Y})$ as the Wasserstein space of order $p$, $W_p(\nu, \nu')$ as the Wasserstein metric on ${\mathcal P}_p({\mathcal Y})$, and $C_p({\mathcal Y})$ as the set of continuous functions with a growth rate of order $p$.

For the product space ${\mathcal X} \times {\mathcal Y}$, we define the corresponding Wasserstein space of order $p$ as
\begin{equation*}
	{\mathcal P}_p({\mathcal X} \times {\mathcal Y}) := \left\{ \pi \in {\mathcal P}({\mathcal X} \times {\mathcal Y}) \Big| \int_{{\mathcal X} \times {\mathcal Y}} d((x, y), (\bar{x}, \bar{y}))^p \pi(dx, dy) < \infty \right\},
\end{equation*}
with the metric $d$ given by
\begin{align*}
	d((x, y), (\bar{x}, \bar{y})) = [d_{{\mathcal X}}(x, \bar{x})^p + d_{{\mathcal Y}}(y, \bar{y})^p]^{1/p}.
\end{align*}
For two probability measures $\pi, \pi' \in {\mathcal P}_p({\mathcal X} \times {\mathcal Y})$, the Wasserstein distance between them is denoted by
\begin{equation*}
	W_p(\pi, \pi') = \left( \inf_{\gamma \in \Pi(\pi, \pi')} \int_{{\mathcal X} \times {\mathcal Y} \times {\mathcal X} \times {\mathcal Y}} d((x, y), (x', y'))^p \gamma (dx, dy, dx', dy') \right)^{1/p}.
\end{equation*}
Similarly, let us introduce
\begin{align*}
	 C_p({\mathcal X}_{1:T} \times {\mathcal Y}_{1:T}) := \Big\{ f \in C({\mathcal X}_{1:T} \times {\mathcal Y}_{1:T}) \Big|  \exists\, & \text{ constant } C > 0, \; |f(x, y)| \leq C(1 + d((x, y), (\bar{x}, \bar{y}))^p ) \\
	& = C \big[ 1 + \sum^T_{t=1} d_{{\mathcal X}_t}(x_t, \bar{x}_t)^p + \sum^T_{t=1} d_{{\mathcal Y}_t}(y_t, \bar{y}_t)^p \big] \Big\}.
\end{align*}
For later use, we denote the sum of Wasserstein distances between $\mu, \mu' \in {\mathcal P}({\mathcal X})$ and  $\nu, \nu' \in {\mathcal P}({\mathcal Y})$ as
\begin{align*}
	W^\oplus_p ((\mu, \nu), (\mu', \nu')) := W_p (\mu, \mu') + W_p (\nu, \nu').
\end{align*}

In this paper, we often consider metrics and spaces with sub time intervals such as $1, ..., t$. With a slight abuse of notation, we interpret $d_{{\mathcal X}}(x_{1:t}, x'_{1:t}) = \left[ \sum^t_{s=1} d_{{\mathcal X}_s}(x_s, x'_s)^p \right]^{1/p}$. Other terms like $d((x_{1:t}, y_{1:t}), (x'_{1:t}, y'_{1:t}))$ and ${\mathcal P}({\mathcal X}_{1:t} \times {\mathcal Y}_{1:t})$ are understood in the same way.

For a Borel probability measure $\mu$ on temporal data, we denote $\mu(dx_{t+1:T}|x_{1:t})$ as a regular conditional probability kernel of $x_{t+1:T}$ given $x_{1:t}$. When there is no confusion, we also write the kernel as $\mu^t$ for simplicity. If $\mu$ is a finite discrete measure, we use $\mu(x_{t+1:T}|x_{1:t})$, which omits the differential, as the conditional probability at the atom $x_{t+1:T}$ given $x_{1:t}$.

\subsection{Bi-causal optimal transport}\label{sec:bic}

There is a vast literature on the OT theory and its applications \citep{villani2009optimal,kuhn2019tutorial,blanchet2021stat}. However, the previous works have been focusing on data without the time dimension. Motivated by the ubiquitous role of temporal data in OR/MS, finance, and economics, a new notion of adapted Wasserstein distance or causal optimal transport (COT) has been proposed to compare distributions on temporal data \citep{lassalle2013causal,backhoff2017causal}. We present a brief review in this subsection as a self-contained article.

Consider two probability measures $\mu \in {\mathcal P}({\mathcal X})$ and $\nu \in {\mathcal P}({\mathcal Y})$. Recall $\Pi(\mu, \nu)$ is the set of all the couplings that admit $\mu$ and $\nu$ as marginals. Suppose transporting one unit of mass from $x$ to $y$ incurs a cost of $c(x,y)$. A generic OT problem is formulated as
\begin{equation*}
	{\mathcal W}(\mu, \nu) := \inf_{ \pi \in \Pi(\mu, \nu)} \int_{{\mathcal X} \times {\mathcal Y}} c(x, y) \pi(dx, dy),
\end{equation*}
with the Wasserstein distance of order $p$ in \eqref{eq:Wp} as a special case.

If the data have a temporal structure as $x = (x_1, ..., x_t, ..., x_T)$ and $y = (y_1, ..., y_t, ..., y_T)$, not all couplings $\pi \in \Pi(\mu, \nu)$ will make sense. A natural requirement of the transport plan $\pi(x, y)$ should be the non-anticipative condition. Informally speaking, if the past of $x$ is given, then the past of $y$ should be independent of the future of $x$ under the measure $\pi$. Mathematically, it means a transport plan $\pi$ should satisfy
\begin{equation}\label{eq:causal}
	\pi (dy_t | dx_{1:T}) = \pi (dy_t | dx_{1:t}), \quad  t = 1, ..., T-1.
\end{equation}
The property \eqref{eq:causal} is known as the causality condition and the transport plan satisfying \eqref{eq:causal} is called {\it causal} by \cite{lassalle2013causal}. If the same condition holds when we exchange the positions of $x$ and $y$, then the transport plan is called bi-causal. Denote $\Pi_{bc}(\mu, \nu)$ as the set of all bi-causal transport plans between $\mu$ and $\nu$. The bi-causal OT problem considers the optimization over $\Pi_{bc}(\mu, \nu)$ only:
\begin{equation}\label{eq:bc}
	{\mathcal W}_{bc}(\mu, \nu) := \inf_{ \pi \in \Pi_{bc}(\mu, \nu)} \int_{{\mathcal X} \times {\mathcal Y}} c(x, y) \pi(dx, dy).
\end{equation}
 For applications of (bi-)causal optimal transport, see \cite{backhoff2020adapted,xu2020cot,acciaio2021cournot,pflug2012distance,acciaio2020causal} for an incomplete list.

 In this paper, we consider two probability measures $\mu$ and $\nu$ with symmetric positions and focus on the bi-causal transport plans only. For a given transport plan $\pi \in \Pi(\mu, \nu)$, we can decompose $\pi$ in terms of successive regular kernels:
\begin{align}
	\pi(dx_{1:T}, dy_{1:T}) = & \bar{\pi} (dx_1, dy_1) \pi(dx_2, dy_2 | x_1, y_1) \pi(dx_3, dy_3 | x_{1:2}, y_{1:2}) ... \nonumber \\
	& \pi(dx_T, dy_T | x_{1:T-1}, y_{1:T-1}). \label{eq:decomp}
\end{align}
By \citet[Proposition 5.1]{backhoff2017causal}, $\pi$ is a bi-causal transport plan if and only if
\begin{itemize}
	\item[(1)] $\bar{\pi} \in \Pi(p^1_* \mu, p^1_* \nu)$, and
	\item[(2)] for each $t = 1, ... , T-1$ and $\pi$-almost every path $(x_{1:t}, y_{1:t})$, the following condition holds:
	\begin{equation*}
		\pi(dx_{t+1}, dy_{t+1}| x_{1:t}, y_{1:t}) \in \Pi( \mu(dx_{t+1}| x_{1:t}), \nu(dy_{t+1}|y_{1:t})).
	\end{equation*}
\end{itemize}
$p^1_* \mu$ (resp. $p^1_* \nu$) is the pushforward of $\mu$ (resp. $\nu$) by the projection $p^1$ onto the first coordinate.

\section{Equilibrium transport}\label{sec:ti}
\subsection{Motivation}
In a classic discrete-time stochastic control problem, the dynamic programming principle (henceforth DPP) is an essential and widely-used tool to break a multi-period optimization problem into a sequence of subproblems that can be solved recursively \citep{bertsekas1978stoch}. The classic problem satisfies the Bellman equation and is time consistent: If a solution is optimal on the time interval $\{t, ..., T-1\}$, then it is also optimal on any subinterval $\{ s, ..., T-1\}$ with $s \geq t$.

Behavioral economics and psychology have found empirical evidence that human preferences may contradict the axiom of rationality and be time-inconsistent \citep{kahneman1979prospect,laibson1997golden}. When the future objective functional is different and DPP fails, an optimal solution on the time interval $\{t, ..., T-1\}$ may be no longer optimal on a future subinterval $\{ s, ..., T-1\}$ with $s > t$. This self-contradiction makes the notion of ``optimality" unclear. Examples of time inconsistency are mentioned in the introduction.

Inspired by the vast literature on time inconsistency, we extend the bi-causal OT problem \eqref{eq:bc} to embrace time-inconsistent objectives. Before we present the general framework, we discuss an illustrative example from regularization as a motivation. In practice, continuous densities are approximated by empirical measures with finite supports. The computational burden of discrete OT can be high, then regularization is adopted in implementation \citep{cuturi2013sinkhorn}. The discrete bi-causal OT with regularization is
\begin{equation}\label{eq:div}
	\inf_{\pi \in \Pi_{bc}(\mu, \nu)} \ang{c(x_{1:T}, y_{1:T}), \pi(x_{1:T}, y_{1:T})}_F + \ang{ f \left( \frac{\pi(x_{1:T}, y_{1:T})}{\mu(x_{1:T}) \otimes \nu(y_{1:T})} \right), \mu(x_{1:T}) \otimes \nu(y_{1:T})}_F
\end{equation}
for a convex function $f$. The last term is known as the $f$-divergence. Here, $\ang{A, B}_F := \sum_{i, j} a_{ij} b_{ij}$ is the Frobenius inner product between two matrices $A= (a_{ij})$ and $B= (b_{ij})$. $\mu \otimes \nu$ is the independent coupling and $\frac{\pi}{\mu \otimes \nu}$ is interpreted as an element-wise division. The sum in the Frobenius inner product is over all path pairs $(x_{1:T}, y_{1:T})$.

A popular choice of $f(\cdot)$ is the Kullback--Leibler (KL) divergence given by $f(x) = x \ln(x)$. In this case, the regularized objective \eqref{eq:div} is still time-consistent. Indeed, the KL divergence is separable in the sense that
\begin{align}
	& \text{KL} (\pi(x_{1:T}, y_{1:T}) \| \mu(x_{1:T}) \otimes \nu(y_{1:T}) ) \nonumber \\

	& \qquad = \ang{ \ln \left( \frac{\pi(x_{1:T}, y_{1:T})}{\mu(x_{1:T}) \otimes \nu(y_{1:T})} \right), \pi(x_{1:T}, y_{1:T})}_F \nonumber  \\
	& \qquad = \Big\langle \Big\langle \ln \left( \frac{\pi(x_T, y_T | x_{1:T-1}, y_{1:T-1})}{\mu(x_T | x_{1:T-1}) \otimes \nu(y_T | y_{1:T-1})} \right), \pi(x_T, y_T | x_{1:T-1}, y_{1:T-1}) \Big\rangle_F \label{eq:div-decom} \\
	& \qquad \qquad + \ln \left( \frac{\pi(x_{1:T-1}, y_{1:T-1})}{\mu(x_{1:T-1}) \otimes \nu(y_{1:T-1})} \right), \pi(x_{1:T-1}, y_{1:T-1})  \Big \rangle_F.  \nonumber
\end{align}
For simplicity, we denote the successive regular kernels of $\pi(x_{1:T}, y_{1:T})$ as $\{ \bar{\pi}, \pi_1, ..., \pi_{T-1} \}$. For example, $\pi_{T-1} = \pi(x_T, y_T | x_{1:T-1}, y_{1:T-1})$. The second equality \eqref{eq:div-decom} shows that the optimization over $\pi_{T-1}$ is only determined by the first term. We can repeat the decomposition on the second term recursively over time. Then the classic DPP is applicable in this setting.  If $\{ \pi^*_t, ..., \pi^*_{T-1} \}$ is optimal on the interval $\{t, ..., T-1\}$, then $\{ \pi^*_s, ..., \pi^*_{T-1} \}$ is also optimal on the subinterval  $\{s, ..., T-1\}$ with $s > t$.

However, when we consider other choices for $f(\cdot)$ as in \cite{bayraktar2022stability}, the objective is no longer separable and time inconsistency arises. Some examples are
\begin{itemize}
	\item[(1)] Squared Hellinger distance: $f(x) = (1 - \sqrt{x})^2$. The $f$-divergence term becomes
	\begin{align*}
		& \ang{f \left( \frac{\pi(x_{1:T}, y_{1:T})}{\mu(x_{1:T}) \otimes \nu(y_{1:T})} \right),  \mu(x_{1:T}) \otimes \nu(y_{1:T})}_F \\
		& \qquad = 2 - 2 \ang{\sqrt{\pi(x_{1:T}, y_{1:T})}, \sqrt{\mu(x_{1:T}) \otimes \nu(y_{1:T})}}_F.
	\end{align*}
	\item[(2)] Le Cam distance: $f(x) = \frac{(x - 1)^2}{2x + 2}$ and
	\begin{align*}
		& \ang{ f \left( \frac{\pi(x_{1:T}, y_{1:T})}{\mu(x_{1:T}) \otimes \nu(y_{1:T})} \right), \mu(x_{1:T}) \otimes \nu(y_{1:T})}_F \\
		& \qquad = \frac{1}{2} \sum_{x_{1:T}, y_{1:T}} \frac{[\pi(x_{1:T}, y_{1:T}) - \mu(x_{1:T}) \otimes \nu(y_{1:T})]^2}{\pi(x_{1:T}, y_{1:T}) + \mu(x_{1:T}) \otimes \nu(y_{1:T})}.
	\end{align*}
	The summation is over all possible paths $x_{1:T}$ and $y_{1:T}$.
	\item[(3)] Jensen-Shannon divergence: $f(x) = x \ln (\frac{2x}{x+1}) + \ln(\frac{2}{x+1})$, and
	\begin{equation*}
		\text{JS}(\mathbb{P}, \mathbb{Q}) = \text{KL} \left(\mathbb{P} \Big\| \frac{\mathbb{P} + \mathbb{Q}}{2} \right) + \text{KL} \left(\mathbb{Q} \Big\| \frac{\mathbb{P} + \mathbb{Q}}{2} \right),
	\end{equation*}
	with $\mathbb{P} = \pi$, $\mathbb{Q} = \mu \otimes \nu$ is the independent coupling.
\end{itemize}
Under these specifications, the objective is no longer separable and the DPP is violated. The cost functional at time $t$ is given by
\begin{align}
	 J(x_{1:t}, y_{1:t}; \pi) =& \ang{c(x_{1:T}, y_{1:T}), \pi(x_{t+1:T}, y_{t+1:T}| x_{1:t}, y_{1:t})}_F \label{eq:div-obj} \\
	& + \ang{ f \left( \frac{\pi(x_{t+1:T}, y_{t+1:T}| x_{1:t}, y_{1:t})}{\mu(x_{t+1:T} | x_{1:t}) \otimes \nu(y_{t+1:T} |y_{1:t})} \right), \mu(x_{t+1:T}|x_{1:t}) \otimes \nu(y_{t+1:T}|y_{1:t})}_F. \nonumber
\end{align}
When the $f$-divergence term is not separable and the iterated expectation relationship (or tower property) fails, a global optimal transport plan on $\{t, ..., T-1\}$ may not be optimal for the subproblems on $\{t+1, ..., T-1\}$. The intuition is that we may sacrifice the optimality on the subinterval $\{t+1, ..., T-1\}$ to attain better solutions on $\{t, ..., T-1\}$.

\subsection{Definition of equilibrium transport}
To handle the time inconsistency and define a reasonable concept of ``optimal" solutions, we should consider consistent plans that the agent can follow, instead of the pre-commitment solutions. The agent who is aware of time inconsistency should keep in mind that the objectives in the future are different. Inspired by the concept of subgame perfect Nash equilibrium, several works \citep{bjork14FS,bjork17FS} propose to reformulate the problem as a non-cooperative game between agents $t+1, ..., T-1$ that are incarnations of the agent at time $t$. An equilibrium solution, denoted as $\{ \pi^*_t, ..., \pi^*_{T-1} \}$ in a discrete-time setting, should satisfy the following property: If ourselves in the future time $t+1, ... , T-1$ stick with $\{ \pi^*_{t+1}, ..., \pi^*_{T-1} \}$, then it is optimal for us at the current time $t$ to adopt $\pi^*_t$. This backward recursive definition guarantees that the agent at each time $t$ will not deviate from $\pi^*_t$.

Inspired by the above discussion, we define the equilibrium transport as follows. We refer to a transport plan $\pi$ by its successive regular kernels $\{ \bar{\pi}, \pi_1, ..., \pi_{T-1} \}$ in \eqref{eq:decomp}. At time $t \in \{0, ..., T-1\}$, denote a generic cost functional as $J(x_{1:t}, y_{1:t}; \pi)$, which implies that the current path is $(x_{1:t}, y_{1:t})$ and the agent uses the transport plan $\pi$ in time $t, ..., T-1$.
\begin{definition}\label{def:equitrans}
	Consider a given bi-causal transport plan $\pi^*$ with successive regular kernels $\{ \bar{\pi}^*, \pi^*_1, ..., \pi^*_{T-1} \}$.
	\begin{itemize}
		\item[1.] For any $t = 1, ..., T-1$ and a path $(x_{1:t}, y_{1:t})$, consider any transport plan $\gamma_t \in \Pi( \mu(dx_{t+1}| x_{1:t}), \nu(dy_{t+1}|y_{1:t}) )$. When $t=0$, interpret $(x_{1:0}, y_{1:0})$ as none and consider $\bar{\gamma} \in \Pi(p^1_* \mu, p^1_* \nu)$.
		\item[2.] Define a perturbed transport plan $\pi^{t, \gamma}$ as
		\begin{align*}
			\pi^{t, \gamma} (dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}) =& \gamma_t( dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t}) \pi^*(dx_{t+2}, dy_{t+2} | x_{1:t+1}, y_{1:t+1}) ... \\
			& \pi^*(dx_T, dy_T | x_{1:T-1}, y_{1:T-1}).
		\end{align*}
	\end{itemize}
	Then $\pi^*$ is a subgame perfect Nash equilibrium bi-causal transport plan if for every $t = 0, ..., T-1$ and $(x_{1:t}, y_{1:t})$, we have
	\begin{equation*}
		\inf_{\gamma_t \in \Pi(\mu(dx_{t+1}|x_{1:t}), \nu(dy_{t+1} | y_{1:t}) )} J(x_{1:t}, y_{1:t}; \pi^{t, \gamma}) = J(x_{1:t}, y_{1:t}; \pi^{t, *}),
	\end{equation*}
	where
	\begin{align*}
		\pi^{t, *} (dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}) =& \pi^*( dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t}) \pi^*(dx_{t+2}, dy_{t+2} | x_{1:t+1}, y_{1:t+1}) ... \\
		& \pi^*(dx_T, dy_T | x_{1:T-1}, y_{1:T-1}).
	\end{align*}
	We call $\pi^*$ an equilibrium transport for simplicity. If $\pi^*$ exists, define the equilibrium value function $V$ as
	\begin{equation*}
		V_t (x_{1:t}, y_{1:t}) = J(x_{1:t}, y_{1:t}; \pi^{t, *}).
	\end{equation*}
\end{definition}

To understand Definition \ref{def:equitrans} with an example, we revisit the regularized objective \eqref{eq:div-obj} and characterize an equilibrium transport as follows. If time moves to $t+1$, the cost functional is given by
\begin{align*}
	 & J(x_{1:t+1}, y_{1:t+1}; \pi) = \ang{c(x_{1:T}, y_{1:T}), \pi(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1})}_F \\
	 & \qquad  + \ang{ f \left( \frac{\pi(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1})}{\mu(x_{t+2:T} | x_{1:t+1}) \otimes \nu(y_{t+2:T} |y_{1:t+1})} \right), \mu(x_{t+2:T}|x_{1:t+1}) \otimes \nu(y_{t+2:T}|y_{1:t+1})}_F.
\end{align*}
If all agents at time $t+1, ..., T-1$ adopt the equilibrium transport $\pi^*$, then
\begin{equation*}
	J(x_{1:t+1}, y_{1:t+1}; \pi^*) = V_{t+1} (x_{1:t+1}, y_{1:t+1}).
\end{equation*}
Next, we argue backwardly and consider the subproblem at time $t$ with a perturbed transport plan $\pi^{t, \gamma} = \pi^*(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})$. We have
\begin{align*}
	& J(x_{1:t}, y_{1:t}; \pi^{t, \gamma})  = \ang{c(x_{1:T}, y_{1:T}), \pi^*(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})}_F \\
	& \qquad + \ang{ f \left( \frac{\pi^*(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})}{\mu(x_{t+1:T} | x_{1:t}) \otimes \nu(y_{t+1:T} |y_{1:t})} \right), \mu(x_{t+1:T}|x_{1:t}) \otimes \nu(y_{t+1:T}|y_{1:t})}_F.
\end{align*}

Therefore, combining the two equations, we obtain the following recursive relationship between $J(x_{1:t}, y_{1:t}; \pi^{t, \gamma})$ and $J(x_{1:t+1}, y_{1:t+1}; \pi^*)$:
\begin{align*}
	&J (x_{1:t}, y_{1:t}; \pi^{t, \gamma}) = \ang{J(x_{1:t+1}, y_{1:t+1}; \pi^*), \gamma_t(x_{t+1}, y_{t+1})}_F \\
	& \quad -  \ang{ f \left( \frac{\pi^*(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1})}{\mu(x_{t+2:T} | x_{1:t+1}) \otimes \nu(y_{t+2:T} |y_{1:t+1})} \right), \mu(x_{t+2:T}|x_{1:t+1}) \otimes \nu(y_{t+2:T}|y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})}_F   \\
	& \quad + \ang{ f \left( \frac{\pi^*(x_{t+2:T}, y_{t+2:T}| x_{1:t+1}, y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})}{\mu(x_{t+1:T} | x_{1:t}) \otimes \nu(y_{t+1:T} |y_{1:t})} \right), \mu(x_{t+1:T}|x_{1:t}) \otimes \nu(y_{t+1:T}|y_{1:t})}_F.
\end{align*}

By Definition \ref{def:equitrans}, if all agents at time $t+1, ..., T-1$ adopt the equilibrium transport $\{\pi^*_{t+1}, ..., \pi^*_{T-1} \}$, then it is also optimal for the agent at time $t$ to use $\pi^*_t$. Therefore, the recursive relationship should lead to the following equation for the value function:
\begin{align*}
	& V_t (x_{1:t}, y_{1:t}) = \inf_{\gamma_t \in \Pi(\mu(x_{t+1}|x_{1:t}), \nu(y_{t+1} | y_{1:t}) )} \Big[ \ang{V_{t+1} (x_{1:t+1}, y_{1:t+1}), \gamma_t(x_{t+1}, y_{t+1})}_F \\
	&\quad  - \ang{f \left( \frac{\pi^*(x_{t+2:T}, y_{t+2:T} | x_{1:t+1}, y_{1:t+1})}{\mu(x_{t+2:T} | x_{1:t+1}) \otimes \nu(y_{t+2:T}| y_{1:t+1})} \right), \mu(x_{t+2:T} | x_{1:t+1}) \otimes \nu(y_{t+2:T}| y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})}_F \\
	&\quad  + \ang{f \left( \frac{\pi^*(x_{t+2:T}, y_{t+2:T} | x_{1:t+1}, y_{1:t+1}) \gamma_t(x_{t+1}, y_{t+1})}{\mu(x_{t+1:T} | x_{1:t}) \otimes \nu(y_{t+1:T}| y_{1:t})} \right), \mu(x_{t+1:T} | x_{1:t}) \otimes \nu(y_{t+1:T}| y_{1:t})}_F  \Big].
\end{align*}
Since the problem is finite and discrete, the infimum over $\gamma_t$ is attained and leads to the conditional kernel $\pi^*_t$ of an equilibrium transport. Moreover, there can be multiple equilibrium transport plans. The value function is continuous automatically under the discrete topology. We follow \cite{bjork14FS,bjork17FS} and call the recursive relationship above of the value function $V_t$ as the extended dynamic programming equation.

As a side remark, if we consider continuous measures $\mu$ and/or $\nu$ instead, the weak topology equipped with Wasserstein distance is not strong enough. Convergence of probabilities under Wasserstein distance does not necessarily lead to convergence under $f$-divergence, including Jensen-Shannon, KL, reverse KL, total variation, and other divergences. See \citet[Example 1]{arjovsky2017wasserstein}. Since our general extended DPP requires the continuity of value function for the recursion, then we can only show the discrete case for the regularized objective \eqref{eq:div-obj} under the current weak topology.

In the next two subsections, we investigate equilibrium transport with nonlinear objectives and state dependence, inspired by \cite{laibson1997golden,basak10,bjork14FS,bjork17FS,bjork14MF}.

\subsection{Equilibrium transport with nonlinear objectives}
Suppose the agent would like to minimize the transport cost given by the following nonlinear objective:
\begin{equation*}
	\inf_{\pi \in \Pi_{bc}(\mu, \nu)} \int c(x_{1:T}, y_{1:T}) \pi(dx_{1:T}, dy_{1:T}) + G \left(\int h(x_{1:T}, y_{1:T}) \pi(dx_{1:T}, dy_{1:T}) \right).
\end{equation*}
An  illustrative example is
\begin{equation*}
	\inf_{\pi \in \Pi_{bc}(\mu, \nu)} \int c(x_{1:T}, y_{1:T}) \pi(dx_{1:T}, dy_{1:T}) + \left(\int h(x_{1:T}, y_{1:T}) \pi(dx_{1:T}, dy_{1:T}) - m \right)^2.
\end{equation*}
The agent wants to minimize the cost represented by $c(\cdot, \cdot)$ and match the expected value of $h(\cdot, \cdot)$ to a given level of $m$. Time inconsistency appears since the tower property fails.

With a slight abuse of notation, for a fixed $t = 0, ... , T-1$, introduce the objective functional as
\begin{align*}
	J(x_{1:t}, y_{1:t}; \pi) :=& \int c(x_{1:T}, y_{1:T}) \pi(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}) \nonumber \\
	& + G \left(\int h(x_{1:T}, y_{1:T}) \pi(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t})  \right)
\end{align*}
with bi-causal $\pi(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t})$. When $t=0$, interpret $(x_{1:0}, y_{1:0})$ as none.

We impose the following assumptions on the objective.
\begin{assumption}\label{a:nlin_obj}
	\begin{itemize}
		\item[(1)] $c(x_{1:T}, y_{1:T}) \in C_p({\mathcal X}_{1:T} \times {\mathcal Y}_{1:T})$.
		\item[(2)] $h(\cdot, \cdot)$ is continuous and
		$$|h(x_{1:T}, y_{1:T})| \leq C(1 + d((x_{1:T}, y_{1:T}), (\bar{x}_{1:T}, \bar{y}_{1:T}))^{1/r})$$
		for some constant $r>0$ satisfying $1/r \leq p$.
		\item[(3)] $G(\cdot)$ is continuous and $|G(h)| \leq C [1 + |h|^{pr}]$.
		\item[(4)] For each given $t = 0, 1, ..., T-1$ and an arbitrary path $(x_{1:t}, y_{1:t}) \in {\mathcal X}_{1:t} \times {\mathcal Y}_{1:t}$, $G(\int_{{\mathcal X}_{t+1:T} \times {\mathcal Y}_{t+1:T}} h(x_{1:T}, y_{1:T}) \pi(d x_{t+1:T}, d y_{t+1:T} |x_{1:t}, y_{1:t}))$ is strictly quasi-convex in any
		\begin{equation*}
			\pi(d x_{t+1:T}, d y_{t+1:T} |x_{1:t}, y_{1:t}) \in \Pi_{bc} (\mu(d x_{t+1:T} |x_{1:t}), \nu(d y_{t+1:T}|y_{1:t})).
		\end{equation*}
	\end{itemize}
\end{assumption}

The assumptions on probability measures are as follows.
\begin{assumption}\label{a:measure}
	\begin{itemize}
		\item[(1)] $x_{1:t} \mapsto \mu(dx_{t+1:T}|x_{1:t})$  is continuous with respect to metric $d_{{\mathcal X}}$ on the domain ${\mathcal X}_{1:t}$ and the $W_p$ metric on the range ${\mathcal P}_p ({\mathcal X}_{1:t})$.
		\item[(2)] For constants $ a = p $ and $a = 1/r$,
		\begin{equation*}
			\int_{{\mathcal X}_{t+1:T}} d_{\mathcal X}( x_{1:T}, \bar{x}_{1:T})^a \mu(dx_{t+1:T}|x_{1:t}) \leq C[1 +  d_{\mathcal X}(x_{1:t}, \bar{x}_{1:t})^a].
		\end{equation*}
	\end{itemize}
	Same assumptions on $\nu$ hold with the metric $d_{\mathcal Y}$.
\end{assumption}

One of the main results of this paper is to characterize the equilibrium transport with the following extended dynamic programming equation of $V_t$:
\begin{align}
	V_t (x_{1:t}, y_{1:t}) =   \inf_{\gamma_t \in \Pi(\mu(dx_{t+1}|x_{1:t}), \nu(dy_{t+1} | y_{1:t}) )} & \Big[ \int V_{t+1} (x_{1:t+1}, y_{1:t+1}) \gamma_t(dx_{t+1}, dy_{t+1}) \label{eq:nlin_dpp} \\
	& - \int G(g_{t+1} (x_{1:t+1}, y_{1:t+1})) \gamma_t(dx_{t+1}, dy_{t+1}) \nonumber \\
	& + G\Big( \int g_{t+1} (x_{1:t+1}, y_{1:t+1}) \gamma_t(dx_{t+1}, dy_{t+1}) \Big) \Big]. \nonumber
\end{align}
In this equation,
\begin{itemize}
	\item[(a)] the boundary condition for $V$ is
	\begin{equation*}
		V_T (x_{1:T}, y_{1:T}) = c(x_{1:T}, y_{1:T}) + G \left( h(x_{1:T}, y_{1:T}) \right);
	\end{equation*}
	\item[(b)] we denote the optimizer for the equation \eqref{eq:nlin_dpp} at time $t$ as $\pi^*(dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t})$. Then we concatenate these kernels to obtain
	\begin{align}
		\pi^*(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}) = & \pi^*(dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t}) ... \nonumber \\
		& \pi^*(dx_T, dy_T | x_{1:T-1}, y_{1:T-1}); \label{eq:pi*_t}
	\end{align}
	\item[(c)] the function sequence $g_{t+1}$ in \eqref{eq:nlin_dpp} is given by
	\begin{align*}
		g_{t+1} (x_{1:t+1}, y_{1:t+1}) := \int h(x_{1:T}, y_{1:T}) \pi^*(dx_{t+2:T}, dy_{t+2:T} | x_{1:t+1}, y_{1:t+1}),
	\end{align*}
	with $\pi^*(dx_{t+2:T}, dy_{t+2:T} | x_{1:t+1}, y_{1:t+1})$ obtained backwardly using \eqref{eq:pi*_t}.
\end{itemize}

In Theorem \ref{thm:nlin_V}, we apply a version of Berge's maximum theorem with strict quasi-concavity, see \citet[Theorem 9.14 and Corollary 9.20]{sundaram1996first}, to prove the existence and uniqueness of the equilibrium transport. The proof details are given in Appendix \ref{sec:proof_nlin}.
\begin{theorem}\label{thm:nlin_V}
	Suppose Assumptions \ref{a:nlin_obj} and \ref{a:measure} hold. Then the optimizer for \eqref{eq:nlin_dpp}, denoted as $\pi^*(dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t})$, is continuous and unique. Moreover, $\pi^*$ defined recursively by \eqref{eq:pi*_t}  is the unique equilibrium bi-causal transport plan in the sense of Definition \ref{def:equitrans}. Additionally, $V_t$ defined in Definition \ref{def:equitrans} satisfies \eqref{eq:nlin_dpp} and $V_t \in C_p({\mathcal X}_{1:t} \times {\mathcal Y}_{1:t})$.
\end{theorem}
\begin{remark}
	 For the special case that $\mu$ and $\nu$ are discrete measures with finite supports, the extended DPP is straightforward. When ${\mathcal X}$ and ${\mathcal Y}$ are equipped with discrete topology, the value function and the correspondence are continuous automatically. The strict quasi-convexity is no longer needed. There exists an equilibrium bi-causal transport plan $\pi^*$ which may not be unique. The value function $V$ satisfies the extended dynamic programming equation
	\begin{align*}
		V_t (x_{1:t}, y_{1:t}) = \inf_{\gamma_t \in \Pi(\mu(x_{t+1}|x_{1:t}), \nu(y_{t+1} | y_{1:t}) )} \Big[ & \ang{V_{t+1} (x_{1:t+1}, y_{1:t+1}), \gamma_t(x_{t+1}, y_{t+1})}_F  \\
		& - \ang{G(g_{t+1} (x_{1:t+1}, y_{1:t+1})), \gamma_t(x_{t+1}, y_{t+1})}_F \\
		& + G\Big( \ang{ g_{t+1} (x_{1:t+1}, y_{1:t+1}), \gamma_t(x_{t+1}, y_{t+1})}_F \Big) \Big].
	\end{align*}
	The boundary condition and $g_{t+1}$ are interpreted accordingly as in Theorem \ref{thm:nlin_V}.
\end{remark}

\subsection{Equilibrium transport with state-dependent objectives}

In this subsection, we consider time inconsistency from state dependence. For a fixed $t = 0, ... , T-1$, consider an objective functional
 \begin{align}\label{eq:state_obj}
 	J(x_{1:t}, y_{1:t}; \pi) =& \int c(x_t, y_t, x_{1:T}, y_{1:T}) \pi(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}).
 \end{align}
With a slight abuse of notation, we still use $c$ while it has four arguments now. We highlight that the first two $x_t$ and $y_t$ play a different role than those in $x_{1:T}$ and $y_{1:T}$. The first two terms are state-dependent: that is, for a different time $s \neq t$, the integrand becomes $c(x_s, y_s, x_{1:T}, y_{1:T})$. Therefore, the integrand is different when the time changes and it leads to time inconsistency.

A classic example is non-exponential discounting \citep{laibson1997golden}. At time $t$, consider
\begin{equation}\label{eq:nonexp}
	\inf_{\pi \in \Pi_{bc}(\mu(dx_{t+1:T}|x_{1:t}), \nu(dy_{t+1:T}|y_{1:t}))} \int  \sum^T_{s = t+1} \varphi(s - t) c(x_s, y_s) \pi(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}),
\end{equation}
where $\varphi(\cdot)$ is a discounting function besides the power function. Note that the form \eqref{eq:nonexp} can be reformulated into \eqref{eq:state_obj} by redefining a new $\tilde x_t = (t, x_t)$. Since transitions of time are deterministic, the reformulation does not change the conditional probability kernels.

We impose the following assumption on $c(x_t, y_t, x_{1:T}, y_{1:T})$.
\begin{assumption}\label{a:state_obj}
	\begin{itemize}
		\item[(1)] For each $t = 0, 1, ..., T$, $c(x_t, y_t, x_{1:T}, y_{1:T}) \in C_p({\mathcal X}_{1:T} \times {\mathcal Y}_{1:T})$.
		\item[(2)] For each given $t = 0, 1, ..., T-1$ and an arbitrary path $(x_{1:t}, y_{1:t}) \in {\mathcal X}_{1:t} \times {\mathcal Y}_{1:t}$, $\int_{{\mathcal X}_{t+1:T} \times {\mathcal Y}_{t+1:T}} c(x_t, y_t, x_{1:T}, y_{1:T}) \pi(d x_{t+1:T}, d y_{t+1:T} |x_{1:t}, y_{1:t})$ is strictly quasi-convex in any
		\begin{equation*}
			\pi(d x_{t+1:T}, d y_{t+1:T} |x_{1:t}, y_{1:t}) \in \Pi_{bc} (\mu(d x_{t+1:T} |x_{1:t}), \nu(d y_{t+1:T}|y_{1:t})).
		\end{equation*}
	\end{itemize}
\end{assumption}

Similar to the nonlinear objective case, we introduce the extended dynamic programming equation as
\begin{align}
	 V_t (x_{1:t}, y_{1:t}) =  \inf_{\gamma_t \in \Pi(\mu(dx_{t+1}|x_{1:t}),  \nu(dy_{t+1} | y_{1:t}) )} \Big[& \int V_{t+1} (x_{1:t+1}, y_{1:t+1}) \gamma_t(dx_{t+1}, dy_{t+1}) \label{eq:state_dpp} \\
	 & - \int b_{t+1} (x_{t+1}, y_{t+1}, x_{1:t+1}, y_{1:t+1}) \gamma_t(dx_{t+1}, dy_{t+1})   \nonumber  \\
	& + \int b_{t+1} (x_t, y_t, x_{1:t+1}, y_{1:t+1}) \gamma_t(dx_{t+1}, dy_{t+1}) \Big]. \nonumber
\end{align}
In this equation,
\begin{itemize}
	\item[(a)] the boundary condition for $V$ is
	\begin{equation*}
		V_T (x_{1:T}, y_{1:T}) = c(x_T, y_T, x_{1:T}, y_{1:T});
	\end{equation*}
	\item[(b)] with a slight abuse of notation, we still denote the optimizer for the equation \eqref{eq:state_dpp} at time $t$ as $\pi^*(dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t})$ and concatenate these kernels to obtain
	\begin{align}
	\pi^*(dx_{t+1:T}, dy_{t+1:T} | x_{1:t}, y_{1:t}) = & \pi^*(dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t}) ... \nonumber \\
	& \pi^*(dx_T, dy_T | x_{1:T-1}, y_{1:T-1}); \label{eq:state_pi*}
	\end{align}
	\item[(c)] the function sequence $b_{t+1}$ is given by
	\begin{align*}
		b_{t+1} (w, v, x_{1:t+1}, y_{1:t+1}) := \int c(w, v, x_{1:T}, y_{1:T}) \pi^*(dx_{t+2:T}, dy_{t+2:T} | x_{1:t+1}, y_{1:t+1}),
	\end{align*}
	with $\pi^*(dx_{t+2:T}, dy_{t+2:T} | x_{1:t+1}, y_{1:t+1})$ obtained backwardly using \eqref{eq:state_pi*}.
\end{itemize}

Similar to Theorem \ref{thm:nlin_V}, we have the following verification theorem for the state-dependent case, with the proof given in the Appendix \ref{sec:proof_state}.
\begin{theorem}\label{thm:state_V}
	Suppose Assumptions \ref{a:measure} and \ref{a:state_obj} hold. Then the optimizer for \eqref{eq:state_dpp}, denoted as $\pi^*(dx_{t+1}, dy_{t+1} | x_{1:t}, y_{1:t})$, is continuous and unique. Moreover, $\pi^*$ defined recursively by \eqref{eq:state_pi*}  is the unique equilibrium bi-causal transport plan in the sense of Definition \ref{def:equitrans}. Additionally, $V_t$ defined in Definition \ref{def:equitrans} satisfies \eqref{eq:state_dpp} and $V_t \in C_p({\mathcal X}_{1:t} \times {\mathcal Y}_{1:t})$.
\end{theorem}
\begin{remark}
	For discrete and finite $\mu$ and $\nu$, strict quasi-convexity is not needed. The value function $V$ satisfies the extended dynamic programming equation
	\begin{align*}
		V_t (x_{1:t}, y_{1:t}) = \inf_{\gamma_t \in \Pi(\mu(x_{t+1}|x_{1:t}), \nu(y_{t+1} | y_{1:t}) )} \Big[ & \ang{V_{t+1} (x_{1:t+1}, y_{1:t+1}), \gamma_t(x_{t+1}, y_{t+1})}_F \\
		& - \ang{b_{t+1} (x_{t+1}, y_{t+1}, x_{1:t+1}, y_{1:t+1}), \gamma_t(x_{t+1}, y_{t+1})}_F \\
		& + \ang{b_{t+1} ( x_t, y_t, x_{1:t+1}, y_{1:t+1}), \gamma_t(x_{t+1}, y_{t+1})}_F \Big].
	\end{align*}
	And there may be multiple equilibrium transport plans.
\end{remark}

\section{A job market model with inertia}\label{sec:model}
Similarly to the classic OT theory, equilibrium transport has many potential interdisciplinary applications. In this work, we examine a job market from an equilibrium transport perspective.

Suppose $\mu$ represents the distribution of certain characteristics $x$ of employers over time. If we consider a $T$-year time horizon, then the time series of characteristics are $(x_1, ..., x_t, ..., x_T) \sim \mu$. It is problem-specific to determine what characteristics to use for employers. If employers prioritize profitability, we propose to rank the net sales of employers. Sales data tend to be more stable than stock market values. If employers are universities, then rankings are readily available. For employees, commonly used characteristics include talent and salary, while talent is not observable. In this paper, we use the ranks of salaries as the characteristics of employees and denote $\nu$ as the distribution of salary ranks.

An approximate relationship between firm sizes and executive wages, known as Roberts' law in \cite{gabaix2008has}, states that CEO compensation is roughly proportional to $(\text{firm size})^{\kappa}$, with a typical empirical exponent $\kappa \simeq 1/3$. A similar trend also exists between faculty wages and university quality (reputation). Nevertheless, empirical data from executive and academic job markets do not always show a perfect employer-employee correspondence. For instance, CEOs may continue to work in the same firm, while they may be underpaid or overpaid compared with peers in the same industry. One well-known example is Warren Buffett, whose average annual compensation from 2012 to 2021 was only \$395,450. Another example is that Nvidia has both larger market values and annual net sales than its competitor AMD.  However, the average annual compensation of Nvidia's CEO, Jen-Hsun Huang, is significantly lower than that of AMD's CEO, Lisa Su, who earned a total of \$55.763 million more than Huang over a ten-year period.

To model the job market, we propose using the bi-causal OT. In the matching process, employers and employees are in symmetric positions when choosing each other, and no information from the future can be used. Thus, bi-causality becomes a natural requirement on the transport plan. First, the Roberts' law motivates a benchmark model called {\it positive assortative matching} (PAM). Suppose the numbers of employers and employees are the same. Then the optimal employer-employee matching should assign the employer ranked number $n$ with the employee ranked number $n$.  Denote $x_t$ as the rank of employers and $y_t$ as the rank of employees at time $t$. The ranks $x_t$ and $y_t$ have the same range since the numbers of employers and employees are equal. A planner aims to minimize the discrepancy between ranks of employers and employees over all bi-causal transport (matching) plan $\pi$:
 \begin{align}\label{eq:PAM}
	\inf_{\pi \in \Pi_{bc}(\mu, \nu)} \int  \sum^T_{t=1} \delta^t |x_t - y_t| \pi(dx_{1:T}, dy_{1:T}),
\end{align}
 where $\delta \in [0, 1]$ is a discount factor.

As mentioned above, some high-ranked firms pay less to executives and some industries are less consistent with the Roberts' law. We conjecture that there are incentives to maintain the original matching, regardless of its efficiency. For example, there may be implicit benefits to continuing to work for the same company, particularly if an employee is the founder. Additionally, employers may want to preserve the previous ranking of wages among their peers. We attempt to quantify this conjecture and refine the benchmark model \eqref{eq:PAM} using the equilibrium transport method. For any time $t > 1$, we introduce a new term to reduce the cost if we keep similar matching and modify the objective \eqref{eq:PAM} to
 \begin{align}\label{eq:incon_match}
	 	\inf_{\pi \in \Pi_{bc}(\mu^t, \nu^t)} \int \Big[ - \alpha e^{- \frac{|x_{t+1} - x_t| + |y_{t+1} - y_t|}{\tau}} + \sum^T_{s=t+1} \delta^s |x_s - y_s| \Big] \pi(dx_{t+1:T}, d y_{t+1:T} |x_{1:t}, y_{1:t}).
	 \end{align}
$\Pi_{bc}(\mu^t, \nu^t)$ uses simplified notations $\mu^t = \mu(dx_{t+1:T}|x_{1:t})$ and $\nu^t = \nu(dy_{t+1:T}|y_{1:t})$. $\alpha$ is a constant to be determined and $\tau > 0$ is a scaling constant. We call the first term the state-dependent preference function. If we interpret $0/0 = 0$ and set $\tau = 0$, the state-dependent function also includes the indicator function $\mathbf{1}_{\{ x_{t+1} = x_t, y_{t+1} = y_t \}}$ as a special case.  The motivation for the state-dependent term can be further explained as follows:
\begin{enumerate}
	\item Federal securities laws require disclosure of compensation paid to high-ranking executives of public companies. These companies are typically required to provide a Summary Compensation Table with executive compensation for the past {\it three} fiscal years{\footnote{\url{https://www.sec.gov/answers/execcomp.htm}}} in their annual financial statement files (usually annual proxy statements). This requirement is intended to help stakeholders better understand year-over-year changes in compensation. As such, the state-dependent term reflects the concerns and practices of stakeholders and regulators.
	\item Commonly, the current firm-wage matching depends on previous ones for financial reasons. For example, executive compensation programs in public companies often have several components, such as an annual base salary, an annual cash incentive, and long-term equity awards. The base salary tends to remain stable over time, while cash incentives and long-term equity awards are often linked to financial performance, such as annual sales, revenues, and market values. These factors can exhibit momentum and depend on past performance. Our state-dependent term captures this potential dependence on the previous firm-wage matching.
	\item From a psychological and behavioral perspective, this specification also aligns with the {\it status quo bias} if $\alpha > 0$. People prefer to keep their current state. The employer may maintain a similar compensation rank among peers in the same industry. The employee may not resign and find a new job even if he/she is underpaid, possibly due to the status quo bias and relocation costs.
\end{enumerate}
Overall, the state-dependent term with $\alpha > 0$ is a form of persistence or inertia to favor matching pairs close to the previous pair $(x_t, y_t)$. In contrast, the PAM favors pairs with $x_{t+1} = y_{t+1}$ regardless of previous states.

The problem \eqref{eq:incon_match} becomes state-dependent and thus is time-inconsistent. Denoting the equilibrium transport plan for a given $\alpha$ as $\pi(\alpha)$, we can calculate the classic Wasserstein distance ${\mathcal W}(\pi(\alpha), \pi_r )$ between $\pi(\alpha)$ and the real transport plan $\pi_r$:
\begin{equation}\label{eq:wass_alpha}
	{\mathcal W}(\pi(\alpha), \pi_r ) := \inf_{\gamma \in \Pi(\pi(\alpha), \pi_r)} \int (|x_t - x'_t| + |y_t - y'_t| ) d \gamma.
\end{equation}
Then varying $\alpha$, we can find $\alpha$ minimizing ${\mathcal W}(\pi(\alpha), \pi_r )$. When the optimal $\alpha$ is negative, it indicates that the data prefer to deviate from the previous matching, while a positive value of $\alpha$ indicates that the data tend to maintain the previous matching.

In the next two sections, we investigate the inertia, measured by state dependence, of two job markets: top-ranking executives and academia (faculty and postdocs). The choice of these markets is mainly due to the availability of data. We examine the inertia measure $\alpha$ and link it with the efficiency of job markets.

\section{Executive job market}\label{sec:exec}

\subsection{Overview}
We focus on a five-year time horizon 2017 -- 2021. This choice is motivated by the following factors. A longer time horizon incurs a heavy computational burden. Moreover, in our final dataset, the median tenure of a CEO is 4 years and the average tenure is 5.97 years. In the literature with different datasets, \citet[Table 1]{taylor2013ceo} reported the average CEO tenure as 7.9 years and the median as 6 years. Overall, a five-year time horizon is close to the CEO's tenure. Besides, the job market can change significantly over a longer period.

To obtain and clean the data, we adopt the following steps:
\begin{enumerate}
	\item We download the net sales data from Compustat and the executive compensation information from Execucomp. There are no restrictions on the net sales or market values of firms at this stage. The key variables are net sales (``sale" in Compustat) and total compensation (``tdc1" in Execucomp). In 2017 -- 2021, we have data on 1998 firms from Execucomp and 11911 firms from Compustat. Execucomp primarily includes salary information for S\&P 1500 components.
	\item We remove firms with missing net sales data for the five-year period. Compensation data are usually available for the CEO, chief financial officer, and the three other most highly compensated executive officers. We use their mean salary as the representative wage paid to top executives by the firm. After merging the sale and wage datasets, there are 1590 firms remaining.
	\item To further filter the firms, we restrict our analysis to companies with investment-grade credit ratings. We download the S\&P Domestic Long-Term Issuer Credit Rating (variable name: ``splticrm") and remove firms with a rating of CCC+ or lower. We also remove firms without ratings.
\end{enumerate}

Our final sample consists of 790 firms, with a total market value of 61.7 trillion U.S. dollars in 2021. 396 firms are S\&P 500 components. Table \ref{tab:data} presents the descriptive statistics for the data. Firms are classified into different industries using the Global Industry Classification Standard (GICS).
 \begin{table}[H]
	\small
	\centering
	\begin{tabular}{cccccc}
		\hline
		Industry & GICS Code & Number of firms & Market value &  Net sales & Compensation \\
		\hline
		Energy & 10  & 48   & 22693.93 & 19713.07  & 4.82  \\
		\hline
		Materials & 15 & 66  &	12889.25 &	 7301.68	& 3.79 \\
		\hline
		Industrials & 20  & 118 & 20781.59	& 12368.6	& 4.14  \\
		\hline
		Consumer   &  \multirow{2}{*}{25}  & \multirow{2}{*}{112} &	\multirow{2}{*}{32048.23}  & \multirow{2}{*}{16761.5}	    &  \multirow{2}{*}{5.80}	   \\
		Discretionary &    &     &    &      & \\
		\hline
		Consumer & \multirow{2}{*}{30}  & \multirow{2}{*}{48} &	\multirow{2}{*}{46617.8}	& \multirow{2}{*}{34468.41}	&  \multirow{2}{*}{5.55} \\
		Staples &       &      &        &             &       \\
		\hline
		Health Care & 35  & 63 & 51570.67 &	33801.81 & 6.07  \\
		\hline
		Financials & 40 & 110 &	24242.18 &	16041.93 & 5.38 \\
		\hline
		Information  & \multirow{2}{*}{45} & \multirow{2}{*}{79} &	\multirow{2}{*}{75539.2}	& \multirow{2}{*}{17750.21}	 &  \multirow{2}{*}{7.25} \\
		Technology & & & & & \\
		\hline
		Communication & \multirow{2}{*}{50} & \multirow{2}{*}{30} &	\multirow{2}{*}{78888.99}	& \multirow{2}{*}{26807.29}	& \multirow{2}{*}{9.02} \\
		Services & & & & & \\
		\hline
		Utilities & 55 & 50 & 18192.87 & 7446.87 & 3.21 \\
		\hline
		Real Estate &  60 & 66 & 11529.68 & 1949.91 & 3.88 \\
		\hline
	\end{tabular}
	\caption{Summary statistics of firms over 2017 -- 2021. Values in the last three columns are measured in millions of USD.}\label{tab:data}
\end{table}

The first question we address is which industries are more aligned with the Roberts' law, such that larger firms pay higher salaries to executives. To answer this question, we calculate the Spearman's \citep{spearman1904proof} and Kendall's \citep{kendall1938new} rank correlations between net sales and compensations in Table \ref{tab:corr}. We find that the Consumer Discretionary (GICS Code 25), Real Estate (GICS Code 60), and Information Technology (GICS Code 45) industries are the top three sectors with weak correlations between net sales and wages. The Utilities industry (GICS Code 55) has the highest correlation.

Next, we examine the potential connection between the sale-wage correlation and inertia. In the following sections, we will calibrate the optimal value of $\alpha$ and find that the inertia is stronger when the relationship between sales and wages is weaker.
\begin{table}
	\scriptsize
	\centering
	\begin{tabular}{cccccccccccc}
		\hline
		Sector & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 & 50 & 55 & 60 \\
		\hline
		Spearman & 0.72 & 0.793 & 0.752 & 0.515 & 0.755 & 0.768 & 0.745 & 0.646 & 0.75 &
		0.871 & 0.559 \\
		\hline
		Kendall & 0.543 & 0.598 & 0.564 & 0.362 & 0.571 & 0.587 & 0.579 & 0.475 & 0.564 & 0.692 & 0.388 \\
		\hline
		\end{tabular}
\caption{Spearman and Kendall rank correlation between net sales and compensations in five years.}\label{tab:corr}
\end{table}

For clarity, we define an efficient job market as follows.
\begin{definition}
	A job market is said to be more efficient if the sale-wage correlation is higher. We call this correlation as the efficiency coefficient.
\end{definition}

\subsection{Number of clusters}\label{sec:num_group}
Ranks are discrete data that are suitable for discrete OT. However, the smallest number of firms in a single industry, as shown in Table \ref{tab:data}, is 30. If we use the original ranks as the variables $x_t$ and $y_t$, a large amount of data would be needed to estimate the transition probability matrices $\mu(dx_{t+1} | x_t)$ and $\nu(dy_{t+1} | y_t)$. This is unrealistic, as the earliest wage data in Execucomp is primarily from 1992. In addition, transitions of original ranks can be noisy and sensitive to fluctuations in the data. There is also a practical concern about the computational burden of a large transport plan matrix.

Therefore, we need to aggregate data into several clusters with orders. There are two questions to consider in this process: the number of clusters and whether to use even-sized or uneven-sized clusters. As for the first question, the number of clusters cannot be too small or too large. A small number of clusters fails to capture the variations in the data. If we use only one cluster, all industries would perfectly match firms with wages. On the other hand, using too many clusters leads to the problems discussed above. There is no theoretical result on the optimal number of clusters. \citet[Definition 1.2]{backhoff2020estimating} suggests using $N^{1/(T+1)}$ clusters, where $N$ is the number of time series. However, this result is only valid for $N \rightarrow \infty$, and would be too small when there are only a few data points. Given these difficulties, we suggest a practical and data-driven rule for determining the number of clusters.

For simplicity, we will use the same number of clusters for each industry to ensure fairness when comparing cross-sectional differences. The main feature we want to preserve after clustering is the relationship between sales and wages. Clusters are ordered and sales/wages in the same cluster are assigned with the same rank. Denote the sale and wage cluster ranks for the firm $i$ at time $t$ as $x^i_t$ and $y^i_t$, respectively. We call the mean absolute value of the difference in wage and sale ranks, that is,
$$\frac{1}{NT}\sum^N_{i=1}\sum^T_{t=1}|x^i_t - y^i_t|,$$
as the {\it sale-wage discrepancy} for this industry. If the number of clusters is chosen properly, the sale-wage discrepancy should be bigger for sectors with lower sale-wage efficiency reported in Table \ref{tab:corr}. In other words, if we calculate the correlation between the sale-wage discrepancy and the sale-wage efficiency, we should prefer the number of clusters that generates strong negative correlations.

Figure \ref{fig:cluster} shows the correlations between the sale-wage discrepancy and efficiency using different numbers of clusters. The clusters are uneven in size in Figure \ref{fig:uneven} and even in size in Figure \ref{fig:even}. The uneven clusters are determined using the Jenks optimization method \citep{jenks} and the original values of wages and log values of sales. One advantage of uneven clusters is that the correlation is roughly monotonic with the number of clusters, as shown in Figure \ref{fig:uneven}. However, a disadvantage is that a large number of clusters is required to achieve a low, negative correlation, which is impractical due to the small amount of data. Even-sized clusters with ranks as the input, shown in Figure \ref{fig:even}, can achieve a low correlation with a relatively small number of clusters. Figure \ref{fig:cluster} reports both Kendall's and Spearman's rank correlations. We choose the first number of clusters that has a low enough correlation. Based on Figure \ref{fig:even}, the suitable number of clusters is 5 or 6, depending on the type of the correlation used. As a robustness check, we also report the results using seven clusters.

\begin{figure}
	\centering
	\begin{minipage}{0.48\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{uneven_split.pdf}
		\subcaption{Uneven-sized clusters}\label{fig:uneven}
	\end{minipage}
	\begin{minipage}{0.48\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{even_split.pdf}
		\subcaption{Even-sized clusters}\label{fig:even}
	\end{minipage}%
	\caption{Selection for the optimal number of clusters.}\label{fig:cluster}
\end{figure}

\subsection{Estimation of transition matrices}
After determining the number of clusters, we estimate the conditional probability kernels which are transition matrices of ranks. We impose the following assumption to further simplify the model specification.
\begin{assumption}
	The conditional kernels $\mu(dx_{t+1}|x_{1:t})$ and $\nu(dy_{t+1}|y_{1:t})$ are Markov and time-homogeneous. That is, $\mu(dx_{t+1}|x_{1:t}) = \mu(dx_{t+1}|x_t)$ and is the same for any $t = 1, ..., T-1$. The same condition holds for $\nu(dy_{t+1}|y_{1:t})$ as well.
\end{assumption}
 To improve the accuracy of these estimates, we use wage-sale data spanning from 2010 to 2021 for the 790 firms considered in Table \ref{tab:data}. We allow the planner to use all available data from 2010 to the final year of 2021 in the estimation process, in order to provide a larger pool of data. The transition matrices are calculated based on the frequencies of observed transitions between ranks over the years.

Figure \ref{fig:exec_trans} gives the estimated transition matrices for sale and wage ranks. Group 1 is the group with the highest wages or sales. To save space, we average the transition matrices across 11 sectors and summarize them in two figures. A common pattern in all sectors is that sale ranks are more stable than wage ranks. Wages are more stable for the highest and the lowest group. For the wage groups in the middle, the probability of a rise or fall in wage ranks is approximately equal.

\begin{figure}
	\centering
	\begin{minipage}{0.45\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{sale_trans.pdf}
		\subcaption{Sales}\label{fig:sale_trans}
	\end{minipage}
	\begin{minipage}{0.45\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{wage_trans.pdf}
		\subcaption{Wages}\label{fig:wage_trans}
	\end{minipage}%
	\caption{Transition matrices of wage and sale ranks. The matrices are averaged across 11 sectors.}\label{fig:exec_trans}
\end{figure}

The transitions between ranks are crucial for our analysis. In the extreme case, $x_{t+1} = x_t$ and $y_{t+1} = y_t$ with probability one, that is, transition matrices are identity matrices. Then the state-dependent term will have no impact on the transport plans, because the marginal constraints have determined the matching already. There is no freedom to alter firm-wage pairs over time when rank transitions are forbidden. However, as long as transitions are not deterministic, we can capture $\alpha$ with enough data.

\subsection{Model validation}\label{sec:valid}
To stabilize the calibration, we adopt two additional methods. First, as our real dataset is of limited size, we use the re-sampling technique (bootstrapping) to generate a larger dataset by drawing with replacement from the real samples. Second, in practice, there are only a few non-zero entries in the one-period transport plans $\pi_t$. However, after multiplying the one-period transport plans as in equation \eqref{eq:decomp} to obtain $\pi(\alpha)$, there are many paths of $(x_{1:T}, y_{1:T})$ with very small probabilities. To speed up the calculation of the Wasserstein distance ${\mathcal W}(\pi(\alpha), \pi_r)$ in \eqref{eq:wass_alpha} and amplify the impact of $\alpha$, we restrict the domain of $\pi(\alpha)$ to a set of paths with an equal size of the bootstrap samples and the highest probabilities.

In our model, several approximations and hyperparameters may cause systematic bias, including the number of clusters, estimation errors in the transition matrices, the choice of the state-dependent preference function, the size of the bootstrap samples, and the restriction of $\pi(\alpha)$ on paths with the highest probabilities. To address concerns about potential bias in our model, we propose to first apply it to synthetic data consisting only of perfectly matched sale-wage pairs. That is, firms ranked number $n$ always pay wages ranked number $n$. In this case, $x_t = y_t$ for all time points $1 \leq t \leq T$. Ideally, we would expect to observe the optimal value of $\alpha$ to be close to zero in this scenario. There are several advantages to this validation procedure. First, it allows us to test whether the framework is functioning correctly in the simplest possible case. Second, it provides a method for selecting appropriate hyperparameters such that the benchmark outcome is close to zero. Finally, we can use the outcome of this analysis to correct any systematic bias in the model.

\subsection{Results}
Consider the number of clusters as 6 and the preference function as the indicator function. We divided the ranks by the number of clusters $n$, such that $x_t$ and $y_t$ are in $[0, 1-1/n]$, with $0$ representing the highest rank. The candidate values for $\alpha$ are discrete and evenly spaced as $\{ -1.5 , -1.44, ..., 0, ..., 1.44, 1.5\}$, with a step size of $0.06$. The discount factor is set to $\delta = 0.9$. We ran the calibration with ten simulation instances. The size of both the synthetic and bootstrap samples is set to 500. We restricted the domain of $\pi(\alpha)$ to 500 paths with the highest probabilities. The first row of Table \ref{tab:6}, named as benchmark $\alpha$, presents the mean values of the optimal $\alpha$ with synthetic perfectly matched data. In most sectors, the optimal benchmark $\alpha$ is close to zero. Moreover, we observe that the model is biased toward the negative values of $\alpha$. Next, we feed the model with bootstrap samples from the real dataset. The row labeled ``raw $\alpha$" in Table \ref{tab:6} reports the optimal $\alpha$ obtained from the resampled real data. Most of these values shift to the positive direction. Sector 50 obtains a smaller value due to a small fluctuation in the Wasserstein distance. The last row in Table \ref{tab:6} calculates the difference between the raw and benchmark $\alpha$. Overall, there is a state-dependent effect or inertia in most industries, regardless of whether the raw $\alpha$ or the adjusted $\alpha$ is used.
\begin{table}[H]
	\scriptsize
	\centering
	\begin{tabular}{cccccccccccc}
		\hline
		Sector & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 & 50 & 55 & 60 \\
		\hline
		Benchmark $\alpha$ & $-0.084$ & $-0.06$  & $-0.06$ & $-0.084$ & $-0.042$ & $-0.114$ & $-0.288$ & $-0.096$ & $-0.258$ & $-0.12$ & $-0.06$ \\
		Raw $\alpha$ &  0.438 & 0.006 &  0.276 &  0.522 &  0.054 & 0.03 & $-0.048$ &  0.66 & $-0.384$ &  0.018 &  0.702 \\
		Adjusted $\alpha$ & 0.522 &  0.066 &  0.336 &  0.606 &  0.096 &  0.144 &  0.24 &  0.756 & $-0.126$ &  0.138 &  0.762 \\
		\hline
	\end{tabular}
	\caption{Mean values of the optimal $\alpha$ in ten simulations. The number of clusters is set as 6. If there are multiple optimal $\alpha$, we choose the one that is closest to zero.}\label{tab:6}
\end{table}

To understand the impact of $\alpha$ on the Wasserstein distance ${\mathcal W}(\pi(\alpha), \pi_r)$ between the model-implied and empirical transport plans, we consider the Materials and the Real Estate sectors as examples. Figures \ref{fig:MTRL} and \ref{fig:RLEST} depict the calibration curves with $\alpha$ as the horizontal axis and ${\mathcal W}(\pi(\alpha), \pi_r)$ in \eqref{eq:wass_alpha} (normalized) as the vertical axis. The Wasserstein distance is divided by ${\mathcal W}(\pi(-1.5), \pi_r)$ such that all curves start from 1 when $\alpha = -1.5$. There are ten curves in each subplot, representing ten simulation instances with independent bootstrap samples. When comparing the results in Figures \ref{fig:MT_zero} and \ref{fig:MT_real} for the Materials sector, we find that the shape of the calibration curves does not change significantly when perfectly matched data are replaced by the resampled real data. This suggests that the Materials sector does not exhibit a significant state-dependent phenomenon or inertia in sale-wage pairs. However, in the Real Estate sector, we find that the inertia is strong. The optimal $\alpha$ shifts significantly toward positive values in Figure \ref{fig:RE_real}, in contrast to the synthetic perfectly matched data in Figure \ref{fig:RE_zero}. Besides, the state-dependent term can reduce the Wasserstein distance on the negative side more effectively, indicating that negative values of $\alpha$ are less plausible.

\begin{figure}[H]
	\centering
	\begin{minipage}{0.45\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{gs15_Mtr_zero.pdf}
		\subcaption{Synthetic data with perfect matching}\label{fig:MT_zero}
	\end{minipage}
	\begin{minipage}{0.45\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{gs15_Mtr_real.pdf}
		\subcaption{Resampled real data}\label{fig:MT_real}
	\end{minipage}%
	\caption{Calibration curves for Materials sector (GICS Code 15). Ten curves in each subplot represent ten independent simulations.}\label{fig:MTRL}
\end{figure}

\begin{figure}[H]
	\centering
	\begin{minipage}{0.45\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{gs60_RealEs_zero.pdf}
		\subcaption{Synthetic data with perfect matching}\label{fig:RE_zero}
	\end{minipage}
	\begin{minipage}{0.45\textwidth}
		\centering
		\includegraphics[width=0.95\textwidth]{gs60_RealEs_real.pdf}
		\subcaption{Resampled real data}\label{fig:RE_real}
	\end{minipage}%
	\caption{Calibration curves for Real Estate sector (GICS Code 60). Ten curves in each subplot represent ten independent simulations.}\label{fig:RLEST}
\end{figure}

Sectors with less efficient job markets tend to exhibit a higher level of inertia. This assertion has been supported by the negative correlations between the optimal $\alpha$ in Table \ref{tab:6} and the sale-wage efficiency coefficient in Table \ref{tab:corr}. Table \ref{tab:exec_test} shows that  the correlation values are generally less than $-0.6$ and statistically significant at the 5\% level. Job markets with stronger inertia are less likely to embrace the PAM and adapt to change. Thus, there is a negative relationship between job market efficiency and inertia within a dynamic framework.

\begin{table}
	\small
	\centering
	\begin{tabular}{ccc}
		\hline
		Correlation & Spearman ($p$-value) & Kendall ($p$-value) \\
		\hline
		Raw $\alpha$ & $-0.645$ $(0.032)$ & $-0.600$ $(0.010)$ \\
		\hline
		Adjusted $\alpha$ & $- 0.773$ $(0.005)$ & $-0.600$ $(0.010)$ \\
		\hline
	\end{tabular}
	\caption{The relation between job market efficiency and the inertia effect. The correlations are between the optimal $\alpha$ in Table \ref{tab:6} and the sale-wage efficiency coefficient in Table \ref{tab:corr}. The number of clusters is 6.}\label{tab:exec_test}
\end{table}

\subsection{Different number of clusters}
As a robustness check, we examine the results with different numbers of clusters. The screening method in Figure \ref{fig:even} suggests that the choice of five clusters is also reasonable. All other specifications are kept the same. The first row of Table \ref{tab:5} shows that the systematic bias moves further toward the negative side as the number of clusters decreases, indicating that the outcomes may be more similar to the one-cluster case. When resampled real data are used, we obtain the same conclusion as in the six-cluster situation. Table \ref{tab:5_test} finds that the correlations are negative and the $p$-values are smaller than 5\%.
\begin{table}
	\scriptsize
	\centering
	\begin{tabular}{cccccccccccc}
		\hline
		Sector & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 & 50 & 55 & 60 \\
		\hline
		Benchmark $\alpha$ & $-0.42$ & $-0.12$ & $-0.318$ & $-0.324$ & $-0.282$ & $-0.366$ & $-0.42$ & $-0.21$ & $-0.252$ & $-0.222$ & $-0.294$ \\
		Raw $\alpha$ & 0.012 & $-0.072$ & $-0.006$ &  0.45 & $-0.132$ & $-0.51$ & $-0.444$ & $-0.066$ & $-0.492$ & $-0.462$ &  0.054 \\
		Adjusted $\alpha$ & 0.432 &  0.048 &  0.312 &  0.774 &  0.15 & $-0.144$ & $-0.024$ &  0.144 & $-0.24$ & $-0.24$ & 0.348 \\
		\hline
	\end{tabular}
	\caption{Mean values of the optimal $\alpha$ using five clusters. Other settings are the same as in the six-cluster counterpart.}\label{tab:5}
\end{table}

\begin{table}
	\small
	\centering
	\begin{tabular}{ccc}
		\hline
		Correlation & Spearman ($p$-value) & Kendall ($p$-value) \\
		\hline
		Raw $\alpha$ & $-0.709$ $(0.015)$ &  $-0.636$ $(0.006)$ \\
		\hline
		Adjusted $\alpha$ & $-0.700$ $(0.016)$ & $-0.636$ $(0.006)$ \\
		\hline
	\end{tabular}
	\caption{Testing the association between job market efficiency and inertia when five clusters are adopted.}\label{tab:5_test}
\end{table}

When seven clusters are used, the systematic bias shifts to the positive side, as shown in the first row of Table \ref{tab:7}. The original and adjusted optimal $\alpha$ tend to be larger than those obtained using five or six clusters. Although the correlations in Table \ref{tab:7_test} are still negative, the power of the test is much lower. One possible explanation is that Figure \ref{fig:even} shows that the choice of seven clusters results in a slightly worse correlation than the five and six clusters options. In addition, estimation errors in the transition matrices can be larger. Unfortunately, using a longer period of empirical data to estimate the transition matrices may not be a viable solution, as the job market may not be stationary over a long time horizon. Our parameter settings were fine-tuned for six clusters, so it is expected that the performance can be worse with a different number of clusters.
\begin{table}
	\scriptsize
	\centering
	\begin{tabular}{cccccccccccc}
		\hline
		Sector & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 & 50 & 55 & 60 \\
		\hline
		Benchmark $\alpha$ & 0.234 & $-0.06$ &  0.33 &  0.324 & $-0.06$ &  0.582 & $-0.096$ &  0.198 & $-0.102$ &  0.066 &  0.078  \\
		Raw $\alpha$ &  0.372 &  $-0.06$ &  0.756 &  0.42 &  0.222 &  0.636 &  0.054 &  0.438 & 0.204 &  0.264 &  0.666 \\
		Adjusted $\alpha$ &0.138 & 0.0  &  0.426 & 0.096 & 0.282 & 0.054 & 0.15 & 0.24 & 0.306 & 0.198 & 0.588 \\
		\hline
	\end{tabular}
	\caption{Mean values of the optimal $\alpha$ using seven clusters. Other settings are the same as in the six-cluster counterpart.}\label{tab:7}
\end{table}

\begin{table}[H]
	\small
	\centering
	\begin{tabular}{ccc}
		\hline
		Correlation & Spearman ($p$-value) & Kendall ($p$-value) \\
		\hline
		Raw $\alpha$ & $-0.345$ $(0.298)$ &  $-0.382$ $(0.121)$ \\
		\hline
		Adjusted $\alpha$ & $-0.245$ $(0.467)$ & $-0.309$ $(0.218)$ \\
		\hline
	\end{tabular}
	\caption{Testing the association between job market efficiency and inertia when seven clusters are adopted.}\label{tab:7_test}
\end{table}

\section{Academic job market}\label{sec:prof}
\subsection{Overview}
In contrast to the executive job market, the academic job market received less attention in the literature. Our analysis fills this gap using the University of California (UC) compensation data from the Government Compensation in California website\footnote{\url{https://publicpay.ca.gov/Reports/RawExport.aspx}} for the period of 2017 -- 2021. The dataset includes total annual wages for various positions. In our analysis, we focus on four specific full-time positions: ``Prof-Ay-B/E/E", ``Assoc Prof-Ay-B/E/E", ``Asst Prof-Ay-B/E/E", and ``Postdoc-Employee". ``AY" is short for Academic Year and ``B/E/E" means Business/Economics/Engineering. In our current dataset, these three areas are reported jointly. The employees under the ``Postdoc-Employee" title are from all departments. These full-time positions may provide a better representation of compensation levels than part-time positions. In the Appendix, Tables \ref{tab:prof} -- \ref{tab:postdoc} present the summary statistics of salaries in 2017 -- 2021 for nine universities in the UC system, except UC San Francisco since it focuses on medical research and does not have faculty in the B/E/E departments. These tables also include the number of employees under each job title for 2021 only.

To quantify the quality of universities, we use the U.S. News rankings. There is a dataset with historical rankings collected by Andrew G. Reiter and available at his website\footnote{\url{https://andyreiter.com/datasets/}}. It is worth noting that faculty salaries are included as a factor in these rankings, albeit with a weight of only 7\% \footnote{ See the article at \url{https://www.usnews.com/education/best-colleges/articles/how-us-news-calculated-the-rankings}}. Therefore, faculty wages and university rankings are not the same automatically. To determine the annual pay levels at each university, we calculate a single value by first ranking all wage observations among a given year and job title, and then selecting the median rank as a measure of compensation levels. This method mitigates the impact of outliers in the data.

Table \ref{tab:uc_corr} presents Spearman's and Kendall's rank correlations between university rankings and wages in 2017 -- 2021. Salaries for B/E/E assistant professors are more consistent with university rankings, compared with executives in most sectors. However, postdoc salaries do not show a strong positive correlation with university rankings. It is the lowest among all job positions and business sectors analyzed in this paper. This pattern motivates us to examine whether the inertia in postdoc salaries is also stronger.
\begin{table}
	\small
	\centering
	\begin{tabular}{ccccc}
		\hline
		Position & Professor & Associate Professor & Assistant Professor & Postdoc \\
		\hline
		Spearman    & 0.78 &  0.747 &  0.868 &  0.403 \\
		Kendall & 0.616 & 0.585 & 0.707 &  0.298 \\
		\hline
	\end{tabular}
	\caption{Correlations between university rankings and wages in 2017--2021.}\label{tab:uc_corr}
\end{table}

\subsection{Number of clusters and estimation of transition matrices}
To determine the number of clusters for our analysis, we consider the U.S. News rankings and domain knowledge, as the small number of universities in our dataset (nine) makes it difficult to use the criterion applied to the executive data. We find that there are usually two universities with closely ranked positions, with the top two universities in the UC system ranking around the 20th place among national universities (and sometimes tied). The other universities tend to rank around the 30th, 40th, 80th, and 90th places. Based on these observations, we set the number of clusters to five, with two universities in each of the top four clusters and one university in the last cluster.

\begin{figure}
	\centering
	\includegraphics[width=0.4\textwidth]{uni_trans.pdf}
	\caption{Transition matrices of university rankings.}\label{fig:uni_trans}
\end{figure}

To estimate the transition matrices of wages and university rankings, we use data from 2013 to 2021, the longest period available for UC employee salaries on the Government Compensation in California website. Figure \ref{fig:uni_trans} shows that university rankings are stable in this time horizon. Most university ranking transitions happen between the second and the third group. For wages, Figure \ref{fig:prof_trans} demonstrates that postdoc wage ranks are relatively unstable over time. The transition matrices become more and more stable when the academic job ranks become higher. The first group with the highest wages is stable among tenure-track and tenured faculty, while it is unstable for postdoc employees.

\begin{figure}[H]
	\centering
	\begin{minipage}{0.4\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{postdoc_trans.pdf}
		\subcaption{Postdoc}\label{fig:postdoc_trans}
	\end{minipage}
	\begin{minipage}{0.4\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{assitprof_trans.pdf}
		\subcaption{Assistant professor}\label{fig:assitprof_trans}
	\end{minipage}
	\begin{minipage}{0.4\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{assocprof_trans.pdf}
		\subcaption{Associate professor}\label{fig:assocprof_trans}
	\end{minipage}
	\begin{minipage}{0.4\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{fullprof_trans.pdf}
		\subcaption{Professor}\label{fig:full_trans}
	\end{minipage}%
	\caption{Transition matrices of wages.}\label{fig:prof_trans}
\end{figure}

\subsection{Results}
Since the number of universities in our dataset is much smaller than the number of firms, we set the size of our bootstrap samples to 200 for each simulation. We use the state-dependent preference function with $\tau = 1$ and $\alpha \in \{ -1.5 , -1.44, ..., 0, ..., 1.44, 1.5\}$ evenly spaced with a step size of $0.06$.

Using the B/E/E Professors data, Figure \ref{fig:prof} plots the Wasserstein distance between the model-implied transport plan $\pi(\alpha)$ and the empirical transport plan $\pi_r$. The Wasserstein distance is divided by ${\mathcal W}(\pi(-1.5), \pi_r)$ such that all curves start from 1 when $\alpha = -1.5$. As a validation test and benchmark, Figure \ref{fig:prof_zero} calibrates $\alpha$ to synthetic data with perfect matching and plots curves from ten simulation instances. It shows that the optimal $\alpha$ is zero and there is no evidence of inertia when synthetic data with perfectly matching are used. In contrast, Figure \ref{fig:prof_real} shows the results using bootstrap samples. The optimal $\alpha$ is still zero across different simulation runs, indicating that the real data for the professor category are very similar to the perfectly matched synthetic data. Table \ref{tab:uc_alph} reports the mean values of optimal $\alpha$ across these ten simulations.

We observe similar curves for associate professors and assistant professors. Their wages are also highly matched with university rankings, with no significant evidence of inertia. We report the mean values of the optimal $\alpha$ in Table \ref{tab:uc_alph} and show the calibration curves in Figures \ref{fig:assoc} -- \ref{fig:asst}, deferred to the Appendix.

\begin{figure}
	\centering
	\begin{minipage}{0.48\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{prof_zero.pdf}
		\subcaption{Synthetic data with perfect matching}\label{fig:prof_zero}
	\end{minipage}
	\begin{minipage}{0.48\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{prof_real.pdf}
		\subcaption{Resampled real data}\label{fig:prof_real}
	\end{minipage}%
	\caption{Calibration curves for B/E/E professors. In each subplot, ten curves represent individual calibrations with ten independent sampled data.}\label{fig:prof}
\end{figure}

In contrast, we observe strong inertia in the salary data for postdocs. Figure \ref{fig:postdoc_zero} demonstrates the effectiveness of our framework by showing that the optimal $\alpha$ is close to zero when using synthetic perfectly matched data. However, when we use the bootstrap samples, the trend in the calibration curves changes significantly. Figure \ref{fig:postdoc_real} shows that the optimal $\alpha$ exceeds 1, indicating strong inertia in the postdoc salary data. This effect is also stronger than in the executive data. The wages of postdocs do not match the rankings of universities. In the past five years from 2017 to 2021, UC Berkeley has paid lower median salaries to postdocs, compared with UCLA and several other universities in the UC system, as shown in Table \ref{tab:postdoc}. However, the difference in postdoc wages is small. Additionally, the change in Wasserstein distance is also small, mainly because the university wage data, with only nine observations, have a simpler structure than the executive data.

\begin{figure}[H]
	\centering
	\begin{minipage}{0.48\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{postdoc_zero.pdf}
		\subcaption{Synthetic data with perfect matching}\label{fig:postdoc_zero}
	\end{minipage}
	\begin{minipage}{0.48\textwidth}
		\centering
		\includegraphics[width=0.9\textwidth]{postdoc_real.pdf}
		\subcaption{Resampled real data}\label{fig:postdoc_real}
	\end{minipage}%
	\caption{Calibration curves for postdocs.}\label{fig:postdoc}
\end{figure}

Figure \ref{fig:postdoc} also shows that the Wasserstein distance curves can be noisy when $\alpha$ varies. To estimate the optimal $\alpha$ robustly, we filter the original curves using the one-dimensional total variation denoising algorithm \citep{condat2013direct}. This algorithm is suitable for signals with piecewise constant behaviors. Table \ref{tab:uc_alph} reports the mean optimal $\alpha$ calculated from the filtered data. There is no significant evidence of inertia in wages for tenured and tenure-track faculty. In contrast, the effect of inertia is very strong in the postdoc sector.
\begin{table}[H]
	\small
	\centering
	\begin{tabular}{ccccc}
		\hline
		Position & Professor & Associate Professor & Assistant Professor & Postdoc \\
		\hline
		Benchmark $\alpha$  & 0.0 & $-0.06$ & $-0.06$ & 0.0 \\
	 	Raw $\alpha$ & 0.0 &  $-0.06$ & $-0.06$ & 1.194 \\
		Adjusted $\alpha$ & 0.0 & 0.0 & 0.0 & 1.194 \\
		\hline
	\end{tabular}
	\caption{Mean values of optimal $\alpha$ after total variation filtering.}\label{tab:uc_alph}
\end{table}
We test the connection between the inertia and university ranking-wage correlations, summarized in Table \ref{tab:uc_finaltest}. Correlations in Table \ref{tab:uc_finaltest} are negative, showing that when the inertia is stronger, wages are less matched with university rankings. However, the power of our tests is low since only four types of jobs are considered.
\begin{table}[H]
	\small
	\centering
	\begin{tabular}{ccc}
		\hline
		Correlation & Spearman ($p$-value) & Kendall ($p$-value) \\
		\hline
		Raw $\alpha$ & $-0.632$ (0.368) &  $-0.548$ (0.279) \\
		\hline
		Adjusted $\alpha$ & $-0.775$ (0.225) & $-0.707$ (0.180) \\
		\hline
	\end{tabular}
	\caption{Correlations between optimal $\alpha$ and the ranking-wage dependence.}\label{tab:uc_finaltest}
\end{table}