restatable[3][#3]

# How Inductive Bias in Machine Learning Aligns with Optimality in Economic
Dynamics

Mahdi Ebrahimi Kahou1 James Yu2 Jesse Perla2 Geoff Pleiss2,3
1Bowdoin College 2University of British Columbia 3Vector Institute
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

This paper examines the alignment of inductive biases in machine learning (ML)
with structural models of economic dynamics. Unlike dynamical systems found in
physical and life sciences, economics models are often specified by
differential equations with a mixture of easy-to-enforce initial conditions
and hard-to-enforce infinite horizon boundary conditions (e.g. transversality
and no-ponzi-scheme conditions). Traditional methods for enforcing these
constraints are computationally expensive and unstable. We investigate
algorithms where those infinite horizon constraints are ignored, simply
training unregularized kernel machines and neural networks to obey the
differential equations. Despite the inherent underspecification of this
approach, our findings reveal that the inductive biases of these ML models
innately enforce the infinite-horizon conditions necessary for the well-
posedness. We theoretically demonstrate that (approximate or exact) min-norm
ML solutions to interpolation problems are sufficient conditions for these
infinite-horizon boundary conditions in a wide class of problems. We then
provide empirical evidence that deep learning and ridgeless kernel methods are
not only theoretically sound with respect to economic assumptions, but may
even dominate classic algorithms in low to medium dimensions. More
importantly, these results give confidence that, despite solving seemingly
ill-posed problems, there are reasons to trust the plethora of black-box ML
algorithms used by economists to solve previously intractable, high-
dimensional dynamical systems—paving the way for future work on estimation of
inverse problems with embedded optimal control problems.

##  1 Introduction

Numerical solutions to dynamical systems are central to many quantitative
fields in economics, including macroeconomics, finance, international trade,
labor, and industrial organization. For example, a model of capital
accumulation and consumption consists of a differential (or difference)
equation enforcing market dynamics, an initial condition defining the starting
capital, and intra-temporal conditions ensuring that supply and demand are
balanced. Statistical models of these economic processes are used as
laboratories for counterfactual scenarios, such as evaluating the impact of
changes in capital gains tax rates on investment. Economists are increasingly
interested in applying ML techniques towards forward and inverse modeling of
these systems, motivated by remarkable successes on dynamical systems in
physical and life sciences [e.g. 10, 13, 39, 27].

#### Challenges of economic dynamical systems.

ML techniques for forward and inverse modeling typically target initial-value
problems (IVP) that are entirely specified by the differential equation and
initial conditions on all variables.111In some cases, mixed with terminal
boundary values at a finite-horizon (e.g., [23], [41], and [35] for theory on
the role of terminal, _finite horizon_ boundary conditions). In contrast,
economic systems often embed infinite-horizon optimal control problems which
include _asymptotic boundary conditions_ (e.g. transversality and no-ponzi-
scheme conditions)—presenting challengers to both classical and ML algorithms.

The asymptotic boundary conditions are necessary to ensure the problem is
_well-posed_ in the sense of Hadamand (i.e., a solution exists, is unique, and
changes smoothly with initial conditions).222The definition of “well-posed”
used here can be traced back to Hadamard’s works [21, 22]. Also, see [43] and
[49] for classic treatments of regularization to solve _ill-posed_ problems,
of which our methods are broadly connected. In particular, without imposing
these asymptotic boundary conditions there will be a continuum of solutions
which fulfill the initial conditions and ODEs. Classical methods to solve
these problems and apply the asymptotic boundary conditions are often
numerically unstable and shooting towards steady states of the system—which
adds numerical instability, and unique challenges when there is _steady-state
multiplicity_ and hysteresis (i.e., requirements to partition the space into
different basins of attraction for each steady state).

Despite these challenges, economists are putting their faith in standard ML
methods for forward and inverse modeling with applications ranging from wealth
inequality [24], financial frictions coupled with inequality [20], the
heterogeneous impacts of climate change [4], portfolio choice problems [2],
heterogeneous agent new Keynesian models [26], human capital accumulation in
the labor market [25], and labor market dynamics in search and matching
[36]—using neural networks [3, 30, 17, 18] or kernel machines [38]. However,
while these methods encode initial values, they _simply ignore the asymptotic
boundary conditions_ and solve the underspecified, _ill-posed_ problems—so
there are no guarantees on which solution the ML model will find.

#### Inductive bias of machine learning.

In this paper, we investigate to what extent the _implicit biases_ of these ML
models will automatically fulfill the asymptotic boundary conditions of
economic systems—despite solving an ill-posed problem. Large neural networks
or ridgeless kernel regressors can learn to satisfy a system of differential
equations on any amount of data/grid points, thus implying that these models
are _overparameterized_. Recent work in optimization suggests that SGD-based
training of these overparameterized models will (approximately) yield the
minimum-norm function capable of interpolating the data [5]. In the case of
dynamical systems modeling, this implies that neural networks and kernel
machines will find the minimum-norm function that (locally) satisfies the
dynamical system and whatever explicit constraints are known. Put differently,
ML models of economic dynamical systems with minimum function norm constraints
may provide sufficient conditions for the asymptotic boundary-value
constraints.

This leads to our key research question: does the minimum norm inductive bias
of machine learning align with asymptotic constraints of economics models?
Outside of general appeals to Occam’s Razor, there is no reason that these two
should align.

#### Contributions.

We provide guidance on when and why we can trust ML solutions for an important
class of dynamical systems used in economics. In particular, we show:

  * •

inductive bias alignment: theoretical and empirical evidence that the minimum
norm implicit bias of large machine learning models aligns with many
asymptotic boundary conditions found in economics problems;

  * •

learning the right set of steady-states: evidence that kernel machines and
neural networks find the optimal steady states of dynamical systems—leading to
very accurate generalization outside of the training data—even without
enforcing asymptotic boundary conditions;

  * •

consistency of ML estimates: guarantees that ML methods converge to the true
minimum-norm solutions (and thus, in many cases, the asymptotic boundary
constrained solutions) in the limit of infinite data; and

  * •

robustness and speed: demonstrations that kernel machines and neural networks
may be competitive in speed and robustness with traditional methods for
modeling economic systems, even on small-scale problems.

Finally, this lays the foundation for future work on solving inverse problems
that will jointly estimate structural parameters and solve dynamic models (in
contrast to approaches like [15], which require solving a dynamical system
through standard techniques as a nested step in forming the likelihood).

##  2 Setup

The economic systems we are concerned with have roots in optimal control
theory, though many will contain additional algebraic constraints not found in
classical optimal control setups. An inherent characteristic of infinite-
horizon optimal control problems, whether deterministic or stochastic, are
asymptotic boundary conditions such as transversality conditions. These
conditions can be formulated sequentially or recursively in a state space, and
in continuous- or discrete-time (see discussions of necessary conditions in
[32, 6, 44]). In this paper, we will focus on an important subset of these
problems—deterministic, continuous-time systems of differential-algebraic
equations (DAE)—which allow for formal proof of results and comparison to
baseline solutions from preexisting algorithms.

#### Problem class.

Within the subset of problems we formally analyze here, classic nonlinear
optimal control decisions are embedded in a large class of economic problems,
resulting in systems of ordinary differential and algebraic equations [1].
Partitioning these into three vectors: _state variables_
${\bm{x}}(t)\in\mathbb{R}^{M_{x}}$ with an initial condition ${\bm{x}}_{0}$;
_jump variables_ ${\bm{y}}(t)\in\mathbb{R}^{M_{y}}$ with accompanying
asymptotic boundary conditions typically from the embedded control problems;
and optional _static variables_ ${\bm{z}}(t)\in\mathbb{R}^{M_{z}}$.333When
required, static variables constrain the solution manifold and are not matched
with boundary or initial values. The connection between the number of _jump
variables_ and stability local to a steady-state is discussed in [7]. Given
this, the dynamical system with primitives
${\bm{F}},{\bm{G}},{\bm{H}},{\bm{B}}$ is

$\displaystyle\dot{{\bm{x}}}(t)$
$\displaystyle={\bm{F}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$ (1)
$\displaystyle\dot{{\bm{y}}}(t)$
$\displaystyle={\bm{G}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$ (2)
$\displaystyle\mathbf{0}$
$\displaystyle={\bm{H}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$ (3) subject to
an initial value ${\bm{x}}(0)={\bm{x}}_{0}$ and boundary conditions
$\displaystyle\mathbf{0}$
$\displaystyle=\lim_{t\to\infty}{\bm{B}}(t,{\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$
(4)

Key necessary conditions to be _well-posed_ are that ${\bm{x}}_{0}$ is given,
there exists a ${\bm{y}}(0)$ uniquely determined using the asymptotic boundary
condition Eq. 4, and—where applicable—there exists a unique ${\bm{z}}(0)$
fulfilling Eq. 3 given ${\bm{y}}(0)$.444We will assume sufficient regularity
conditions to make this possible, such as Lipshitz ${\bm{F}},{\bm{G}}$ and
${\bm{H}}$, and that $\det(\nabla{\bm{H}})\neq 0$.

#### Algorithmic challenges with classic, non-ML methods.

Since this problem looks deceptively standard, it is worth emphasizing the
challenges before explaining how ML algorithms can tackle them. If an initial
condition for ${\bm{y}}(0)$ was given, then Eqs. 1 and 2 is an initial value
problem (IVP) ODE which is well-posed and can be solved with high dimensional
${\bm{x}}(t)$ and ${\bm{y}}(t)$ using a variety of classical methods.555The
optional DAE term with ${\bm{z}}(t)$ and Eq. 3 adds additional, but
surmountable, challenges for IVPs. The central challenge to this class of
problems is that applying Eq. 4 numerically is not directly possible, as it is
asymptotic, and approximations to learn ${\bm{y}}(0)$ are numerically unstable
since the solution path lies on a saddle-path (see [1]). These problems are
compounded when there is steady-state multiplicity (i.e., regions of initial
conditions $x_{0}$ which converge to separated $x_{\infty}$). Central to our
algorithm in Sec. 3 is the insight that the least-explosive solutions to the
under-determined system Eqs. 1, 2 and 3 given ${\bm{x}}(0)={\bm{x}}_{0}$
initial condition will be those that fulfill Eq. 4.

##  3 Method

In this paper we use a machine learning approach for forward modeling of the
dynamical system described in Sec. 2. Specifically, we will train large
machine learning models to learn functions ${\hat{\bm{x}}}$, ${\hat{\bm{y}}}$,
${\hat{\bm{z}}}$ that minimize the following objective function:

$\displaystyle\min_{{\hat{\bm{x}}},{\hat{\bm{y}}},{\hat{\bm{z}}}}\sum_{t_{i}\in\mathcal{D}}$
$\displaystyle\left[\eta_{1}\left\|{\hat{\dot{\bm{x}}}}(t_{i})-{\bm{F}}({\hat{\bm{x}}}(t_{i}),{\hat{\bm{y}}}(t_{i}),{\hat{\bm{z}}}(t_{i}))\right\|_{2}^{2}+\eta_{2}\left\|{\hat{\dot{\bm{y}}}}(t_{i})-{\bm{G}}({\hat{\bm{x}}}(t_{i}),{\hat{\bm{y}}}(t_{i}),{\hat{\bm{z}}}(t_{i}))\right\|_{2}^{2}\right.$
$\displaystyle\left.+\eta_{3}\left\|{\bm{H}}({\hat{\bm{x}}}(t_{i}),{\hat{\bm{y}}}(t_{i}),{\hat{\bm{z}}}(t_{i}))\right\|_{2}^{2}+\eta_{4}\left\|{\hat{\bm{x}}}(0)-{\bm{x}}_{0}\right\|_{2}^{2}\right],$
(5)

where $\mathcal{D}$ is a set of time points from some fixed interval $[0,T]$.
At first glance, there are several potential problems with this approach. Most
critically, Eq. 5 completely ignores the asymptotic boundary condition in Eq.
4. The resulting learning problem is thus ill-posed. Despite this
underspecification, we will show that many overparameterized ML models trained
with gradient descent tend towards economically-optimal solutions—i.e. those
that (approximately) satisfy Eq. 4. These results justify the use of ML models
for solving economic models with embedded optimal control problems and also
hint at scalability and numerical benefits over classic algorithms.

#### Intuition: overparameterized ML models tend towards minimum norm
solutions.

Overparameterized models by definition have more parameters than training data
points, and thus there are many settings of parameters capable of
“interpolating” (perfectly fitting) the training data. Recent work has
suggested that overparameterized models in ML tend towards the interpolating
solution with the smallest norm in the hypothesis class [e.g. 5]. This
property is exactly true for linear models like kernel machines, and many SGD-
trained neural networks are implicitly regularized towards small-norm
solutions [29]. We will see that—with careful design—this minimum-norm
inductive bias can be used as an approximate stand-in for the infinite-horizon
boundary condition in Eq. 4.

#### Goal: minimum Sobolev seminorm solutions.

In particular, we will use overparameterized ML models for Eq. 5 that are
biased towards minimum _Sobolev seminorm_ interpolating solutions. Given a
function $w:[0,\infty)\to\mathbb{R}$ from
$\mathcal{W}^{P+1,2}(\mathbb{R}^{+})$, the Sobolev-$P\\!+\\!1,2$ seminorm of
$w$ is defined as the Sobolev-$P,2$ norm of its weak derivative:
$\left|w\right|_{\mathcal{W}^{P+1,2}(\mathbb{R}^{+})}:=\left\|\dot{w}\right\|_{\mathcal{W}^{P,2}(\mathbb{R}^{+})}.$
The minimal structure required for many economic models is that will have a
single derivative almost everywhere, and hence we can assume solutions to our
problem classes lie in Sobolev spaces. From a theoretical perspective, there
is an established connection between infinite-horizon optimal solutions of
economic growth models and Sobolev spaces [44] and [14]. (In Sec. 4.1 we
explore an economic model where the minimum Sobolev seminorm solution provably
coincides with infinite-horizon boundary conditions.) More intuitively, a
Sobolev seminorm bias is a natural choice to resolve the instability issues of
classic algorithms. Solutions to Eqs. 1, 2 and 3 either diverge, in which case
$\dot{{\bm{x}}}(t)$ or $\dot{{\bm{y}}}(t)$ are large in absolute value, or
they convergence and fulfil Eq. 4. Hence,
$\|\dot{\bm{x}}\|_{L_{2}(\mathbb{R}^{+})}=|{\bm{x}}|_{W^{1,2}(\mathbb{R}^{+})}$
would be large for instable/incorrect solutions and small for stable/correct
solutions.

#### The inductive bias of ridgeless kernel regression.

Our first ML model class that we consider is ridgeless kernel regression. This
model class has the advantage that we can formally characterize the function
norms and underlying inductive bias. Specifically, we represent
${\hat{\bm{x}}},{\hat{\bm{y}}},{\hat{\bm{z}}}$ as

$\begin{gathered}{\hat{\bm{x}}}(t)={\bm{x}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{x}}}}(\tau)d\tau,\qquad{\hat{\bm{y}}}(t)={\hat{\bm{y}}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{y}}}}(\tau)d\tau,\qquad{\hat{\bm{z}}}(t)={\hat{\bm{z}}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{z}}}}(\tau)d\tau,\qquad\\\
\dot{\hat{\bm{x}}}(t)=\operatorname{{\textstyle\sum}}_{j=1}^{N}\bm{\alpha}^{x}_{j}K(t,t_{j}),\qquad\dot{\hat{\bm{y}}}(t)=\operatorname{{\textstyle\sum}}_{j=1}^{N}\bm{\alpha}^{y}_{j}K(t,t_{j}),\qquad\dot{\hat{\bm{z}}}(t)=\operatorname{{\textstyle\sum}}_{j=1}^{N}\bm{\alpha}^{z}_{j}K(t,t_{j}),\end{gathered}$
(6)

where
$\bm{\alpha}^{x}_{j},\bm{\alpha}^{y}_{j},\bm{\alpha}^{z}_{j},{\hat{\bm{y}}}_{0},{\hat{\bm{z}}}_{0}$
are learned parameters and $K(\cdot,\cdot)$ is the Matérn kernel with
smoothness $\nu$ and lengthscale $\ell$. As we will demonstrate shortly, it is
necessary to model
${\hat{\dot{\bm{x}}}},{\hat{\dot{\bm{y}}}},{\hat{\dot{\bm{z}}}}$ rather than
${\hat{\bm{x}}},{\hat{\bm{y}}},{\hat{\bm{z}}}$ to target the Sobolev seminorm
rather than the Sobolev norm.

If we add an infinitesimal
$\sum_{i,j=1}^{N}\bm{\alpha}_{i}^{x\top}\bm{\alpha}^{x}_{j}K(t_{i},t_{j})$
regularization term to the objective function in Eq. 5 (and analogous terms
for $y$ and $z$), then the solution to the optimization problem will be
infinitesimally close to the minimum-norm solution:

$\displaystyle\min_{\begin{subarray}{c}{\hat{\bm{x}}}\in\mathcal{H}^{M_{x}},{\hat{\bm{y}}}\in\mathcal{H}^{M_{y}},{\hat{\bm{z}}}\in\mathcal{H}^{M_{z}}\\\
{\hat{\bm{y}}}_{0}\in\mathbb{R}^{M_{y}},{\hat{\bm{z}}}_{0}\in\mathbb{R}^{M_{z}}\end{subarray}}\,$
$\displaystyle\left(\operatorname{{\textstyle\sum}}_{m=1}^{N_{x}}\,\|\dot{\hat{{x}}}^{(m)}\|^{2}_{\mathcal{H}}+\operatorname{{\textstyle\sum}}_{m=1}^{N_{y}}\,\|\dot{\hat{{y}}}^{(m)}\|^{2}_{\mathcal{H}}+\operatorname{{\textstyle\sum}}_{m=1}^{N_{z}}\,\|\dot{\hat{{z}}}^{(m)}\|^{2}_{\mathcal{H}}\right)$
(7) s.t.
$\displaystyle\dot{{\hat{\bm{x}}}}(t_{i})={\bm{F}}({\hat{\bm{x}}}(t_{i}),{\hat{\bm{y}}}(t_{i}),{\hat{\bm{z}}}(t_{i}),\quad\text{for
all }t_{i}\in\mathcal{D}$
$\displaystyle\dot{{\hat{\bm{y}}}}(t_{i})={\bm{G}}({\hat{\bm{x}}}(t_{i}),{\hat{\bm{y}}}(t_{i}),{\hat{\bm{z}}}(t_{i})),\quad\text{for
all }t_{i}\in\mathcal{D}$
$\displaystyle\mathbf{0}={\bm{H}}({\hat{\bm{x}}}(t_{i}),{\hat{\bm{y}}}(t_{i}),{\hat{\bm{z}}}(t_{i})),\quad\text{for
all }t_{i}\in\mathcal{D}$

where $\|\cdot\|_{\mathcal{H}}$ denotes the RKHS norm associated with the
Matérn kernel. For functions defined over compact domains, it is well
established that the Matérn RKHS with $\nu=(P-1/2)$ and a specific value of
$\ell$ is exactly equal to the $\mathcal{W}^{P,2}([0,T])$ norm for all $P\geq
1$ (see Appx. A). Thus, the Matérn kernel with $\nu=1/2$ (and appropriate
lengthscale) produces the minimum $\mathcal{W}^{2,2}([0,T])$-seminorm
solution. While this does not perfectly align with the
$\mathcal{W}^{1,2}(\mathbb{R}^{+})$ seminorm of interest,666 Note that we
cannot target the Sobolev-$1,2$ seminorm with Eq. 6, as
$|w|_{\mathcal{W}^{1,2}}=\|\dot{w}\|_{L_{2}}$ and $L_{2}$ is not an RKHS.
Moreover, we cannot target Sobolev spaces of $[0,\infty)$ functions as the
Matérn kernel only provably aligns with Sobolev spaces over functions with
Lipschitz boundaries (see Appx. A).  we experimentally demonstrate in Sec. 4
that minimum $\mathcal{W}^{2,2}([0,T])$ solutions closely align with those
that satisfy Eq. 4.

Finally, we prove that ${\hat{\bm{x}}},{\hat{\bm{y}}},{\hat{\bm{z}}}$ are
consistent estimators of the true minimum $\mathcal{W}^{2,2}([0,T])$-seminorm
solution to Eqs. 1, 2 and 3. (See Appx. B for the proof.)
theoremthm:consistency Given some $0<K<\infty$, let ${\mathbb{S}}$ be the set
of functions $({\bm{x}},{\bm{y}},{\bm{z}})$ that satisfy Eqs. 1, 2 and 3 with
${\bm{x}}(0)={\bm{x}}_{0}$ and
$\|{\bm{y}}(0)\|_{\infty},\|{\bm{z}}(0)\|_{\infty}\leq K$. If

  * •

$t\in\mathcal{D}$ are drawn uniformly i.i.d. from $[0,T]$;

  * •

${\bm{F}}$, ${\bm{G}}$, and ${\bm{H}}$ in Eqs. 1, 2 and 3 are Lipschitz with
respect to $\|\cdot\|_{\infty}$; and

  * •

${\bm{F}}$ and ${\bm{G}}$ have Lipschitz first derivatives and ${\bm{H}}$ has
Lipschitz first and second derivatives,

then then the solutions
${\hat{\bm{x}}}_{N},{\hat{\bm{y}}}_{N},{\hat{\bm{z}}}_{N}$ from Eq. 7 with the
Matern-$1/2$ kernel satisfies Eqs. 1, 2 and 3 almost everywhere in the limit
as $N\to\infty$ and

$\displaystyle\textstyle{\lim_{N\to\infty}}$
$\displaystyle\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,|{\hat{{x}}}^{(m)}_{N}|^{2}_{\mathcal{W}^{2,2}([0,T])}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,|{\hat{{y}}}^{(m)}_{N}|^{2}_{\mathcal{W}^{2,2}([0,T])}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,|{\hat{{z}}}^{(m)}_{N}|^{2}_{\mathcal{W}^{2,2}([0,T])}$
$\displaystyle\overset{\mathrm{a.s.}}{=}\>\textstyle{\inf_{({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}}}$
$\displaystyle\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,|{\hat{{x}}}^{(m)}|^{2}_{\mathcal{W}^{2,2}([0,T])}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,|{\hat{{y}}}^{(m)}|^{2}_{\mathcal{W}^{2,2}([0,T])}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,|{\hat{{z}}}^{(m)}|^{2}_{\mathcal{W}^{2,2}([0,T])}.$

#### The inductive bias of neural networks.

Alternatively, we consider solutions to Eq. 5 where
${\hat{\bm{x}}},{\hat{\bm{y}}},{\hat{\bm{z}}}$ are represented by
overparameterized neural networks trained with stochastic gradient descent.
While the inductive bias of neural networks is difficult to formally
characterize, recent work suggests that they are approximately biased towards
minimum Sobolev seminorm solutions. For infinite width neural networks, the
inductive bias is towards the minimum norm solution in the RKHS of the neural
tangent kernel [5]. For finite width neural networks, the implicit bias of SGD
regularizes the Sobolev seminorm of the solution [29]. We empirically
demonstrate in Sec. 4 that neural networks, like ridgeless kernel regression,
recover solutions that satisfy Eq. 4.

##  4 Results

We solve two standard baselines in dynamic economics: a model of risk-neutral
asset pricing (Sec. 4.1) and the neoclassical growth model (Sec. 4.2). These
problems are chosen because they are standard examples in textbooks [e.g. 28],
they admit reference solutions from classical methods, and they have
established results regarding the set of solutions to the ill-posed versions
without the asymptotic boundary conditions.777For asset pricing models, see
[8, 16], which characterizes the set of solutions not fulfilling
transversality and connects them to economic bubbles. [1] covers saddle-path
analysis of the neoclassical growth model as a way to analyze failures of the
transversality condition. In Sec. 4.3, we show a case where multiple steady
states exist—a setting where our methods are especially helpful—and summarize
other experiments in Sec. 4.4.888In our baseline examples, the models have
been simplified to eliminate static variables (i.e., no ${\hat{\bm{z}}}(t)$),
but examples of the broader problem class are discussed in Sec. 4.4.

For both baseline examples, we compare two ML methods trained on
$\mathcal{D}:=\\{0,1,2,\ldots,30\\}$:

#### Neural networks.

We model ${\hat{\bm{x}}}(t)$ and ${\hat{\bm{y}}}(t)$ (and ${\hat{\bm{z}}}(t)$
when applicable) with overparameterized neural networks and minimize Eq. 5
directly. Crucially, there is no explicit regularization and Eq. 3 is not
imposed directly, so this model relies exclusively on the inductive bias of
the optimization process to align with the economically-optimal solution.

#### Matérn kernels.

We solve with the kernel approximation in Eq. 7, which directly minimizes the
function norm in the RKHS. We use the Matérn kernel with $\nu=1/2$ to align
with the $\mathcal{W}^{2,2}([0,T])$ seminorm, which is “close” to the
$\mathcal{W}^{1,2}(\mathbb{R}^{+})$ seminorm of interest.

###  4.1 Asset Pricing

Models of asset pricing and “rational bubbles” are relatively simple and
benefit from closed-form solutions. These models have traditionally served a
pedagogical role in demonstrating the necessity of transversality conditions;
for instance, see [28]. The model values a sequence of dividends where, in
cases of deviation from the fundamentals, the asset price reflects a “bubble.”
Let ${\bm{x}}(t)\in\mathbb{R}$ be the flow payoffs from a claim to an asset,
and ${\bm{y}}(t)\in\mathbb{R}$ be the price of a claim to that asset. For
simplicity, we assume that the flow payoffs, ${\bm{x}}$, follow a
deterministic process. For a given $x_{0}$, the key equations are999Our
baseline parameters are $x_{0}=1.0,c=0.02,g=-0.2$, and $r=0.1$.

$\displaystyle\dot{{\bm{x}}}(t)$
$\displaystyle=c+g{\bm{x}}(t):={\bm{F}}({\bm{x}}(t),{\bm{y}}(t))$ (8)
$\displaystyle\dot{{\bm{y}}}(t)$
$\displaystyle=r{\bm{y}}(t)-{\bm{x}}(t):={\bm{G}}({\bm{x}}(t),{\bm{y}}(t))$
(9) $\displaystyle 0$
$\displaystyle=\lim_{t\rightarrow\infty}e^{-rt}{\bm{y}}(t):={\bm{B}}({\bm{x}}(t),{\bm{y}}(t))$
(10)

where constants $c,g$ govern the process for dividends and $r>0$ is a risk-
neutral investor’s discount rate. Eq. 10 is called a “no-bubble” condition in
this case. The set of solutions to Eqs. 9 and 8, without imposing Eq. 10, can
be found analytically (see [28, 16]):

$\displaystyle{\bm{y}}(t)$
$\displaystyle=\operatorname{{\textstyle\int}}_{0}^{\infty}e^{-r\tau}{\bm{x}}(t+\tau)\mathrm{d}\tau={\bm{y}}_{f}(t)+\zeta
e^{rt},$ (11) $\displaystyle{\bm{y}}_{f}(t)$
$\displaystyle:=\tfrac{c}{r-g}+\left({\bm{x}}_{0}-\tfrac{c}{r-g}\right)e^{(r-g)t}$
(12)

where ${\bm{y}}_{f}(t)$ is interpreted as the “fundamental” price of the
asset, and $\zeta\geq 0$ is indeterminate. However, when the “no-bubble”
condition (i.e., Eq. 10) is imposed, this problem is well-posed, with a unique
solution of ${\bm{y}}(t)={\bm{y}}_{f}(t)$ (i.e., $\zeta=0$).

#### Sufficiency of the min-norm solutions to select the “no-bubble” solution.

A key property of this example is that Eq. 11 characterizes the full set of
solutions to Eqs. 8 and 9 when not imposing Eq. 10. For a given function norm,
apply the triangle inequality to the set of solutions from Eq. 11 to yield
$||{\bm{y}}_{f}||_{\mathcal{H}}\leq||{\bm{y}}||_{\mathcal{H}}\leq||{\bm{y}}_{f}||_{\mathcal{H}}+\zeta||e^{rt}||$.
This shows that among the possible solutions to the ill-posed problem, the
norm $||{\bm{y}}||_{\mathcal{H}}$ is minimized when $\zeta=0$, at which point
${\bm{y}}(t)={\bm{y}}_{f}(t)$ and the “no-bubble” condition in Eq. 10 is
fulfilled.

#### Results.

Fig. 1 shows the results of the kernel machine and neural network
approximations.101010The neural network maps $t$ to ${\bm{y}}(t)$. We use 4
hidden layers with 128 nodes and the Tanh activation function. In the final
layer, the Softplus is used to ensure positive values for prices. For
optimization, we use Adam with a learning rate of $10^{-3}$. Despite not
applying the “no-bubble” condition in Eq. 10, both the kernel interpolator and
the neural network recover the unique equilibrium of economic interest almost
perfectly.

Figure 1: Linear asset pricing model (Eqs. 9 and 8) solved without imposing
the “no bubble” condition (Eq. 10). The left panel shows the kernel machine
and neural network approximations relative to the benchmark. The right panel
displays the absolute errors relative to the benchmark on a test set of 100
points from $t=0$ to $t=40$. In all cases, the ML methods recover the
benchmark and even extrapolate well beyond the training data, despite solving
an indeterminate problem and not directly calculating the steady state.

###  4.2 Neoclassical Growth Model

The neoclassical growth model (also known as the Ramsey–Cass–Koopman model) is
the benchmark problem in macroeconomics to explore optimal control in
dynamical systems (see [1] for its many variations). A planner optimizes paths
of capital, ${\bm{x}}(t)\in\mathbb{R}$, and consumption
${\bm{y}}(t)\in\mathbb{R}$, given an initial condition $x(0)=x_{0}$. Key
parameters are a depreciation rate $\delta>0$, discount rate $r>0$, and a
production function $f(x):=x^{a}$ where $a$ is the capital share of
output.111111Our baseline parameters are $x_{0}=1.0,\delta=0.1,r=0.11,$ and
$a=\frac{1}{3}$.

$\displaystyle\dot{{\bm{x}}}(t)$
$\displaystyle=f({\bm{x}}(t))-{\bm{y}}(t)-\delta{\bm{x}}(t):={\bm{F}}({\bm{x}}(t),{\bm{y}}(t))$
(13) $\displaystyle\dot{{\bm{y}}}(t)$
$\displaystyle={\bm{y}}(t)\left[f^{\prime}({\bm{y}}(t))-\delta-r\right]:={\bm{G}}({\bm{x}}(t),{\bm{y}}(t))$
(14) $\displaystyle 0$
$\displaystyle=\lim_{t\rightarrow\infty}e^{-rt}{\bm{x}}(t)/{\bm{y}}(t):={\bm{B}}({\bm{x}}(t),{\bm{y}}(t))$
(15)

Figure 2: Neoclassical growth model (Eqs. 14 and 13) solved without imposing
the transversality condition (Eq. 15). The top panels show the results for
capital, ${\bm{x}}(t)$, and the bottom panels show the results for
consumption, ${\bm{y}}(t)$. The left panels show the kernel machine and neural
network approximations relative to the benchmark. The right panels display the
absolute errors relative to the benchmark on a test set of 100 points from
$t=0$ to $t=40$. In all cases, the ML methods recover the benchmark and even
extrapolate well beyond the training data, despite solving an indeterminate
problem and not directly calculating the steady state.

#### Results.

Fig. 2 shows the results for both experiments.121212The neural network maps
$t$ to $\begin{bmatrix}{\bm{x}}(t)&{\bm{y}}(t)\end{bmatrix}$. We use 4 hidden
layers with 128 nodes and the Tanh activation function. In the final layer,
the Softplus is used to ensure positive values for consumption and capital.
For optimization, we use Adam with a learning rate of $10^{-3}$. For the
weights in Eq. 5 we use $\eta_{1}=0.4,\eta_{2}=0.4,\eta_{4}=0.2$. Both the
kernel interpolator and the neural network recover the unique equilibrium of
economic interest almost perfectly, fulfilling Eq. 15 despite not being
provided the steady-state as a boundary condition.

#### Robustness.

While LABEL:thm:consistency demonstrates consistency, it is helpful to check
our methods’ sensitivity to hyperparameters and features of $\mathcal{D}$:
Sec. C.1 shows that our methods perform well with a much sparser set of time
points; Sec. C.2 indicates low sensitivity to different kernel
hyperparameters; and Sec. C.3 demonstrates that the approximation remains
effective in the short- to medium-term even if $\mathcal{D}$ does not approach
the steady state.

###  4.3 Neoclassical Growth Model with Multiple Steady-States

We now turn to a more complex version of the neoclassical growth model (Eqs.
14, 13 and 15) where $f(x):=A\max\\{x^{a},b_{1}x^{a}-b_{2}\\}$, as in
[42].131313This function has a kink in $f^{\prime}(x)$ at
$\bar{x}=\left(\frac{b_{2}}{b_{1}-1}\right)^{\frac{1}{a}}$. Our baseline
parameters are $A=0.5,b_{1}=3.0,b_{2}=2.5,k_{0}=1.0,\delta=0.1,r=0.11$. This
problem is still well-posed (as long as Eq. 15 is included) and hence has a
unique trajectory for a given initial condition—which would lead towards one
of the two steady states. This provides a significant challenge for classical
algorithms, such as shooting methods, which rely on finding the steady state
and using it as a finite-horizon boundary condition. All steady states need to
be calculated, and the space of initial conditions must be sharply partitioned
into different domains of attraction to know which steady state to use as the
boundary condition for a particular initial condition.

Figure 3:  Neoclassical growth model (Eqs. 14 and 13) with multiple steady
states due to the the concave-convex production function in Sec. 4.3. The left
panel shows the solution trajectories for capital, ${\bm{x}}(t)$, using Matérn
Kernels for $70$ different initial conditions ${\bm{x}}_{0}\in[0.5,4]$. The
right panel shows the solution trajectories for consumption, ${\bm{y}}(t)$,
for the same initial conditions. In all cases the ML method gives the correct
trajectory for a given initial condition, despite not directly calculating the
steady states. Furthermore, the left panel shows that the method is effective
at learning the basins of attraction for the different steady states.

#### Results.

Fig. 3 shows the results using Matérn kernels for 70 different initial
conditions across the basins of attraction for the two steady states. While
the error relative to the benchmark solution is not shown, it is of comparable
magnitude to that in Fig. 2.

It is worth noting that the kernel regression solutions find the “correct”
steady state; i.e. the steady state that correctly corresponds with the
supplied initial condition ${\bm{x}}_{0}$. Intuitively, this behavior is a
consequence of the minimum Sobolev seminorm inductive bias. Consider the two
possible trajectories ${\bm{x}}_{0}$ to each of the two steady states: the
trajectory with smaller gradients and less steep dynamics will have a smaller
norm. In contrast, classical methods—which do not have a minimum norm
inductive bias—require a mapping between initial conditions and the correct
steady states.

###  4.4 Other Examples

In Sec. C.4 we examine a variation of the DAE formulation of the neoclassical
growth problem of Sec. 4.2; i.e. a dynamical system where the ${\bm{z}}(t)$
static variables are not eliminated. In addition, we test our methods on a
standard model of human capital and economic growth (see Sec. D.1) as well as
a model of optimal advertising used in the marketing literature [40] (see Sec.
D.2).

##  5 Related Work

#### Alternative methods for linearized systems.

Perturbation solutions to models of economic dynamics use classic stability
analysis of linear quadratic (LQ) control to find solutions that fulfill
asymptotic boundary conditions (e.g., [7]). These methods use the linearity of
time-invariant policies in LQ control and enforce stability by spectral
methods (i.e., remove explosive roots that violate Eq. 4 by solving Ricatti
Equations). Our paper takes inspiration from this broader approach: algorithms
which select non-explosive roots will be on the solution manifold and
automatically fulfill Eq. 4.

#### ML for optimal control and economic dynamical systems.

Recent literature on optimal control and robotics explores how to achieve
stable solutions with deep neural networks [33, 34, 12, 11]. More directly
within the economics literature, [18] discusses the intuitive connection
between the inductive bias of deep neural networks and turnpikes [31] in
dynamic economics, but does not provide theory. Our paper contributes to this
literature, providing a formal argument on why the inductive bias of ML
algorithms promotes stability in infinite-horizon control.

##  6 Discussion

While the class of problems in Sec. 2 focused on deterministic, continuous-
time dynamical systems, these challenges exist in all dynamic models with
forward-looking behavior or optimal control. Many questions in this research
agenda remain:

  * •

How to model inequality and complementarity constraints? Can Secs. 2 and 3 be
extended to include inequality constraints (e.g., an additional primitive
${\bm{C}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))\geq\mathbf{0}$)? This extension
would allow us to incorporate a larger class of important models (e.g.,
lifecycle models of consumption, macroeconomic models of financial frictions,
and Stackelberg games).

  * •

What is the role of inductive bias of ML algorithms for state-space
representations and stochastic models? In Sec. 2, we characterized
deterministic, infinite-horizon models written in a sequence space, but many
dynamic economic models are stochastic and more naturally written in a
recursive, state-space representation. Does the inductive bias of ML
algorithms applied to state-space formulations (e.g., approximating policy or
value functions of the state) fulfill Eq. 4?141414Transversality conditions
remain necessary conditions for problems to be well-posed in state-space
representations of dynamic models. Even verifying them is difficult in
practice because they must hold at every point in the state space and are
written as expected discounted values of random trajectories. For more
details, see [28] and [19].

  * •

When are kernel methods and neural networks preferable to classic algorithms?
Beyond its purpose in formalizing the analysis, we have shown that kernel
methods are fast and robust in practice, even in low dimensions. Classic
algorithms such as shooting methods will likely dominate kernels when
carefully tuned and where steady states are easy to solve, but ML methods
might be directly useful in this problem domain.

## References

  * Acemoglu [2008] D. Acemoglu.  _Introduction to modern economic growth_.  Princeton university press, 2008.
  * Azinovic and Žemlička [2023] M. Azinovic and J. Žemlička.  Economics-inspired neural networks with stabilizing homotopies.  _arXiv preprint arXiv:2303.14802_ , 2023.
  * Azinovic et al. [2022] M. Azinovic, L. Gaegauf, and S. Scheidegger.  Deep equilibrium nets.  _International Economic Review_ , n/a(n/a), 2022.  doi: https://doi.org/10.1111/iere.12575.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/iere.12575.
  * Barnett et al. [2023] M. Barnett, W. Brock, L. P. Hansen, R. Hu, and J. Huang.  A deep learning analysis of climate change, innovation, and uncertainty, 2023.
  * Belkin [2021] M. Belkin.  Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation.  _Acta Numerica_ , 30:203–248, 2021.
  * Benveniste and Scheinkman [1982] L. Benveniste and J. Scheinkman.  Duality theory for dynamic optimization models of economics: The continuous time case.  _Journal of Economic Theory_ , 27(1):1–19, 1982\.  ISSN 0022-0531.  doi: https://doi.org/10.1016/0022-0531(82)90012-6.  URL https://www.sciencedirect.com/science/article/pii/0022053182900126.
  * Blanchard and Kahn [1980] O. J. Blanchard and C. M. Kahn.  The solution of linear difference models under rational expectations.  _Econometrica: Journal of the Econometric Society_ , pages 1305–1311, 1980.
  * Blanchard and Watson [1982] O. J. Blanchard and M. W. Watson.  Bubbles, rational expectations and financial markets, 1982.
  * Brezis and Brézis [2011] H. Brezis and H. Brézis.  _Functional analysis, Sobolev spaces and partial differential equations_ , volume 2.  Springer, 2011.
  * Calderhead et al. [2008] B. Calderhead, M. Girolami, and N. Lawrence.  Accelerating Bayesian inference over nonlinear differential equations with Gaussian processes.  _Advances in neural information processing systems_ , 21, 2008.
  * Chang et al. [2019] Y.-C. Chang, N. Roohi, and S. Gao.  Neural lyapunov control.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/2647c1dba23bc0e0f9cdf75339e120d2-Paper.pdf.
  * Chen [2023] G. Chen.  Deep neural network approximations for the stable manifolds of the hamilton-jacobi equations, 2023.
  * Chen et al. [2018] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31, 2018.
  * Chichilnisky [1977] G. Chichilnisky.  Nonlinear functional analysis and optimal economic growth.  _Journal of Mathematical Analysis and Applications_ , 61(2):504–520, 1977.  ISSN 0022-247X.  doi: https://doi.org/10.1016/0022-247X(77)90134-2.  URL https://www.sciencedirect.com/science/article/pii/0022247X77901342.
  * Childers et al. [2022] D. Childers, J. Fernández-Villaverde, J. Perla, C. Rackauckas, and P. Wu.  Differentiable state space models and hamiltonian monte carlo estimation.  _Working Paper_ , 2022.
  * Diba and Grossman [1988] B. T. Diba and H. I. Grossman.  The theory of rational bubbles in stock prices.  _The Economic Journal_ , 98(392):746–754, 1988\.  ISSN 00130133, 14680297.  URL http://www.jstor.org/stable/2233912.
  * Ebrahimi Kahou et al. [2021] M. Ebrahimi Kahou, J. Fernández-Villaverde, J. Perla, and A. Sood.  Exploiting symmetry in high-dimensional dynamic programming.  Working Paper 28981, National Bureau of Economic Research, July 2021.  URL http://www.nber.org/papers/w28981.
  * Ebrahimi Kahou et al. [2024] M. Ebrahimi Kahou, J. Fernández-Villaverde, S. Gomez, J. Perla, and J. Rosa.  Spooky boundaries at a distance: Exploring transversality and stability with deep learning.  _Working Paper_ , 2024.
  * Fernández-Villaverde et al. [2016] J. Fernández-Villaverde, J. F. Rubio-Ramírez, and F. Schorfheide.  Solution and estimation methods for dsge models.  In _Handbook of macroeconomics_ , volume 2, pages 527–724. Elsevier, 2016.
  * Fernández-Villaverde et al. [2023] J. Fernández-Villaverde, S. Hurtado, and G. Nuño.  Financial frictions and the wealth distribution.  _Econometrica_ , 91(3):869–901, 2023.  doi: https://doi.org/10.3982/ECTA18180.  URL https://onlinelibrary.wiley.com/doi/abs/10.3982/ECTA18180.
  * Hadamard [1902] J. Hadamard.  Sur les problèmes aux dérivées partielles et leur signification physique.  _Princeton university bulletin_ , pages 49–52, 1902.
  * Hadamard and Hadamard [1932] J. Hadamard and J. Hadamard.  Le problème de cauchy et les équations aux dérivées partielles linéaires hyperboliques: leçons professées à l’université yale.  1932\.
  * Han et al. [2018] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, Aug. 2018.  ISSN 1091-6490.  doi: 10.1073/pnas.1718942115.  URL http://dx.doi.org/10.1073/pnas.1718942115.
  * Han et al. [2022] J. Han, Y. Yang, and W. E.  Deepham: A global solution method for heterogeneous agent models with aggregate shocks, 2022.
  * Jungerman [2023] W. Jungerman.  Dynamic monopsony and human capital.  Working Paper, 2023.
  * Kase et al. [2022] H. Kase, L. Melosi, and M. Rottner.  _Estimating nonlinear heterogeneous agents models with neural networks_.  Centre for Economic Policy Research, 2022.
  * Kovachki et al. [2023] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Neural operator: Learning maps between function spaces with applications to pdes.  _Journal of Machine Learning Research_ , 24(89):1–97, 2023.
  * Ljungqvist and Sargent [2018] L. Ljungqvist and T. Sargent.  _Recursive Macroeconomic Theory, fourth edition_.  MIT Press, 2018.  ISBN 9780262038669.  URL https://books.google.ca/books?id=IG1qDwAAQBAJ.
  * Ma and Ying [2021] C. Ma and L. Ying.  On linear stability of sgd and input-smoothness of neural networks.  _Advances in Neural Information Processing Systems_ , 34:16805–16817, 2021.
  * Maliar et al. [2021] L. Maliar, S. Maliar, and P. Winant.  Deep learning for solving dynamic economic models.  _Journal of Monetary Economics_ , 122:76–101, 2021.  ISSN 0304-3932.  doi: https://doi.org/10.1016/j.jmoneco.2021.07.004.  URL https://www.sciencedirect.com/science/article/pii/S0304393221000799.
  * McKenzie [1976] L. W. McKenzie.  Turnpike theory.  _Econometrica_ , 44(5):841–865, 1976.  ISSN 00129682, 14680262.  URL http://www.jstor.org/stable/1911532.
  * Michel [1982] P. Michel.  On the transversality condition in infinite horizon optimal problems.  _Econometrica: Journal of the Econometric Society_ , pages 975–985, 1982.
  * Nakamura-Zimmerer et al. [2022a] T. Nakamura-Zimmerer, Q. Gong, and W. Kang.  Neural network optimal feedback control with guaranteed local stability.  _IEEE Open Journal of Control Systems_ , 1:210–222, 2022a.  doi: 10.1109/OJCSYS.2022.3205863.
  * Nakamura-Zimmerer et al. [2022b] T. Nakamura-Zimmerer, Q. Gong, and W. Kang.  Neural network optimal feedback control with enhanced closed loop stability.  In _2022 American Control Conference (ACC)_ , pages 2373–2378. IEEE, 2022b.
  * Pardoux and Tang [1999] E. Pardoux and S. Tang.  Forward-backward stochastic differential equations and quasilinear parabolic pdes.  _Probability theory and related fields_ , 114:123–150, 1999\.
  * Payne et al. [2024] J. Payne, A. Revei, and Y. Yang.  Deep learning for search and matching models, 2024.  URL https://dx.doi.org/10.2139/ssrn.4768566.
  * Rasmussen and Williams [2006] C. E. Rasmussen and C. Williams.  _Gaussian processes for machine learning_ , volume 1.  MIT Press, 2006.
  * Scheidegger and Bilionis [2019] S. Scheidegger and I. Bilionis.  Machine learning for high-dimensional dynamic stochastic economies.  _Journal of Computational Science_ , 33:68–82, 2019.  ISSN 1877-7503.  doi: https://doi.org/10.1016/j.jocs.2019.03.004.  URL https://www.sciencedirect.com/science/article/pii/S1877750318306161.
  * Schober et al. [2014] M. Schober, D. K. Duvenaud, and P. Hennig.  Probabilistic ode solvers with runge-kutta means.  _Advances in neural information processing systems_ , 27, 2014.
  * Sethi [1973] S. P. Sethi.  Optimal control of the vidale-wolfe advertising model.  _Operations research_ , 21(4):998–1013, 1973\.
  * Sirignano and Spiliopoulos [2018] J. Sirignano and K. Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  _Journal of Computational Physics_ , 375:1339–1364, Dec. 2018.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2018.08.029.  URL http://dx.doi.org/10.1016/j.jcp.2018.08.029.
  * Skiba [1978] A. K. Skiba.  Optimal growth with a convex-concave production function.  _Econometrica_ , 46:527–539, 1978.  ISSN 00129682, 14680262.  URL http://www.jstor.org/stable/1914229.
  * Tikhonov [1963] A. N. Tikhonov.  On the solution of ill-posed problems and the method of regularization.  In _Doklady akademii nauk_ , volume 151, pages 501–504. Russian Academy of Sciences, 1963.
  * Van et al. [2007] C. L. Van, R. Boucekkine, and C. Saglam.  Optimal control in infinite horizon problems: a sobolev space approach.  _Economic Theory_ , 32(3):497–509, 2007.
  * Vidale and Wolfe [1957] M. L. Vidale and H. Wolfe.  An operations-research study of sales response to advertising.  _Operations research_ , 5(3):370–381, 1957.
  * Wainwright [2019] M. J. Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge university press, 2019.
  * Weber [2006] T. A. Weber.  An infinite-horizon maximum principle with bounds on the adjoint variable.  _Journal of Economic Dynamics and Control_ , 30(2):229–241, 2006.
  * Wendland [2004] H. Wendland.  Scattered data approximation.  2004\.
  * Willoughby [1979] R. A. Willoughby.  Solutions of ill-posed problems (a. n. tikhonov and v. y. arsenin).  _SIAM Review_ , 21(2):266–267, 1979.  doi: 10.1137/1021044.  URL https://doi.org/10.1137/1021044.

##  Appendix A Connection Between Sobolev-${P,2}$ Spaces and Matérn RKHS

Let $k_{\nu,\ell}(\cdot,\cdot):\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ denote
the Matérn covariance function with smoothness $\nu$ and lengthscale $\ell$.
For the purposes of this paper, we will define $k_{\nu,\ell}$ as:

$k_{\nu,\ell}(t,t^{\prime})=\kappa(|t-t^{\prime}|),\qquad\widehat{\kappa}(\omega):=\left(1+\ell^{2}\omega^{2}\right)^{-\nu-1/2},$
(16)

where $\widehat{(\cdot)}$ corresponds to the Fourier transform. Note that Eq.
16 corresponds to the standard Matérn kernel definition [37] after
appropriately scaling the inputs and outputs.
$\nabla^{(r)}k_{\nu,\ell}(\cdot,\cdot)$ will denote the $r^{\mathrm{th}}$
derivative of $k_{\nu,\ell}$ with respect to its its first argument.

Given some interval $[0,T]\subseteq\mathbb{R}$,
$\mathcal{H}^{\nu,\ell}([0,T])$ denotes the RKHS of $[0,T]\to\mathbb{R}$
functions where the reproducing kernel is equal to $k_{\nu,\ell}$.
$\mathcal{W}^{P,2}([0,T])$ denotes to the Sobolev-$P,2$ space of
$[0,T]\to\mathbb{R}$ functions. Whenever possible we will drop the
superscripts for $\mathcal{H}^{\nu,\ell}([0,T])$.

######  Proposition 1 (Equivalence of $\mathcal{W}^{P,2}$ and
$\mathcal{H}^{\nu,\ell}$).

For any positive integer $P$, any $\ell>0$, and any $\nu=P-1/2$ there exists
positive constants $C_{1}(\nu,\ell)$ and $C_{2}(\nu,\ell)$ so that

$C_{1}(\nu,\ell)\|w\|_{\mathcal{W}^{P,2}([0,T])}\leq\|w\|_{\mathcal{H}^{\nu,\ell}([0,T])}\leq
C_{2}(\nu,\ell)\|w\|^{\mathcal{W}^{P,2}([0,T])}.$

for all $w\in\mathcal{W}^{P,2}$. In other words, the Sobolev norm
$\|\cdot\|_{\mathcal{W}^{P,2}([0,T])}$ and the RKHS norm
$\|\cdot\|_{\mathcal{H}^{\nu,\ell}([0,T])}$ are equivalent, and thus
$\mathcal{W}^{P,2}([0,T])=\mathcal{H}^{\nu,\ell}([0,T])$.

###### Proof.

We begin by first establishing an equivalence between
$\mathcal{W}^{P,2}(\mathbb{R})$ and $\mathcal{H}^{\ell,\nu}(\mathbb{R})$. In
the Fourier domain, the Sobolev norm for $\mathbb{R}\to\mathbb{R}$ functions
is given by

$\|w\|_{\mathcal{W}^{P,2}(\mathbb{R})}=\left\|\widehat{w}(\cdot)\left(1+(\cdot)^{2}\right)^{P/2}\right\|_{L_{2}(\mathbb{R})},$

and, for any RKHS $\mathcal{H}(\mathbb{R})$ with stationary reproducing
kernels, the RKHS norm for $\mathbb{R}\to\mathbb{R}$ functions is given by

$\|w\|_{\mathcal{H}(\mathbb{R})}=\left\|\widehat{w}(\cdot)\left(\widehat{\kappa}(\cdot)\right)^{-1/2}\right\|_{L_{2}(\mathbb{R})}$

The Matérn kernel, as defined in Eq. 16, has a Fourier transform that decays
at at a rate of $(1+|\cdot|^{2})^{-P}$ and thus
$\|\cdot\|_{\mathcal{H}^{\nu,\ell}(\mathbb{R})}$ is bounded above and below by
a constant multiple of $\|\cdot\|_{\mathcal{W}^{P,2}(\mathbb{R})}$ where the
constant only depends on $\ell$. This argument can be generalized to prove an
equivalence between $\mathcal{W}^{P,2}([0,T])$ and
$\mathcal{H}^{\ell,\nu}([0,T])$, as the domain $[0,T]$ trivially has a
Lipschitz boundary. (See [48, Corollary 10.48] for details.) ∎

######  Corollary 0.1 (Equality of $\mathcal{W}^{P,2}$ and
$\mathcal{H}^{\nu,\ell}$ with specific $\nu,\ell$ values).

For any $P$, there exists a value of $\ell$ so that, for all
$w\in\mathcal{W}^{P,2}([0,T])$:

$\|w\|_{\mathcal{W}^{P,2}([0,T])}=\|w\|_{\mathcal{H}^{P-1/2,\ell}([0,T])}.$

##  Appendix B Proof of Theorem LABEL:thm:consistency

Throughout this proof we will drop the domain $[0,T]$ from the $\mathcal{W}$
and $\mathcal{H}$ notation for brevity. We break up this proof into a series
of lemmas.

######  Lemma 1.

For every $({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}$, the functions
${x}^{(1)},\ldots,{x}^{(M_{x})}$, ${y}^{(1)},\ldots,{y}^{(M_{y})}$, and
${z}^{(1)},\ldots,{z}^{(M_{z})}$, are all elements of $\mathcal{W}^{P+1,2}$.

###### Proof.

By Lipschitz continuity of ${\bm{F}}$, ${\bm{G}}$, ${\bm{H}}$ and the Cauchy-
Lipschitz-Picard theorem [e.g. 9, Thm. 7.3], we have that every
$({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}$ are continuous and continuously
differentiable over $[0,\infty)$. This fact implies that
${\bm{x}},{\bm{y}},{\bm{z}}$ and $\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}}$ are
bounded over $[0,T]$. Assuming that $\dot{\bm{x}}$, $\ldots$,
$\overset{(k)}{{\bm{x}}}$, $\dot{\bm{y}}$, $\ldots$,
$\overset{(k)}{{\bm{y}}}$, $\dot{\bm{z}}$, $\ldots$, $\overset{(k)}{{\bm{z}}}$
are also bounded, consider $\overset{(k+1)}{{\bm{x}}}$,
$\overset{(k+1)}{{\bm{y}}}$, $\overset{(k+1)}{{\bm{z}}}$ for $k\leq P$:

$\displaystyle\overset{(k+1)}{{\bm{x}}}(t)$
$\displaystyle=\overset{(k)}{{\bm{F}}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$
$\displaystyle\overset{(k+1)}{{\bm{y}}}(t)$
$\displaystyle=\overset{(k)}{{\bm{G}}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$
$\displaystyle\mathbf{0}$
$\displaystyle=\overset{(k)}{{\bm{H}}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$

By Faà di Bruno’s formula,
$\overset{(k)}{{\bm{F}}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$ is a polynomial
involving the terms

  * •

${\bm{F}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$, $\ldots$,
$\overset{(k+1)}{{\bm{F}}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$,

  * •

${\bm{x}}(t)$, $\ldots$, $\overset{(k-1)}{{\bm{x}}}(t)$,

  * •

${\bm{y}}(t)$, $\ldots$, $\overset{(k-1)}{{\bm{y}}}(t)$, and

  * •

${\bm{z}}(t)$, $\ldots$, $\overset{(k-1)}{{\bm{z}}}(t)$.

By inductive assumption, the latter terms are bounded over $[0,T]$. Moreover,
the range of ${\bm{x}}(t)$, ${\bm{y}}(t)$, and ${\bm{z}}(t)$ will be a compact
interval, and thus ${\bm{F}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$, $\ldots$,
$\overset{(k+1)}{{\bm{F}}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$, will be
bounded by Lipschitz continuity. Therefore, ${\bm{x}}$, $\ldots$,
$\overset{(P+1)}{{\bm{x}}}$ are all bounded over $[0,T]$. Analogous results
hold for ${\bm{y}}$, $\ldots$, $\overset{(P+1)}{{\bm{y}}}$ and ${\bm{z}}$,
$\ldots$, $\overset{(P+1)}{{\bm{z}}}$. Since bounded functions are trivially
square integrable, we have that

  * •

$\|{x}^{(1)}\|_{\mathcal{W}^{P+1,2}},\ldots,\|{x}^{(M_{x})}\|_{\mathcal{W}^{P+1,2}}<\infty$,

  * •

$\|{y}^{(1)}\|_{\mathcal{W}^{P+1,2}},\ldots,\|{y}^{(M_{y})}\|_{\mathcal{W}^{P+1,2}}<\infty$,
and

  * •

$\|{z}^{(1)}\|_{\mathcal{W}^{P+1,2}},\ldots,\|{z}^{(M_{z})}\|_{\mathcal{W}^{P+1,2}}<\infty$.

∎

######  Lemma 2.

Define $M=M_{x}+M_{y}+M_{z}$, and define $\mathcal{W}^{P+1,2}_{M}$ as the
vector-valued reproducing kernel Hilbert space of functions
$({\bm{x}}(\cdot),{\bm{y}}(\cdot),{\bm{z}}(\cdot)):[0,T]\to\mathbb{R}^{M}$
with the property that

${\bm{w}}\in\mathcal{W}^{P+1,2}_{M}\Leftrightarrow
w^{(m)}\in\mathcal{W}^{P+1,2}\>\>\forall m\in[1,M]$

Then the set ${\mathbb{S}}$ is closed with respect to
$\|\cdot\|_{\mathcal{W}^{P+1,2}_{M}}$.

###### Proof.

By standard theory of vector-valued reproducing kernel Hilbert spaces, the
norm associated with $\mathcal{W}^{P+1,2}_{M}$ is

$\|({\bm{x}},{\bm{y}},{\bm{z}})\|_{\mathcal{W}^{P+1,2}_{M}}:=\sqrt{\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,\|{\hat{{x}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P+1,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,\|{\hat{{y}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P+1,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,\|{\hat{{z}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P+1,2}}}.$

Consider a Cauchy sequence
$({\bm{x}}_{n},{\bm{y}}_{n},{\bm{z}}_{n})\in{\mathbb{S}}$ w.r.t.
$\|\cdot\|_{\mathcal{W}^{P+1,2}_{M}}$. By Lemma 1, we have that
${x}_{n}^{(m)}\in\mathcal{W}^{P+1,2}$ for all $m\in[1,M_{x}]$. By completeness
of $\mathcal{W}^{P+1,2}$, we have that ${x}_{n}^{(m)}\to{x}^{(m)}$ for some
${x}^{(m)}\in\mathcal{W}^{P+1,2}$; i.e. for every $\epsilon$, there exists
some $n^{\prime}$ such that
$\|{y}_{n^{\prime}}^{(m)}-{y}^{(m)}\|_{\mathcal{W}^{P+1,2}}<\epsilon$.

$\displaystyle\sup_{t\in[0,T]}\|{x}^{(m)}_{n^{\prime}}(t)-{x}^{(m)}(t)\|_{\infty}$
$\displaystyle=\sup_{t\in[0,T]}\left\langle
k(t,\cdot),{x}^{(m)}_{n^{\prime}}(t)-{x}^{(m)}(t)\right\rangle_{\mathcal{W}^{P+1,2}}$
$\displaystyle\leq\sup_{t\in[0,T]}\|k(t,\cdot)\|_{\mathcal{W}^{P+1,2}}\>\|{x}^{(m)}_{n^{\prime}}(t)-{x}^{(m)}(t)\|_{\mathcal{W}^{P+1,2}}$
$\displaystyle\leq
C\|{x}^{(m)}_{n^{\prime}}(t)-{x}^{(m)}(t)\|_{\mathcal{W}^{P+1,2}}<C\epsilon,$

where $k(\cdot,\cdot)$ is the reproducing kernel associated with
$\mathcal{W}^{P+1,2}$ and $C$ is some universal constant. The penultimate
inequality comes from fact that $\mathcal{W}^{P+1,2}$ is equivalent to a
Matérn RKHS which has a bounded-everywhere reproducing kernel. Thus,
${x}^{(m)}_{n^{\prime}}$ converges uniformly to ${x}^{(m)}$, and so

$\displaystyle\dot{\bm{x}}(t)=\lim_{n\to\infty}\dot{\bm{x}}_{n}(t)$
$\displaystyle=\lim_{n\to\infty}{\bm{F}}({\bm{x}}_{n}(t),{\bm{y}}_{n}(t),{\bm{z}}_{n}(t))$
$\displaystyle={\bm{F}}(\lim_{n\to\infty}{\bm{x}}_{n}(t),\lim_{n\to\infty}{\bm{y}}_{n}(t),\lim_{n\to\infty}{\bm{z}}_{n}(t))={\bm{F}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t))$

where the penultimate equality comes from the continuity of ${\bm{F}}$.
Analogous results hold for ${\bm{y}}_{n}$ and ${\bm{z}}_{n}$. Moreover, by
uniform convergence we have that $\|{\bm{y}}(0)-{\bm{y}}_{n}(0)\|_{\infty}$
and $\|{\bm{z}}(0)-{\bm{z}}_{n}(0)\|_{\infty}$ are arbitrarily small and thus
$\|{\bm{y}}(0)\|_{\infty},\|{\bm{z}}(0)\|_{\infty}<K$. Therefore
$({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}$. ∎

######  Lemma 3.

Denote $M$ and $\mathcal{W}^{P+1,2}_{M}$ as in Lemma 2. Let

$\displaystyle B:$
$\displaystyle={\inf_{({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}}}\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,|{x}^{(m)}_{N}|^{2}_{\mathcal{W}^{P+1,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,|{y}^{(m)}_{N}|^{2}_{\mathcal{W}^{P+1,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,|{z}^{(m)}_{N}|^{2}_{\mathcal{W}^{P+1,2}}$
$\displaystyle={\inf_{({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}}}\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,\|\dot{x}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,\|\dot{y}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,\|\dot{z}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P,2}}$
$\displaystyle={\inf_{({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}}}\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}$

Then there exists some
$({\bm{x}}^{*},{\bm{y}}^{*},{\bm{z}}^{*})\in{\mathbb{S}}$ that achieves this
infimum.

###### Proof.

Define the operator $D:\mathcal{W}^{P+1,2}_{M}\to\mathcal{W}^{P+1,2}_{M}$ as
$D({\bm{x}},{\bm{y}},{\bm{z}})=(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})$,
where here $\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}}$ denote weak derivatives.
Note that $D$ is a surjective and bounded linear operator between two Hilbert
spaces:

$\displaystyle\|({\bm{x}},{\bm{y}},{\bm{z}})\|_{\mathcal{W}^{P+1,2}_{M}}^{2}$
$\displaystyle=\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,\|{\hat{{x}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P+1,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,\|{\hat{{y}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P+1,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,\|{\hat{{z}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P+1,2}}$
$\displaystyle=\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,\|\dot{\hat{{x}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,\|\dot{\hat{{y}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P,2}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,\|\dot{\hat{{z}}}^{(m)}_{N}\|^{2}_{\mathcal{W}^{P,2}}$
$\displaystyle=\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}^{2}$

Furthermore, note that the nullspace of $D$ (i.e. the set of constant
functions) is a closed set. Therefore, since
${\mathbb{S}}\subset\mathcal{W}^{P+1,2}_{M}$ is a closed subset of a Hilbert
space (Lemma 2),
$\dot{\mathbb{S}}:=D({\mathbb{S}})\subset\mathcal{W}^{P+1,2}_{M}$ is also
closed subset of a Hilbert space [see e.g. 9, Exercise 2.10]. By the existence
portion of the Hilbert projection theorem, there exists a (potentially non-
unique) minimum $\mathcal{W}^{P,2}_{M}$-norm element of $\dot{\mathbb{S}}$
(i.e. there exists some
$(\dot{\bm{x}}^{*},\dot{\bm{y}}^{*},\dot{\bm{z}}^{*})\in\dot{\mathbb{S}}$ such
that
$\|(\dot{\bm{x}}^{*},\dot{\bm{y}}^{*},\dot{\bm{z}}^{*})\|_{\mathcal{W}^{P,2}_{M}}=\inf_{({\bm{x}},{\bm{y}},{\bm{z}})\in{\mathbb{S}}^{\prime}}\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}.$
We conclude the proof by setting $({\bm{x}}^{*},{\bm{y}}^{*},{\bm{z}}^{*})$ to
be some element in ${\mathbb{S}}$ such that
$(\dot{\bm{x}}^{*\prime},\dot{\bm{y}}^{*\prime},\dot{\bm{z}}^{*\prime})=D({\bm{x}}^{*},{\bm{y}}^{*},{\bm{z}}^{*})$.
∎

######  Lemma 4.

For any $0<C<\infty$, define the sets

$\begin{gathered}{\mathbb{F}}^{C}:=\\{w\in\mathcal{W}^{P,2}\>:\>\|w\|_{\mathcal{W}^{P,2}}\leq
C\\}\\\
\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}:=\\{\operatorname{{\textstyle\int}}_{0}^{(\cdot)}w(\tau)d\tau\>:\>w\in{\mathbb{F}}\\}.\end{gathered}$

Denoting $\widehat{\mathcal{R}}_{N}$ as the empirical Rademacher complexity
for some dataset $t_{1},\ldots,t_{N}\in[0,T]$, we have that

$\widehat{\mathcal{R}}_{N}({\mathbb{F}}^{C})\lesssim
CN^{-1/2},\qquad\widehat{\mathcal{R}}_{N}(\operatorname{{\textstyle\int}}{\mathbb{F}}^{C})\lesssim
TCN^{-1/2}.$

###### Proof.

The Rademacher complexity $\widehat{\mathcal{R}}_{N}({\mathbb{F}}^{C})\lesssim
CN^{-1/2}$ follows a standard result for reproducing kernel Hilbert spaces,
using the fact that $\mathcal{W}^{P,2}$ is equivalent to the Matérn RKHS which
has a bounded-everywhere reproducing kernel. Bounding the Rademacher
complexity of $\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}$ mirrors the
standard proof of the $\widehat{\mathcal{R}}_{N}({\mathbb{F}}^{C})$ bound:

$\displaystyle\widehat{\mathcal{R}}(\operatorname{{\textstyle\int}}{\mathbb{F}}^{C})$
$\displaystyle:=\mathbb{E}_{\epsilon_{i}}\left[\sup_{w\in\mathcal{W}^{P,2}}\frac{1}{N}\sum_{i=1}^{N}\epsilon_{i}\int_{0}^{t_{i}}w(\tau)d\tau\right]$
($\epsilon_{i}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Rad}$)
$\displaystyle=\mathbb{E}_{\epsilon_{i}}\left[\sup_{w\in\mathcal{W}^{P,2}}\left\langle\frac{1}{N}\sum_{i=1}^{N}\epsilon_{i}\int_{0}^{t_{i}}k(\tau,\cdot)d\tau,\>w(\cdot)\right\rangle_{\mathcal{W}^{P,2}}\right]$
$\displaystyle\leq\mathbb{E}_{\epsilon_{i}}\left[\left\|\frac{C}{N}\sum_{i=1}^{N}\epsilon_{i}\int_{0}^{t_{i}}k(\tau,\cdot)d\tau\right\|_{\mathcal{W}^{P,2}}\right]$
(Cauchy-Schwarz inequality)
$\displaystyle\leq\sqrt{\mathbb{E}_{\epsilon_{i}}\left[\left\|\frac{C}{N}\sum_{i=1}^{N}\epsilon_{i}\int_{0}^{t_{i}}k(\tau,\cdot)d\tau\right\|_{\mathcal{W}^{P,2}}^{2}\right]}$
(Jensen inequality)
$\displaystyle=\sqrt{\mathbb{E}_{\epsilon_{i}}\left[\frac{C^{2}}{N^{2}}\sum_{i,j=1}^{N}\epsilon_{i}\epsilon_{j}\left\langle\int_{0}^{t_{i}}k(\tau,\cdot)d\tau,\>\int_{0}^{t_{j}}k(\tau,\cdot)d\tau\right\rangle_{\mathcal{W}^{P,2}}\right]}$
$\displaystyle=\sqrt{\frac{C^{2}}{N^{2}}\sum_{i}^{N}\left\|\int_{0}^{t_{i}}k(\tau,\cdot)d\tau\right\|_{\mathcal{W}^{P,2}}^{2}}$
($\epsilon_{i}$ are uncorrelated)
$\displaystyle\leq\sqrt{\frac{C^{2}}{N^{2}}\sum_{i}^{N}\left(\int_{0}^{t_{i}}\left\|k(\tau,\cdot)\right\|_{\mathcal{W}^{P,2}}d\tau\right)^{2}}$
(triangle inequality)
$\displaystyle\leq\sqrt{\frac{C^{2}}{N^{2}}\sum_{i}^{N}\left(\int_{0}^{T}\sup_{t\in[0,T]}\left\|k(t,\cdot)\right\|_{\mathcal{W}^{P,2}}d\tau\right)^{2}}$
$\displaystyle=TCN^{-1/2}\sup_{t\in[0,T]}\left\|k(t,\cdot)\right\|_{\mathcal{W}^{P,2}}.$

Recognizing that $k(t,\cdot)$ is equivalent to a bounded-everywhere
reproducing kernel completes the proof. ∎

######  Lemma 5.

Define $M$ and $\mathcal{W}^{P,2}_{D}$ as in Lemma 2. For any $0<C<\infty$,
define the sets

$\begin{gathered}{\mathbb{F}}^{C}_{M}:=\left\\{(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\in\mathcal{W}^{P,2}_{M}\>:\>\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{D}}\leq
C\right\\}\\\
\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}_{M}:=\left\\{\operatorname{{\textstyle\int}}_{0}^{(\cdot)}(\dot{\bm{x}}(\tau),\dot{\bm{y}}(\tau),\dot{\bm{z}}(\tau))d\tau\>:\>(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\in{\mathbb{F}}\right\\}\end{gathered}$

Then $\widehat{\mathcal{R}}_{N}({\mathbb{F}}^{C}_{M})\lesssim CMN^{-1/2}$ and
$\widehat{\mathcal{R}}_{N}(\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}_{M})\lesssim
TMCN^{-1/2}$.

###### Proof.

The proof follows a standard summation argument for Rademacher complexity:

$\displaystyle\widehat{\mathcal{R}}({\mathbb{F}}^{C}_{M})$ $\displaystyle:=$
$\displaystyle\mathbb{E}\left[\sup_{(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\in\mathcal{W}^{P,2}_{M}}\frac{1}{N}\sum_{i=1}^{N}\left(\sum_{j_{x}=1}^{M_{x}}\epsilon_{ij_{x}}\dot{x}^{(j_{x})}(t_{i})+\sum_{j_{y}=1}^{M_{y}}\epsilon_{ij_{y}}\dot{y}^{(j_{y})}(t_{i})+\sum_{j_{z}=1}^{M_{y}}\epsilon_{ij_{z}}\dot{z}^{(j_{z})}(t_{i})\right)\right]$
($\epsilon_{ij_{x}},\epsilon_{ij_{y}},\epsilon_{ij_{z}}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Rad}$)
$\displaystyle\leq$
$\displaystyle\sum_{j_{x}=1}^{M_{x}}\mathbb{E}\left[\sup_{\dot{x}^{(j_{x})}\in\mathcal{W}^{P,2}}\frac{1}{N}\sum_{i=1}^{N}\epsilon_{ij_{x}}\dot{x}^{(j_{x})}(t_{i})\right]+\sum_{j_{y}=1}^{M_{y}}\mathbb{E}\left[\sup_{\dot{y}^{(j_{y})}\in\mathcal{W}^{P,2}}\frac{1}{N}\sum_{i=1}^{N}\epsilon_{ij_{y}}\dot{y}^{(j_{y})}(t_{i})\right]$
$\displaystyle+\sum_{j_{z}=1}^{M_{z}}\mathbb{E}\left[\sup_{\dot{z}^{(j_{z})}\in\mathcal{W}^{P,2}}\frac{1}{N}\sum_{i=1}^{N}\epsilon_{ij_{z}}\dot{z}^{(j_{z})}(t_{i})\right]$
$\displaystyle\lesssim$ $\displaystyle\>CMN^{-1/2},$ (Lemma 4)

where the last inequality comes from the fact that
$\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}\leq C$
implies that $\|\dot{\bm{x}}^{(m)}\|_{\mathcal{W}^{P,2}}$,
$\|\dot{\bm{y}}^{(m)}\|_{\mathcal{W}^{P,2}}$,
$\|\dot{\bm{z}}^{(m)}\|_{\mathcal{W}^{P,2}}\leq C$ for all $m$. An analogous
proof holds for
$\widehat{\mathcal{R}}(\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}_{M})$.
∎

######  Lemma 6.

Denote $M$ and $\mathcal{W}^{P,2}_{M}$ as in Lemma 2. For any
$(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\in\mathcal{W}^{P,2}_{M}$,
$\hat{\bm{y}}_{0}\in\mathbb{R}^{M_{y}}$,
$\hat{\bm{z}}_{0}\in\mathbb{R}^{M_{z}}$, define the differential equation
error function

$\begin{gathered}e_{\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(\cdot):=\left\|\begin{bmatrix}\dot{\bm{x}}(\cdot)-{\bm{F}}\left({\bm{x}},{\bm{y}},{\bm{z}}\right)\\\
\dot{\bm{y}}(\cdot)-{\bm{G}}\left({\bm{x}},{\bm{y}},{\bm{z}}\right)\\\
{\bm{H}}\left({\bm{x}},{\bm{y}},{\bm{z}}\right)\end{bmatrix}\right\|_{\infty}\\\
{\bm{x}}(t):={\bm{x}}_{0}+\int_{0}^{t}\dot{\bm{x}}(\tau)d\tau,\quad{\bm{y}}(t):=\hat{\bm{y}}_{0}+\int_{0}^{t}\dot{\bm{y}}(\tau)d\tau,\quad{\bm{z}}(t):=\hat{\bm{z}}_{0}+\int_{0}^{t}\dot{\bm{z}}(\tau)d\tau,\end{gathered}$
(17)

For any $0<C<\infty$, define ${\mathbb{G}}^{C,K}$ as the set of error
functions

$\left\\{e_{\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(\cdot)\>:\>\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}\leq
C,\|\hat{\bm{y}}_{0}\|_{\infty}\leq K,\|\hat{\bm{z}}_{0}\|_{\infty}\leq
K\right\\}.$

Then every error function in ${\mathbb{G}}^{C,K}$ is bounded by some constant
$\tilde{C}$ and $\widehat{\mathcal{R}}({\mathbb{G}}^{C,K})\lesssim CN^{-1/2}$.

###### Proof.

Note that each error function in ${\mathbb{G}}^{C,K}$ is a Lipschitz function
($\|\cdot\|_{\infty}$), each of which is applied to the summation of two sub-
functions:

  1. 1.

a vector-valued RKHS function
($(\dot{\bm{x}}(\cdot),\dot{\bm{y}}(\cdot),\dot{\bm{z}}(\cdot))$) with norm
less than $C$ (i.e. an element of ${\mathbb{F}}^{C}_{M}$, as defined in Lemma
5), and

  2. 2.

lipschitz functions (${\bm{F}}$, ${\bm{G}}$, ${\bm{H}}$) applied to the
integral of a vector-valued RKHS function
($\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}}$) with norm less than $C$ (i.e. a
Lipschitz function applied to an element of
$\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}_{M}$, as defined in Lemma 5).

Boundedness of the error functions falls from the fact
$({\bm{x}},{\bm{y}},{\bm{z}})$ are bounded,
$\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}$ are bounded, and ${\bm{F}}$, ${\bm{G}}$,
${\bm{H}}$ are continuous (and thus bounded over $[0,T]$). The Rademacher
complexity falls from standard Lipschitz and summation rules:
$\widehat{\mathcal{R}}({\mathbb{G}}^{C,K})\lesssim\widehat{\mathcal{R}}({\mathbb{F}}^{C}_{M})+\widehat{\mathcal{R}}(\operatorname{{\textstyle\int}}{\mathbb{F}}^{C}_{M})\lesssim
CMN^{-1/2}.$ ∎

Now we are ready to prove LABEL:thm:consistency.

###### Proof of LABEL:thm:consistency.

We begin by noting that

$\displaystyle B:=$
$\displaystyle\inf_{({\bm{x}},{\bm{y}},{\bm{z}}\in{\mathbb{S}}}\operatorname{{\textstyle\sum}}_{m=1}^{M_{x}}\,|{\hat{{x}}}^{(m)}|^{2}_{\mathcal{H}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{y}}\,|{\hat{{y}}}^{(m)}|^{2}_{\mathcal{H}}+\operatorname{{\textstyle\sum}}_{m=1}^{M_{z}}\,|{\hat{{z}}}^{(m)}|^{2}_{\mathcal{H}}$
$\displaystyle\asymp$
$\displaystyle\inf_{({\bm{x}},{\bm{y}},{\bm{z}}\in{\mathbb{S}}}\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}<\infty$

is implied by Lemma 1. (The equivalence in norms is given by combining the
definition of $\mathcal{W}^{P,2}_{M}$ in Lemma 2 with the equivalence between
Matérn and Sobolev spaces.) Let $({\bm{x}}^{*},{\bm{y}}^{*},{\bm{z}}^{*})$ be
some element in ${\mathbb{S}}$ that achieves this infimum (the existence of
which is guaranteed by Lemma 3). We know that
$\|(\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N})\|_{\mathcal{W}^{P,2}_{M}}\leq\|(\dot{\bm{x}}^{*},\dot{\bm{y}}^{*},\dot{\bm{z}}^{*})\|_{\mathcal{W}^{P,2}_{M}}=B$—since
$({\bm{x}}^{*},{\bm{y}}^{*},{\bm{z}}^{*})$ satisfies the constraints for Eq.
7—and thus
$(\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N})\in{\mathbb{F}}^{B}_{M}$
(as defined by Lemma 5).

Defining
$e_{\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(\cdot)$
as in Lemma 6, we have that
$e_{\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(t_{i})=0$
for each $t_{i}$ in $\mathcal{D}$. Applying a standard uniform large law
argument [e.g. 46, Thm. 4.2] we have that, for any $\delta>0$,

$\displaystyle\int_{0}^{T}e_{\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(\tau)d\tau$
$\displaystyle=\left|\frac{1}{N}\sum_{i=1}^{N}e_{\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(t_{i})-\int_{0}^{T}e_{\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(\tau)d\tau\right|$
$\displaystyle\leq
2\mathbb{E}_{t_{i}}\left[\widehat{\mathcal{R}}({\mathbb{G}}^{B,K})\right]+\delta$

with probability $1-2\exp(-\frac{N\delta^{2}}{8\tilde{C}^{2}})$ (where
$\tilde{C}$ is the constant defined in Lemma 6). Since
$\widehat{\mathcal{R}}({\mathbb{G}}^{B,K})\lesssim BMN^{-1/2}$ (Lemma 6), we
have that
$\int_{0}^{T}e_{\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N},\hat{\bm{y}}_{0},\hat{\bm{z}}_{0}}(\tau)d\tau\overset{\mathrm{a.s.}}{\longrightarrow}0$,
which implies that
$\lim_{N\to\infty}({\hat{\bm{x}}}_{N},{\hat{\bm{y}}}_{N},{\hat{\bm{z}}}_{N})$
satisfies the differential equation almost everywhere.

Define $(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})$ as the continuously-
differentiable representative of
$\lim_{N\to\infty}(\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N})$
(see [e.g. 9, Thm. 8.2]). Since $(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})$ is
continuously-differentiable and satisfies the differential equation
everywhere, it must also be an element of ${\mathbb{S}}$. All together, this
implies that
$\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}\geq B$,
and so
$\|(\dot{\bm{x}},\dot{\bm{y}},\dot{\bm{z}})\|_{\mathcal{W}^{P,2}_{M}}=\|\lim_{N\to\infty}(\dot{\hat{\bm{x}}}_{N},\dot{\hat{\bm{y}}}_{N},\dot{\hat{\bm{z}}}_{N})\|_{\mathcal{W}^{P,2}_{M}}=B$.
∎

##  Appendix C Robustness Checks and Variations

###  C.1 Sparse training data and data efficiency

Figure 4: The Neoclassical Growth model (i.e., Eqs. 14 and 13) solved using
the Matérn kernel for a sparse training data. The top-left panel shows the
approximate and benchmark capital paths. The top-right panel shows the
relative errors between the approximate and benchmark solutions for capital.
The bottom-left panel shows the approximate and benchmark consumption paths.
The bottom-right panel shows the relative errors between the approximate and
benchmark solutions for consumption. Accurate solutions can be obtained even
with sparse training data.

The cardinality of the training set, $\mathcal{D}$, becomes an impediment in
higher dimensions, where we have many state and jump variables. Therefore, it
is crucial that we can obtain accurate approximate solutions with a sparse
training data.

Fig. 4 shows the result of the Neoclassical Growth model (i.e., Eqs. 14 and
13) solved with the Matérn kernel for sparse training data,
$\mathcal{D}:={1,5,10,15,20,25,30}$. The top-left panel shows the approximate
and benchmark capital paths, denoted by $\hat{x}(t)$ and $x(t)$, respectively.
The top-right panel shows the relative errors between the approximate and
benchmark solutions for capital, denoted by $\varepsilon_{x}(t)$. The bottom-
left panel shows the approximate and benchmark consumption paths, denoted by
$\hat{y}(t)$ and $y(t)$, respectively. The bottom-right panel shows the
relative errors between the approximate and benchmark solutions for
consumption, denoted by $\varepsilon_{y}(t)$.

These results show that one can obtain very accurate approximate solution,
even with a very sparse training data.

###  C.2 Robustness to the choice of the kernel and kernel parameters

Sec. C.2 shows the result of the approximate solution of the Neoclassical
Growth model, described in Eqs. 14 and 13), for different Matérn kernels and
kernel parameters.

The first three rows show the performance of the approximate solutions for
three different kernels. We report the maximum and minimum of the absolute
value of relative errors for both the capital path, $\hat{x}(t)$, and the
consumption path, $\hat{y}(t)$. The first row presents the baseline solution,
Matérn kernel with $\nu=\frac{1}{2}$. The second row shows results for Matérn
kernel with $\nu=\frac{3}{2}$, defined as

$\displaystyle
K(t_{i},t_{j})=C_{\frac{3}{2}}(t_{i},t_{j})=\sigma^{2}\left(1+\frac{\sqrt{3}|t_{i}-t_{j}|}{\ell}\right)\exp\left(-\frac{\sqrt{3}|t_{i}-t_{j}|}{\ell}\right),$

and the third row for $\nu=\frac{5}{2}$, defined as

$\displaystyle
K(t_{i},t_{j})=C_{\frac{5}{2}}(t_{i},t_{j})=\sigma^{2}\left(1+\frac{\sqrt{5}|t_{i}-t_{j}|}{\ell}+\frac{5|t_{i}-t_{j}|^{2}}{3\ell^{2}}\right)\exp\left(-\frac{\sqrt{5}|t_{i}-t_{j}|}{\ell}\right).$

Table 1: The robustness of the approximate solutions of the Neoclassical
Growth model (i.e., Eqs. 14 and 13) is tested using different Matérn kernels,
$\nu=\frac{1}{2},\frac{3}{2},\frac{5}{2}$, and length scales $\ell=2,20$. We
report the maximum and minimum of the absolute value of relative errors for
both the capital path, $\hat{x}(t)$, and the consumption path, $\hat{y}(t)$.
Our method is not sensitive to the choice of Matérn kernels and length scales.

The last two rows show the performance of the approximate solutions for two
different length scales, $\ell=2$ and $\ell=20$.

Throughout these experiments, we achieve highly accurate approximate
solutions. Therefore, the results demonstrate insensitivity to the selection
of Matérn kernels and the length scales.

###  C.3 Smaller time horizons: accurate short-run dynamics

One might suspect that achieving an accurate optimal solution, which does not
violate the transversality condition, Eq. 15, is only possible if one uses a
large time horizon in the training data. For instance, we use
$\mathcal{D}={0,1,2,\ldots,30}$ to obtain the results depicted in Fig. 2. In
this experiment, we establish that we can still achieve accurate short-run
dynamics by using a smaller time horizon.

Figure 5: The Neoclassical Growth model (i.e., Eqs. 14 and 13) solved with
Matérn kernel for training data with a smaller time horizon. The top-left
panel shows the approximate and benchmark capital paths. The top-right panel
shows the relative errors between the approximate and benchmark solutions for
capital. The bottom-left panel shows the approximate and benchmark consumption
paths. The bottom-right panel shows the relative errors between the
approximate and benchmark solutions for consumption. Very accurate short-run
dynamics can be obtained with training data using a smaller time horizon.

Fig. 5 shows the approximate solutions for the Neoclassical Growth model
(i.e., Eqs. 14 and 13) for training data with smaller time horizon, defined as
$\mathcal{D}:=\\{0,1,2,\cdots,10\\}$. The top-left panel shows the approximate
and benchmark capital paths, denoted by $\hat{x}(t)$ and $x(t)$, respectively.
The top-right panel shows the relative errors between the approximate and
benchmark solutions for capital, denoted by $\varepsilon_{x}(t)$. The bottom-
left panel shows the approximate and benchmark consumption paths, denoted by
$\hat{y}(t)$ and $y(t)$, respectively. The bottom-right panel shows the
relative errors between the approximate and benchmark solutions for
consumption, denoted by $\varepsilon_{y}(t)$

###  C.4 Neoclassical Growth model: differential-algebraic equations

In this section, we solve a different formulation of the problem introduced in
section 4.2. Introducing a static variable, ${\bm{z}}(t)\in\mathbb{R}$,
defined as

$\displaystyle{\bm{z}}(t):=f^{\prime}({\bm{x}}(t))-\delta,$

we can reformulate Eqs. 13 and 14 in the following way

$\displaystyle\dot{{\bm{x}}}(t)$
$\displaystyle=f({\bm{x}}(t))-{\bm{y}}(t)-\delta{\bm{x}}(t),$ (18)
$\displaystyle\dot{{\bm{y}}}(t)$
$\displaystyle={\bm{y}}(t)\left[{\bm{z}}(t)-r\right],$ (19) $\displaystyle 0$
$\displaystyle={\bm{z}}(t)-f^{\prime}({\bm{x}}(t))+\delta,$ (20)

where $f(x)=x^{a}$, $f^{\prime}(x)=ax^{a-1}$, and the transversality
condition, Eq. 15 remains unchanged. The static variable, ${\bm{z}}(t)$,
represents the real interest rate in this economy.

What distinguishes this example from the Neoclassical Growth model is the
presence of an algebraic equation, Eq. 20, in the equations. Formally, using
the notation introduced in Eq. 3

$\displaystyle{\bm{H}}({\bm{x}}(t),{\bm{y}}(t),{\bm{z}}(t)):={\bm{z}}(t)-f^{\prime}({\bm{x}}(t))+\delta.$

In this example, we aim to investigate whether our method can handle
differential-algebraic equations.

Using
${\hat{\bm{x}}}(t)={\bm{x}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{x}}}}(\tau)d\tau$,
${\hat{\bm{y}}}(t)={\hat{\bm{y}}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{y}}}}(\tau)d\tau$,
and
${\hat{\bm{z}}}(t)={\hat{\bm{z}}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{z}}}}(\tau)d\tau$,
we solve the following optimization

$\displaystyle\min_{\hat{{\bm{x}}},\hat{{\bm{y}}},\hat{{\bm{z}}},{\hat{\bm{y}}}_{0},{\hat{\bm{z}}}_{0}}$
$\displaystyle\,\|\dot{{\hat{\bm{x}}}}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{y}}}}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{z}}}}\|^{2}_{\mathcal{H}}$
s.t.
$\displaystyle\dot{{\hat{\bm{x}}}}(t_{i})=f\left(\hat{{\bm{x}}}(t_{i})\right)-\hat{{\bm{y}}}(t_{i})-\delta\hat{{\bm{x}}}(t_{i}),\quad\text{for
all }t_{i}\in\mathcal{D},$
$\displaystyle\dot{{\hat{\bm{y}}}}(t_{i})=\hat{{\bm{y}}}(t_{i})\left[\hat{{\bm{z}}}(t_{i})-r\right],\quad\text{for
all }t_{i}\in\mathcal{D},$ $\displaystyle
0=\hat{{\bm{z}}}(t_{i})-f^{\prime}\left(\hat{{\bm{x}}}(t_{i})\right)+\delta,\quad\text{for
all }t_{i}\in\mathcal{D},$

to obtain approximate solutions for $x,y$, and $z$.

Here,
$\|\dot{{\hat{\bm{x}}}}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{x\top}K\bm{\alpha}^{x}$,
$\|\dot{{\hat{\bm{y}}}}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{y\top}K\bm{\alpha}^{y}$,
and
$\|\dot{{\hat{\bm{z}}}}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{z\top}K\bm{\alpha}^{z}$.
We use the same numerical values introduced in section 4.2.

Figure 6: DAE formulation of the Neoclassical Growth model (i.e., Eqs. 18 and
20) solved with Matérn kernel. The top-left panel shows the approximate and
benchmark capital paths. The top-right panel shows the relative errors between
the approximate and benchmark solutions for capital. The bottom-left panel
shows the approximate and benchmark consumption paths. The bottom-right panel
shows the relative errors between the approximate and benchmark solutions for
consumption. Accurate approximate solutions can be obtained for the DAE
formulation.

Fig. 6 shows the approximate solutions of the DAE formulation of the
Neoclassical Growth model for the training data
$\mathcal{D}=\\{0,1,\cdots,30\\}$. The top-left panel shows the approximate
and benchmark capital paths, denoted by $\hat{x}(t)$ and $x(t)$, respectively.
The top-right panel shows the relative errors between the approximate and
benchmark solutions for capital, denoted by $\varepsilon_{x}(t)$. The bottom-
left panel shows the approximate and benchmark consumption paths, denoted by
$\hat{y}(t)$ and $y(t)$, respectively. The bottom-right panel shows the
relative errors between the approximate and benchmark solutions for
consumption, denoted by $\varepsilon_{y}(t)$

This result establishes that our method can find accurate approximate
solutions for the DAE formulation of the Neoclassical Growth model.

##  Appendix D More Applications

In this section, we provide more applications in economics and business. Each
application differs fundamentally from the problems introduced in Sec. 4. In
Sec. D.1, we solve a model that incorporates the time evolution of human
capital and its interaction with physical capital in economic growth. This
model has five variables: two state variables and three jump variables. In
Sec. D.2, we solve an optimal advertising model. This model is considered a
building block of modern dynamic advertising in business economics and
operations management.

###  D.1 Human capital and growth

In this example, we solve the Neoclassical Growth model with human and
physical capital as illustrated in [1]. The optimal paths
${\bm{x}}_{k}(t),{\bm{x}}_{h}(t),{\bm{y}}_{c}(t),{\bm{y}}_{k}(t),{\bm{y}}_{h}(t)$,
solve

$\displaystyle\dot{{\bm{x}}}_{k}(t)$
$\displaystyle={\bm{y}}_{k}(t)-\delta_{k}{\bm{x}}_{k}(t),$ (21)
$\displaystyle\dot{{\bm{x}}}_{h}(t)$
$\displaystyle={\bm{y}}_{h}(t)-\delta_{k}{\bm{x}}_{h}(t),$ (22)
$\displaystyle\dot{{\bm{y}}}_{c}(t)$
$\displaystyle={\bm{y}}_{c}(t)\left[f_{1}\left({\bm{x}}_{k}(t),{\bm{x}}_{h}(t)\right)-\delta_{k}-r\right],$
(23) $\displaystyle 0$
$\displaystyle=f\left({\bm{x}}_{k}(t),{\bm{x}}_{h}(t)\right)-{\bm{y}}_{c}(t)-{\bm{y}}_{k}(t)-{\bm{y}}_{h}(t),$
(24) $\displaystyle 0$
$\displaystyle=f_{2}\left({\bm{x}}_{k}(t),{\bm{x}}_{h}(t)\right)-f_{1}\left({\bm{x}}_{k}(t),{\bm{x}}_{h}(t)\right)+\delta_{k}-\delta_{h},$
(25)

for given initial conditions ${\bm{x}}_{k}(0)={\bm{x}}_{k_{0}}$,
${\bm{x}}_{h}(0)={\bm{x}}_{h_{0}}$, and two transversality conditions

$\displaystyle 0$
$\displaystyle=\lim_{t\rightarrow\infty}e^{-rt}\frac{{\bm{x}}_{k}(t)}{{\bm{y}}_{c}(t)},$
(26) $\displaystyle 0$
$\displaystyle=\lim_{t\rightarrow\infty}e^{-rt}\frac{{\bm{x}}_{h}(t)}{{\bm{y}}_{c}(t)}.$
(27)

The production function is defined as
$f\left({\bm{x}}_{k}(t),{\bm{x}}_{h}(t)\right)={\bm{x}}_{k}(t)^{a_{k}}{\bm{x}}_{h}(t)^{a_{h}}$.
Here, $f_{1}\left(\cdot,\cdot\right)$ is the derivative with respect to the
first input and $f_{2}\left(\cdot,\cdot\right)$ is the derivative with respect
to the second input. The two constants in the production function, $a_{k}$ and
$a_{h}$, are positive numbers, such that $a_{k}+a_{h}<1$. Additionally,
$\delta_{k}>0$, $\delta_{h}>0$, and $r>0$.

Human capital is denoted by ${\bm{x}}_{h}(t)$, physical capital by
${\bm{x}}_{k}(t)$, consumption by ${\bm{y}}_{c}(t)$, investment in human
capital by ${\bm{y}}_{h}(t)$, and investment in physical capital by
${\bm{y}}_{k}(t)$. Here ${\bm{y}}_{c}(t),{\bm{y}}_{h}(t)$, and
${\bm{y}}_{k}(t)$ are the jump variables. The state variables are
${\bm{x}}_{h}(t)$ and ${\bm{x}}_{k}(t)$ .

This problem is more challenging than the Neoclassical Growth model introduced
in Sec. 4 because it is high-dimensional and involves two algebraic equations,
Eqs. 24 and 25.

Using
${\hat{\bm{x}}}_{k}(t)={\bm{x}}_{k_{0}}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{x}}}}_{k}(\tau)d\tau$,
${\hat{\bm{x}}}_{h}(t)={\bm{x}}_{h_{0}}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{x}}}}_{h}(\tau)d\tau$,
${\hat{\bm{y}}}_{c}(t)={\hat{\bm{y}}}_{c_{0}}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{y}}}}_{c}(\tau)d\tau$,
${\hat{\bm{y}}}_{h}(t)={\hat{\bm{y}}}_{h_{0}}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{y}}}}_{h}(\tau)d\tau$,
${\hat{\bm{y}}}_{k}(t)={\hat{\bm{y}}}_{k_{0}}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{y}}}}_{k}(\tau)d\tau$,
we solve the following optimization

$\displaystyle\min_{\hat{{\bm{x}}}_{k},\hat{{\bm{x}}}_{h},\hat{{\bm{y}}}_{c},\hat{{\bm{y}}}_{h},\hat{{\bm{y}}}_{k},{\hat{\bm{y}}}_{c_{0}},{\hat{\bm{y}}}_{h_{0}},{\hat{\bm{y}}}_{k_{0}}}$
$\displaystyle\,\|\dot{{\hat{\bm{x}}}}_{k}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{x}}}}_{h}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{y}}}}_{c}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{y}}}}_{h}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{y}}}}_{k}\|^{2}_{\mathcal{H}}$
s.t.
$\displaystyle\dot{{\hat{\bm{x}}}}_{k}(t_{i})=\hat{{\bm{y}}}_{k}(t_{i})-\delta_{k}\hat{{\bm{x}}}_{k}(t_{i}),\quad\text{for
all }t_{i}\in\mathcal{D},$
$\displaystyle\dot{{\hat{\bm{x}}}}_{h}(t_{i})=\hat{{\bm{y}}}_{h}(t_{i})-\delta_{h}\hat{{\bm{x}}}_{h}(t_{i}),\quad\text{for
all }t_{i}\in\mathcal{D},$
$\displaystyle\dot{{\hat{\bm{y}}}}_{c}(t_{i})=\hat{{\bm{y}}}_{c}(t_{i})\left[f_{1}\left(\hat{{\bm{x}}}_{k}(t_{i}),\hat{{\bm{x}}}_{h}(t_{i})\right)-\delta_{k}-r\right],\quad\text{for
all }t_{i}\in\mathcal{D},$ $\displaystyle
0=f\left(\hat{{\bm{x}}}_{k}(t_{i}),\hat{{\bm{x}}}_{h}(t_{i})\right)-\hat{{\bm{y}}}_{c}(t_{i})-\hat{{\bm{y}}}_{k}(t_{i})-\hat{{\bm{y}}}_{h}(t_{i}),\quad\text{for
all }t_{i}\in\mathcal{D},$ $\displaystyle
0=f_{2}\left(\hat{{\bm{x}}}_{k}(t_{i}),\hat{{\bm{x}}}_{h}(t_{i})\right)-f_{1}\left(\hat{{\bm{x}}}_{k}(t_{i}),\hat{{\bm{x}}}_{h}(t_{i})\right)+\delta_{k}-\delta_{h},\quad\text{for
all }t_{i}\in\mathcal{D},$

to obtain approximate solutions for
${\bm{x}}_{k}(t),{\bm{x}}_{h}(t),{\bm{y}}_{c}(t),{\bm{y}}_{k}(t),{\bm{y}}_{h}(t)$.

Here,
$\|\dot{{\bm{x}}}_{k}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{x_{k}\top}K\bm{\alpha}^{x_{k}}$,
$\|\dot{{\hat{\bm{x}}}}_{k}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{x_{h}\top}K\bm{\alpha}^{x_{h}}$,
$\|\dot{{\hat{\bm{y}}}}_{c}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{y_{c}\top}K\bm{\alpha}^{y_{c}}$,
$\|\dot{{\hat{\bm{y}}}}_{h}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{y_{h}\top}K\bm{\alpha}^{y_{h}}$,
and
$\|\dot{{\hat{\bm{y}}}}_{k}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{y_{k}\top}K\bm{\alpha}^{y_{k}}$.

In this experiment we, use $\delta_{k}=0.1$, $\delta_{h}=0.05$,
$\alpha_{k}=\frac{1}{3}$, $\alpha_{h}=\frac{1}{4}$, $r=0.11$, $x_{k_{0}}=1.5$,
and $x_{h_{0}}=1.37$ as the numerical values for the economic parameters. We
use $\mathcal{D}=\\{0,1,\cdots,80\\}$ as the training data.

Figure 7: The Neoclassical Growth model with human capital (i.e., Eqs. 21, 22,
23, 24 and 25) is solved using the Matérn kernel for the training data
$\mathcal{D}=\\{0,1,\cdots,80\\}$. The top-left panel shows the approximate
physical capital path. The top-right panel shows the approximate human capital
path. The middle-left panel shows the approximate consumption paths. The
middle-right panel shows the approximate investment in physical capital path.
The bottom-left panel shows the approximate investment in human capital path.
The approximate solutions demonstrate convergence to the correct set of steady
states. Therefore, they satisfy the transversality conditions.

Fig. 7 shows the approximate physical capital path, denoted by
$\hat{x}_{k}(t)$, the approximate human capital path, denoted by
$\hat{x}_{h}(t)$, the approximate consumption path, denoted by
$\hat{y}_{c}(t)$, the approximate investment in physical capital path, denoted
by $\hat{y}_{k}(t)$, and the approximate investment in human capital path,
denoted by $\hat{y}_{h}(t)$. The vertical dashed line shows the boundary
between the interpolation and extrapolation regions.

This result establishes that even in higher dimensions with more complex
algebraic equations, the approximate solutions converge to the correct set of
steady states. Therefore, the solutions satisfy the transversality conditions.

#### The correct set of steady states.

How do we know the approximate solutions converge to the correct set of steady
states? A set of solutions that violates the transversality conditions is
characterized by a consumption path that converges to zero, and human and
physical capital paths that converge to finite but very large values. Hence,
in this case,
$\lim_{t\rightarrow\infty}e^{-rt}\frac{{\bm{x}}_{k}(t)}{{\bm{y}}_{c}(t)}=\lim_{t\rightarrow\infty}e^{-rt}\frac{{\bm{x}}_{h}(t)}{{\bm{y}}_{c}(t)}=\infty$.
Since the approximate solutions shown in Fig. 7 do not exhibit this behavior,
we can be confident that they converge to the correct set of steady states.

###  D.2 Optimal Advertising

In this example we solve an optimal advertising model based on the classical
model of expenditure on advertising introduced in [45]. The optimal paths
${\bm{x}}(t)$ and ${\bm{y}}(t)$, solve

$\displaystyle\dot{{\bm{x}}}(t)$
$\displaystyle=\left[1-{\bm{x}}(t)\right]^{\frac{1}{1-\kappa}}\left[\kappa{\bm{y}}(t)\right]^{\frac{\kappa}{1-\kappa}}-\beta{\bm{x}}(t),$
(28) $\displaystyle\dot{{\bm{y}}}(t)$
$\displaystyle=-\gamma+(r+\beta){\bm{y}}(t)+{\bm{y}}(t)^{\frac{1}{1-\kappa}}\bigg{(}\kappa\big{[}1-{\bm{x}}(t)\big{]}\bigg{)}^{\frac{\kappa}{1-\kappa}}$
(29)

for a given initial condition ${\bm{x}}(0)={\bm{x}}_{0}$, and a transversality
condition

$\displaystyle\lim_{t\rightarrow\infty}e^{-rt}{\bm{y}}(t)=0,$ (30)

Here, ${\bm{x}}$ is the market share of the company, and ${\bm{y}}$ is the
costate variable.

The parameter $\kappa$ is a constant between $0$ and $1$, $\beta$ is strictly
positive, $r$ is the discount rate, the constant $\gamma$ is defined as
$\gamma\equiv\frac{\beta+r}{c}$, and $c$ is the cost of advertising. See [47,
40] for a detailed treatment of this problem.

Similar to the linear asset pricing model, the optimal advertising problem has
a unique steady state, which is a saddle path. Consequently, any sub-optimal
path diverges. However, unlike the linear asset pricing model, this is a two-
dimensional problem, which makes finding the solution a more challenging task.

Using
${\hat{\bm{x}}}(t)={\bm{x}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{x}}}}(\tau)d\tau$
and
${\hat{\bm{y}}}(t)={\hat{\bm{y}}}_{0}+\operatorname{{\textstyle\int}}_{0}^{t}{\hat{\dot{\bm{y}}}}(\tau)d\tau$
we solve the following optimization

$\displaystyle\min_{\hat{{\bm{x}}},\hat{{\bm{y}}},{\hat{\bm{y}}}_{0}}$
$\displaystyle\,\|\dot{{\hat{\bm{x}}}}\|^{2}_{\mathcal{H}}+\|\dot{{\hat{\bm{y}}}}\|^{2}_{\mathcal{H}}$
s.t.
$\displaystyle\dot{{\hat{\bm{x}}}}(t_{i})=\left[1-\hat{{\bm{x}}}(t_{i})\right]^{\frac{1}{1-\kappa}}\left[\kappa\hat{{\bm{y}}}(t_{i})\right]^{\frac{\kappa}{1-\kappa}}-\beta\hat{{\bm{x}}}(t_{i}),\quad\text{for
all }t_{i}\in\mathcal{D},$
$\displaystyle\dot{{\hat{\bm{y}}}}(t_{i})=-\gamma+(r+\beta)\hat{{\bm{y}}}(t_{i})+\hat{{\bm{y}}}(t_{i})^{\frac{1}{1-\kappa}}\bigg{(}\kappa\left[1-\hat{{\bm{x}}}(t_{i})\right]\bigg{)}^{\frac{\kappa}{1-\kappa}},\quad\text{for
all }t_{i}\in\mathcal{D},$

to obtain approximate solutions for ${\bm{x}}(t)$ and ${\bm{y}}(t)$. Here,
$\|\dot{{\hat{\bm{x}}}}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{x\top}K\bm{\alpha}^{x}$,
and
$\|\dot{{\hat{\bm{y}}}}\|^{2}_{\mathcal{H}}=\bm{\alpha}^{y\top}K\bm{\alpha}^{y}$.

Figure 8: The Optimal Advertising model (i.e., Eqs. 28, 28 and 29) solved
using the Matérn kernel for the training data
$\mathcal{D}=\\{0,1,\cdots,30\\}$. The left panel shows the approximate market
share and the right panel shows the approximate costate variable. The
approximate solutions demonstrate convergence to the correct set of steady
states. Therefore, they satisfy the transversality condition.

In this example we use $x_{0}=0.4$, $r=0.11$, $c=0.5$, $\beta=0.05$,
$\kappa=0.5$, and $\gamma=\frac{\beta+r}{c}=0.32$ as the numerical values for
the parameters. We use $\mathcal{D}=\\{0,1,\cdots,30\\}$ as the training data.

Fig. 8 shows the approximate market share, denoted by $\hat{x}(t)$ and the
approximate costate variable, denoted by $\hat{y}(t)$. The vertical dashed
line shows the boundary between the interpolation and extrapolation regions.

This result establishes that the approximate solutions demonstrate convergence
to the correct set of steady states. Therefore, they satisfy the
transversality condition.

#### The correct set of steady states.

How do we know the approximate solutions converge to the correct set of steady
states? As shown in Fig. 8, the approximate costate variable approaches a
finite number. Therefore,
$\lim_{t\rightarrow\infty}e^{-rt}\hat{{\bm{y}}}(t)=0$.