Modeling Financial System with Interbank Flows, Borrowing, and Investing

In our model, private actors with interbank cash flows similar to, but nore general than (Carmona, Fouque, Sun, 2013) borrow from the outside economy at a certain interest rate, controlled by the central bank, and invest in risky assets. Each private actor aims to maximize its expected terminal logarithmic utility. The central bank, in turn, aims to control the overall economy by means of an exponential utility function. We solve all stochastic optimal control problems explicitly. We are able to recreate occasions such as liquidity trap. We study distribution of the number of defaults (net worth of a private actor going below a certain threshold).


Introduction
Modeling financial and banking systems is a focus of much recent research. Of particular interest is the concept of systemic risk, which can be informally described as the probability of a large number of banks defaulting or getting into financial trouble. We are interested in probability of this undesirable event; of the mechanism of such failure; and of the financial contagion, when failure of a few banks leads to many more failures. We refer the reader to the handbook [13] containing many different approaches to systemic risk. Our work is inspired by the model introduced in [7] and also described in [4,Section 5.5].
If X i (t) is the wealth and Y i (t) := log X i (t), authors in [7] model the banking system as a system of N continuous-time stochastic processes Y 1 , . . . , Y N , with multidimensional Ornstein-Uhlenbeck dynamics and the SDE given by: The constant a is referred to as the interbank flow rate. In [7], these mean-reverting drifts are generated by the decisions of banks to borrow money from one another. Their decisions are done by minimizing a certain cost functional, which measures, roughly speaking, the preference of a bank to borrow from other banks, as opposed to borrowing from the central bank. Authors discuss both the finite player solution and the mean-field limit of the problem in the context of systemic risk. In this paper, we further explore the individual decision-making of the private banks and extend the role of the central bank. Furthermore, we analyze how this decisionmaking affects the stability of the system. More specifically, we extend the model by assuming that each private bank invests in a risky portfolio of assets, borrows money from the general economy to invest in this portfolio (with interest rate controlled by the central bank), pockets the profit, and pays back the interest. Unlike in [7] , where the decision making of the central bank is not analyzed, in our model the central bank uses interest rate as a policy tool in order to govern the behavior of the private banks. It is derived as a solution of the control problem solved by central bank.
We assume two kinds of players in our model: private banks and the central bank. Private banks want to maximize the terminal logarithmic wealth: (1.3) sup E [log X i (T )] , i = 1, . . . , N, through borrowing and investing in the portfolio of risk assets. For simplicity, we assume the portfolios of private bank to be correlated geometric Brownian motions. The choice of logarithmic utility function allows us to solve the corresponding Hamilton-Jacobi-Bellman (HJB) equation explicitly.
On the other hand, the central bank wants to control the overall size of the financial system, using interest rate r > 0 as a monetary policy instrument, which measures how attractive it is for banks to borrow from the general economy and invest, as opposed to sitting on cash. As we see in Section 4, this interest rate r controls the overall size of the system, measured by The central bank chooses r to maximize terminal expected CARA (exponential) utility (1.4) E − exp(−λY (T )) for some λ > 0.
This corresponds to the central bank being even more risk-averse than private banks: Private banks have utility function which is linear in Y i (t), and the central bank has utility function which is concave in these variables. Given the choice of the interest rate r, the problem becomes a stochastic N -player game for private banks. We can solve the stochastic optimal control problem explicitly for each player: private banks and the central bank. This is due to the special choice of logarithmic utility function for private banks in (1.3), and for the central bank in (1.4). Therefore, there is no need to use mean-field game theory. For other choices of utility function, it is probably impossible to solve the stochastic game explicitly. Then one could try to use mean-field limits, as in [8,9,19]. This topic is left for future research.
This setup somewhat resembles the principal-agent problem: the principal (now the central bank) allows private banks to borrow from the economy, and private banks (agents) maximize their expected logarithmic terminal utility (their contract).
Besides incorporating the optimal strategy of the central bank, we also generalize (1.1) by allowing interbank flow rates from bank i to bank j to depend on the banks i, j, and on time t. We denote such rate by c ij (t): (1.5) dY i (t) = 1 N to be correlated: that is, we assume W = (W 1 , . . . , W N ) is an N -dimensional Brownian motion with drift vector µ and covariance matrix A.
As mentioned before, systemic risk is the prospect of multiple defaults. In order to quantify this, we define default of the ith bank as having Y i (t) < D for some t ∈ [0, T ], where T is a given time horizon, and D < 0 is a certain threshold. That is, a bank is in default if its capital at some point dropped below the level e D . We also consider the number of defaults: where 1(·) is the indicator random variable. We compute numerically the distribution of the (random) number of defaults. Authors in [7] show that as long as banks (at least, most of them) stay above the threshold, their interaction through Ornstein-Uhlenbeck drift terms keeps each of them from default. But as soon as one or more of the banks are in default, the other banks are pushed towards default too. We derive this in subsection 3.4. This conclusion stays true for the current model. To summarize, the topic of this paper is optimal decisions of private banks, and interaction of these individual decisions with each other, as well as with that of the central bank. Systemic risk is a consequence of the optimal decision making of both the central bank and the private banks in the economy.

Contributions.
Our work here is inspired by Carmona et. al. [7], however, unlike them, we view the systemic risk through the lens of utility maximization for both the central bank and the private banks. We first solve for the optimal investment of the private banks in their portfolios of assets, given the interest rate set by the central bank. Next, we turn around and find the optimal interest rate for the central bank, given the strategy of the private banks and its own risk aversion. To the best of our knowledge, we are the first to present systemic risk using a joint framework of central bank and private banks in a game to maximize their expected utility. Unlike in [7], where the interbank flow rate was constant, here we generalize the flow rates to be different for each pair of private banks.
1.2. Organization of the paper. In Section 2, we describe the model in terms of stochastic control problem for a system of stochastic differential equations. In Section 3, we solve the stochastic control problem for each private bank, and in Section 4, for the central bank (given optimal control for each private bank). Section 5 contains results on long-term stability of the system: the fact that the capitals of banks tend to stay close, as opposed to splitting into two or more groups. Section 6 is devoted to concluding remarks and suggestions for future research. The Appendix contains the statements of technical lemmata, used in the proof of the long-term behavior of this system Y (t) in Theorem 5.1.

1.3.
Review of related models. The fundamental wealth dynamics (1.1) of this model has been studied under different settings. For example, in [6], a similar system was studied with time delay. Large deviations were studied in [14]. Without attempting to give an exhaustive survey, let us mention the following papers, which use stochastic differential equations and interacting Brownian particles to model dynamics of capitals of banks or other financial agents. Fouque and Ichiba [12] use a system of stochastic differential equations with Bessel-type diffusion coefficients to model simultaneous defaults (in this model, a default is when the capital reaches zero). Authors in [3,26] combine such Bessel-type diffusion coefficients with a meanfield-type drift term, with [3] having an additional jump term (therefore the processes there are jump-diffusions).
In the paper [23], the financial system is modeled by independent geometric Brownian motions, with defaults happening at hitting times of some lower threshold. Once a bank defaults, other banks see their capital decrease by a certain amount, possibly triggering a cascade of defaults. Carmona et. al. [5], introduce mean-field game of timing. The term game of timing refers to a game where each player chooses an optimal stopping time. Mean-field game of timing is when a player competes against a "crowd" of other agents, instead of individual competitors; this can be informally viewed as a limit of games of timing as the number of players tends to infinity. This models a bank run, continuing the research in a celebrated paper [10] by Diamond & Dyvbig. Let us also mention the paper [20], which models the dynamics by geometric Brownian motions. There, they study mean-field relative performance criteria. That is, an agent again competes against a "crowd", and we maximize the agent's performance compared to the performance of the "crowd".
1.4. Notation. For a vector or a matrix a, its transpose is denoted by a . We usually think of vectors as column-vectors. The dot product of two vectors a and b is denoted by a · b. The term standard Brownian motion stands for a one-dimensional Brownian motion with drift coefficient 0 and diffusion coefficient 1. For V ≡ 1, this is called the total variation norm. Fix a dimension N ≥ 2. Then e ∈ R N is a vector (1, . . . , 1) with unit components, and we define the following hyperplane in R N : Define the (closed) ball of radius r on Π centered at the origin: The (N − 1)-dimensional Lebesgue measure on Π is denoted by mes Π (·). As mentioned above, the symbol 1(A) or 1 A stands for the indicator function of an event A.
2. Description of the model 2.1. Formal description. Consider a system of N banks which continuously lend money to each other, borrow from the outside economy, pay back the interest, and invest in some risky portfolios. We operate on filtered probability space (Ω, F t , (F t ) t≥0 , P) with the filteration satisfying the usual conditions. All the process which we consider are adapted to the (F t ) t≥0 . Let X i (t) > 0 be the monetary reserves of the ith bank at time t, for i = 1, . . . , N . Let Z i (t) be the amount borrowed at the moment t by the ith private bank from the outside economy. Assume the interest rate for such borrowing is r(t) ≥ 0, controlled by the central bank. Then during the time interval [t, t + dt], the ith bank pays back interest r(t)Z i (t) dt. At time t, the ith bank has at its disposal the amount X i (t) + Z i (t): its own capital plus borrowed amount. This amount Z i (t) ≥ 0 is controlled by the ith bank.
Alternatively, the ith bank might decide to not borrow anything, and instead to even put aside some of its own money in cash (which does not earn any interest). This happens if the investment is not very profitable, or, more precisely, if the return does not outweigh the risk. In this case, we let Z i (t) < 0, and define −Z i (t) to be the quantity of cash put aside. The amount invested is still X i (t) + Z i (t), but the bank does not pay or receive any interest.
We combine these two cases: the ith bank invests the amount X i (t) + Z i (t) at time t into a risky portfolio, and pays interest r(t)(Z i (t)) + dt during the time interval [t, t + dt].
At time t, the i th bank invests in a portfolio of risky assets with value S i (t). The ith bank buys (X i (t) + Z i (t))/S i (t) units of this portfolio. Net profit for the time interval [t, t + dt] is Combining all of the above, we get the following system of equations: Next, we make some assumptions on S i , the dynamics of the portfolio processes. A separate question is how banks construct these portfolios out of stocks and other risky assets. This question is separate from the topic of this paper, and we shall not study it here. Instead, we assume that these are geometric Brownian motions. This assumption is very simplifying, but we believe it captures to some extent the features of portfolios. The processes form an N -dimensional Brownian motion with drift vector µ = (µ 1 , . . . , µ N ) and covariance matrix A = (a ij ) i,j=1,...,N . In particular, each M i , i = 1, . . . , N , is a Brownian motion with drift coefficient µ i and diffusion coefficient σ 2 i := a ii , so it can be represented as where W i is a one-dimensional standard Brownian motion. Although the portfolio process (2.2) is driven by only one Brownian motion, a more general representation: where (B 1 , . . . , B m ) are brownian motions, can also be considered in our framework. Since where W i is a brownian motion, we fall back to the original portfolio process (2.2).
The covariance structrure of Brownian motions (W 1 , . . . , W N ) can be modeled in various ways. We consider four cases. Case 1. All W 1 , . . . , W N , are independent. Then the matrix A is diagonal: . This means that the portfolios of banks are independent.
Case 3. An intermediate case: for some i.i.d. Brownian motionsW i , i = 0, . . . , N , and some coefficients ρ 0 ,ρ 0 with ρ 2 0 +ρ 2 0 = 1 we have: Case 4. One can also split N banks into subsets and construct dependence as in Case 3 for each subset; portfolio processes corresponding to different subsets are assumed to be independent.
Apply Itô's formula to find the dynamics of Y i (t) := log X i (t): i (t) Superimposing (2.1) and (2.6), we get: where we define the relative investment ratio: and the following quantity: Finally, the ith bank also borrows and/or lends from other banks. We superimpose on top of (2.7) the Ornstein-Uhlenbeck-type drift term from (1.5), similar to the one in [7], but more general. The model (1.5) gives us more flexibility, compared to (1.1). Some possible cases: Then there are no cash flows between banks.
For a constant c, this is the model from [7].
After superimposing these Ornstein-Uhlenbeck-type drifts on top of (2.7), it takes the form Equation (2.10) resembles the model from [7]. However, it also has significant differences: the volatility in (2.10) can be controlled, unlike in [7]; and the drift coefficient in (2.10) is a bit more complicated. For homogeneous rates: c ij (t) ≡ c(t), the equation (2.10) takes the form where D is a given threshold, stipulated by the central bank. The central bank would like to stimulate the activity of banks by persuading them to take risks, but not too much, lest they may become bankrupt. Equation (2.7) means that the central bank can use interest rate r(t) as a monetary policy tool to alter the behaviour of the private-banks.
Assume that banks start borrowing too much money and investing them in risky assets (leveraging). By doing this, they increase their probability of default. Then the central bank can raise this interest rate to discourage private banks from excessive borrowing. Conversely, if banks are too cautious in borrowing against future profits and risk-taking, then the central bank can stimulate them by lowering the interest rate. As we see later, this interest rate affects the overall state of the system.
The parameter r(t) is determined by the central bank and given to all private banks. These private banks then determine the investment rates α i (t), independently of each other. In light of decision-making of the banks, the central bank needs to determine optimal values of these parameters. This setup is similar to the principal-agent problem, but with many agents.

Optimal behavior of private banks
3.1. Statement of the problem. We assume that the ith bank takes as given the capital of other banks: X j (t), j = i (or, equivalently, Y j (t) := log X j (t), j = i), as well as the interest rate r(t) (the instrument of the monetary policy). The bank is trying to choose the relative investment ratio α i (t), or, equivalently, the amount borrowed Z i (t), to maximize its expected terminal logarithmic wealth: where the supremum in (3.1) is taken over all bounded adapted controls Assume that the interest rate r(t) is already set by the central bank. Then this stochastic control problem is, in fact, a stochastic game between N players: the private banks. In this section, we solve this stochastic control problem explicitly. This corresponds to the agent's problem in the principal-agent framework. In the next section, we discuss the optimal policy choices of the central bank (the principal).

3.2.
Solution of the problem. This specific choice of the utility function, which is linear in Y i , and the interbank flows, which are also linear in Y i in (2.7), makes this optimization problem tractable: We can solve this explicitly.
Theorem 3.1. For the optimization problem stated in (3.1), where α i is bounded, adapted on [0, T ], the following value of α i is optimal for the ith private bank: Remark 3.1. In particular, if µ i ≤ σ 2 i , that is, the return on the investment does not outweigh its risks, then the ith bank does not borrow anything to invest. On the contrary, this bank sets aside money as cash. If µ i ≥ σ 2 i , the investment is attractive for borrowing, but a high enough interest rate: i can preclude the ith bank from borrowing; then this bank will invest only its own money into the portfolio. Only if the interest rate is low enough: r(t) < µ i − σ 2 i , the ith bank borrows money to invest.
Remark 3.2. Note that the optimal strategies (3.2) do not depend on the flow rates c ij , because of the special choice of logarithmic utility function, which is linear in Y i . Although logarithmic utility function leads to myopic agents, this assumption is important for mathematical tractability of the results. CRRA utility function is another popular form used in the literature, however we were unable to evaluate the optimal control for it even in the mean field case.
Proof. The dynamic programming principle tells us that the function where we take the supremum over all α i which are bounded and adapted on [t, T ], satisfies the Hamilton-Jacobi-Bellman (HJB) equation: with terminal condition Φ i (T, y) = y i . We assume all α j , j = i, are already chosen. Try the following Anzats, linear in y j : Because it is linear, the second-order derivatives in (3.3) turn out to be zero. Therefore, the only term in (3.3) which needs to be maximized is h(α i , r(t)). The solution to this maximization problem is given by the value α * i from (3.2). This is a simple algebraic exercise; detailed calculations are given in Lemma 6.3 in the Appendix.
The maximal value of h i (α, r(t)) is This means that the ith bank chooses the control value α i := α * i . This value is independent of terminal time T , and of the values of Y j , j = 1, . . . , N . This corresponds to the classical solution of the Merton problem. If r is constant (independent of t), then α * i and h * i are also constant. Comparing (3.4) with the terminal condition, we have: Next, plug the anzats (3.4) into (3.3). Note that all second-order derivatives of the anzats (3.4) are equal to zero, and first-order derivatives are In addition, the time derivative of this value function Φ from (3.4) is Combining (3.2), (3.5), (3.7), (3.8), we get that the HJB equation (3.3) takes the form Comparing coefficients in (3.9) at each y j , we see that The free terms in (3.9) sum up to Together with terminal conditions (3.6), this system (3.10) and (3.11) of N + 1 linear ODEs has a unique solution g i0 , . . . , g iN . This solves the HJB equation.
To complete the proof, let us do the verification argument. Take a bounded adapted control (3.12) Using the boundedness of α j = (α j (t), 0 ≤ t ≤ T ), we get that the stochastic integral term in (3.12) has expectation zero. Combining (3.3) with (3.12), we get that (Φ i (t, Y (t)), t ≥ 0) is a supermartingale for all admissible (adapted bounded) controls α i , but a martingale for the control , with equality for the control α * i . From here it immediately follows that α * i is indeed the optimal control.
3.3. The dynamics of banks under optimal control. Under the optimal control (3.2), the processes Y i , i = 1, . . . , N (we shall denote them by Y * i ) satisfy the following system of stochastic differential equations: . If r = const, then M * is an N -dimensional Brownian motion with drift vector and covariance matrix given by (3.14) µ * = (µ * 1 , . . . , µ * N ), µ * i := h i (α * i , r). The dynamics (3.13) is similar to that in [7]. If r(t) does not depend on t, then M * = (M * 1 , . . . , M * N ) , like (M 1 , . . . , M N ) , is an N -dimensional Brownian motion, but with different drift vector and covariance matrix. As in (1.2), we define Averaging equations (3.13) and using the symmetry property c ij = c ji , we have: The interest rate r controls the overall size of the system, measured by Y . Express (3.16) as: where W is a standard Brownian motion, and the coefficients g(·) and ρ(·) are defined as: for all i; that is, all portfolios are profitable to invest, at least for zero interest rate r = 0. First, in Figure 1 we assume (2.4), that is, the portfolio processes S 1 , . . . , S N , are independent. We also assume that there are no flows: We take three interest rates r: 0%, 12%, and 20% respectively. As expected, increasing the interest rate forces the banks to borrow less and thus optimal investment ratio α * becomes 0 in Figure 3.3)c) while it varied between 2 to 12 in Figure 3.3)a).
We observe that the banks with significant flows tend to have wealth dynamics which are more "tied" together. Moreover, this adds to the stability to the system as we observe lesser defaults for i = 1, . . . , 10 compared to j = 11, . . . , 30. Finally, in Figure 3 we assume that the portfolio processes are correlated, as in (2.5), with ρ 0 := 0.5, flows are given by (3.20) and interest rate r is 8%. Compared to Figure 2, the impact of correlation on the dynamics of the banks is clearly visible through movement of the wealth dynamics strongly tied together.    Figure 5 we present empirical cumulative distribution function (CDF) for the number of defaulted banks assuming correlated portfolios and no interbank flows for different interest rates. The corresponding histogram is presented in figure 6. As stated in previous studies, increase in correlation increases the probability of large defaults and at the same time reducing the small default probabilities, similar to flocking behaviour in various biological studies. So, in Figure 7 we present empirical estimates of P(D > 60) and P(D < 5), for N = 100 banks, as a function of the correlation between their portfolio process at different interest rates. As expected, we observe the increase in probability of large and small default as the correlation increases. However, increasing the interest rate reduces the probability of large default at the expense of small default probability. Finally, in Figure 8 we present the empirical CDF for the number of defaulted banks assuming correlated portfolio process and constant interbank flows c ij = a for i, j = 1, . . . , N , where a ∈ {0, 0.5, 1}. We observe that interbank flows help stabilize the system and reduce the probability of default.

Optimal central bank policy
In this section, we assume that the central bank has to choose the interest rate r in an optimal way, so that, after banks make their choice as in the previous section, optimal policy choice is achieved. We assume that banks make optimal (for them) choices and we omit all asterisks from notation of processes. This can be thought of as a principal's problem within the principal-agent problem framework. Let us now revisit the description of policymaking by the central bank.
Its tool is the interest rate r, which the central bank uses to control the overall amount of capital in the system, measured by the Y * from (3.17). If the interest rate is low, the growth rate g(r) from (3.18) and the volatility ρ 2 (r) from (3.19) are large. A more risk-averse central bank can choose therefore a larger r. One can apply a concave utility function to Y * (t), and solve the stochastic control problem for this r. In this section, we apply the exponential (CARA) utility function to Y * (t). The private banks wish to maximize their expected logarithmic net worth Y i (t). In other words, they have logarithmic utility U (x) := log x, which shows their aversion to risk. However, in terms of logarithmic capital, their utility function is linear. Now, if the central bank was as risk-averse as private banks, she would also try to maximize E(Y 1 (T ) + . . . + Y N (T )), or, alternatively, EY (T ), for a time horizon T > 0. Below, we show that the central bank would then choose zero interest rate r = 0, because this would produce the same result as the private banks were aiming for. Now, suppose the central bank is even more risk-averse than private banks. This should manifest itself in the utility function being concave (rather than linear) even in logarithmic terms. Consider, for example, a commonly used exponential (CARA) utility function: Assume the central bank maximizes expected terminal utility: where the supremum in (4.2) is chosen over all bounded adapted controls r. We can alternatively choose instead of (4.1) the utility function as There is no difference between (4.1) and (4.3) when we try to maximize (4.2), but writing (4.3) highlights the risk-aversion of the central bank. As λ ↓ 0, the function U λ from (4.3) satisfies: The commonly used absolute risk aversion is calculated for (4.3) as follows: In other words, λ > 0 is the coefficient of risk aversion (of the central bank relative to private banks). For λ = 0 the central bank is not risk-averse at all (at least not more than private banks).
Theorem 4.1. An optimal interest rate r(t) for the problem (4.2) is given by a constant r = r * which maximizes the following expression: Remark 4.1. It is interesting to note that the optimal interest rate r does not depend on the flow rates c ij . This is because we measure the size of the system by the stochastic process Y (t). This process satisfies a stochastic differential equation with coefficients independent of c ij . These coefficients do depend on the optimal controls α * i . However, as we mentioned in Remark 3.2, these optimal controls α * i , in turn, do not depend on the flow rates, because of our special choice of logarithmic utility function.
Proof. The HJB equation for the function where the supremum is taken over all bounded adapted controls r = (r(t), 0 ≤ t ≤ T ), takes the form with terminal condition Φ(T, y) = U λ (y). Try the following form: From (4.6), we can calculate derivatives with respect to t and y: Plug (4.7) into (4.5). Because Φ < 0, we can rewrite (4.5) as This, in turn, is equivalent to Since we have Φ(T, y) < 0 and U λ (y) < 0, for compatibility we need to show that f (t) > 0 for all t. From the terminal condition Φ(T, y) = U λ (y) combined with (4.6), we have: f (T ) = 1. The equation (4.8) can be written as f (t) = λk 0 f (t), which gives us f (t) = exp (λk 0 (t − T )). Therefore, f (t) is positive. Finally, let us do the verification argument to complete the proof. The idea is similar to the verification argument in Theorem 3.1. Assume r * = (r * (t), 0 ≤ t ≤ T ) is our constant control from (4.4), found from (4.5), and r = (r(t), 0 ≤ t ≤ T ) is some other admissible (adapted bounded) control. Apply the function Φ(t, ·) to the process Y . By Itô's formula, we get: (4.9) Comparing (4.5) with (4.9), we get that Φ(t, Y (t)) is a supermartingale for the control r, but a martingale for the control r * . Indeed, by boundedness of r(t), the expectation of the stochastic integral in (4.9) is zero. Since Φ(T, y) = U λ (y), we get: , with equality for the control r * . The result immediately follows from here.
Let us find the r which corresponds to the maximum in the right-hand side of (4.8). This depends on the structure of the vector g and the matrix A.
If µ i ≤ σ 2 i for all i = 1, . . . , N , then all investments are too unprofitable to borrow money for them. Then the interest rate policy cannot influence the behavior of private banks. This corresponds to the case of the liquidity trap, when conventional monetary policy no longer works. From now on until the end of this section, let us assume that all investments are attractive: . . , N. Case 1. Assume S 1 = . . . = S N : all investments are the same. Then we have: . . = g n =: g, and σ 1 = . . . = σ N =: σ; and at any r ≥ µ − σ 2 for λ > λ * . This has the following meaning: the case λ < λ * corresponds to less risk-averse central bank, and in order to increase the total quantity of capital in the system, it wishes to slash the interest rate to zero. For the case λ > λ * , however, the central bank is very risk-averse, and it increases the interest rate to prevent excessive borrowing and overheating of the financial system. Case 2. Independent portfolio process: a ij = σ 2 i δ ij . Then This function attains maximum: In the general case, we do not have an explicit form for the optimal r in case λ ∈ [λ min , λ max ]. But in case µ 1 = . . . = µ N = µ and σ 1 = . . . = σ N = σ, we have λ min = λ max . Note that here the central bank chooses expansionary monetary policy (zero interest rate r = 0) for larger values of λ than in Case 1. This has the following interpretation: If the portfolios of banks are independent, then this creates diversification in the system and reduces risk. Therefore, even a relatively risk-averse central bank (large λ) can pursue aggressive expansionary monetary policy. Case 3. Correlated portfolio process with same growth rates µ = µ i and volatilities σ 2 = σ 2 i . Assume the driving Brownian motions of these portfolio process are correlated as in (2.5). After calculation , we get: Then we can find optimal r: this is Note that for ρ 0 = 1 we get Case 1, and for ρ 0 = 0 we get Case 2. Case 3 is intermediate: there is diversification between portfolios of private banks, but this diversification is not complete. Therefore, a risk-averse central bank can pursue more expansionary monetary policy than in Case 1, but less so than in Case 2.
To illustrate the impact of risk aversion λ on the optimal interest rate, we simulate three scenarios. First, in Figure 9 we assume uncorrelated portfolio process, each with same mean and volatility µ i = σ i = 0.1 for i = 1, . . . , N . This is Case 1, which is discussed above in this section.
Next, in Figure 10 we assume independent portfolio process but with mean and standard deviation µ i , σ i , i = 1, . . . , N i.i.d uniform on [0.1, 0.2]. This is Case 2, which is discussed above in this section. However, to our surprise, we observe the optimal interest rate to have only one jump as we increase the risk aversion parameter λ.
Finally, in Figure 11 we assume correlated portfolio process with ρ 0 = 0.8 and mean and standard deviation µ i , σ i , i = 1, . . . , N drawn from i.i.d uniform [0.1, 0.2]. This is a generalized version of Case 3 discussed above. We observe that due to correlation in the portfolio process, even a relatively less risk averse central bank is forced to raise the interest rate.

Long-term stability
In this section, we analyze the long-term behavior of the centered process: It takes values in the hyperplane Π := {y ∈ R N | y 1 + . . . + y N = 0}. In other words, we are trying to find whether log capitals of banks stay together as time t goes to infinity, or they split into two or more "clouds". A similar problem was posed in [1] and solved in [2] for rankbased models of financial markets. In that paper, log capitalizations of stocks are modeled as Brownian particles with drift and diffusion coefficients depending only on the current rank of this particle relative to other particles. This class of systems, known also as first-order models or competing Brownian particles, was a subject of much recent research. See also a related paper [25], with Lévy instead of Brownian particles.
The key parameters in our case are rates c ij of interbank cash flows. It turns out that, under certain fairly general conditions on these rates, the process (5.1) is ergodic: It has a unique statonary distribution; and for any initial conditions, it converges to this distribution as t → ∞. This section has two results. Theorem 5.1 deals with the case of flow rates c ij being time-independent: c ij (t) ≡ c ij . Lemma 5.3 covers the general case.
Assume the central bank has already chosen the interest rate r = r * , as above. Then Y is a Brownian motion with drift coefficient g(r * ) and diffusion coefficient ρ 2 (r * ). We have: Here, the following process: is an N -dimensional Brownian motion with drift vector and covariance matrix µ * = (µ * 1 , . . . , µ * N ), A * = (a * ij ) i,j=1,...,N from (3.14) and (3.15). The centered process (5.1) satisfies the SDE (a)Ỹ has a unique stationary distribution π on Π, which is a multivariate normal distribution; (b) the transition function satisfies for some constants c, λ, k > 0: (c) for any bounded measurable function f : Π → R we have, a.s.: Proof. From the properties of solutions of SDE and nondegeneracy of the covariance matrixÃ * of M , we have the following positivity property: The generator ofỸ for all twice continuously differentiable functions f : Π → R is given by: Here, M = (m ij ) i,j=1,...,N is the following matrix: Now, plug this function V from (5.6) for a suitable λ into the generator (5.8). Then Combining (5.10) with (5.8), we get: Using Lemma 5.2 below, we get: Using (5.12), the estimate x Ã * x ≤ a 0 x 2 for all x ∈ Π and a constant a 0 > 0, we get from (5.11): Choose λ := c 0 /a 0 , then (5.13) takes the form (5.14) Note that, as x → ∞, we have: K(x) → −∞. Therefore, for some constants c 1 , c 2 > 0, Recall the definition of the ball B(c 2 ) in (1.7). Since LV and V are continuous, we have: Combining (5.14) with (5.15) and (5.16), we get: Finally, combine (5.17) with the Feller property ofỸ (i.e. for a bounded continuous function f, the map x → P t f (x) is also bounded and continuous for the transition function P ofỸ ), and with the positivity property (5.7). Apply Lemma 6.1 from Appendix to Lebesgue reference measure ψ and the function V from (5.6). This completes the proof of (a) (the uniqueness of a stationary distribution), as well as of (b). The fact that this stationary distribution π is multivariate normal follows from the observation thatỸ is a multidimensional Ornstein-Uhlenbeck process on the hyperplane Π. To Proof. Note that M is a generator matrix for a continuous-time Markov chain Q = (Q(t), t ≥ 0) on {1, . . . , N }. This Markov chain can be viewed as a biased random walk on the graph G: As it wants to jump out of a state i, it chooses one of its nearest neighbors j, such that i and j are connected, only not with uniform probability. This graph G is connected. Therefore, this Markov chain is irreducible. Since it is finite, it is positive recurrent. From the standard results on continuous-time Markov chains, see for example [16,Theorem 2.7.15], this Markov chain Q has a unique stationary distribution This stationary distribution satisfies π Q M = 0. But the columns of the matrix M sum up to zero. Therefore, e M = 0, and e/N can serve as a stationary distribution. By uniqueness, π Q = e/N . Next, let λ 1 , . . . λ N and v 1 , . . . , v N be the eigenvalues and eigenvectors of the matrix M. In other words, we have: The eigenvectors of the matrix M are all real, because M is symmetric. Next, nonzero eigenvalues are negative: This follows from [15,Exercise 8.1]. Each zero eigenvalue λ i has eigenvector v i which satisfies v i M = 0, that is, v i is proportional to a stationary distribution. But the stationary distribution is unique, so we have (without loss of generality): Now, take an x ∈ R N . Assume v 1 , . . . , v N are normalized: v i = 1, i = 1, . . . , N . Because M is symmetric, v 1 , . . . , v N form an orthonormal basis in R N . Therefore, we can decompose For a vector x ∈ Π, we have: x · e = 0, and therefore x · v 1 = 0. Thus, (5.20) takes the form Apply the matrix M to this vector in (5.21) and use (5.19). We have: Take the dot product of (5.21) and (5.22). Because v 1 , . . . , v N form an orthonormal basis of R N , we have: In addition, multiplying (5.21) by itself, we get: Let c(M) := min (|λ 2 |, . . . , |λ N |) > 0. Comparing (5.23) and (5.24), we get (5.18).
Lemma 5.3. Assume the flow rates are given by and c ij are real numbers as in Theorem 5.1. Then the conclusion of Theorem 5.1 is the same, minus the conclusion that π is multivariate normal.
Proof. Similar to Theorem 5.1, but with the following changes: Instead of (5.11), we have: There exist c 4 , c 5 > 0 such that f (x) ≥ c 4 for x ∈ Π, x ≥ c 5 . Therefore, for such x, the estimate (5.14) is preserved with c 0 changed to c 0 c 4 . The rest of the proof is similar to that of Theorem 5.1.
Note, however, that if the graph G is disconnected, then this stability breaks down. Indeed, assume G has connected components G 1 and G 2 (only two for sake of notational simplicity; analysis is the same for more than two connected components), and the flow rates c ij are positive constants if i and j are adjacent, c ij = 0 if not. By Theorem 5.1, we get: are ergodic, that is, they satisfy an inequality similar to (5.6). Here, But these averages from (5.25) are, in fact, Brownian motions with certain drift and diffusion coefficients, which are easy to calculate from (5.2). They are correlated, but not perfectly. Therefore, Y 1 − Y 2 is not ergodic, and the processỸ defined in (5.1) is also not ergodic. Private banks are separated into two groups, which "drift" from each other.

Concluding remarks
We studied a model of N private banks exchanging money through interbank flows, borrowing from the general economy under an interest rate set by the central bank, and investing in portfolios consisting of risky assets; these portfolios are modeled by correlated geometric Brownian motions. This represents an enhancement of the model (1.1), which is obtained in [7] as a result of banks borrowing from each other. We generalize the interbank flows from [7], making them heterogeneous.
Each private bank maximizes its expected terminal logarithmic utility. The central bank maximizes CARA (exponential) utility function of the total size of the system. We are able to solve the control problems for each private banks and the central bank because of this special choice of utility functions. The resulting dynamics looks a bit like (1.1), except that each private bank has its own growth rate and volatility in the driving Brownian motion, and the flow rates c ij depend on i and j.
For future research, one can consider the case when some but not all portfolios S i satisfy µ i ≤ σ 2 i (and are therefore unprofitable). It might be interesting to consider different utility functions for private banks, for example power utility. Since the corresponding HJB equations are likely to be intractable, the resulting problem might be analyzed using mean-field formulation, when each bank is competing against the "mass of banks". Lemma 6.1. Take a Feller continuous strong Markov process X = (X(t), t ≥ 0) on the metric state space E, with transition function P t (x, ·), and generator L. Denote by P x the probability measure under which X(0) = x. Assume for some positive reference measure ψ and a function V : E → [1, ∞) in the domain D(L) of the generator L, we have: (a) for some compact subset C ⊆ E, we have ψ(C) > 0; (b) for all ψ-positive subsets A ⊆ E, x ∈ E, t > 0, we have: P t (x, A) > 0.
The following Strong Law of Large Numbers is taken from [17, Theorem 4.1, Theorem 4.2]. It holds under an assumption which can be called uniform positive recurrence, and which can be deduced from existence of Lyapunov functions. Assume that E = R d above, and X is the solution of an SDE with a certain drift vector, and the covariance matrix A(·). Let τ C := inf{t ≥ 0 | X(t) ∈ C} be the hitting time of a subset C ⊆ R d . Assume there exists a unique stationary distribution π. Then P x -a.s. for every x ∈ R d and bounded measurable function f : R d → R, we have: Lemma 6.3. Fix µ ∈ R and σ > 0. Take a function h : R → R, defined as Its global maximum is reached at the point x * and is equal to h * = h(x * ), where: Proof. We can write h(x) = µx − σ 2 2 x 2 − r(x − 1), x ≥ 1; µx − σ 2 2 x 2 , x ≤ 1.