Economics-Inspired Neural Networks with Stabilizing Homotopies
Contemporary deep learning based solution methods used to compute approximate
equilibria of high-dimensional dynamic stochastic economic models are often
faced with two pain points. The first problem is that the loss function
typically encodes a diverse set of equilibrium conditions, such as market
clearing and households' or firms' optimality conditions. Hence the training
algorithm trades off errors between those -- potentially very different --
equilibrium conditions. This renders the interpretation of the remaining errors
challenging. The second problem is that portfolio choice in models with
multiple assets is only pinned down for low errors in the corresponding
equilibrium conditions. In the beginning of training, this can lead to
fluctuating policies for different assets, which hampers the training process.
To alleviate these issues, we propose two complementary innovations. First, we
introduce Market Clearing Layers, a neural network architecture that
automatically enforces all the market clearing conditions and borrowing
constraints in the economy. Encoding economic constraints into the neural
network architecture reduces the number of terms in the loss function and
enhances the interpretability of the remaining equilibrium errors. Furthermore,
we present a homotopy algorithm for solving portfolio choice problems with
multiple assets, which ameliorates numerical instabilities arising in the
context of deep learning. To illustrate our method we solve an overlapping
generations model with two permanent risk aversion types, three distinct
assets, and aggregate shocks.
Marlon Azinovic, Jan Žemlička
# Economics-Inspired Neural Networks with Stabilizing Homotopies+

###### Abstract

Contemporary deep learning based solution methods used to compute approximate equilibria of high-dimensional dynamic stochastic economic models are often faced with two pain points. The first problem is that the loss function typically encodes a diverse set of equilibrium conditions, such as market clearing and households' or firms' optimality conditions. Hence the training algorithm trades off errors between those -- potentially very different -- equilibrium conditions. This renders the interpretation of the remaining errors challenging. The second problem is that portfolio choice in models with multiple assets is only pinned down for low errors in the corresponding equilibrium conditions. In the beginning of training, this can lead to fluctuating policies for different assets, which hampers the training process. To alleviate these issues, we propose two complementary innovations. First, we introduce _Market Clearing Layers_, a neural network architecture that automatically enforces all the market clearing conditions and borrowing constraints in the economy. Encoding economic constraints into the neural network architecture reduces the number of terms in the loss function and enhances the interpretability of the remaining equilibrium errors. Furthermore, we present a homotopy algorithm for solving portfolio choice problems with multiple assets, which ameliorates numerical instabilities arising in the context of deep learning. To illustrate our method we solve an overlapping generations model with two permanent risk aversion types, three distinct assets, and aggregate shocks.

_JEL classification_: C61, C63, C68, E21.

_Keywords_: computational economics, deep learning, deep neural networks, implicit layers, market clearing, global solution method, life-cycle, occasionally binding constraints, overlapping generations.
Introduction

A rapidly growing literature (see, _e.g._Duarte, 2018; Valaitis and Villa, 2021; Azinovic et al., 2022; Ebrahimi Kahou et al., 2021; Han et al., 2021; Maliar et al., 2021; Kase et al., 2022; Duarte et al., 2021; Folini et al., 2021; Bretscher et al., 2022) uses deep neural networks to approximate price, policy and value functions in economic models. To find good approximations, the learning algorithm often aims directly at minimizing the errors in the equilibrium conditions, such as Euler equations, market clearing conditions, and budget constraints, on simulated paths of the economy.1 Since the equilibrium conditions will not be fulfilled exactly, the researcher has to, implicitly, trade-off errors between different equilibrium conditions. This can be challenging, especially for equilibrium conditions of different types and units.

Footnote 1: Examples of deep learning algorithm which directly minimize errors in the equilibrium conditions are Azinovic et al. (2022); Maliar et al. (2021); Ebrahimi Kahou et al. (2021); Kase et al. (2022).

We propose a new architecture for deep neural networks, which we call market clearing layers, that, by construction, enforces exact market clearing conditions. Encoding known properties, coming from the economic model, directly into the architecture has the advantage that it reduces the search space in the space of neural network approximators, to only the subset consistent with market clearing. Hence the neural network does not need to learn a property, we know ex-ante. Furthermore, the number of different error terms in the loss function is reduced, and hence a given loss value is easier to interpret. Next to a simple and easily implementable way to ensure market clearing, we also introduce an implicit layer, which additionally ensures that borrowing constraints are always satisfied. To do so, the last layer of the network solves a quadratic program to obtain adjustments to the layer inputs, with computationally desirable properties. We are not the first to leverage economic insights for the design of neural network architectures. Ebrahimi Kahou et al. (2021) and Han et al. (2021) encode a symmetry property into the neural network architecture. Han et al. (2021) and Azinovic et al. (2022) use suitable activation functions in the output layers to ensure borrowing constraints are satisfied.

Additionally, we propose a homotopy method to stabilize the training progress of deep learning based solution methods for economic models with multiple assets. Homotopy methods are a class of methods, which start by solving a simple problem and then gradually transform the simple problem into the harder problem of interest. In our case, the simple problem to start from is an economic model with a single asset. The hard problem of interest is a model with multiple assets. Multiple assets pose a problem for deep learning based solution methods since a precise solution is required for agents to correctly distinguish between multiple assets. Especially during the beginning of training, when the policies are not yet precisely approximated, this can lead to erratic policies with spurious bang-bang allocations to different assets, which may lead to a failure of the training procedure. To address this, we propose to firstly solve a single-asset model, by setting the supply and short-sale constraints of the remaining assets to zero. Once the single asset model is solved, we keep training the model and slowly increase the supply of the second asset. We proceed this way until the full model is solved. The idea to use homotopy methods to solve for equilibria in economic models has been used, for example, in Schmedders (1998), where a homotopy is used to gradually lift bounds on asset sales. The idea to introduce small amounts of an additional asset is reminiscent of the idea in Judd and Guu (2000), who first solve a model with deterministic asset returns and then uses perturbation methods to introduce small amounts of uncertainty into the asset returns.

The two contributions of this paper go hand in hand. Since the homotopy method we are proposing, slowly and one after the other, increases the supply of additional assets in the economy, it naturally
leverages the market clearing neural network architectures.

To the best of our knowledge, we are the first to study market clearing architectures for neural networks and to lay out a homotopy algorithm to step-wise introduce multiple assets into a deep learning based solution method.

## 2 Market Clearing Layers

To introduce market clearing layers, we first consider a neural network approximating policy and price functions for a single asset and a finite number of agents, abstracting from the rest of the economy, such as the households' objective functions. In section 2.5 we show how market clearing layers can be applied in a concrete economic model.

Consider an economic model, with a set of finitely many agents, indexed with \(i\in\mathcal{I}:=\{1,\ldots,H\}\) and associated population weights \(\mu_{i}\). In an overlapping generations model, for example, \(i\) could index a specific age group and \(\mu_{i}\) could denote its size.2 Let \(\mathbf{x}\in\mathbb{R}^{N}\) denote the state of the economy. Assume that agents can save in a single asset \(b\) with price \(p^{b}\) and aggregate supply \(B\). Let \(b^{i}(\mathbf{x})\) denote the policy function of agent \(i\). Market clearing requires that the total asset demand by households adds up to the total asset supply:

Footnote 2: In an economy with a continuum of agents, \(i\) could index a specific combination of shocks and asset holdings, and \(\mu_{i}\) would denote the mass of agents at a specific histogram point.

\[\sum_{i=1}^{H}\mu_{i}b^{i}(\mathbf{x})=B. \tag{1}\]

We are interested in using a neural network to approximate the households' asset policies. Let

\[\mathcal{N}_{\rho}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{H}\times\mathbb{R}_{> 0},\mathcal{N}_{\rho}(\mathbf{x})=\begin{pmatrix}\tilde{b}^{1}(\mathbf{x})\\ \tilde{b}^{2}(\mathbf{x})\\ \ldots\\ \tilde{b}^{H}(\mathbf{x})\\ \hat{p}^{b}(\mathbf{x})\end{pmatrix} \tag{2}\]

denote a neural network with parameters \(\rho\) that maps the state of the economy to an approximation of the household savings choices and the price of the asset. The goal of this paper is to devise a neural network architecture that, by design, ensures that market clearing is satisfied while also retaining desirable computational properties.

### Simple Adjustment for Pure Market Clearing

Let the aggregate asset demand, which is implied by the households' policies \(\tilde{b}^{i}(\mathbf{x})\), be denoted with

\[B^{H}(\mathbf{x}):=\left(\sum_{i=1}^{H}\mu_{i}\tilde{b}^{i}(\mathbf{x})\right). \tag{3}\]
The corresponding excess demand is given by

\[\Delta B(\mathbf{x})=B^{H}(\mathbf{x})-B. \tag{4}\]

We would like to adjust the policies \(\tilde{b}^{i}(\mathbf{x})\), such that they are consistent with market clearing, _i.e._ such that the excess demand is zero. Let \(b^{i}(\mathbf{x})\) denote the adjusted policy functions. There are multiple ways to adjust the policies \(\tilde{b}^{i}\) to make up for the excess demand \(\Delta B(\mathbf{x})\). For example, one way would be to rescale the policies to obtain new policies

\[\forall i\in\mathcal{I}:\ b^{i}(\mathbf{x})=\tilde{b}^{i}(\mathbf{x})\frac{B} {B^{H}(\mathbf{x})}. \tag{5}\]

However, in general this algorithm would not be convenient, for example in settings when \(B^{H}(\mathbf{x})=0\). Another way would be to adjust the policy of only a single household, for example

\[j\in\mathcal{I}:b^{j}(\mathbf{x}) =\tilde{b}^{j}(\mathbf{x})-\frac{\Delta B(\mathbf{x})}{\mu^{j}} \tag{6}\] \[\forall i\in\mathcal{I},i\neq j:b^{i}(\mathbf{x}) =\tilde{b}^{i}(\mathbf{x}). \tag{7}\]

This way would load all the excess demand onto a single agent and cause an asymmetry in the accuracy of the policies across different households.

We propose the solution to the following problem as a desirable adjustment mechanism:

\[\{b^{i}(\mathbf{x})\}_{i\in\mathcal{I}} =\arg\min\sum_{i\in\mathcal{I}}\frac{1}{2}\mu_{i}(b^{i}(\mathbf{x })-\tilde{b}^{i}(\mathbf{x}))^{2}\] \[\text{subject to}: \sum_{i\in\mathcal{I}}\mu_{i}b^{i}(\mathbf{x}) =B. \tag{8}\]

The problem has the simple solution that all households policy is adjusted by an equal amount

\[\forall i\in\mathcal{I}:b^{i}(\mathbf{x})=\tilde{b}^{i}(\mathbf{x})-\Delta B (\mathbf{x}). \tag{9}\]

Next to avoiding a distortion in the symmetry between the Euler equations across households this formulation has the advantage that it also works for assets in zero net supply as well as for assets where the net supply is state-dependent.3

Footnote 3: The latter is true for the previously mentioned to formulations as well.

Despite its computational simplicity, a disadvantage of this approach is that, while the adjusted policies are consistent with market clearing, they are not necessarily consistent with borrowing constraints. A softplus activation function, for example, could ensure that \(\tilde{b}^{i}(\mathbf{x})\) are non-negative, however this would not imply that the adjusted, market clearing, \(b^{i}(\mathbf{x})\) are non-negative as well. Ideally, we would like the adjustment to be able to ensure that both, market clearing and borrowing constraints, are always satisfied. We address this point in the next section.
### Implicit Layer to Encode Market Clearing and Borrowing Constraints

The market clearing adjustment described in the previous section does not ensure that borrowing constraints are satisfied. One way to deal with this would be to separately predict Karush-Kuhn-Tucker (KKT) multipliers and to include the KKT conditions into the loss function. An alternative would be to use the Fischer-Burmeister equation (see Jiang, 1999), as for in example Maliar et al. (2021), to encode the Euler equation error and the violation of the borrowing constraint into a single error term. We propose a modification to the market clearing mechanism described in the previously in section 2.1, which simultaneously ensures market clearing and that the borrowing constraints are satisfied, by adding the borrowing constraint to the above minimization problem:

\[\{b^{i}(\mathbf{x})\}_{i\in\mathcal{I}} =\arg\min\sum_{i\in\mathcal{I}}\frac{1}{2}\mu_{i}(b^{i}(\mathbf{x })-\tilde{b}^{i}(\mathbf{x}))^{2}\] \[\text{subject to}:\sum_{i\in\mathcal{I}}\mu_{i}b^{i}(\mathbf{x}) =B\] \[b^{i}(\mathbf{x}) \geq\underline{b}^{i}(\mathbf{x}) \tag{10}\]

where \(\underline{b}^{i}(\mathbf{x})\) denotes the borrowing limit of agent \(i\). The solution to this problem can be obtained with solvers for box-constrained Quadratic Programs, which are meanwhile implemented in modern deep learning libraries, such as JAX, and provide efficiently computed derivatives for the use in backpropagation algorithms.4 The draw-back of these algorithms is that they substantially slow down the training process and are often programmed for a general class of Quadratic Programs and do not exploit symmetries, which are often present in economic models.

Footnote 4: We used the BoxQSQP solver from the jaxopt library: [https://jaxopt.github.io/stable/quadratic_programming.html](https://jaxopt.github.io/stable/quadratic_programming.html).

Consider, for example, a borrowing constraint at zero for all agents: \(\forall i:\ b^{i}(\mathbf{x})\geq 0\). In this case we can ensure that \(\tilde{b}^{i}\geq 0\) by using a softplus activation function in the output layer of the neural network. If the excess demand is negative, _i.e._\(\Delta B(\mathbf{x})<0\), the solution to the constrained quadratic program coincides with the simple solution described in section 2.1. The adjusted asset demand of all agents is given by \(b^{i}(\mathbf{x})=\tilde{b}^{i}(\mathbf{x})-\Delta B(\mathbf{x})\) and no constrained is violated. When the excess demand is positive, all asset choices \(\tilde{b}^{i}\) below a threshold value \(\tilde{b}^{\text{threshold}}\) are adjusted to lie on the borrowing constraint and all remaining asset holdings are adjusted by the same amount such that markets clear. Let \(\mathcal{J}\) denote the set of households \(j\) with \(\tilde{b}^{j}<\tilde{b}^{\text{threshold}}\). The solution to the constrained quadratic program is given by

\[b^{i}=\begin{cases}0&\text{for }i\in\mathcal{J}\\ \tilde{b}^{i}-\underbrace{(\Delta B(\mathbf{x})-\sum_{j\in\mathcal{J}}\mu_{j}( \tilde{b}^{j}-0))}_{\text{remaining excess demand}}&\text{for }i\in\mathcal{I} \setminus\mathcal{J}\end{cases} \tag{11}\]

Solving the constrained quadratic program hence simplifies to finding a single number \(\tilde{b}^{\text{threshold}}\). Since \(\tilde{b}^{\text{threshold}}\in[\min_{i}\tilde{b}^{i},\max_{i}\tilde{b}^{i}+\epsilon]\), for \(\epsilon>0\), it can be efficiently computed with a few bisection steps. For our examples, this turned out to be faster than using the general solver for constrained quadratic programs.
#### 2.2.1 Policy updating for falsely constrained agents

One problem, when ensuring the borrowing constraint as described above, is that for constrained households, _i.e._ households with \(\tilde{b}^{i}<\tilde{b}^{\text{threshold}}\) and \(b^{i}=0\), an infinitesimal change in the neural network output \(\tilde{b}^{i}\) does not influence the adjusted policy \(b^{i}\). For a loss function purely based on the post-adjustment prediction \(b^{i}\), this is a problem if an agent is predicted to be constrained, though it should not be. The derivative of the loss function with respect to the pre-adjustment prediction \(\tilde{b}^{i}\) would be zero and the gradient would hence not help to improve the corresponding parameters in the neural network.5 To address this issue, we provide a way to propagate the signal from the loss function to the neural network, so that the gradient descent step will lead to improved neural network parameters, which increase the prediction for \(\tilde{b}^{i}\) in such cases. For clarity, we detail it in context of the simple example below in section 2.5. While our solution allows us to retain the signal while enforcing the borrowing constraints by design, we admit that the solution is based on additional terms in the loss function, which is unfortunate given the very goal of this paper. Another solution, of course, is to follow the simple adjustment described above in section 2.1. Using reformulations like the Fischer-Burmeister equation (see Jiang, 1999),6 or the Garcia Zangwill-trick (see Zangwill and Garcia, 1981; Judd et al., 2000), the KKT conditions can be summarized by a single equation and no additional terms in the loss function are required.

Footnote 5: A similar problem could arise in models where the last layer has a softplus activation when the pre-activated value is very negative. While the softplus activation would still guarantee a positive derivative, the derivative could be vanishingly small.

Footnote 6: The Fischer-Burmeister function is \(\psi(x,y):=x+y-\sqrt{x^{2}+y^{2}}\). While the Fischer-Burmeister function is differentiable, its roots satisfy \(x\geq 0\), \(y\geq 0\) and \(xy=0\).

### Choosing Between the Two Market-Clearing Algorithms

In sections 2.1 and 2.2 we laid out two ways to ensure that market clearing is satisfied exactly. The simple algorithm described in section 2.1 has the advantage that we can obtain the solution to problem (8) in closed form, as given in equation (9). The market clearing policies can hence be computed with minimal computational complexity. The disadvantage is that the adjusted policies may violate the borrowing constraints. Consequently, the borrowing constraints must be encoded in the loss function, for example by using the Fischer-Burmeister equation (see Jiang, 1999).

The solver-based adjustment described in section 2.2 has the advantage that the solution to the problem (10) ensures that, additionally to market clearing, all borrowing constraints are always satisfied. The disadvantage is that a closed form solution for constrained quadratic programs is not available and we therefore need to invoke a solver, which slows down the training. Furthermore, the solver-based adjustment faces the problem that we need to add auxiliary terms to the loss-function in order to ensure that the neural network weights are updated if an agent is falsely predicted to be constrained, as described in section 2.2.1.

Which of the two market clearing algorithm is better may depend on the model at hand and on whether a strict enforcement of the borrowing constraints adds stability to the training algorithm. When the borrowing constraints can be conveniently encoded using the Fisher-Burmeister equation, we found that the advantage of ensuring the borrowing constraints by design is not large enough to make up for the loss in speed resulting from the invocation of a quadratic solver each time the neural network is evaluated.
### Learning algorithm

We use a version of the Deep Equilibrium Nets method developed by Azinovic et al. (2022). Following the Deep Equilibrium Nets algorithm, we proceed in three steps. First, we parameterize the equilibrium functions of the economy (i.e. value, policies, prices, etc) using a deep neural network. Second, we write down a loss function, which is defined as a mean squared residual of model equations, which characterize the equilibrium functions, over a batch of states. Third, we update network weights by minimizing the loss function along simulated paths of the economy.

Because the loss function includes all the equilibrium conditions of the economy, achieving a sufficiently low level of loss over some region of the state space indicates that the neural network learned how to construct an approximate model solution over that region. By minimizing loss function over states obtained from a long simulation of the economy we make sure that the neural network learns to approximate the model solution on the approximate ergodic set of the economy. As discussed in Judd et al. (2011), drawing states from the approximate ergodic set instead of sampling from some of its enclosing hypercubes allows for a drastic reduction of the computational complexity of obtaining an approximate solution.7

Footnote 7: Sparse grids (see Krueger and Kubler, 2004) and adaptive sparse grids (see Brumm and Scheidegger, 2017), for example, require a hyper-cubic domain.

Relative to the original Deep Equilibrium Nets algorithm, our algorithm differs along two dimensions. First, our neural networks feature Market Clearing Layers, hence their outputs automatically satisfy market clearing conditions. Therefore, market clearing conditions do not have to be included in the loss function. Moreover, the states we obtain from the simulation are always consistent with market clearing. Second, we use a modified simulation algorithm. Instead of simulating several paths of the economy for many periods before retraining the neural network, we simulate a large number of states for only one period forward before updating the network weights by training on those newly simulated states.8

Footnote 8: Results shown by Zemlicka (2023) suggest, that this approach improves the stability of the training process.

Specifically, we simulate a number of parallel trajectories for one period forward. Let \(N^{\text{trajectories}}\) denote the number of parallel trajectories. To simulate the model forward, we use the policies as approximated by the neural network and draw random exogenous shocks according to their transition probabilities. In each simulation step, we hence obtain \(N^{\text{trajectories}}\) new states. These newly simulated state constitute our dataset to update the approximated policies by training the neural network on the dataset to minimize the loss function. Let \(N^{\text{epochs}}\) denote the number of epochs we train on each simulated dataset, and let \(N^{\text{mini-batch}}\) denotes the mini-batch size. An epoch refers to a single training pass through the entire dataset.9 The mini-batch size refers to the number of data-point used for each step of stochastic gradient descent.10 Hence, each simulated dataset of \(N^{\text{trajectories}}\) new states is used for \(\frac{N^{\text{trajectories}}N^{\text{epochs}}}{N^{\text{mini-batch}}}\) steps of stochastic gradient descent. After training the neural network, we simulate a new dataset. We refer to one simulation of new training data, together with subsequent training of the neural network on that data, as one episode. Let \(N^{\text{episodes}}\) denote the total number of episodes used to train the neural network. Our simulation procedure starts out with \(N^{\text{trajectories}}\) copies of a single feasible state of the economy. As the policies improve and as we simulate the model forward, the set of states, generated by simulation, converges toward an approximation of the ergodic set of states of the economy.

Footnote 10: See Goodfellow et al. (2016) for a general introduction to deep learning.

Footnote 10: We use the Adam optimizer (Kingma and Ba, 2014), which is a variant of stochastic gradiant descent and the de-facto standard in the deep learning literature.

In order to evaluate the accuracy of our approximate solution, we simulate the model forward without retraining the neural network parameters and evaluate the equilibrium conditions on the newly generated
set of states.

### Simple Example

We illustrate the market clearing mechanism with an overlapping generations economy, with a single risk free asset and borrowing constraints.

#### 2.5.1 Model

Time is discrete. Every period, a shock \(z_{t}\) realizes, which determines the aggregate labor income. We assume that \(\log(z_{t})\) follows an AR(1) process

\[\log(z_{t})=\rho\log(z_{t-1})+\sigma\epsilon_{t}, \tag{12}\]

where \(\rho\) denotes the persistence and \(\sigma\) the volatility, and where the innovation \(\epsilon_{t}\) is drawn form a standard normal distribution.

Every period a new cohort of measure one enters the economy and lives deterministically for \(H<\infty\) periods. Households of age \(h\in\{1\dots H\}\) receive an exogenous income \(y_{t}^{h}=z_{t}y^{h}\) which is given by an age-specific share \(y^{h}\) of aggregate labor income \(z_{t}\). Households rent housing services with fixed supply \(H^{r}\), consume and can trade a risk free one-period bond with fixed supply \(B\) subject to borrowing constraints. Let \(b_{t}^{h}\) denote the bond holdings of age group \(h\) in the beginning of period \(t\) and let \(\underline{b}\) denote the exogenous borrowing constraint: \(b_{t+1}^{h+1}\geq\underline{b}\). We assume that households are born and die without assets, \(b_{t}^{1}=b_{t}^{H+1}=0\). To adjust their bond holdings, households have to pay an adjustment cost \(\zeta^{b}\frac{1}{2}(b_{t+1}^{h+1}-b_{t}^{h})^{2}\). Households of age \(h\) maximize their remaining time-separable discounted expected lifetime utility given by

\[\sum_{i=0}^{H-h}\mathrm{E}_{t}\left[\beta^{i}\left(u(c_{t+i}^{h+i} )+\psi^{h}v(h_{t+i}^{r,h+i})\right)\right], \tag{13}\]

where \(c_{t}^{h}\) denotes the consumption of age group \(h\) in period t, \(h_{t}^{r,h}\) denotes the consumption of housing services of age group \(h\) in period \(t\), and \(\beta\) denotes the patience. We assume constant relative risk aversion utility for both types of consumption with coefficient of relative risk aversion \(\gamma\), \(u(c)=\frac{c^{1-\gamma}}{1-\gamma}\) and \(v(h^{r})=\frac{(h+h^{r})^{1-\gamma}}{1-\gamma}\). The weight for housing services, \(\psi^{h}\), is age dependent. Households budget constrained is given by

\[c_{t}^{h}=z_{t}y^{h}+b_{t}^{h}-p_{t}^{b}b_{t+1}^{h+1}-p_{t}^{b} \zeta^{b}\frac{1}{2}(b_{t+1}^{h+1}-b_{t}^{h})^{2}-p_{t}^{r}h^{r}h_{t+1}^{h+1}, \tag{14}\]

where \(p_{t}^{b}\) denotes the price of the bond and \(p_{t}^{r}\) denotes the price for renting housing services. The Bellman equation corresponding to the households' problem is given by

\[V^{h}(z_{t},b_{t}^{h}) =\max_{b_{t+1}^{h+1},h_{t}^{r,h}}\left\{u(c_{t}^{h})+\psi^{h}v(h_ {t}^{r,h})+\beta\mathrm{E}\left[V^{h+1}(z_{t+1},b_{t+1}^{h+1}\right]\right\}\] (15) subject to : \[c_{t}^{h} =z_{t}y^{h}+b_{t}^{h}-b_{t+1}^{h+1}p_{t}^{b}-h_{t}^{r,h}p_{t}^{r} -\zeta^{b}\frac{1}{2}\left(b_{t+1}^{h+1}-b_{t}^{h}\right)^{2} \tag{16}\] \[0 \leq b_{t+1}^{h+1}-\underline{b} \tag{17}\]
The Karush-Kuhn-Tucker conditions for the savings choice of and households' savings choices are given by

\[(p_{t}^{b}+p_{t}^{b}\zeta^{b}(b_{t+1}^{h+1}-b_{t}^{h}))u^{\prime}(c_{t }^{h}) =\beta\mathrm{E}_{t}\left[u^{\prime}(c_{t+1}^{h+1})(1+p_{t+1}^{b} \zeta^{b}(b_{t+2}^{h+2}-b_{t+1}^{h+1}))\right]+\lambda_{t}^{h} \tag{18}\] \[(b_{t+1}^{h+1}-\underline{b}) \geq 0\] (19) \[\lambda_{t}^{h} \geq 0\] (20) \[\lambda_{t}^{h}(b_{t+1}^{h+1}-\underline{b}) =0. \tag{21}\]

Similarly, optimality conditions for the intratemporal renting choice are given by

\[p_{t}^{r}u^{\prime}(c_{t}^{h})=\psi^{h}v^{\prime}(h_{t}^{r,h}). \tag{22}\]

We provide details on parameter choices in the appendix A.

Functional Rational Expectations EquilibriumLet \(\mathbf{b}:=[b^{1},\ldots,b^{H}]\in\mathbb{R}^{H}\) denote the distribution of asset holdings across age groups. Let \(\mathbf{x}:=[z,\mathbf{b}]\in\mathbb{R}^{H+1}\) denote the state of the economy. Following Spear (1988) and Krueger and Kubler (2004) a functional rational expectations equilibrium is given by a set of \(H-1+H\) policy functions, which map the state of the economy the bond purchases and the consumption of rental services of each age group, together with two price functions, which map the state of the economy to the prices of the bond and rental services, such that the policies and prices are consistent with market clearing as well as the households' optimality conditions.

Approximating the equilibrium functions with deep neural networksFollowing Azinovic et al. (2022), we aim to solve the model by approximating the equilibrium functions with a deep neural network. Let \(\mathcal{N}_{\rho}:\mathbb{R}^{H+1}\rightarrow\mathbb{R}^{2H+1}\) denote a neural network we use to approximate the price and policy functions, such that

\[\mathcal{N}_{\rho}(\mathbf{x}_{t})=\begin{bmatrix}\tilde{b}_{t}^{1}\\ \ldots\\ \tilde{b}_{t}^{H-1}\\ \tilde{b}_{t}^{r,1}\\ \ldots\\ \tilde{h}_{t}^{r,H-1}\\ \hat{p}_{t}^{b}\\ \hat{p}_{t}^{r}\end{bmatrix} \tag{23}\]

In difference to Azinovic et al. (2022) we now apply the market clearing layers described in section 2 to obtain approximations for the households' policies \(\hat{b}_{t}^{1},\ldots,\hat{b}_{t}^{H-1}\) and \(\hat{h}_{t}^{r,1},\ldots,\hat{h}_{t}^{r,H-1}\), which are already consistent with market clearing. We illustrate both ways, the simple way described in section 2.1 and the adjustment described in section 2.2, which additionally ensures that the borrowing constraints are satisfied. Due to the Inada property of the utility from housing services, there is no borrowing constraint for renting services and the two adjustments coincide for renting.

Loss functionSince our algorithm ensures that market clearing holds exactly for all predicted policy functions, the remaining equilibrium conditions are the households' first order conditions. As in Azinovic
et al. (2022) we construct a loss function as the means squared error in the remaining equilibrium conditions and the train the neural network to minimize those errors on simulated paths of the economy. Let \(\hat{c}_{t}^{h}\) denote the consumption of age group \(h\), which the neural network, together with the market clearing layer, predicts for state \(\mathbf{x}\). We use the Fischer-Burmeister function, \(\psi(x,y):=x+y-\sqrt{x^{2}+y^{2}}\) to encode the KKT conditions for each age group into a single equation. Following Judd (1998) we rewrite the equations, such that deviations from fullfilling them exactly are interpretable as relative consumption errors. The resulting errors in the equilibrium conditions for age group \(h\) are given by

\[\mathrm{err}_{\rho}^{c,h}(\mathbf{x}_{t}) =\psi\left(\frac{(u^{\prime})^{-1}\left(\frac{\beta\mathrm{E}_{ \mathrm{g}}\left[u^{\prime}(\hat{c}_{t+1}^{h+1})(1+\hat{p}_{t+1}^{h}\hat{c}^{b} (\hat{b}_{t+2}^{h+2}-\hat{b}_{t+1}^{h+1}))\right]}{\hat{p}_{t}^{h}+\hat{p}_{t} ^{h}\hat{c}_{t}^{h}(\hat{b}_{t+1}^{h+1}-\hat{b}_{t}^{h})}\right)}{\hat{c}_{t}^ {h}}-1,\frac{\hat{b}_{t+1}^{h+1}-\hat{b}}{\hat{c}_{t}^{h}}\right) \tag{24}\] \[\mathrm{err}_{\rho}^{r,h}(\mathbf{x}_{t}) =\frac{(u^{\prime})^{-1}\left(\frac{\psi^{h}v^{\prime}(\hat{p}_{ t}^{r,*})}{\hat{p}_{t}^{r}}\right)}{\hat{c}_{t}^{h}}-1. \tag{25}\]

Equation (24) encodes the errors in the optimality conditions for the choice of saving in the bond, and equation (25) encodes the errors in the optimality conditions for the choice of renting housing services. Let \(\mathcal{D}\) denote a set \(\left|\mathcal{D}\right|\) states. Our loss function, which encodes all remaining equilibrium conditions, is given by

\[l_{\rho}(\mathcal{D}):=\frac{1}{\left|\mathcal{D}\right|}\sum_{ \mathbf{x}\in\mathcal{D}}\left(\frac{1}{H-1}\sum_{h=1}^{H-1}\left(\mathrm{err }_{\rho}^{c,h}(\mathbf{x})\right)^{2}+\frac{1}{H}\sum_{h=1}^{H}\left(\mathrm{ err}_{\rho}^{r,h}(\mathbf{x})\right)^{2}\right). \tag{26}\]

Updating the policies of falsely constrained householdsWhen using our solver based method, there would be no feedback from the loss function to a neural network prediction \(\tilde{b}^{h}\) when the solver sets the household onto the borrowing constraint. This is problematic when the household is predicted to be constrained, but should not be. In such cases, \(\frac{\hat{b}_{t+1}^{h+1}-\hat{b}}{\hat{c}_{t}^{h}}=0\), while \(\frac{(u^{\prime})^{-1}\left(\frac{\beta\mathrm{E}_{\mathrm{g}}\left[u^{\prime }(\hat{c}_{t+1}^{h+1})(1+\hat{c}_{t}^{b}(\hat{b}_{t+2}^{h+2}-\hat{b}_{t+1}^{h+1 }))\right]}{\hat{p}_{t}^{h}+\hat{c}_{t}^{h}(\hat{b}_{t+1}^{h+1}-\hat{b}_{t}^{h })}\right)}{\hat{c}_{t}^{h}}-1<0\). This would correctly lead to a positive term in the loss function, \(\left(\mathrm{err}_{\rho}^{c,h}(\mathbf{x}_{t})\right)^{2}>0\), but to no feedback to the neural network parameters to increase the prediction \(\tilde{b}^{h}\). To restore this feedback, we add the square of the following term to the loss function

\[\mathrm{err}_{\rho}^{h,\text{false binding}}(\mathbf{x}_{t}):=\frac{1}{1+ \tilde{b}_{t+1}^{h+1}}\times e^{-\frac{\left(\hat{b}_{t+1}^{h+1}-\hat{b} \right)^{2}}{10^{-5}}}\times\max\left\{-\frac{(u^{\prime})^{-1}\left(\frac{ \beta\mathrm{E}_{\mathrm{g}}\left[u^{\prime}(\hat{c}_{t+1}^{h+1})(1+\hat{c}_{t }^{b}(\hat{b}_{t+2}^{h+2}-\hat{b}_{t+1}^{h+1}))\right]}{\hat{p}_{t}^{h}+\hat{c }_{t}^{b}(\hat{b}_{t+1}^{h+1}-\hat{b}_{t}^{h})}\right)}{\hat{c}_{t}^{h}}+1,0 \right\}. \tag{27}\]

The intuition behind these terms is as follows. The last term, \(\max\left\{-\frac{(u^{\prime})^{-1}\left(\frac{\beta\mathrm{E}_{\mathrm{g}} \left[u^{\prime}(\hat{c}_{t+1}^{h+1})(1+\hat{c}_{t}^{b}(\hat{b}_{t+2}^{h+2}- \hat{b}_{t+1}^{h+1}))\right]}{\hat{p}_{t}+\hat{c}_{t}^{b}(\hat{b}_{t+1}^{h+2}- \hat{b}_{t}^{h+1})}\right)}{\hat{p}_{t}^{h}}+1,0\right\}\), is always non-negative and only different from zero, when the household is currently saving too little, in which case it should not be on the borrowing constraint. The second term, \(e^{-\frac{\left(\hat{b}_{t+1}^{h+1}-\hat{b}\right)^{2}}{10^{-5}}}\), is a differentiable approximation to a function, which is always zero, except when \(\hat{b}_{t+1}^{h+1}-\hat{b}\), in which case it is equal to one. Hence, it is always non-negative and only positive if the household is predicted to be constrained.
Taken together, the last two terms in equation (27) ensure that the term is only positive if the household is falsely predicted to be constrained. Finally, the first term, \(\frac{1}{1+\hat{b}_{t+1}^{n+1}}\), ensures that the overall error term is reduced by increasing \(\hat{b}_{t+1}^{h+1}\) and is thus restoring a pass through from the loss function to the neural network parameters governing the prediction for the household's policy.

IntegrationAn evaluation of the loss function implied by the model requires the computation of conditional expectations in the Euler equations, which characterize the households' behavior. Despite a large number of state variables, this economy features only one source of uncertainty: a stationary AR(1) process with Gaussian innovations. This structure allows us to approximate conditional expectations using standard Gauss-Hermite quadrature (for more details on numerical integration in the context of economic modeling, see Judd, 1998). We denote the number of quadrature nodes with \(N^{\text{integration}}\) and choose \(N^{\text{integration}}=8\).

#### 2.5.2 Training

To parameterize equilibrium objects of the economy (policies and prices), we use a densely connected feed-forward neural network with two hidden layers. The dimensionality of its input is \(N^{\text{input}}=21\) since the state vector includes both, the aggregate productivity shock, and the asset distribution across all twenty age groups in the economy. For the two hidden layers, we choose \(N^{\text{hidden 1}}=N^{\text{hidden 2}}=400\) neurons with relu activation functions.11 Finally, the network output is a vector of \(N^{\text{output}}=41\) elements consisting of bond savings, housing services consumption policy functions, and two price functions for the price of the bond and the price for renting housing services. Outputs corresponding to prices are transformed using the softplus activation in order to ensure that the predicted prices are always positive.12

Footnote 11: The relu activation function is given by \(\text{relu}(x):=\max\{0,x\}\). Relu is short for the rectified linear unit.

Footnote 12: The softplus activation function is given by \(\text{softplus}(x):=\log(1+\exp(x)\).

To train the network to approximate equilibrium price and policy functions, we use the training procedure described in section 2.4. To obtain training states, we simulate \(N^{\text{trajectories}}=8192\) parallel state trajectories for \(N^{\text{episodes}}=3584\) episodes. For each episode (_i.e._ each simulation step), we train the neural network for \(N^{\text{epochs}}=10\) epochs. In each epoch, we split the simulated dataset of \(N^{\text{trajectories}}\) states into batches of \(N^{\text{minibatch}}=128\) states. Each batch is used to compute one stochastic gradient descent update of the neural network weights. To calculate those updates, we use the Adam optimizer (Kingma and Ba, 2014), a version of the stochastic gradient descent algorithm, with a learning rate of \(\alpha^{\text{learn}}=1\times 10^{-5}\).13

Footnote 13: To improve robustness of our learning algorithm, we transformed loss function gradients using the zero_nans transformation before using them as an input for the Adam optimizer. Specifically, we use adam and zero_nans provided by the optax library: [https://optax.readthedocs.io/en/latest/](https://optax.readthedocs.io/en/latest/). The zero_nans transformation replaces NaI,when they occur, with 0, avoiding a failure of the training. Furthermore, we found it useful to re-start internal state of Adam optimizer to zero at the beginning of training episode.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline Parameters & \(N^{\text{input}}\) & \(N^{\text{hidden 1}}\) & \(N^{\text{hidden 2}}\) & \(N^{\text{output}}\) \\  & & Activations & Activations & Activations \\ \hline \multirow{2}{*}{Values} & \multirow{2}{*}{21} & 400 & 400 & 41 \\  & & relu & relu & see text \\ \hline \hline \end{tabular}
\end{table}
Table 1: Network Architecture chosen for the single asset model.
#### 2.5.3 Accuracy

To evaluate the quality of our solution, we evaluate statistics of the errors in the equilibrium conditions along simulated paths of the economy. We simulate a total of \(2^{13}=8192\) new states without retraining the neural network. Figure 1 reports statistics on the remaining errors in the KKT conditions for each age group, expressed as relative consumption errors. As the figure shows, the models are solved precisely, with the 99th percentile of errors in the optimality conditions well below \(1\%\) for each age group and both policies. The left panel in figure 1 shows the simple market clearing architecture described in section 2.1, the right panel shows the solver-based architecture, which additionally guarantees that the borrowing constraints are always satisfied. Indeed we can see that the equilibrium errors are fulfilled exactly for the young agents, for which the borrowing constraint is always binding.

## 3 Homotopy Algorithm for Portfolio Choice

While deep learning based solution methods can handle high-dimensional state spaces comparatively well, portfolio choice problems pose a challenge for this approach. Portfolio choice is pinned down only at low levels of errors in the households' optimality conditions. In the beginning of training the errors are still high and it is therefore hard to learn the portfolio choice correctly. As a consequence, the training process is prone to get stuck at spurious solutions featuring a bang-bang type of portfolio allocations. In general, those oscillations are difficult to overcome.

In simulation-based solution algorithms, this issue is amplified, since the intermediate "solution" is used to simulate the dynamics of the economy to produce new state values, on which the network is subsequently trained. Bad policies might hence generate economically non-sensible states, which are far from the domain on which the neural network was trained previously. Those oscillations can hence create a feedback loop, which might hamper the convergence of the whole algorithm.

To ameliorate the propensity of portfolio choice models to generate zig-zag oscillations, we propose a homotopy procedure. Homotopy methods (see Schmedders, 1998) are a general class of methods. The idea is to start by solving a simple problem and then slowly transform the problem into the harder problem of interest. If the transformations of the problem are gradual enough, the previously obtained solution to the nearby problem provides a good starting point to solve the slightly harder, transformed problem. In order to leverage this concept for solving portfolio choice models, we start by solving a single-asset economy, which is then gradually enriched with new assets. By starting with solving a simple problem up to a high degree of precision, and then creating a sequence of perturbed problems, where every two consecutive problems are very similar to each other, we are able to guide the solution procedure and avoid spurious zig-zag solutions, which may otherwise prevent the solution procedure from converging to a high-quality solution.

For a practical implementation of the homotopy procedure, we utilize a natural nesting structure present

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline Parameters & \(N^{\text{episodes}}\) & \(N^{\text{trajectories}}\) & \(N^{\text{epochs}}\) & \(N^{\text{minibatch}}\) & \(N^{\text{integration}}\) & \(\alpha^{\text{learn}}\) \\ \hline Values & 3584 & 8192 & 10 & 128 & 8 & \(10^{-5}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Hyperparameters chosen for the single asset model.
Figure 1: Statistics on the violation of the equilibrium conditions by age group on a simulated path of 8192 states by age group. The left panel shows the statistics for the simple market clearing architecture and the right panel shows the results for the solver-based architecture, which additionally guarantees that the borrowing constraints are always satisfied. The top row shows the errors in the KKT conditions for the bond and the bottom row shows the errors in the Euler equation for the renting of housing services. The different lines, from the lowest to the highest, show the minimum, the 10th percentile, the average, the 90th percentile and the maximum.
in portfolio choice models. A multi-asset economy can be transformed into a single-asset one by setting the net supply of all but one asset to zero and imposing a strict no short-sale constraint. Given a zero net supply and a strict no short-sale constraint, equilibrium allocations of that asset are uniformly zero. Moreover, the loss function used to train the network can be parameterized with weights on Euler equations of different assets. By setting the weights on Euler equations for some assets to zero and and imposing zero net supply coupled with a strict no short-sale constraint, the loss function constructed to represent the multiple-asset economy characterizes the single-asset economy. In addition to parameters governing asset supply and the weights in the loss function, we also employ a mask parameter that multiplies all the asset policies. We initially set the mask to zero for all assets, which are not yet included into the economy. The mask parameter is increased to small positive values when the corresponding asset is added to the economy.

This nesting structure allows us to implement our homotopy procedure using a simple loop that implements the deep learning procedure for a particular parameterization of the economy and the loss function. When adding the N-th asset, we first use the \(N-1\) asset economy and the N-th asset Euler equations to construct a price function for the N-th asset. Next, we introduce a small amount of the N-th asset and take the policy and the price functions constructed in the previous step as an initial guess. Then we run the deep learning procedure to solve for an equilibrium of this perturbed economy.

Our homotopy procedure can be summarized in the following steps:

* Parameterize all price and policy functions using a neural network equipped with market clearing layers. This network takes as an input a full state vector of the multi-asset economy and outputs all prices and portfolio allocations. Initialize its parameters.
* For all but one asset set the net supply to zero, and impose strict no short-sale constraints on the remaining assets. Set the weights in the loss function on Euler equations associated with the assets in zero net supply to zero. Set the mask parameter for the corresponding policy functions to zero.
* Select initial state values for the training procedure.
* Starting from initial weights and state values, run the training procedure to solve for an equilibrium of the single asset economy.
* Set the weights in the loss function on the Euler equations associated with the second asset to its full value (e.g. 1). Keep the zero net supply and the no short-sale constraint.
* Use the Euler equations associated with the second asset to solve for its price function. Since the weight on this Euler equations is non-zero, the network will be able to learn a price function implied by zero allocations of that asset. Recall that this second asset is in zero net supply, and agents are bound by strict no short-sale when trading it.
* Set the mask parameter for the second asset to a small, but non-zero value (_e.g._ 0.01). Now the network can predict non-zero allocations of the second asset. The choice of a small mask value helps to avoid instability associated with a potentially large and discontinuous jump.
* Set the net supply of the second asset at a small, but non-zero value. Use network weights and state values from the previous step as an initial guess for solving this economy.
* Increase the net supply of the second asset by another step. Re-solve the economy starting from weights and states from the previous step. Repeat this procedure until the supply target for the second asset is reached.
* If required, perform the same procedure with no-short sale constraint (_i.e._ relax it slowly).
* Set the weight in the loss function on the Euler equations associated with the third asset to its full value (_e.g._ 1).
* Since the weight on the Euler equations is non-zero, the network will be able to learn a price function implied by zero allocations of that asset. Recall that this third asset is in zero net supply, and agents are bound by strict no short-sale when trading it.
* Set the mask parameter for the third asset to a small, but non-zero value (_e.g._ 0.01). Now the network can predict non-zero allocations of the third asset. The choice of a small mask value helps to avoid instability associated with a potentially large and discontinuous jump.
* Set the net supply of the third asset at some small, but non-zero level. Solve this economy starting from network weights from the previous step.
* Increase the net supply of the third asset by another step. Re-solve the economy starting from weights and states from the previous step. Repeat this procedure until the supply target for the third asset is reached.
* If required, perform the same procedure with no-short sale constraint (_i.e._ relax it slowly).
* In case of adding the fourth and further asset, perform the analogous procedure.

### Example with Three Assets

#### 3.1.1 Model

To illustrate the homotopy method, together with our market-clearing neural network architecture, we enrich the overlapping generations model introduced in section 2.5.1 with two additional assets. Next to a one-period risk-free bond, households can also trade claims to a Lucas tree as well as claims to the stream of rent payments for housing services. We think of the former as a stock and of the latter as house ownership. The households hence have to allocate their total savings between three different assets. This portfolio choice problem is challenging for deep learning based solution methods, rendering the model a suitable testing ground for our homotopy method. The remaining components, are the same as in the model we previously described in 2.5.1.

Asset structureAs before, households can purchase a risk-free one-period bond in net-supply \(B\) at equilibrium prices \(p_{t}^{b}\). Additionally, households can trade claims to a Lucas tree at equilibrium prices \(p_{t}^{s}\), which in net-supply \(S\).14 Every period, the Lucas Tree pays an exogenous amount of dividends \(d_{t}=dz_{t}\), which is perfectly correlated with aggregate labor income. Next to deriving utility from consumption, household derive age-specific utility from consuming housing services, which they rent every period at equilibrium rent \(p_{t}^{r}\). Finally households can trade claims to the aggregate housing stock at equilibrium
price \(p_{t}^{o}\), which is in net-supply \(H^{o}\). The total housing stock in the economy, which can be rented by household is the sum of the locally owned housing supply \(H^{o}\), which are the houses owned by households, and an housing supply \(H^{ex}\), with owners outside the model. In this model a claim to the housing stock is a claim to the streams of rents. Owning a house is hence similar to owning a Lucas tree, where the dividends are replaced by the rent payments the house owners receive. None of the assets can be sold short, and all of the assets are traded subject to asset-specific quadratic adjustment costs \(p_{t}^{x}\zeta^{x}\frac{1}{2}(x_{t+1}^{h+1}-x_{t}^{h})^{2}\), for \(x\in\{b,s,o\}\). We model the bond as the most liquid asset and house ownership as the least liquid asset, _i.e._\(\zeta^{h}>\zeta^{s}>\zeta^{b}\). The Bellman equation corresponding to the households' problem is given by

\[V^{h}(z_{t},b_{t}^{h},s_{t}^{h},h_{t}^{o,h}) =\max_{h_{r,t}^{h},b_{t+1}^{h+1},s_{t+1}^{h+1},h_{t+1}^{o,h+1}} \left\{u(c_{t}^{h})+\psi^{h}v(h_{t}^{r,h})+\beta\mathrm{E}\left[V^{h+1}(z_{t+1 },b_{t+1}^{h+1},s_{t+1}^{h+1},h_{t+1}^{o,h+1})\right]\right\} \tag{28}\] \[c_{t}^{h} =z_{t}y^{h}+b_{t}^{h}+s_{t}^{h}(p_{t}^{s}+d_{t}^{s})+h_{t}^{o,h}( p_{t}^{o}+p_{t}^{r})\] \[-b_{t+1}^{h+1}p_{t}^{h}-s_{t+1}^{h+1}p_{t+1}^{s}-h_{t+1}^{o,h+1}p_ {t}^{o}-h_{t}^{r,h}p_{t}^{r}\] \[-p_{t}^{b}\zeta^{b}\frac{1}{2}\left(b_{t+1}^{h+1}-b_{t}^{h} \right)^{2}-p_{t}^{s}\zeta^{s}\frac{1}{2}\left(s_{t+1}^{h+1}-s_{t}^{h}\right) ^{2}-p_{t}^{o}\zeta^{h}\frac{1}{2}\left(h_{t+1}^{o,h+1}-h_{t}^{o,h}\right)^{2}\] (29) subject to : \[0 \leq b_{t+1}^{h+1}-\underline{b}\] (30) \[0 \leq s_{t+1}^{h+1}-\underline{s}\] (31) \[0 \leq h_{t+1}^{o,h+1}-\underline{b}^{o} \tag{32}\]

The corresponding optimality conditions, which, together with market clearing, characterize the functional rational expectations equilibrium, are given in Appendix B. We provide details on parameter choices in Appendix C.

Heterogeneity in risk aversionWe continue to model an overlapping generations economy, where households live deterministically for \(H\) periods and receive an age specific share of aggregate labor income. Additionally, we now assume that in each cohort, there are two types of households of equal mass. One type has a lower risk aversion, the other type as a higher risk aversion.

State of the economyThe state of the economy is given by the exogenous shock, as well as the distribution of assets accross households and risk-aversion types. Let the superscript 1 denote the low risk aversion households and superscript 2 denote the high risk aversion households. The state of the economy is given by

Footnote 15: The \((H-1)\) stems from the fact that agents are born without assets and their asset holding is hence always constant and does not need to be included in the state.

\[\mathbf{x}_{t}=[z_{t},\underbrace{b_{t}^{1,1},\ldots,b_{t}^{H,1},b_{t}^{1,2},\ldots,b_{t}^{H,2}}_{\text{bond holdings}},\underbrace{s_{t}^{1,1},\ldots,s _{t}^{H,1},s_{t}^{1,2},\ldots,s_{t}^{H,2}}_{\text{stock holdings}},\underbrace{h_{t}^{ o,1,1},\ldots,h_{t}^{o,H,1},h_{t}^{o,1,2},\ldots,h_{t}^{o,H,2}}_{\text{ housing ownership}} \tag{33}\]

The state space is hence \(1+2\times 3\times(H-1)\)-dimensional.16 All policy and price functions are a function of the whole state, rendering this model a formidable testing ground for our algorithm.

Footnote 16: The \((H-1)\) stems from the fact that agents are born without assets and their asset holding is hence always constant and does not need to be included in the state.

Equilibrium functions to approximateWe need to approximate \(2\times(H-1)\times 4+2\) household policies: for each risk aversion type and each age group, we need the bond policy, the stock policy, the renting policy
as well as the house ownership policy. While we only need to approximate \(H-1\) policies per asset and risk aversion type, we need to approximate \(H\) policies for the intra-temporal choice of renting housing services. Furthermore, we need to approximate four equilibrium price functions: for the bond price, the tree price, the price for renting, and the price for house owner ship.

We approximate the mapping from the \(1+2\times 3\times(H-1)\) dimensional state to the \(2\times(H-1)\times 4+2+4\) endogenous variables by a deep neural network.

\[\mathcal{N}_{\rho}:\mathbb{R}^{1+2\times 3\times H}\rightarrow\mathbb{R }^{2\times(H-1)\times 4+2+4}\] \[\mathcal{N}_{\rho}(\mathbf{x})=[\tilde{b}_{t+1}^{2,1},\dots, \tilde{b}_{t+1}^{H,1},\tilde{b}_{t+1}^{2,2},\dots,\tilde{b}_{t+1}^{H,2},\] \[\tilde{s}_{t+1}^{2,1},\dots,\tilde{s}_{t+1}^{H,1},\tilde{s}_{t+1}^ {2,2},\dots,\tilde{s}_{t+1}^{H,2},\] \[\tilde{h}_{t+1}^{o,2,1},\dots,\tilde{h}_{t}^{o,H,1},\tilde{h}_{t+ 1}^{o,2,2},\dots,\tilde{h}_{t+1}^{o,H,2},\] \[\tilde{h}_{t}^{r,1,1},\dots,\tilde{h}_{t}^{r,H,1},\tilde{h}_{t}^ {r,2,1},\dots,\tilde{h}_{t}^{r,H,2},\] \[\hat{p}_{t}^{b},\hat{p}_{t}^{s},\hat{p}_{t}^{o},\hat{p}_{t}^{r}] \tag{34}\]

We use the simple adjustment described in section 2.1 to ensure that market clearing conditions are always satisfied and use our homotopy algorithm to solve the model step-wise by slowly increasing the supply and number of assets.

Loss functionSince our neural network architecture already enforces market clearing, the only remaining equilibrium conditions are the households' optimality conditions for each of the four choices. Analogously to equation (24), we use the Fischer Burmeister equation to summarize each set of KKT conditions into a single equation. Let \(\left\{\{\mathrm{err}_{\rho}^{b,i,h}(\mathbf{x}_{t}),\mathrm{err}_{\rho}^{s,i,h}(\mathbf{x}_{t}),\mathrm{err}_{\rho}^{o,i,h}(\mathbf{x}_{t})\}_{h=1}^{H-1}, \{\mathrm{err}_{\rho}^{h^{r},i,h}(\mathbf{x}_{t})\}_{h=1}^{H}\right\}_{i=1}^{2}\) denote the corresponding errors in the equilibrium conditions. Our loss function is given by

\[l_{\rho}(\mathcal{D}) :=\frac{1}{2}\left(\right.\] \[w_{b}\frac{1}{|\mathcal{D}|}\frac{1}{H-1}\left(\sum_{\mathbf{x} \in\mathcal{D}}\sum_{h=1}^{H-1}\left(\mathrm{err}_{\rho}^{b,1,h}(\mathbf{x}_{ t})\right)^{2}+\left(\mathrm{err}_{\rho}^{b,2,h}(\mathbf{x}_{t})\right)^{2}\right)\] \[+w_{s}\frac{1}{|\mathcal{D}|}\frac{1}{H-1}\left(\sum_{\mathbf{x} \in\mathcal{D}}\sum_{h=1}^{H-1}\left(\mathrm{err}_{\rho}^{s,1,h}(\mathbf{x}_{ t})\right)^{2}+\left(\mathrm{err}_{\rho}^{s,2,h}(\mathbf{x}_{t})\right)^{2}\right)\] \[+w_{o}\frac{1}{|\mathcal{D}|}\frac{1}{H-1}\left(\sum_{\mathbf{x} \in\mathcal{D}}\sum_{h=1}^{H-1}\left(\mathrm{err}_{\rho}^{o,1,h}(\mathbf{x}_{t}) \right)^{2}+\left(\mathrm{err}_{\rho}^{o,2,h}(\mathbf{x}_{t})\right)^{2}\right)\] \[+w_{h^{r}}\frac{1}{|\mathcal{D}|}\frac{1}{H}\left(\sum_{\mathbf{x} \in\mathcal{D}}\sum_{h=1}^{H}\left(\mathrm{err}_{\rho}^{r,1,h}(\mathbf{x}_{t}) \right)^{2}+\left(\mathrm{err}_{\rho}^{r,2,h}(\mathbf{x}_{t})\right)^{2}\right)\] \[\left.\right) \tag{35}\]

The weights \(w_{x}\) allow to vary the weight on the equilibrium conditions. The loss function hence nests a single and two-asset model, by allowing the corresponding weight in the loss function to be set to zero, while also ensuring that the associated policies are equal to zero.
#### 3.1.2 Training

To parameterize equilibrium objects of the economy (_i.e._ policies and prices), we use a densely connected feed-forward neural network with two hidden layers. The dimensionality of its input is \(N^{\text{input}}=161\) since the state vector includes aggregate productivity shock and asset distribution across all twenty age groups for both types in the economy, and the distribution of dividend income, which we choose as an auxiliary state variable.16 For the two hidden layers, we choose \(N^{\text{hidden 1}}=N^{\text{hidden 2}}=400\) neurons with relu activation functions. Finally, the network output is a vector of \(N^{\text{output}}=158\) elements consisting of the policies for the three assets, the policies for renting housing services, and four price functions for the price of the bond, the price of stocks, the price for house ownership, and the price for renting housing services.17 Outputs corresponding to prices are transformed using the softplus activation in order to ensure that the predicted prices are always positive.

Footnote 16: The minimal state vector includes 121 state variables, that includes distribution of three assets across 20 age groups for each of two risk aversion groups and the aggregate shock. Following Azinovic et al. (2022), we augmented the input of the neural network by a vector of auxiliary statistics. In our case, we choose a distribution of dividend income as an auxiliary statistic.

Footnote 17: For each asset category, we approximate the 19 savings function, exploiting the fact that the oldest age group does not save in our model. Furthermore, we approximate the 20 housing service consumption function, since housing services are demanded by all age groups.

To ameliorate numerical instabilities associated with portfolio choice, we employ the homotopy algorithm described at the beginning of section 3. Specifically, we firstly solve a bond-only economy. To do so, we set the \(S^{\text{initial}}=0\), \(H^{\text{o,initial}}=0\), and \(H^{\text{ex,initial}}=1\) and run the training procedure described in section 2.4.18 Then, we approximate the stock pricing function implied by the consumption dynamics of the bond-only economy.19 Using the price and policy functions of the bond-only economy as an initial guess, we compute equilibrium of the two asset economy featuring a small supply of stocks. Starting from this solution, we gradually increase stock supply and retrain the neural network after each increase. We do so until we reach the stock supply target. In total, we perform \(S^{\text{steps}}=10\) homotopy steps to increase the stock supply from \(S^{\text{initial}}=0\) to \(S^{\text{final}}=1\). We proceed analogously for the case of housing ownership. In that case, we perform \(H^{\text{steps}}=20\) housing supply steps. We simultaneously increase the supply of private housing \(H^{o}\) from \(H^{\text{o, initial}}=0\) to \(H^{\text{o, final}}=1\) and decrease the external housing supply \(H^{\text{ex}}\) from \(H^{\text{ex, initial}}=1\) to \(H^{\text{ex, final}}=0\), such that the total housing supply \(H^{o}+H^{ex}\) remains constant and equal to one.

Footnote 18: We use the same training hyperparameters as in subsequent homotopy steps, with exception of \(N^{\text{episodes}}\), which we set to 512, in order to obtain a high-precision solution to start from. The value for the remaining hyperparameter are provided in table 4.

Footnote 19: In this step, we use the same training hyperparameters as in subsequent homotopy steps. See table 4.

In each step of the homotopy algorithm we continue to train the neural network to approximate the equilibrium functions for the economy with the corresponding values for the aggregate supply of assets in the economy, _i.e._ for the intermediate values for \((B,S,H^{o},H^{ex})\). Again, we use the learning procedure described in section 2.4. To obtain training states, we simulate \(N^{\text{trajectories}}=8192\) parallel state trajectories for \(N^{\text{episodes}}=256\) episodes.20 For each episode (_i.e._ each simulation step), we train the neural network for \(N^{\text{epochs}}=10\) epochs with a learning rate of \(\alpha^{\text{learn}}=1\times 10^{-6}\).21

Footnote 20: In each epoch, we split the simulated dataset of \(N^{\text{trajectories}}\) states into batches of \(N^{\text{minibatch}}=128\) states. Each batch is used to compute one stochastic gradient descent update of the neural network weights. As for the single asset model, we transformed loss function gradients using the zero_nams transformation before using them as an input for the Adam optimizer and re-start internal state of Adam optimizer to zero at the beginning of training episode.

Figure 2 shows the asset policies during different stages of our homotopy training algorithm. Figure 3
shows the corresponding policies for consumption and the renting of housing services. The top panel in figure 2 shows the trained bond policy for each risk aversion type and age group, when the policies for stock and house ownership are still masked and the corresponding terms in the loss function are weighted with zero. The second panel shows the policies, when a small amount of stock is introduced to the model, the policies for house ownership are still masked and equal to zero. The third panel shows the policies, when a small amount of house ownership is introduced to the model. Finally, the last panel shows the household policies, when the full amount of house ownership is introduced to the economy. As the figures illustrate, the policies remain stable throughout the training.

Figure 3 the corresponding profiles for consumption and the renting of housing services over the life-cycle during the stages of the homotopy algorithm. Increasing the supply of assets allows for a smoother consumption profile over the life-cycle. In the top panel, which shows the economy with only a risk-free bond, the consumption profile over the life-cycle is slightly decreasing for unconstrained agents due to the low interest rate. As more assets are added to the economy, the bond price drops and the interest rate rises. In the bottom panel, the higher interest rate leads to an increasing consumption profile over the life-cycle.

Comparing the households with high risk aversion to the households with low risk aversion, figure 3 shows that households with high risk aversion increase their consumption and save less early in life, at the cost of less consumption later in life. For time separable expected utility with constant relative risk aversion, the risk aversion also pins down the intertemporal elasticity of substitution. The households with higher risk aversion are hence less willing to substitute consumption between periods and choose a more stable consumption sequence over their life-cycle.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Parameters & \(N^{\text{episodes}}\) & \(N^{\text{trajectories}}\) & \(N^{\text{epochs}}\) & \(N^{\text{minibatch}}\) & \(N^{\text{integration}}\) & \(\alpha^{\text{learn}}\) \\ \hline Values & 256 & 8192 & 10 & 128 & 8 & \(10^{-6}\) \\ \hline \hline \end{tabular}
\end{table}
Table 4: Hyperparameters for training steps within homotopy loop.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Parameters & \(N^{\text{input}}\) & \(N^{\text{hidden 1}}\) & \(N^{\text{hidden 2}}\) & \(N^{\text{output}}\) \\  & & Activations & Activations & Activations \\ \hline \multirow{2}{*}{Values} & \multirow{2}{*}{161} & 400 & 400 & 158 \\  & & relu & relu & see text \\ \hline \hline \end{tabular}
\end{table}
Table 3: Network Architecture chosen for the three asset model with risk aversion heterogeneity.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline Parameters & \(S^{\text{initial}}\) & \(S^{\text{final}}\) & \(H^{\text{o,initial}}\) & \(H^{\text{o,final}}\) & \(H^{\text{ex,initial}}\) & \(H^{\text{ex,final}}\) & \(S^{\text{steps}}\) & \(H^{\text{steps}}\) \\ \hline Values & 0 & 1 & 0 & 1 & 1 & 0 & 10 & 20 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Homotopy loop hyperparameters
Figure 2: Policy functions for the different assets at different stages of trainig of the homotopy algorithm.
Figure 3: Policy functions for consumption and the renting of housing services at different stages of training of the homotopy algorithm.
#### 3.1.3 Accuracy

To assess the accuracy of our solution we compute the errors in the equilibrium conditions on an approximate ergodic set of the economy. To obtain the approximate ergodic set, we simulate the dataset of \(N^{\text{trajectories}}=8192\) states, generated by the last step of our training procedure, forward without retraining the neural network.22 As shown in the figure 4, the remaining errors are low, with the 99th percentile well below 1% for each age group and optimality condition.

Footnote 22: We simulate the final dataset 256 periods forward without further training to test the “out-of-sample” performance of our approximate solution.

## 4 Conclusion

We make two distinct but complementary contibutions in the context of deep learning based solution methods. First, we show how market clearing and borrowing constraints can be encoded directly into the architecture of the neural network. Second, we present a homotopy algorithm for solving portfolio choice problems. The idea of our homotopy algorithm is to start solving a single asset and then add additional assets one by one and to slowly increase their supply in the economy.
Figure 4: Statistics on the errors in the remaining equilibrium conditions after training, expressed in % and as relative consumption errors. The black solid line shows the mean error by age group, the blue dashed line shows the 90th percentile, the yellow dash-dotted line shows the 99th percentile and the dotted red line shows the maximum error by age group. The left column shows the errors for the households with low risk aversion, the right panel shows the errors for the households with high risk aversion. The top row shows the errors in the optimality conditions for bond purchases, the second row shows the errors in the equilibrium conditions for the choice of housing services, the third row shows the errors for the equilibrium conditions for house ownership, and the last row shows the errors in the equilibrium conditions for the stock choice.