Numerical simulations for Ample Liquidity¶

This notebook compares FCEs to Ample Reserve equilibria via numerical simulations. See Quintin (2026) for the details of the story.

In [1]:
# Usual suspects, whether we need them or not

import numpy as np 
import matplotlib.pyplot as plt 
# import scipy.stats as stats
# import datetime

# from scipy.linalg import expm 
# from scipy.stats import lognorm
from scipy.optimize import fsolve
#from scipy.integrate import quad
# import scipy.optimize as opt
# from scipy.stats import gaussian_kde

# import pandas as pd 
import math
In [2]:
# Declare parameters

alpha=0.5 # curvature of production function
A=10 # scale parameter
delta=1;

plow=0.3
# plow=0.2
phigh=0.7 # range for aggregate shocks
# assume uniform to start

w=4 # outside option

K=10 # max capacity

eta=.1 # fraction of agents that get a signal

# Supply of funds is isoelastic S=s0 r^theta

s0=2
theta=1 # start with linear schedule, why not

# assume phi is uniform over [0,1] to start, that's f
# lest's show fhat for p=0.5
In [3]:
# start computing equilibria

# get a pgrid going

pgrid=np.linspace(plow,phigh, 500)
# pgrid=pgrid/(pgrid+(1-pgrid)*(1-eta))
kgrid=np.linspace(0,K, 50)
know=(A*pgrid*alpha)**(1/(1-alpha)) # will need to eventually check that this is less that K

# now we need to find the market clearing interest rate

rguess=np.zeros(np.size(pgrid))

# Define unconstrained output given k and p

def Hnow(k,p):

    kstar=(A*p*alpha/delta)**(1/(1-alpha)) # optimal scale given delta
    
    if k> kstar:
        return A*p*kstar**(alpha) + delta*(k-kstar)
    else:
        return A*p*k**alpha

# Define the first derivative of unconstrained output given k and p

def H1now(k,p):

    kstar=(A*p*alpha/delta)**(1/(1-alpha)) # optimal scale given delta
    
    if k> kstar:
        return delta
    else:
        return A*p*alpha*k**(alpha-1)
        
# compute the optimal saving (should combine with above, but this is cheap)

def Snow(k,p):

    kstar=(A*p*alpha/delta)**(1/(1-alpha)) # optimal scale given delta
    
    if k> kstar:
        return k-kstar
    else:
        return 0

fig, axs  = plt.subplots(1,2,figsize=(12, 8))

axs[0].plot(np.linspace(0,K, 50),[Hnow(k,0.5) for k in kgrid], label=r"$H(k,0.5)$")
# axs[0].plot(np.linspace(0,K, 50),[H1now(k,0.5) for k in kgrid], label=r"$H_1(k,0.5)$")
axs[0].plot(np.linspace(0,K, 50),[Snow(k,0.5) for k in kgrid], label=r"$s(k,0.5)$")
axs[0].set_xlabel(r'$k$')
axs[0].legend()
axs[0]

axs[1].plot(np.linspace(plow,phigh, 500),[Hnow(5,p) for p in pgrid], label=r"$H(5,p)$")
# axs[1].plot(np.linspace(plow,phigh, 500),[H1now(5,p) for p in pgrid], label=r"$H_1(k,0.5)$")
axs[1].plot(np.linspace(plow,phigh, 500),[Snow(5,p) for p in pgrid], label=r"$s(5,p)$")
axs[1].set_xlabel(r'$p$')
axs[1].legend()
Out[3]:
<matplotlib.legend.Legend at 0x230326cfed0>
No description has been provided for this image
In [4]:
# start searching for FCEs

# record output and divestment per p

# Step1: guess lambda (eventually update so that participation constraint binds)

lam=4.082

# step 2: start looking for the average r that clears markets
# dichotomy is your always safe best friend here, slow though it may be

rbarguesslow=0 # I'm looking for solutions for average rates in the [0,0.2] range, this can always be widened if binding
rbarguesshigh=0.2

itbar=0;

while itbar<20: # precision of 2^-20 is quite enough 
    
    rbar=(rbarguesslow+rbarguesshigh)/2


# Step 2: find k
# Here we will assume that the distributions of phi and p are uniform, so that all integrals are means
# we will vary that in different experiments

# we look for k via dichotomy since the problem is well behaved
    
    klow=0
    khigh=K # here the bound is not arbitrary but a fundamental part of the model
    # this allow for k=K as a solution, which must in fact happen for A high enough

    it=1
    while it<20:

        kguess=(klow+khigh)/2
    
        devnow=np.array([H1now(kguess,p) for p in pgrid])
        intnow=devnow.mean() # this is the marginal product of capital at t=0

        if intnow>1+rbar:
            klow=kguess
        else:
            khigh=kguess
    
        it+=1

# next we need to compute the rhat schedule which, given lambda and k, clears markets
# again we dichotomize because it's safe, no numerical snafus possible here

    rhat=np.zeros(np.size(pgrid))

    pcount=0
    for p in pgrid:

        philownow=eta*(1-p) 
        phibarnow=(1+philownow)/2 # this is the expected arrival rate given eta and p 

        rlow=0
        rhigh=4*rbar # arbitrary bound, we can loosen if needed
    
        it=1
        while it<20:

            rguess=(rlow+rhigh)/2       
    
            ednow=phibarnow*lam/(1+rguess)+kguess-K-s0*rguess**theta # this excess demand for liquidity
            # we need that to be zero for all p in FCE

            if ednow>0:
                rlow=rguess
            else:
                rhigh=rguess
    
            it+=1

        rhat[pcount]=rguess
        pcount+=1

    if rhat.mean()>rbar: # the average rate must eventually match the initial guess
        rbarguesslow=rbar
    else:
        rbarguesshigh=rbar

    itbar+=1

yhat=np.array([Hnow(kguess,p) for p in pgrid])
shat=np.array([Snow(kguess,p) for p in pgrid]) # optimal divestment policy given p in equilibrium

kFCE=kguess # capital formation in FCE
In [5]:
# Now we go after an ample liquidity equilibrium 
# First we set r* to the average rate in FCE
# we don't need to do that but that's certainly a natural r* for central bankers to consider

rstar=rhat.mean()

# we find capital formation given r*
# we look for k via dichotomy since the problem is well behaved
    
klow=0
khigh=K # bounds imposed by the model

it=1
while it<20:

    kguess=(klow+khigh)/2
    
    devnow=H1now(kguess,pgrid.mean()) # this is the marginal product of capital given the average p

    if devnow>1+rstar:
        klow=kguess
    else:
        khigh=kguess
    
    it+=1

kample=kguess # capital formation at the ample solution

# now we compute output in FCE

HFCE=np.zeros(np.size(pgrid))

pcount=0

while pcount<np.size(pgrid):
    
    pnow=pgrid[pcount]
    HFCE[pcount]=A*pnow*kample**alpha
    pcount+=1
In [6]:
import random

# Now we back out the w that makes all of it work
# via MC so as not to have to integrate

ndraws=500
draws = [random.randint(0, np.size(pgrid)-1) for _ in range(ndraws)]

# f i uniform on [0,1]
# lower bound on phihat
philow=[eta*(1-p) for p in pgrid] 
meanphi=[(1+phil)/2 for phil in philow]


Utils=np.array([meanphi[p]*math.log(lam/(1+rhat[p]))+ (1-meanphi[p])*math.log(lam) for p in draws])

w=math.exp(Utils.mean()) # we need this w to match the parameter choice
# update lambda above to get this to what it needs to be
# could embed this in the equilibrium loop
# but a few iterations suffice so manually does it fast too
In [7]:
w # we need w to be near 4, then we know that the lambda above is correct
Out[7]:
3.9993265595507324
In [8]:
# bayesian updating in Ample Liquidity equilibrium

phigrid=np.linspace((1-phigh)*eta,1,500)
# pgrid=np.linspace(0,1,50)
gpgivenphi=np.zeros((np.size(phigrid),np.size(pgrid)))

phicount=0

while phicount<np.size(phigrid):

    pcount=0
    while pcount<np.size(pgrid):

        if phigrid[phicount]>eta*(1-plow):
            kappa=eta/(np.log(1-eta*(1-phigh))-np.log(1-eta*(1-plow)))
        else:
            kappa=eta/(np.log(1-eta*(1-phigh))-np.log(1-phigrid[phicount]))
        

        if pgrid[pcount]>1-phigrid[phicount]/eta:
            gpgivenphi[phicount,pcount]=kappa/(1-eta*(1-pgrid[pcount]))
        else:
            gpgivenphi[phicount,pcount]=0

        pcount+=1

    phicount+=1

# now we can calculate pbar given phi


pbarphi=np.array([np.dot(pgrid,gpgivenphi[phicount,:])/np.sum(gpgivenphi[phicount,:]) for phicount in range(np.size(phigrid))])

C:\Users\quintin\AppData\Local\Temp\ipykernel_31624\2726994044.py:17: RuntimeWarning: divide by zero encountered in scalar divide
  kappa=eta/(np.log(1-eta*(1-phigh))-np.log(1-phigrid[phicount]))
C:\Users\quintin\AppData\Local\Temp\ipykernel_31624\2726994044.py:32: RuntimeWarning: invalid value encountered in scalar divide
  pbarphi=np.array([np.dot(pgrid,gpgivenphi[phicount,:])/np.sum(gpgivenphi[phicount,:]) for phicount in range(np.size(phigrid))])
In [17]:
fig, axs  = plt.subplots(2,1,figsize=(12, 8))

axs[0].set_title("Posterior densities given arrival rates")
axs[0].plot(pgrid,gpgivenphi[9,:], label=r'$\hat{\phi}=0.02$')
axs[0].plot(pgrid,gpgivenphi[18,:], label=r'$\hat{\phi}=0.04$')
axs[0].plot(pgrid,gpgivenphi[90,:], label=r'$\hat{\phi}>\eta (1-p_L)$')
axs[0].legend()
axs[1].set_title("Expected aggregate shock given arrival rates")
axs[1].plot(phigrid,pbarphi)
Out[17]:
[<matplotlib.lines.Line2D at 0x23037711f90>]
No description has been provided for this image
In [10]:
# now we can compute the ample liquidity solution

rstar=rhat.mean()

# we find capital formation given r*
# we look for k via dichotomy since the problem is well behaved
    
klow=0
khigh=K # bounds imposed by the model

it=1
while it<20:

    kguess=(klow+khigh)/2
    
    # devnow=np.array([H1now(kguess,p) for p in pbarphi]).mean() # this is the marginal produce of capital given the average p
    devnow=np.nanmean(np.array([H1now(kguess,p) for p in pbarphi]))

    if devnow>1+rstar:
        klow=kguess
    else:
        khigh=kguess
    
    it+=1

kample=kguess # capital formation at the ample solution

# now we compute output in ample liquidity equilibrium

Hample=np.zeros((np.size(pgrid),np.size(phigrid)))
Sample=np.zeros((np.size(pgrid),np.size(phigrid)))

pcount=0
while pcount<np.size(pgrid):
    phicount=0
    while phicount<np.size(phigrid):
        snow=Snow(kample,pbarphi[phicount])
        pnow=pgrid[pcount]
        Hample[pcount,phicount]=A*pnow*(kample-snow)**alpha + snow*delta
        Sample[pcount,phicount]=snow
        phicount+=1
    pcount+=1

Hamplemean=[Hample[pcount,:].mean() for pcount in range(np.size(pgrid))] # average output across locations given a shock
Samplemean=[Sample[pcount,:].mean() for pcount in range(np.size(pgrid))] # average divestment across locations given a shock
In [11]:
kample
Out[11]:
5.842342376708984
In [12]:
# start showing some time series

import random
ndraws=100 

draws = [random.randint(0, np.size(pgrid)-1) for _ in range(ndraws)] # draw 50 aggregate shocks, iid draws, with replacement

fig, axs  = plt.subplots(2,2,figsize=(12, 8))

axs[0][0].plot(range(ndraws),[pgrid[p] for p in draws])
axs[0][0].set_title(r"Aggregate shock ($p$)")
axs[0][1].plot(range(ndraws),[100*rhat[p] for p in draws], label='FCE')
axs[0][1].plot(range(ndraws),[100*rstar for p in draws],label='Ample liquidity')
axs[0][1].set_title(r"Interest rate $(\hat{r})$, in percent")
axs[0][1].set_ylabel(r"$\%$")
axs[0][1].legend()
axs[1][0].plot(range(ndraws),[yhat[p] for p in draws], label='FCE')
axs[1][0].plot(range(ndraws),[Hamplemean[p] for p in draws], label='Ample liquidity')

axs[1][0].legend()
axs[1][0].set_title(r"Average output across locations ($H$)")
axs[1][1].plot(range(ndraws),[shat[p] for p in draws], label='FCE')
axs[1][1].plot(range(ndraws),[Samplemean[p] for p in draws], label='Ample liquidity')
axs[1][1].set_title("Average divestment across locations ($s$)")
axs[1][1].legend()
Out[12]:
<matplotlib.legend.Legend at 0x23036d347d0>
No description has been provided for this image
In [13]:
# plot the rate schedule

plt.plot(pgrid,100*rhat,label="Rate schedule in FCE")
plt.plot(pgrid,[100*rstar for p in pgrid],label="Rate schedule in Ample Liquidity equilibrium")
plt.ylabel(r"%")
plt.xlabel(r"$p$")
plt.legend()
Out[13]:
<matplotlib.legend.Legend at 0x2303723ec10>
No description has been provided for this image
In [14]:
kample
Out[14]:
5.842342376708984
In [15]:
kFCE
Out[15]:
8.017864227294922
In [16]:
kample/kFCE
Out[16]:
0.7286656659537976
In [17]:
np.array([yhat[p] for p in draws]).std()
Out[17]:
np.float64(2.815959042605213)
In [18]:
np.array([Hamplemean[p] for p in draws]).std()/np.array([yhat[p] for p in draws]).std()
Out[18]:
np.float64(1.0009732341002573)
In [ ]: