Original paper found here: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3983535
Last updated: 3/13/2024
Code requires Julia version 1.9.2 and the latest version of all packages to function.
using LinearAlgebra
using SymPy
using Pkg
using Plots
using Distributions
using DataFrames
using CSV
using Combinatorics
Creates the complete, ring, gamma-convex, and delta-connected networks
# Create different networks in matrix form
function complete_network(Dweight,n)
# identity matrix n x n
eye = Diagonal(ones(n,n))-zeros(n,n);
# base debt weight matrix
D = ones(n,n)-eye;
# complete network
Dc = D.*Dweight/(n-1);
return Dc
end
function ring_network(Dweight,n)
# eye matrix
eye = Diagonal(ones(n,n))-zeros(n,n);
# ring network
lc = eye[[n],:];
mc = eye[1:n-1,:];
Dr = Dweight*vcat(lc,mc);
return Dr
end
function gamma_network(Dweight,n,gamma)
Dc = complete_network(Dweight,n)
Dr = ring_network(Dweight,n)
# gamma-convex of ring and complete
Dg = (gamma*Dc)+((1-gamma)*Dr);
end
# delta network is assumed to be 2 complete network components
function delta_network(Dweight,n,delta)
r = Int(n/2);
weight = delta*Dweight
main_weight = Dweight-weight*r
Dp = complete_network(main_weight,r)
Z = weight*ones(r,r);
Dd = [Dp Z ; Z Dp]
return Dd
end
function gamma_delta_network(Dweight,n,delta,gamma)
r = Int(n/2);
weight = delta*Dweight
main_weight = Dweight-weight*r
Dp = gamma_network(main_weight,r,gamma)
Z = weight*ones(r,r);
Dd = [Dp Z ; Z Dp]
return Dd
end
function hetero_delta_network(Dweight,n,r,delta,gamma)
weight = delta*Dweight # given to other cluster
main_weight = Dweight-weight*(n-r)
main_weight2 = Dweight-weight*r
# big cluster
m1 = gamma_network(main_weight,r,1)
m2 = ones(n-r,r)*weight
m3 = ones(r,n-r)*weight
# little cluster
m4 = gamma_network(main_weight2,n-r,1)
Dd = [m1 m3 ; m2 m4]
return Dd
end
# helper function for making weights according to the weighting rule
function weights(D,C,p)
# first create an nxn matrix of excess collateral values
EC = D-C.*D.*p;
# replace values with 0 if value is negative
EC[EC.<0].=0;
# create total collateral per agent matrix values for denominator
den = sum(EC,dims=1);
#den[den.<0].=0;
# create weights
Q = EC./den;
# replace undefined values with 1
Q[(isnan.(Q)).|(Q.==Inf).|(Q.==-Inf)] .= 1;
return Q
end
# solves for equilibrium x payment values
function x_solve(d_hat,z_hat,Q,xi,p,s,zeta)
x_hat = d_hat;
check=true;
while check
#calculate l_hat
if p >=s*zeta
l_hat = max.(zeros(n,1),min.(1/zeta*(d_hat-Q*x_hat-z_hat),xi.*ones(n,1)));
else
l_hat = xi.*ones(n,1);
end
x_hat2 = max.(zeros(n,1),min.(Q*x_hat+z_hat+zeta*l_hat,d_hat)); # calculate x_hat directly, which is the sum of all agents' payments
#check if x_hat = x_hat2 with some margin, if not then we continue
#loop and reassign x_hat to the minimum statement
test = sum((x_hat-x_hat2).^2);
if test <= 1e-10
check = false;
end
x_hat = x_hat2;
end
return x_hat
end
# calculates c underbar
function c_underbar(Dweight,n,ei,s,h0_w,k)
output = (Dweight-((n-k)/k)*ei)/(Dweight*s)+h0_w/Dweight;
return output
end
# calculates c upperbar
function c_upperbar(s,n,k,ei,h0_w)
output = (min(s,(n-k)*ei/(k*h0_w)))^-1;
return output
end
c_upperbar (generic function with 1 method)
Definition 1: For given $ (N,C,D,e,h, s, \omega) $, if liquidation decisions $ \{l_j(p)\} $ satisfy the liquidation rule, payments $ \{x_{ij}(p)\} $ satisfy the payment rule, $ \{m_j(p)\} $ is determined by net wealth equation, fire sale amount $\{\phi_j(p)\}$ determined by the fire-sale equation, and price $ p $ clears the market according to the market clearing condition then $ (\{x_{ij}\},\{l_j\},\{m_j\}, \{\phi_j\},p) $ is a full equilibrium.
This function identifies the full equilibrium as defined in definition 1 and returns $p$, total liquidations, and the number of defaults. Agent $j$ defaults when $m_j \leq 0$.
function equ(s,n,e0,epsilon,xi,zeta,C,D,omega,h0)
# check for while loop to end when equilibrium is found
look = true;
# start with p equalling fair value s
p=s;
# assign increment for when loop decreases p
inc = 0.01;
# stop after 1 iteration
stop = 1
while look
# Collateral matrix in monetary terms (n x n)
cdp = C.*D.*p;
# Collateral matrix in non-monetary terms (n x n)
cd = C.*D;
# Debt-Collateral matrix (owed net amount) (n x n)
d_cdp = D-cdp;
# Weighting matrix, pro-rata (n x n)
Q = weights(D,C,p);
# total debt owed per agent w/ collateral netting (n x 1)
d_cdp[d_cdp.<0].=0;
d_hat = d_cdp;
d_hat = transpose(sum(d_hat,dims=1));
# collateral matrix in monetary terms (n x n)
cp = C*p;
# establish indicator matrix denoting which collateral are greater than 1 or less than 1 (n x n)
b1 = cp.>1;
b2 = cp.<=1;
# construct cash, asset and z_hat vectors (n x 1)
e_hat = e0 + (D.*b1)*ones(n,1) - (transpose(D.*b1))*ones(n,1);
h_hat = h0 + (cd.*b2)*ones(n,1) - (transpose(cd.*b2))*ones(n,1);
z_hat = e_hat + h_hat*p - omega*epsilon;
# solve for x payments
x = x_solve(d_hat,z_hat,Q,xi,p,s,zeta)
# indifference threshold between holding asset and liquidating
thres_p = s*zeta
# plug in x to solve for l_hat in numeric terms
if p >=thres_p
A = 1/zeta*(d_hat-Q*x-z_hat);
B = xi.*ones(n,1);
term = min.(A,B);
l_hat = max.(term,zeros(n,1)); # maximum term
else
l_hat = xi.*ones(n,1);
end
# use threshold (tol = 1.1102e-16) to ensure cases like boolean 18.6 == 18.600000 are
# recognized to equal 1 and not 0
tol = eps(0.5);
d_solvent = (abs.(x-d_hat).<tol);
d_default = ones(n,1) - d_solvent;
# full payment matrix of all agents, sum of defaulting and solvent agents' payments (n x n)
all_X = d_cdp.*(transpose(d_solvent))+ Q.*(transpose(x).*transpose(d_default));
# agents' wealth (n x 1)
w = e_hat + h_hat*p - omega*epsilon + zeta*l_hat + all_X*ones(n,1) - x;
# net wealth
total_wealth = sum(max.(zeros(n,1),w))
# total # of liquidations
l = length(l_hat[l_hat.>0]);
# total value of liquidations
tot_liq = round(sum(l_hat),digits=2);
# sum of non-defaulting agents' assets (s)
nd_assets = 0;
# of defaulting agents
def = 0;
# second term for the numerator in price equation
num = 0;
# count the number of defaults
for i= 1:n
if w[i,1] > 1e-10
nd_assets = nd_assets+h0[i,1];
num = num+max(0,h0[i,1]*p-w[i,1]);
else
def = def+1;
end
end
# check if price can be supported at the threshold otherwise continue to regular market clearing condition
if p == thres_p
# calculate the remaining supply of available liquidations
full_l = xi*ones(n,1)
remainder_l = sum((full_l-l_hat))
remainder_l_zeta =zeta*remainder_l
# calculate wealth discrepancy to support price
wealth_needed = thres_p*sum(h0)
remainder_w = wealth_needed-total_wealth
# if agents can liquidate a little more to support price at the threshold, return the threshold
if remainder_l_zeta >= remainder_w
add_l = remainder_w/zeta
final_liquidation = add_l+tot_liq
# calculate social surplus involving payoffs in time t=2 which includes assets being at full value (s)
ss = sum(e0)+s*sum(h0)+n*xi-(1-zeta)*final_liquidation
return round(final_liquidation,digits=2), def,thres_p,ss
end
end
# fire sales
temp = min.(max.(h0*p-w,zeros(n,1)),h0*p);
phi = sum(temp);
# check if there's enough cash in the market to buy fire sales, if
# so then price will be fair value; else calculate liquidity
# constrained price
condition = sum(max.(zeros(n,1),w.-h0*s));
if p == s && condition>=phi
p_new = s;
else
if nd_assets == 0
p_new = 0; # everyone defaulting meaning the asset is worthless
else
p_new = total_wealth/sum(h0); # liquidity constrained price
end
end
# calculate difference between inputted p and re-calculated p
d=round(abs(p-p_new),digits=3);
# calculate social surplus involving payoffs in time t=2 which includes assets being at full value (s)
ss = sum(e0)+s*sum(h0)+n*xi-(1-zeta)*tot_liq
# check if market clearing p is the guessed p
if d <= 0.01||p_new > p
p = p_new;
# stop looking and return p
look = false;
# return liquidations, defaults, price, social surplus
return tot_liq,def,p,ss
else
# try a different p value
p = round(p-inc,digits=2);
end
# market crash if price is 0
if p < 0.0001
p = 0;
return tot_liq,def,p,ss
end
end
end
equ (generic function with 1 method)
function multi_simulator(n,e0,h0,h0_w,s,xi,zeta,omega,ei,D,Dweight,e_star)
# calculate upperbound, a little larger than epsilon star
ub = 1.5*e_star;
# set range: iterate in increments of 0.1 from 0 to a little bit above collateral upperbound
range = 0:0.1:ub;
num_x = length(range);
c_ub = c_upperbar(s,n,k,ei,h0_w)*1.1;
c_range = Array(0.01:0.01:c_ub);
total_obs = num_x*length(c_range);
liquidations = zeros(total_obs,1);
prices = zeros(total_obs,1);
social_surplus = zeros(total_obs,1);
ind = 1
# iterate through the cases of liquidity shocks and collateral ratios
for i = 1:num_x
m = range[i];
epsilon = ei*m;
for j in 1:length(c_range)
c = c_range[j];
C = ones(n,n)*c
# calculate liquidations, number of defaults, collateral price, and social surplus
tot_liq,def,p,ss = equ(s,n,e0,epsilon,xi,zeta,C,D,omega,h0);
liquidations[ind] = tot_liq;
prices[ind] = p;
social_surplus[ind] = ss;
# add to indicator
ind+=1
end
end
return Array(liquidations),Array(prices),Array(social_surplus)
end
multi_simulator (generic function with 1 method)
# denote whether liquidations are nonzero
nonzero = true
if nonzero
zeta = 0.1 # liquidation inefficiency
else
zeta = 1e-10
end
n = 20 # 20 agents
s = 1 # fair value of 1
ei = 1 # cash unit of 1
e0 = ei*ones(n,1) # vector of cash values
xi = 10 # long term project return
# set shock and debt thresholds
if nonzero
epsilon_star = ei*n+zeta*n*xi
dstar = (n-1)*(ei+zeta*xi)
else
epsilon_star = ei*n # epsilon star value
dstar = (n-1)*ei
end
# total agent liability
Dweight = 2*dstar
# construct networks
Dc = complete_network(Dweight,n)
Dr = ring_network(Dweight,n)
Dg = gamma_network(Dweight,n,0.5)
Dd = delta_network(Dweight,n,0) # delta=0 network
# hit only one agent with liqudity shock
k = 1;
omega = zeros(n,1);
for i = 1:k
omega[i,1]=1;
end
# set asset number
h0_w = 4;
h0 = h0_w*ones(n,1);
cub = c_underbar(Dweight,n,ei,s,h0_w,k)
println("c underbar:",cub)
# Epsilon ranges
println("Regular Epsilon:",ei*n)
println("Lower Epsilon:",ei*n+xi)
println("Upper Epsilon:",epsilon_star)
epsilon = 1.5*epsilon_star;
c underbar:0.8026315789473684 Regular Epsilon:20 Lower Epsilon:30 Upper Epsilon:40.0
# run simulations for each type of network structure and produce results
liquidations_r, prices_r,social_surplus_r=multi_simulator(n,e0,h0,h0_w,s,xi,zeta,omega,ei,Dr,Dweight,epsilon_star);
liquidations_c, prices_c,social_surplus_c=multi_simulator(n,e0,h0,h0_w,s,xi,zeta,omega,ei,Dc,Dweight,epsilon_star);
liquidations_g, prices_g,social_surplus_g=multi_simulator(n,e0,h0,h0_w,s,xi,zeta,omega,ei,Dg,Dweight,epsilon_star);
liquidations_d,prices_d,social_surplus_d=multi_simulator(n,e0,h0,h0_w,s,xi,zeta,omega,ei,Dd,Dweight,epsilon_star);
# can write results into csvs
mytable = DataFrame(Dict("Ring"=>vec(liquidations_r),"Complete"=> vec(liquidations_c),"Gamma"=> vec(liquidations_g),"Delta"=>vec(liquidations_d)))
CSV.write("simulation_results_liquidations.csv",mytable)
mytable = DataFrame(Dict("Ring"=>vec(prices_r),"Complete"=> vec(prices_c),"Gamma"=> vec(prices_g),"Delta"=>vec(prices_d)))
CSV.write("simulation_results_prices.csv",mytable)
mytable = DataFrame(Dict("Ring"=>vec(social_surplus_r),"Complete"=> vec(social_surplus_c),"Gamma"=> vec(social_surplus_g),"Delta"=>vec(social_surplus_d)))
CSV.write("simulation_results_social_surplus.csv",mytable)
"simulation_results_social_surplus.csv"
# read results in without having to rerun code
sim_results = CSV.read("simulation_results_liquidations.csv",DataFrame);
liquidations_c = sim_results.Complete;
liquidations_r = sim_results.Ring;
liquidations_g = sim_results.Gamma;
liquidations_d = sim_results.Delta;
sim_results2 = CSV.read("simulation_results_prices.csv",DataFrame);
prices_c = sim_results2.Complete;
prices_r = sim_results2.Ring;
prices_g = sim_results2.Gamma;
prices_d = sim_results2.Delta;
sim_results3 = CSV.read("simulation_results_social_surplus.csv",DataFrame);
social_surplus_c = sim_results3.Complete;
social_surplus_r = sim_results3.Ring;
social_surplus_g = sim_results3.Gamma;
social_surplus_d = sim_results3.Delta;
# create x (collateral ratio) and y axes (liquidity shock)
e_range = 0:0.1:epsilon;
c_ub = c_upperbar(s,n,1,ei,h0_w)*1.1;
c_range = Array(0.01:0.01:c_ub);
# set special color gradient for price
col2 = cgrad(:inferno, scale = :exp);
Plots.contourf(c_range,e_range,vec(social_surplus_c),lw=0,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Complete Network")
Plots.contourf(c_range,e_range,vec(social_surplus_r),lw=0,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Ring Network")
Plots.contourf(c_range,e_range,vec(social_surplus_g),lw=0,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Gamma-Convex Network")
# add extra strip of values to rescale results from 125 to 300
social_surplus_d2 = vec(copy(social_surplus_d));
ub = epsilon+0.1
e_range2 = 0:0.1:ub;
for i in 1:length(c_range)
push!(social_surplus_d2,125)
end
Plots.contourf(c_range,e_range2,vec(social_surplus_d2),lw=0,ylim = (0,60),xlabel = "Collateral Ratio",ylabel = "Liquidity Shock",title="Delta-Connected Network")
# save figures for paper
Plots.contourf(c_range,e_range,vec(social_surplus_c),lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12, xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi=1000)
savefig("nonzero_zeta/v2/complete_ss.png")
Plots.contourf(c_range,e_range,vec(social_surplus_r),lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12, xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi=1000)
savefig("nonzero_zeta/v2/ring_ss.png")
Plots.contourf(c_range,e_range,vec(social_surplus_g),lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12, xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi=1000)
savefig("nonzero_zeta/v2/gamma_ss.png")
Plots.contourf(c_range,e_range2,vec(social_surplus_d2),ylim = (0,60),lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi = 1000)
savefig("nonzero_zeta/v2/delta_ss.png")
"/Users/gracechuan/Documents/network contagion/nonzero_zeta/v2/delta_ss.png"
# Complete
Plots.contourf(c_range,e_range,vec(prices_c),c=col2,lw=0,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Complete Network")
# Ring
Plots.contourf(c_range,e_range,vec(prices_r), c=col2,lw=0,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Ring Network")
# Gamma-convex
Plots.contourf(c_range,e_range,vec(prices_g),c=col2,lw=0,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Gamma-Convex Network")
# add extra strip of values to rescale results from 0 to 1
prices_d2 = vec(copy(prices_d));
ub = epsilon+0.1
e_range2 = 0:0.1:ub;
for i in 1:length(c_range)
push!(prices_d2,0)
end
# Delta-connected
Plots.contourf(c_range,e_range2,vec(prices_d2),lw=0,c=col2,xlabel="Collateral Ratio",ylabel="Liquidity Shock",title="Delta-Connected Network")
# save price figures for paper
Plots.contourf(c_range,e_range,vec(prices_c),c=col2,lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12, xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi=1000)
savefig("nonzero_zeta/v2/complete_price.png")
Plots.contourf(c_range,e_range,vec(prices_r),c=col2,lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12, xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi=1000)
savefig("nonzero_zeta/v2/ring_price.png")
Plots.contourf(c_range,e_range,vec(prices_g),c=col2,lw=0,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12, xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi=1000)
savefig("nonzero_zeta/v2/gamma_price.png")
Plots.contourf(c_range,e_range2,vec(prices_d2),lw=0,c=col2,cbar = false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,xlabel="Collateral Ratio",ylabel="Liquidity Shock",dpi = 1000)
savefig("nonzero_zeta/v2/delta_price.png")
"/Users/gracechuan/Documents/network contagion/nonzero_zeta/v2/delta_price.png"
# get the axes for plotting
function axes_3d_plots(n,e0,h0,h0_w,s,xi,zeta,omega,ei,Dweight,epsilon)
# set range: iterate in increments of 0.1 from 0.9 to upperbound
range = 0:0.1:epsilon;
num_x = length(range);
c_ub = c_upperbar(s,n,k,ei,h0_w)*1.1;
c_range = Array(0.01:0.01:c_ub);
total_obs = num_x*length(c_range);
epsilons = zeros(total_obs,1);
collateral_ratios = zeros(total_obs,1);
ind = 1
for i = 1:num_x
m = range[i];
epsilon = ei*m;
for j in 1:length(c_range)
c = c_range[j];
C = ones(n,n)*c
# add epsilon and c to arrays
epsilons[ind] = epsilon
collateral_ratios[ind] = c
# add to indicator
ind+=1
end
end
return epsilons, collateral_ratios
end
epsilons, collateral_ratios = axes_3d_plots(n,e0,h0,h0_w,s,xi,zeta,omega,ei,Dweight,epsilon);
# 3d versions of plots
# Social surplus
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(social_surplus_c),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Social Surplus")
savefig("nonzero_zeta/v3/complete_ss_3d.png")
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(social_surplus_r),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Social Surplus")
savefig("nonzero_zeta/v3/ring_ss_3d.png")
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(social_surplus_g),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Social Surplus")
savefig("nonzero_zeta/v3/gamma_ss_3d.png")
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(social_surplus_d),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Social Surplus")
savefig("nonzero_zeta/v3/delta_ss_3d.png")
# Prices
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(prices_c),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Price",c=col2)
savefig("nonzero_zeta/v2/complete_price_3d.png")
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(prices_r),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Price",c=col2)
savefig("nonzero_zeta/v2/ring_price_3d.png")
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(prices_d),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Price",c=col2,zlim=(0,1))
savefig("nonzero_zeta/v2/delta_price_3d.png")
Plots.surface(vec(collateral_ratios),vec(epsilons),vec(prices_g),camera = (28,8),cbar=false,xlabel="Collateral Ratio",ylabel="Liquidity Shock",zlabel="Price",c=col2)
savefig("nonzero_zeta/v2/gamma_price_3d.png")
"/Users/gracechuan/Documents/network contagion/nonzero_zeta/v3/delta_ss_3d.png"
# denote whether liquidations are nonzero
nonzero = false
if nonzero
zeta = 0.1 # liquidation inefficiency
else
zeta = 1e-10
end
n = 20 # 20 agents
s = 1 # fair value of 1
ei = 1 # cash unit of 1
e0 = ei*ones(n,1) # vector of cash values
xi = 10 # long term project return
r = 5; # number of agents in smaller component
# set shock and debt thresholds
if nonzero
epsilon_star = ei*n+zeta*n*xi
dstar = (n-1)*(ei+zeta*xi)
else
epsilon_star = ei*n # epsilon star value
dstar = (n-1)*ei
end
# total agent liability
Dweight = 2*dstar
# construct networks
Dc = complete_network(Dweight,n)
Dr = ring_network(Dweight,n)
Dg = gamma_network(Dweight,n,0.5)
Dd = delta_network(Dweight,n,0) # delta=0 network
# hit only one agent with liqudity shock
k = 1;
omega = zeros(n,1);
for i = 1:k
omega[i,1]=1;
end
# set assets number
h0_w = 2;
h0 = h0_w*ones(n,1);
cub = c_underbar(Dweight,n,ei,s,h0_w,k)
println("c underbar:",cub)
# Epsilon ranges
println("Regular Epsilon:",ei*n)
println("Lower Epsilon:",ei*n+xi)
println("Upper Epsilon:",epsilon_star)
epsilon = 1.5*epsilon_star;
c underbar:0.5526315789473684 Regular Epsilon:20 Lower Epsilon:30 Upper Epsilon:20
# for a given delta-connected network of 2 gamma-convex (of complete and ring networks) components, find delta star (network collapse threshold)
# note: when gamma = 1, the components are simply complete, when gamma =0, components are ring
function find_delta_star(n,r,Dweight,cub,s,e0,epsilon_star,xi,zeta,omega,h0,gamma=1)
C = ones(n,n)*(cub-0.01)
epsilon = 2*epsilon_star
delta_ub = 1/(n-1)
delta_range = 0:0.0001:delta_ub
len_dr = length(delta_range)
for ind in 1:len_dr
delta = delta_range[ind]
Dld = hetero_delta_network(Dweight,n,r,delta,gamma)
tot_liq,def, p,ss = equ(s,n,e0,epsilon,xi,zeta,C,Dld,omega,h0)
if tot_liq == n*xi
return delta
end
end
return "no delta star"
end
# runs multiple simulations to plot results in figures
function graph_agg_delta(n,r,Dweight,s,e0,ei,epsilon_star,xi,zeta,omega,h0,k,gamma=1,delta_ub=0.01)
epsilon = 2*epsilon_star
c_bound = c_underbar(Dweight,n,ei,s,h0_w,k)
collateral_range = 0:0.01:(c_bound+0.15)
delta_range = 0:0.0001:delta_ub
len_dr = length(delta_range)
liquidations = []
social_surplus = []
prices = []
for c in collateral_range
C = ones(n,n)*c
for ind in 1:len_dr
delta = delta_range[ind]
Dld = hetero_delta_network(Dweight,n,r,delta,gamma)
tot_liq,def, p,ss = equ(s,n,e0,epsilon,xi,zeta,C,Dld,omega,h0)
push!(liquidations,tot_liq)
push!(prices,p)
push!(social_surplus,ss)
end
end
return social_surplus, liquidations, prices
end
graph_agg_delta (generic function with 3 methods)
# create x and y axes of contour plot
delta_range = 0:0.01:1
c_bound = cub
collateral_range = 0:0.01:(c_bound+0.15)
# create delta and collateral values used for each simulation
delta_values = []
collateral_values = []
for c in collateral_range
for ind in 1:length(delta_range)
delta = delta_range[ind] #in percentage
push!(delta_values,delta)
push!(collateral_values,c)
end
end
# create x and y axes of extended contour plot
delta_ub = 6
delta_range_extended = 0:0.01:delta_ub
# create delta and collateral values used for each simulation
delta_values_extended = []
collateral_values_extended = []
for c in collateral_range
for ind in 1:length(delta_range_extended)
delta = delta_range_extended[ind] #in percentage
push!(delta_values_extended,delta)
push!(collateral_values_extended,c)
end
end
# Run simulations
# homogeneous delta-connected network
social_surplus, liquidations, prices=graph_agg_delta(n,10,Dweight,s,e0,ei,epsilon_star,xi,zeta,omega,h0,k,1);
# heterogeneous delta-connected networks
social_surplus2, liquidations2, prices2=graph_agg_delta(n,5,Dweight,s,e0,ei,epsilon_star,xi,zeta,omega,h0,k,1);
social_surplus3, liquidations3, prices3=graph_agg_delta(n,15,Dweight,s,e0,ei,epsilon_star,xi,zeta,omega,h0,k,1);
# heterogeneous delta-connected network extended for 5 agents
social_surplus2e, liquidations2e, prices2e=graph_agg_delta(n,5,Dweight,s,e0,ei,epsilon_star,xi,zeta,omega,h0,k,1,0.06);
# if we wanted to write results into csv
mytable = DataFrame(Dict("delta"=> vec(delta_values), "collateral"=> vec(collateral_values),"S"=>social_surplus,"L"=>vec(liquidations),"price"=> vec(prices)))
CSV.write("simulation_results_homogenous_delta.csv",mytable)
mytable = DataFrame(Dict("delta"=> vec(delta_values), "collateral"=> vec(collateral_values),"S5" => social_surplus2,"S15"=> social_surplus3,"L5"=>vec(liquidations2),"price5"=> vec(prices2),"L15"=> vec(liquidations3),"price15"=>vec(prices3)))
CSV.write("simulation_results_heterogenous_delta.csv",mytable)
mytable = DataFrame(Dict("delta"=> vec(delta_values_extended), "collateral"=> vec(collateral_values_extended),"S5"=> social_surplus2e,"L5"=>vec(liquidations2e),"price5"=> vec(prices2e)))
CSV.write("simulation_results_extended_heterogenous_delta.csv",mytable)
"simulation_results_extended_heterogenous_delta.csv"
# read results in without having to rerun code
sim_results = CSV.read("simulation_results_homogenous_delta.csv",DataFrame);
social_surplus =sim_results.S;
liquidations = sim_results.L;
prices = sim_results.price;
sim_results2 = CSV.read("simulation_results_heterogenous_delta.csv",DataFrame);
social_surplus2 =sim_results2.S5;
liquidations2 = sim_results2.L5;
prices2 = sim_results2.price5;
social_surplus3 =sim_results2.S15;
liquidations3 = sim_results2.L15;
prices3 = sim_results2.price15;
sim_results2e = CSV.read("simulation_results_extended_heterogenous_delta.csv",DataFrame);
social_surplus2e =sim_results2e.S;
liquidations2e = sim_results2e.L5;
prices2e = sim_results2e.price5;
# Plot social surplus homogeneous case
homo_ss = Plots.contourf(delta_range, collateral_range, social_surplus, lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",title="Homogeneous Delta-Network Social Surplus")
# Plot prices homogeneous case
homo_price = Plots.contourf(delta_range, collateral_range, prices, c= col2, lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",title="Homogeneous Delta-Network Price")
hetero5_ss = Plots.contourf(delta_range_extended, collateral_range, social_surplus2e,lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",title="Shocked 5-agent component: Social Surplus")
hetero5_p = Plots.contourf(delta_range_extended, collateral_range, prices2e, c=col2,lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",title="Shocked 5-agent component: Price")
hetero15_ss = Plots.contourf(delta_range, collateral_range, social_surplus3,lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",title="Shocked 15-agent component: Social Surplus")
hetero15_p = Plots.contourf(delta_range, collateral_range, prices3, c=col2,lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",title="Shocked 15-agent component: Price")
# Save graphs: ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12
Plots.contourf(delta_range, collateral_range, social_surplus, lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",cbar=false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/homogeneous_case_ss.png")
Plots.contourf(delta_range, collateral_range, prices, c=col2,lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",cbar=false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/homogeneous_case_price.png")
Plots.contourf(delta_range_extended, collateral_range,social_surplus2e, lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",cbar=false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous_case_ss5.png")
Plots.contourf(delta_range_extended, collateral_range, prices2e,c=col2, lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",cbar=false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous_case_price5.png")
Plots.contourf(delta_range, collateral_range, social_surplus3, lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",cbar=false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous_case_ss15.png")
Plots.contourf(delta_range, collateral_range, prices3, c=col2,lw=0,xlabel = "Intercomponent Exposure (%)", ylabel = "Collateral Ratio",cbar=false,ylabelfontsize = 14, xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous_case_price15.png")
"/Users/gracechuan/Documents/network contagion/delta_connected/heterogeneous_case_price15.png"
Note: $c_1$, $c_2$, and $c_3$ are defined differently here in the code than in the paper. In the code, $c_1$ refers to the collateral ratio consistent across loans in the 5-agent component, $c_2$ refers to the collateral ratio consistent across loans within the 15-agent component, and $c_3$ refers to the collateral ratio consistent across loans between the components.
function hetero_coll_experiment(r,c1,c2,c3,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega,h0)
# create collateral network
C = hetero_collat_network(r,c1,c2,c3);
# create Debt network
D = hetero_delta_network(Dweight,n,r,delta,1);
# calculate equilibrium
tot_liq,def,p, ss = equ(s,n,e0,epsilon,xi,zeta,C,D,omega,h0);
return tot_liq,ss, p
end
function hetero_collat_network(r,c1,c2,c3)
# little cluster
m1 = ones(r,r)*c1
# intercomponent loans
m2 = ones(n-r,r)*c3
m3 = ones(r,n-r)*c3
# big cluster
m4 = ones(n-r,n-r)*c2
C = [m1 m3 ; m2 m4]
return C
end
hetero_collat_network (generic function with 1 method)
c_values = [0.2,0.3,0.4,0.5, 0.55,0.6, 0.7, 0.8,0.9,1]
num_c = length(c_values)
num_cells = num_c^2
100
# shock the 5 agent component
price_matrix = ones(num_c,num_c);
ss_matrix = ones(num_c,num_c)
delta = find_delta_star(n,r,Dweight,0.55,s,e0,epsilon_star,xi,zeta,omega,h0)
ind = 1
for i in 1:length(c_values)
for j in 1:length(c_values)
c1 = c_values[i]
c3 = c_values[j]
# calculate utilities and liquidations
tot_liq,ss,p =hetero_coll_experiment(r,c1,0.55,c3,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega,h0)
ss_matrix[i,j] = ss;
price_matrix[i,j] = p;
# add to indicator
ind+=1
end
end
# shock the 15 agent component
price_matrix2 = ones(num_c,num_c);
ss_matrix2 = ones(num_c,num_c)
omega2 = zeros(n,1)
omega2[n] = 1
delta = find_delta_star(n,r,Dweight,0.55,s,e0,epsilon_star,xi,zeta,omega2,h0)
ind = 1
for i in 1:length(c_values)
for j in 1:length(c_values)
c2 = c_values[i]
c3 = c_values[j]
# calculate utilities and liquidations
tot_liq,ss,p =hetero_coll_experiment(r,0.4,c2,c3,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega2,h0)
ss_matrix2[i,j] = ss;
price_matrix2[i,j] = p;
# add to indicator
ind+=1
end
end
Plots.heatmap(c_values,c_values,ss_matrix2,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",title = "Shocked 15-agent component: Social Surplus")
Plots.heatmap(c_values,c_values,price_matrix2,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",title = "Shocked 15-agent component: Price")
Plots.heatmap(c_values,c_values,ss_matrix,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",title = "Shocked 5-agent component: Social Surplus")
Plots.heatmap(c_values,c_values,price_matrix,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",title = "Shocked 5-agent component: Price")
Plots.heatmap(c_values,c_values,ss_matrix2,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",legend=false,ylabelfontsize = 14,xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous/heterogeneous_col_ss15.png")
Plots.heatmap(c_values,c_values,price_matrix2,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",legend=false,ylabelfontsize = 14,xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous/heterogeneous_col_p15.png")
Plots.heatmap(c_values,c_values,ss_matrix,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",legend=false,ylabelfontsize = 14,xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous/heterogeneous_col_ss5.png")
Plots.heatmap(c_values,c_values,price_matrix,xlabel = "Intercomponent Collateral Ratio", ylabel = "Within Collateral Ratio",legend=false,ylabelfontsize = 14,xlabelfontsize = 14,tickfontsize=12,dpi=1000)
savefig("delta_connected/heterogeneous/heterogeneous_col_p5.png")
"/Users/gracechuan/Documents/network contagion/delta_connected/heterogeneous/heterogeneous_col_p5.png"
function comparative_statics(r,c_set,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega,h0)
c1_values = []
c2_values = []
c3_values = []
results = []
results2 = []
results3 = []
for c in c_set
c1 = c[1]
c2 = c[2]
c3 = c[3]
# run experiment
liq,ss,p = hetero_coll_experiment(r,c1,c2,c3,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega,h0)
# add to lists
push!(c1_values,c1)
push!(c2_values,c2)
push!(c3_values,c3)
push!(results,liq)
push!(results2,ss)
push!(results3,p)
end
return DataFrame(Dict("5-agent Component"=>c1_values, "15-agent Component"=> c2_values,"Intercomponent"=>c3_values,"liqd"=> results,"social_surplus"=> results2,"price" => results3))
end
comparative_statics (generic function with 1 method)
# 5 agent shocked case
c_combos = samplewithreplacement(c_values,3)
delta = find_delta_star(n,r,Dweight,0.55,s,e0,epsilon_star,xi,zeta,omega,h0)
table = comparative_statics(r,c_combos,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega,h0)
CSV.write("simulation_results_heterogenous_collateral5.csv",table)
"simulation_results_heterogenous_collateral5.csv"
# 15 agent shocked case
omega2 = zeros(n,1)
omega2[n] = 1
delta = find_delta_star(n,r,Dweight,0.55,s,e0,epsilon_star,xi,zeta,omega2,h0)
table2 = comparative_statics(r,c_combos,Dweight,delta,s,n,e0,epsilon,xi,zeta,omega2,h0)
CSV.write("simulation_results_heterogenous_collateral15.csv",table2)
"simulation_results_heterogenous_collateral15.csv"